A major question regarding limb development concerns how the number and positioning of skeletal elements along the proximodistal (P-D) and anteroposterior (A-P) axes is determined [4]. There is a general P-D increase in the number of skeletal elements which occurs even in cases, such as the chicken forelimb (Fig. 1), in which the A-P width remains essentially constant while the skeletal pattern is being laid out. In fish, non-tetrapod vertebrates, the fin endoskeleton is a mixture of bars and nodules which have no discernable P-D numerical trend [5].

Cartilage differentiation, or chondrogenesis, is preceded by “condensation” of the precartilage mesenchyme, in which cell density increases and the cells enter into broad, transient contact with one another [6]. The precartilage cells are embedded in a dilute extracellular matrix (ECM) and condensation is accompanied by and dependent on local accumulation of the ECM molecule fibronectin [7], with markers of prospective condensation appearing earlier than ECM and morphological changes [8]. Molecules secreted by the dorsal and ventral ectoderm, including FGFs and Wnt, inhibit chondrogenesis [8], [9], thereby confining the developing one-bar proximal cartilage primordium (stylopod, i.e., humerus, femur) to a central planar sector of the paddle-shaped limb bud [4]. As development proceeds, the skeleton remains confined to this plane but expands laterally in more distal regions as the stylopod gives way to the two-bar (zeugopod, i.e., radius and ulna, tibia and fibula) and multiple-bar (autopod, i.e., digits) primordia of the mid and terminal regions of the limb (Fig. 1C). This brings the developing skeletal elements increasingly closer to the anterior and posterior edges of the limb bud, reflecting attenuation of the peripheral inhibitory effect. Attenuation of inhibition can also be seen in the proximity of the more distal elements to the dorsal and ventral surfaces as the limb bud tapers towards its tip, and to the apical boundary, as the potency of the AER wanes [10].

There is a large literature that treats limb pattern formation as a question of establishing spatiotemporal informational fields that cells subsequently interpret according to their genetic programs [11], [12]. Other approaches to this problem (and those of other developing systems), take into account the conditional, self-organizational properties of interacting cells and tissues (see refs. [13] and [14] for reviews). This latter approach has the advantage over the “positional information” perspective in its ability to formulate hypotheses for why the limbs (or other embryonic structures) assume the forms they actually do.

A mechanism that has gained increasing attention as a generator of spot- and stripe-like patterns in developmental systems is related to the chemical reaction-diffusion process described by Turing in 1952 [15], [16]. This has been schematized in a biological “local autoactivation-lateral inhibition” (LALI) framework by Meinhardt and Gierer [17]. LALI systems, while formally similar to reaction-diffusion systems, are more suitable to biological applications, since they include cases where the activator and inhibitor terms are mediated by cellular “reactors” rather than simple chemical reactions [18], and spatial transport can be mediated by mechanisms in addition to simple diffusion [19]. LALI systems are self-organizing: that is, with appropriate choices of activator-inhibitor and transport parameters, in domains of appropriate size, shape and boundary conditions, a spatially uniform distribution of the morphogens becomes dynamically unstable, giving way to nonuniform distributions, typically periodic arrangements of stripes and spots in two dimensions, and their three-dimensional analogs (bars, nodules) [13], [20].

The developing limb bud contains all the components of a skeletogenic LALI mechanism [21], [22], [23]. The locus of activity of this process is a population of multipotent skeletogenic progenitor cells (termed “precartilage” cells in what follows) [24], [25], [26], which are maintained in a non-terminally differentiated state in the distal-most 0.3 mm of the avian and mammalian limb bud by the action of the apical ectodermal ridge (AER), a narrow A-P-oriented ectodermal thickening, along the limb bud tip [27]. Whereas the precartilage cells can give rise to cartilage, connective tissue, and bone, and can undergo programmed cell death, the limb musculature is formed by a separate population of cells [26]. The non-myogenic limb mesenchyme thus forms the skeletal pattern by being caused to choose among its potential fates of cartilage, soft connective tissue fibroblasts and apoptotic cells, in a spatiotemporally regulated fashion.

