The results presented here confirm that a reaction-diffusion process operating in a reshaping domain in the context of a distal suppressive gradient can reproduce major aspects of skeletal patterning in vertebrate limbs. In particular, the proximodistal temporal sequence of development of skeletal elements (), the proximodistally increasing number of elements in fully developed limbs () and the dependence of development on the suppressive FGF gradient (i.e., the AER), with distal truncation occurring upon its removal (), are all readily accounted for by this mechanism. The “bare bones” patterning mechanism also produces increased numbers of elements when the apical domain is expanded in an anteroposterior direction during development, as occurs with ectopic grafts of the ZPA, with double knockouts of Shh and Gli3, development of the chicken talpid2 mutant or the dogfish ().
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 
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
. Since limb bud shape is controlled by a different set of molecular determinants from the ones regulating the initiation of chondrogenesis (reviewed in ref. 
), a more complete model would incorporate an independent system of equations for limb bud shaping 
involving the growth 
and viscoelastic 
properties of the limb bud tissue regulated by molecules such as Wnt, Shh, Gli3, and FGFs 
In our model, an increase in the number of parallel elements occurs when either the P-D length of the LALI zone decreases (
; first figure of File S3
) or its A-P width increases (). While certain other conditions may lead to peak splitting in activator-inhibitor systems 
, 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 
, 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 
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., ; Movie S1
). Oscillations are well known to occur in reaction-diffusion system in certain parameter regimes 
, 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 
, 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 
. 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
in a P-D level-specific fashion; that is, different values were used for the stylopod, zeugopod and autopod (see Files S3
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 
, the variation of
(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.
were defined by a detailed mathematical analysis 
in which factors relating to cell movement and extracellular matrix production were folded into functions (U
) governing the activating and inhibitory morphogens operating within a system encompassing a fuller range of biological phenomena 
. This simplification was necessitated by limitations on the complexity of systems that can be simulated by available finite element algorithms 
. 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
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
. 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 
, components of the activator subnetwork (), whereas a transcriptional global control region (GCR) in this domain activates the limb-specific expression of the BMP antagonist Gremlin 
, 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 
, 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 () 
. Hoxa9, abundantly expressed in more proximal regions, modulates TGF-β superfamily signaling 
. 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 
. 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 
, 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 () and a bare bones pattern generating scheme (), 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
), 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 
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 
. 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.
While our model provides a plausible account of the general form of the limb skeletal pattern (a capability absent in competing models), its main role at present is as a framework for further experimental tests. In particular, it will be important to gain knowledge of the earliest acting activators and inhibitors of precartilage condensation, the details of their transmission and interaction, and the role of FGF or other signals from the AER in defining the morphogenetically active region of the developing limb.