Precartilage condensation and subsequent chondrogenesis are promoted by morphogens of the TGF-β superfamily [28], [29] that form an “activator subnetwork” (Fig. 1B). In brief, one or more TGF-βs and activins set off a train of events whereby BMPs (one or more of BMP2, 4 and 7), acting via the receptor BMPR-Ia, induce regions of high BMP signaling activity, marked by phosphorylated Smads [30]. Other classes of molecules, such as galactose-binding lectins (“galectins”) [31], may also be involved in this initiation step. Signaling by the activator subnetwork induces extracellular matrix, matricellular and cell adhesion molecules, promoting mesenchymal condensation and subsequent chondrogenesis [32], [33] (Fig. 1B).

The spatial profile of BMP action is dictated not by its receptor distribution, which is uniform, but by a “prepattern” (i.e., a molecular pattern prefiguring the cellular one) of the diffusible morphogen itself, which in the autopod takes the form of narrow crescents at the tips of the forming digits [30]. This prepattern arises from the spatiotemporal dynamics of the activator network morphogens interacting with inhibitory factors. The latter include the extracellular antagonistic BMP binding partners Noggin, Gremlin and Ventroptin (CHL2), the BMP receptor antagonist BAMBI and the antagonistic intracellular co-receptor Smad6, the expression of which are all induced by the activator subnetwork [30], [34]. Consistent with this, the prospective digits and interdigits are uncommitted as to their fates late into autopod formation [30], and this plasticity likely obtains in the fields of precartilage cells that form the stylopod and zeugopod at earlier stages. In particular, exogenous TGF-β can cause noncondensing limb bud precartilage mesenchyme to condense and form ectopic cartilage in vitro [28] and in vivo [34].

The FGF pathway [8] and the Notch pathway [35], [36] are involved in restricting the expansion of precartilage condensations as they form. Along with the inhibitory factors described above, these perinodular chondrogenesis-restricting components constitute what we refer to as the “inhibitor subnetwork” (I) (Fig. 1B).

The activator (A) and inhibitor (I) subnetworks together constitute a core network for the regulation of chondrogenesis (Fig. 1B). All cells of the precartilage mesenchyme appear to be equally capable of producing all components of A and I. Subnetwork A functions in a paracrine and (in part via its TGF-β component [29]), a positively autoregulatory fashion. It induces the production of cell-to-cell transmissible components of subnetwork I, which in turn limit subnetwork A's radius of action (Fig. 1B).

Limb bud mesenchyme can in fact self-organize into cartilage patterns with limb-like features. Dissociated and reaggregated limb mesoderm cells, packed into an ectodermal hull, form limb-like arrays of rod-like and nodular skeletal elements [37]. In high-density culture these cells form spots and stripes of cartilage with dimensions and spacing similar to the in vivo patterns [38]. When the physical properties of the culture microenvironment [39] or putative LALI components are manipulated in vitro and vivo, the pattern changes are predictable from LALI models [8], [28], [29], [30].

In most amniote species, emergence of the cartilage primordia of the limb skeleton occurs first in proximal regions of the limb bud and then in increasingly more distal regions [40] (Figs. 1; ​;2,2, left). This process is dependent on an intact AER, the removal of which leads to truncations of the skeleton with progressively more distal elements forming the longer the AER remains in place [40]). The growth of the skeletal structures (and portions thereof) that are spared after AER removal is unaffected, with these components attaining nearly their normal size and shapes, except for terminal deletions, as development proceeds [40]. The limbs of certain amphibians seem to represent partial exceptions to the proximal-to-distal rule [41].

Although the idea that a LALI mechanism could account for major features of limb development in its normal, genetically and experimentally perturbed, and evolutionary aspects is several decades old [23], what has been lacking for most of the succeeding period has been a fast computational modeling environment that would enable simulation of a LALI system as a continuous (“real-time”) process, in domains with natural contours and changing geometries. There has been progress in this direction more recently, using finite element algorithms, producing modeling environments suitable for investigating time- and shape-dependent changes in patterns generated by reaction-diffusion systems of biological interest [42], [43]. Here we use one such framework to test the capability of a LALI limb model to generate realistic patterns under normal and biologically altered conditions.

All of these simulations employed model limbs with curved apical contours, which are truer to the living limb shapes than the simple geometries generally used in reaction-diffusion simulations. The ability to perform such simulations in domains with natural contours and continuously changing nonstandard shapes was enabled by our previous development of new finite element computational methods which are applicable to any biological system in which fields of chemicals or other mobile agents are produced by activator-inhibitor interactions on domains of changing size and shape [42].

By our main hypothesis, limb skeletal patterning is established by a LALI process within a reshaping zone of tissue in the presence of a distal suppressive gradient. While the skeletal patterns in all the limbs we have simulated are mediated by the same self-organizing system, the changing limb bud shapes within which this system operates were imposed arbitrarily according to the schemes described in Files S3, S4, S5. Since limb bud shape is controlled by a different set of molecular determinants from the ones regulating the initiation of chondrogenesis (reviewed in ref. [2]), a more complete model would incorporate an independent system of equations for limb bud shaping [64], [65] involving the growth [66] and viscoelastic [67] properties of the limb bud tissue regulated by molecules such as Wnt, Shh, Gli3, and FGFs [68], [69].

In our model, an increase in the number of parallel elements occurs when either the P-D length of the LALI zone decreases ([45]; first figure of File S3) or its A-P width increases (Fig. 6). While certain other conditions may lead to peak splitting in activator-inhibitor systems [22], [70], we have seen only occasional examples of intercalation of new peaks between preexisting ones with parameter change (e.g., second figure of File S3). Since such patterning modes characterize certain variant limb types [22], this may indicate a limitation of our representation that would improve with the introduction of additional modulatory parameters corresponding to the molecular complexity of the biological system [3], [4], [11].

Joints can form in our model as a result of discontinuities in the pattern of elements along the AP axis as different spatial solutions become stable as the LALI zone changes in size and shape (e.g., Fig. 2; Movie S1). Oscillations are well known to occur in reaction-diffusion system in certain parameter regimes [71], including our morphostatic system (H.G.E. Hentschel and S.A. Newman, unpublished data). Although such oscillations would not generally serve to segment stripe-like LALI pattern elements [70], the elongated elements in our system are typically formed by temporal persistence of spot-like patterns rather than as de novo stripes (e.g., Brachypterygius simulation in File S5), making oscillatory modes plausible factors in joint formation. In embryonic limbs the generation of joints depends on members of the BMP superfamily beyond those specified in our core mechanism, such as GDF5 [72]. These factors have activating and inhibitory effects on chondrogenesis that would result in more complex spatiotemporal waves than those seen in the basic model.

In all the simulations shown we have varied the values of the reaction kinetic parameters and in a P-D level-specific fashion; that is, different values were used for the stylopod, zeugopod and autopod (see Files S3, S4, S5 for details). Although we determined previously that the contraction of the P-D length of the LALI zone (a consequence of the attenuation of the AER suppressive gradient), was by itself sufficient to increase in the number of parallel skeletal elements [45], the variation of and (which does not change the network topology of the core patterning system nor impart any element-specific positional information to the model limb), fine tuned the auto- and cross-regulatory interactions between the morphogens and led, with appropriate choices to relatively authentic skeletal patterns.

The parameters and were defined by a detailed mathematical analysis [45] in which factors relating to cell movement and extracellular matrix production were folded into functions (U and V) governing the activating and inhibitory morphogens operating within a system encompassing a fuller range of biological phenomena [18]. This simplification was necessitated by limitations on the complexity of systems that can be simulated by available finite element algorithms [42]. While attribution of molecular functions to the parameters in the original eight-equation system is more straightforward, it is nonetheless still possible to discern the roles of and in the morphostatic system (1) and thereby attempt to put a “molecular face” on these parameters.

Both parameters appear in the production rates of the activator and inhibitor morphogens, and is also the association rate constant of the activator and the inhibitor. The parameter denotes the activator morphogen concentration at which there is a transition from a linear response to a saturation phase. In particular, the Turing bifurcation (a transition in the behavioral characteristics of the system) that enables pattern formation only occurs in a narrow range of values [45]. This constraint is in evidence in all the simulations shown here (with displacement of ±2% around its mean value significantly affecting the pattern), and can be considered a required relationship for pattern formation among the system's activator and inhibitory components under our biological assumptions and their mathematical representations.

The limb deformity (ld) locus in the mouse encodes several functionalities corresponding to the role of in our model. Formin1, disruption of which leads to the absence of the fibula and reduction in digit number, probably regulates expression of BMPs [73], components of the activator subnetwork (Fig. 1b), whereas a transcriptional global control region (GCR) in this domain activates the limb-specific expression of the BMP antagonist Gremlin [74], a component of the inhibitor subnetwork. Not included in the present model, however, are the buffering systems accumulated over 400 million years of evolution that protect the biological equivalent of from deviating from its prescribed range of values in present-day limbs.

In contrast to , the value of can vary extensively without compromising the pattern-forming capacity of the system. Variations in lead to extensive changes in the number and arrangement of skeletal elements, particularly in the autopod, the most variable region of normal, mutant and fossil limbs. As suggested above, the behavior of the system under parameter variation may provide insight into the elusive functions of Hox transcription factors in the developing limb. Different members of the Hoxa family vary in abundance in the limb tip in a stage-dependent fashion [57], and like the parameter , influence activator dynamics. Hoxa13, for example, which is the most highly expressed member in the prospective autopod, alters the expression in the developing limb of two morphogens of the TGF-β superfamily, members of which comprise the activator subnetwork (Fig. 1b) [75]. Hoxa9, abundantly expressed in more proximal regions, modulates TGF-β superfamily signaling [76]. We can thus tentatively identify the role of in our scheme with certain functions of the Hoxa gene products.

A different set of Hox gene products, the Hoxd class, vary along the limb bud A-P axis during development [57]. A more detailed model could potentially account for the morphological differences between the radius and ulna, or the various fingers, by parameter variations in this direction as well. As with the dynamical network for limb bud shaping which is not yet included in our framework (see above), equations for setting parameter values in a biological fashion, ultimately based on the Hox gene regulatory network [58], [77], would fill an important gap in the model. It should be noted, however, that such shape and parameter regulatory systems presuppose a core chondrogenic mechanism (Fig. 1B) and a bare bones pattern generating scheme (Fig. 1C), as described above. They do not by themselves provide positional information to the cells of the developing limb.

Despite the extensive sampling of parameter space represented by the simulation results shown, and scores of additional simulations with other parameter combinations, some in which the parameter varied by more than 10-fold (Files S3, S4, S5), it was difficult to obtain an outcome that did not resemble a limb skeleton. This in no way means that any arbitrary pattern could be produced: the simulated skeletons are all composed of spot- and stripe-like elements. The LALI system (1), in the context in which we have simulated it, is therefore inherently “skeletogenic.” It appears, moreover, that the topology of the core network, rather than the specific identity of the relevant gene products (of which our knowledge is incomplete, and for which there is variation between tetrapod species), may be decisive for this process.

A recent study has provided evidence that Tiktaalik [78] and other vertebrates known from the fossil record to have noncanonical bony limb skeletons previously thought to be transitional to definitive tetrapod limbs actually coexisted with early tetrapods [79]. The drastically narrowed the time span during which evolution of all known limb morphotypes must have occurred calls for a skeletogenic mechanism with a propensity to generate a profusion of patterns due to small genetic changes affecting limb bud shaping and the rates and strengths of core interactions. The mechanism we propose is the only empirically based one currently under consideration that has these properties.

File S1

File S2

File S3

File S4

File S5

Movie S1

Movie S2

Movie S3