{'fraction_non_alphanumeric': 0.024767984867293254, 'fraction_numerical': 0.019270556245197138, 'mean_word_length': 4.635980011104942, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Digit patterning during limb development as a result of the BMP-receptor interaction

Turing models have been proposed to explain the emergence of digits during limb development. However, so far the molecular components that would give rise to Turing patterns are elusive. We have recently shown that a particular type of receptor-ligand interaction can give rise to Schnakenberg-type Turing patterns, which reproduce patterning during lung and kidney branching morphogenesis. Recent knockout experiments have identified Smad4 as a key protein in digit patterning. We show here that the BMP-receptor interaction meets the conditions for a Schnakenberg-type Turing pattern, and that the resulting model reproduces available wildtype and mutant data on the expression patterns of BMP, its receptor, and Fgfs in the apical ectodermal ridge (AER) when solved on a realistic 2D domain that we extracted from limb bud images of E11.5 mouse embryos. We propose that receptor-ligand-based mechanisms serve as a molecular basis for the emergence of Turing patterns in many developing tissues.


Introduction
The limb has long served as a paradigm to study organ development. Its easy access enabled developmental biologists to use tissue grafting to define organizing centers, and advanced genetics has led to the identification of the key regulatory circuits 1,2 . Early grafting experiments showed that tissue from the posterior limb bud (referred to as zone of polarizing activity, ZPA) induced ectopic digits when implanted in the anterior limb bud 3 . These extra digits formed a mirror image of the normal forelimb, i.e. 1234554321 in mice and 234432 in chicken. Here numbers refer to digits with the lowest number referring to the most anterior digit in the wildtype (i.e. the thumb in humans). The effect was shown to be concentration-dependent in that a smaller number of grafted cells induced fewer digits and failed to induce the most posterior identities in the anterior limb bud 3 . Wolpert's French Flag model explained these experiments by suggesting that digit number and identity were determined by the local concentration of a chemical compound 4 . This compound was predicted to be produced in the ZPA and to diffuse across the limb bud, thereby creating a concentration gradient. Decades after the grafting experiments had suggested that patterning was induced by a diffusible compound, Sonic hedgehog (SHH) was isolated and shown to affect patterning in a concentration-dependent manner 5,6 . In agreement with the morphogen model, SHH is expressed in cells of the ZPA and diffuses across the limb bud 7 . In mice that do not express SHH only digit 1 develops 8 , and in the polydactylous mouse mutant extra toes (Xt) SHH is expressed also in the anterior region of the limb bud 9 .
However, several other experiments question this simple morphogen model and suggest that a more complex mechanism determines digit patterning during limb bud development. For one, polydactyly is still observed even when SHH expression is abolished in the extra toe (Xt) mutant 10 . Moreover, an experiment in which SHH expression is stopped prematurely at different stages during limb development results in an alternating anterior-posterior sequence of digit formation that cannot be reconciled with the standard morphogen model 11 . Thus the shortest pulse that yielded an additional digit induced formation of digit 4, followed by digit 2, and then digits 3 and 5. Based on the morphogen model digits would be expected to emerge simultaneously or in an anterior-posterior sequence.
Experiments in which SHH-expressing cells were marked revealed that descendants of the SHH-producing cells could be found far beyond the small posterior zone in which SHH expression can be detected 7,12 . Descendants of the SHH-producing cells filled almost half of the limb bud and encompassed all cells in the two most posterior digits and also contributed to the middle digit 7,12 . Based on this observation it was suggested that it is the length for which cells express and secrete SHH (and therefore experience a high SHH concentration) that determines digit specification 7,12 . This model, however, still cannot explain the observed sequence of digit formation.
A number of mathematical models have been developed to explain limb bud development and digit patterning. Based on experimental observations, the earliest model defined some 40 years ago a set of simple rules to recapitulate the observed changes in the limb buds' size and shape 13 . The modelling effort highlighted the importance of oriented cell division to obtain the elongated shape that is characteristic of the limb bud. Recent 3D imaging data in combination with a Navier-stokes description of the tissue as an viscous fluid indeed confirmed that isotropic proliferation alone is insufficient to explain the shape and size changes during limb bud development 14 . Local gain rates much higher than the measured proliferation rates as well as loss rates had to be included; these may reflect oriented cell division and directed cell migration. The description of the limb tissue as a viscous fluid was first introduced by Dillon and Othmer who studied a model where the local growth rate (the source term in the Navier stokes equation) depended on the concentration of two signaling factors, SHH and FGF 15 . A more recent model for limb bud growth explores the impact of differential tissue elasticity 16 . Other models have focused on particular regulatory interactions that pattern the limb expression domains, i.e. the AER-ZPA regulatory interaction 17 , the establishment of the proximal and distal signaling centers 18 , and the diffusivity of SHH 19 .
Most models, however, have focused on the skeletal patterns that emerge during limb bud development [20][21][22][23][24][25][26][27][28][29] . A number of different model types have been studied ranging from continuous to cell-based models and from rule-based stochastic models to reaction-diffusion models; for reviews see 30,31 . Turing mechanisms are a particular type of reaction-diffusion models 32 and have been frequently invoked to explain the emergence of pattern during development 33 ; they have been studied in the context of limb patterning since 1979 34 . In Turing models the number of modes (patterns) increases with larger domain size. Newman and co-workers showed that the different number of skeletal elements in the stylopod, zeugopod, and autopod in anteriorposterior direction can in principle be achieved as a result of the increase in the proximal-distal length; here they have to assume that patterns emerge sequential in stylopod, zeugopod, and autopod and are fixed (frozen) upon their emergence 23 . Newman and co-workers further showed that parameter values can be identified that give realistic steady state pattern if the equations are solved on static domains whose shapes are based on the limb shapes of various fossils 35 . However, there is no evidence that the patterning domain in the embryo has the same shape as the adult fossil and patterning is well known to occur on a growing rather than static domain. Miura and co-workers explored how the speed of pattern emergence could be affected 24 and studied Turing patterning on growing limb domains 22 . In particular, Miura and co-workers show that a Turing model can explain the increased number of digits in the doublefoot mutant with the observed increased domain size 22 .
The many computational studies show that Turing models, in principle, are sufficiently flexible to account for the experimentally observed pattern in wildtype and mutant limb buds. However, in spite of the demonstrated sensitivity of Turing models to domain size and particular reaction types and kinetics the models have so far not been solved on realistic domains with biochemically realistic reaction kinetics. Recent efforts in this direction have been made by the Newman group who propose TGF-beta family ligands as the activator and FGFs, Noggin, Notch and/or CHL2 as the inhibitor in their Turing models with Schnakenberg kinetics 21,23,35 . However, mutant mice with defects in TGF-beta 36 , FGF 37-39 , Noggin 40 , and Notch 41 signaling all have digits, while CHL2 expression has been found to be restricted to chondrocytes of various developing joint cartilage surfaces and connective tissues in reproductive organs 42 . While it is possible that the apparent robustness of the patterning mechanism to these mutations is a consequence of redundant regulatory interactions, we wondered whether the reported regulatory interactions would allow us to construct a model that would be consistent with available experimental data and that would not have to rest on missing data and redundancy arguments.
Advances in experimental techniques now provide us with detailed information regarding the domain geometries 14 , the timing of the processes, and the regulatory interactions 1,2,43 . Limb buds first form around embryonic day (E)9.5 while expression of Sox9, an early marker for tissue condensations, is first observed around E10.5 in the  central mesenchyme of the limb bud (Fig. 1A 44 ). Around E11 the Sox9 expression pattern splits and a region in the stalk patterns the stylopod, two patches mark the zeugopod, and a region of strong expression marks the autopod (Fig. 1B,F, 44,45 ). Around E11.5 the autopod expression pattern splits into distinct patches that mark the different digits (Fig. 1C,G, 44,45 ). Over the next day more defined and elongated Sox9 expression pattern emerge that resemble the final digit structure (Fig. 1D,E,H 44,45 ). Sox9 expression has been shown to be stimulated by BMP-2 46 and the BMP-2 dependent increase in Sox9 then triggers increased expression of Noggin 47 which acts as an antagonist of BMP signaling by sequestering BMP in inactive complexes 48 . Noggin null mice form digits, but the digit condensations are much wider than in wildtype 40 . Much as Sox9, Noggin is expressed in digit condensations and serves as a marker of endochondral differentiation 40 .
Even though BMPs also induce the expression of the BMP antagonist Gremlin 49 , Gremlin (unlike Noggin) is expressed in the interdigital space rather than in the digit condensations 50 . Another study notes a concentration-dependent regulation of Gremlin expression in that low concentrations of BMP2 upregulate Gremlin while high concentrations of BMP2 downregulate Gremlin in limb mesenchyme cultures 51 . It is thus possible that stronger BMP signaling in the digit condensations induces Sox9 while weaker BMP signaling in the interdigital space induces Gremlin. Similarly, it has been suggested that low BMP concentrations support expression of FGF8 while high BMP concentrations repress FGF4 and FGF8 expression in the apical ectodermal ridge (AER) 52,53 . Two different BMP receptors are expressed in the limb mesenchyme, BMPR1a (ALK3) and BMPR1b (ALK6) that bind the different BMP ligands with different affinities. BMPR1a binds mainly BMP4 while BMPR1b binds BMP4, BMP2 and less well to BMP7 54, 55 . The BMP4-BMPR1b affinity has been established in a Scatchard plot as Kd 5 530 pM 55 . Activin receptors (ActRI/ALK2 and ActRII) are also expressed in the limb mesenchyme 56 but appear to bind their BMP ligands with lower affinity 57 . BMP7 binds ActRII with 3-fold higher affinity compared to BMP2 (Kd 5 1.7 mM versus 5.4 mM) and binds ActRI/ALK2 only with low affinity (Kd 5 143 mM) 57 . It should be noted that the values are based on Biacore measurements and thus do not adequately reflect the cooperative binding effects.
In both mouse and chicken expression of BMPR1b is restricted to precartilaginous condensations 58,59 while Bmpr1a is expressed at low levels throughout the limb bud mesenchyme 59,60 . In the chicken Bmpr1a expression is elevated at the border between precartilaginous condensations and mesenchyme, but low or absent in the precartilaginous condensations. Accordingly, overexpression of dominant negative (DN)-BRK-2 (the chicken receptor that corresponds to the mouse BMPR1B receptor) blocks chondrogenesis, but not overexpression of (DN)-BRK-1 (which corresponds to mouse BMPR1A) or DN-ActR1 60,61 . BMPR1B thus appears to be the major transducer of BMP signals in chicken limb condensations.
In the chicken the two receptors have been shown to be engaged in different regulatory interactions. Thus BMP beads induce noggin, tgfb2, and bmpr1b expression only when implanted into the tip of the chicken digit, but not when implanted in the interdigital space. Also bmpr1a is not expressed when BMP beads are implanted in the interdigital space while expression of the bmpR-1b gene was enhanced in the adjacent digits with the spatial distribution appearing displaced toward the margin of the digits adjacent to the bead 62 . It therefore seems that in the chicken only BMPR1b (but not BMPR1a) enhances its own expression in response to ligand binding. TGF-beta beads also induce the expression of noggin in the interdigital regions but only with a time delay: noggin expression is detected after 30 hours, while bmpr1b expression becomes detectable in the interdigital region already after 10 hours. TGF-beta thus most likely induces noggin expression indirectly by inducing bmpr1b expression 62 .
There are important differences between mouse and chicken limb development as there are differences between forelimb and hindlimb development. We will therefore mainly focus on mouse forelimb bud development and comment on other limb buds as appropriate. In the mouse both Bmpr1a and BmpR1b expression is elevated in the precartilaginous condensations and Bmpr1a and Bmpr1b can compensate for each other 59,60 . BMPR1b null mice show normal digit patterning and no differences in Sox9 expression are observed in wild-type and BmprIB2/2 limbs up to and including E12.5, indicating that the prechondrogenic limb mesenchyme is specified and is able to form condensations up to this time 60,63,64 . BmprIA2/2 embryos are not viable and the severity of the BmprIA2/2 conditional knock-out depends on the stage at which BmprIA2/2 is removed. Cre-mediated excision in Col2-Cre mice seems to occur only in forming cartilaginous condensations, while in Prx1-Cre mice, it acts in limb mesenchyme at much earlier stages (95% active by E.10.5) 63 . If BMP receptor type IA is removed with Col2-Cre digit condensations can still be observed while in Prx1-Cre mice no digits form in the forelimb 63 . The mutants have shortened limbs and show almost complete agenesis of the autopod because of reduced cell proliferation 63 . Moreover, the expression of BMP target genes Msx1/2 and Grem are severely reduced. BMPR1b appears to rescue Msx1/2 expression. However, the expression is then colocalized with the patchy BMPR1B expression 63 . In the mouse BMPR1a thus appears to be the major transducer of BMP signals in limb condensations.
In the limb there are three important BMP ligands: BMP2, BMP4, and BMP7 65 . BMP2 and BMP4 appear to be the most important BMP ligands, but also in the BMP2: BMP4 double conditional knock-out some digits still form 66 . At E12.5 BMP2 is expressed predominantly in the interdigital webbing, BMP4 is expressed along the whole AER, and BMP7 expression is ubiquitous in the limb bud 65 . Bmp2 conditional mutants form five digits, conditional inactivation of a gene for BMP4 results in a polydactylous phenotype 67 , but in BMP2:BMP4 conditional mutants two posterior digits are missing even though the limb field is broader 66 . Bmp7 null mice show an occasional anterior polydactyly 65 . Importantly, even though BMP4 expression is reduced in the conditional BMP4 knock-out BMP signaling activity is increased 68 . To avoid the redundant role of multiple BMP ligands and receptors Zeller and co-workers conditionally removed the downstream signaling protein Smad4 (Co-Smad). The mutants still express Sox9 in the autopod but the pattern does not split into digits 69 . Smad4-dependent signaling is therefore necessary for digit formation.
BMP signaling is embedded in a larger network that most importantly comprises FGFs and SHH in the limb. BMP signaling interferes with FGF-dependent signaling, at least in part by repressing expression of its receptor FGFR1 70 . FGF-loaded beads in turn repress BMPR1b and Noggin expression in chicken limb buds 62 . In mouse limb buds FGFs repress the expression of the BMP antagonists Gremlin 71 and enhance the expression of SHH in the ZPA, which in turn enhances Gremlin expression. SHH induces the expression of BMP-2 and BMP-7 in chicken limb buds, and with a delay that of BRK2 (BMPR1b) 61 . Without SHH digits cannot form unless expression of Gli3 is removed as well 10 . In the absence of SHH signaling Gli3 forms the Gli3 repressor which prevents the expression of Gremlin and many other genes. However, in the absence of both SHH and Gli3 more digits appear than in the wildtype. Neither SHH nor Gli3 are therefore necessary for digit formation. Both Fgf4 and SHH expression terminate around E12 43 but Fgf8 continues to be expressed in the AER.
Given the complex regulatory interactions it is difficult to understand and predict the regulatory outcome by verbal reasoning. We therefore sought to build and analyse a computational model of these interactions. In light of the presence of digits without Shh/Gli3, and the absence of digits in BMPR1A and SMAD4 conditional mutants www.nature.com/scientificreports SCIENTIFIC REPORTS | 2 : 991 | DOI: 10.1038/srep00991 we focus on BMP-dependent signaling to explain the observed patterning of Sox9 in the limb bud. We find that the interactions between BMP and BMP receptor give rise to a Schnakenberg-type Turing mechanism. When we solve the model on a realistic domain we obtain patterns that are similar to those observed in experiments with wildtype and mutant mice. We propose that the observed patterning during limb bud development may be the consequence of a Turing type patterning mechanism that arises from the interactions between BMPs and their receptor.

Results
Model. The core signaling pathways that regulate limb patterning have been well characterized and comprise, SHH, GLI, Gremlin, BMP, FGF, and their receptors 1,72 . Digit patterning is still observed in the absence of Gremlin 73 , Noggin 40 , SHH and GLI3 10 , but not in the absence of BMPR1A 63 and Smad4, an essential protein in BMP signaling 69 . We therefore neglect all dispensable network elements and focus on BMP signaling to explain digit patterning. We will include FGF as an important modulator and growth factor. Shhdependent signaling is also an important modulator of BMP signaling in the limb bud and is clearly important to establish the anterior-posterior polarity. Turing patterns have been shown to persist in the presence of such a modulating gradient 74 . Accordingly, we focus on the mechanism that enables the emergence of digits in the autopod and we will neglect the detailed effects of Shh and how these enable the establishment of anterior-posterior polarity that results in the formation of a thumb and a pinky.
BMP (which we denote by B) and FGF (denoted by F) are secreted proteins and diffuse much faster than the BMP receptor (denoted by R) which resides in the plasma membrane. Experiments on the BMP Drosophila homologue Dpp suggest that the majority of ligandbound receptors, C, are internalised rapidly and reside mainly inside the cell 75 . We will therefore ignore diffusion of the receptor-ligand complex. We write D j (j 5 B, R, F) for the diffusion coefficients with 76,77 . We write D: D : ½ for the diffusion fluxes where D denotes the Laplacian operator in Cartesian coordinates, and concentration. The characteristic length of gradients depends both on the speed of diffusion and their removal. In the absence of contrary experimental evidence we will assume the simplest relation, linear decay, at rates d j [j] for all components (i.e. j 5 B, R, C, F).
BMPs are dimers and one BMP molecule can therefore bind two receptors 78 . The rate of BMP-receptor binding is therefore proportional to R 2 B, and we use k on and k off as the binding and dissociation rate constants. BMP2 signaling has been shown to reduce BMP2 expression in the limb bud 79 and we therefore make the rate of BMP production negatively dependent on the BMP-receptor complex, C, i.e. we write for the BMP production rate P B K B K B z C ½ and thus obtain for the BMP and BMP-receptor dynamics Signaling of BMP-bound receptors positively regulates receptor production 62 , and we therefore must make receptor production dependent on the concentration of C. In the absence of contrary data we will use the simplest possible relation for the receptor production rate and write P R ([C]) 5 p R 1 p C [C] where p R and p C are constants, We can simplify this set of equations if we assume that the dynamics of the receptor-ligand complex are fast such that the receptor-ligand complex assumes a quasi steady-state. The concentration of bound receptor, C, is then proportional to R 2 B, i.e.
Equations (1,3) are sufficient for pattern formation and would reduce to the classical Turing model of Schnakenberg-type if we were to set p C 5 2d C , and d B 5 0. BMP expression and activity is modulated by FGF and Gremlin 68 . BMP induces Gremlin expression and Gremlin then binds and sequesters BMP in an inactive complex 51,80 . Gremlin can in principle be secreted, but experimental evidence suggests that BMP4 activation and secretion are also negatively regulated by an intracellular Gremlin-BMP4 interaction 81 . We therefore do not consider Gremlin diffusion. Since experiments suggest that also Gremlin is a dimer 82  and Hill coefficient n, and assume that dissociation of the complex (at rate k ' off ) is much faster than its degradation then the previous system of equations would need to be expanded to include Experiments show that Gremlin is induced by BMP already early during limb development 68 . If we assume rapid formation of the BMP-Gremlin complex in comparison with other reactions we can use a quasi-steady state approximation for complex concentration: The terms for complex formation thus disappear from the equations for Gremlin and BMP and the equation for Gremlin uncouples from the other equations. Gremlin seems thus not to affect the patterning module directly. We therefore ignore the Gremlin-BMP interaction in this model, and we will show that the observed oligodactily in the Gremlin knock-outs can be explained with the reduced limb bud size. Given the importance of FGFs for limb patterning we extended the model to also include FGF signaling. Much as BMP, FGFs can diffuse rapidly 83 . BMP expression is induced by FGF and the parameter P B in (1) thus becomes a function of FGF, F, i.e.
where K BF and n and the Hill constant and Hill coefficient respectively. Signaling of BMP-bound receptors, C , R 2 B (4), has been found to both stimulate and inhibit FGF-dependent processes 53,84 . Such inconsistent observations may be the result of a concentration-dependency of signaling, and it has been suggested that high levels of BMP signaling interfere with FGFdependent processes while lower levels may positively regulate FGF expression in the AER 85 . Hence FGF activity is best described as , where K F1 = K F2 are the Hill constants for the activating and the inhibitory impacts of BMP signaling and n is the Hill coefficient. Here we note that FGF expression is restricted to the apical ectodermal ridge (AER), and p F is therefore non-zero only on the distal boundary of the domain (dashed part of Fig. 2C; simulated expression patterns are shown in Fig. S1, S3G-J). We then write for the FGF dynamics FGFs repress the expression of the BMP antagonist Gremlin, and enhance SHH expression in the ZPA, which in turn enhances Gremlin expression. However, as discussed above the Gremlin dynamics uncouple from the other equations, while the expression of SHH terminates before the expression pattern of Sox9 has assumed the characteristic digit pattern [43][44][45] and digit patterning can still be observed in limb buds when expression of Shh (and Gli3) is removed in knock-outs 10 . We thus neither include Gremlin nor SHH in the model.

Summary.
We propose that the observed patterning in BMP signaling is the result of regulatory interactions between BMP, B, its receptor R, and FGFs, F as summarized graphically in Fig. 2A.
The effective interactions in this model are graphically summarized in Fig. 2B. As final set of equations we have We note that FGFs are not part of the core patterning mechanism. Parameter values. The measurement of the parameter values (i.e. diffusion coefficients, production and degradation rates) in vivo is complicated and has only been carried out in few model systems, but not in the limb [86][87][88] LR The dimensionless model contains five parameters less and the patterning mechanism no longer depends on absolute diffusion and decay constants, but only on the relative diffusion coefficients and the relative decay rates. Similarly, the absolute protein concentrations do not matter but only the relative concentrations (as result from the relative expression and decay rates) relative to the Hill constant and the effective binding constant K C .
Domain and boundary conditions. The equations were solved on a domain whose shape was extracted from limb bud images at E12.5 as well as on a growing domain. The idealized domain is shown in Fig. 2C. We used zero-flux boundary conditions for B and R as the developing limb does not exchange with the surrounding except at the flank which is ignored here. FGF production in the AER was implemented as a flux boundary condition, i.e.
whereñ refers to the unit normal vector. We also considered alternative implementations with production in a thin layer (Fig.  S1B-E) or with Eq. 11 defined on the boundary (Fig. S1F). A thin layer with flux boundary conditions (Fig. S1C) and with Eq. 11 defined on the boundary (Fig. S1C) give similar patterns to those observed with flux boundary condition on the outer boundary (Fig.  S1A). A thin layer with either constant (Fig. S1E) or BMP-dependent FGF production (Fig. S1D) yielded overall the same pattern (same number of spots) but the shape of the spots was different.  Modelling domain growth. Proliferation is not uniform within the limb domain but the velocity field has not been published 14 . We therefore measured the proximal-distal axis of limb buds over developmental time and find that the growth rate is approximately linear (Fig. 3A) and we are therefore using a linear growth law with growth rate v g . We can then write for the proximal-distal length of the limb bud L t ð Þ~ L 0 z v g t where L 0 denotes the initial proximaldistal length of the limb bud, and t denotes time.
To conserve mass (rather than concentration) the reactiondiffusion equations for the soluble factors B and F, equations 9 and 11 must be expanded to include the advection and dilution terms 89 where u denotes the growth speed. For homogeneous growth at rate v g in proximal-distal direction y we then have u~_ To model local growth at the limb boundary (Fig. 4) the direction of growth was set to be normal to the boundary and proportional to the local FGF concentration, i.e.ũ~v g F mñ where vg is the growth speed andñ is the unit normal vector on the boundary. The exponent m was set to four, even though we note that similar results can be obtained with m 5 2.
Patterning on the limb domain. Since the model is an example of a Schnakenberg Turing-type model 90 we expected to find parameter ranges for which we would observe the emergence of patterns. To judge whether such mechanism could yield realistic patterns in a realistic time frame we solved the model on a domain that had been extracted from limb bud images taken at E11.5. Sox9 expression is induced in response to BMP signaling 91 , and we were therefore interested in the pattern of the BMP-receptor complex (C 5 R 2 B) which will mark the regions of active BMP signaling. This pattern is necessarily distinct from the pattern of BMP expression which remains uniform. In spite of uniform expression of BMP we observe the emergence of pattern in the BMP-receptor signaling domains that resemble the Sox9 pattern observed in experiments (Fig. 1J-L). Thus initially (at t 5 150) we observe strong BMP signaling (BMP-receptor binding) in the stalk (Fig. 1J) (marking the stylopod (s) in the limb bud images in Fig. 1B or the humerus (h) in Fig. 1F), two distinct spots in the transition zone from stalk to handplate (marking the zeugopod (z) in the limb bud images in Fig. 1B or radia (r) and ulna (u) in Fig. 1F) and a homogenous expression in the handplate (marking the autopod (a) in the limb bud images in Fig. 1B; digit 4 can be discerned already in Fig. 1F). Within Dt 5 200 the pattern in the handplate splits in the simulation (Fig. 1K) and distinct condensations can be discerned that mark digits 2-5 in the limb bud images (Fig. 1C,G). For the horseshoe pattern to split the BMP expression rate at the proximal end of the domain (rate r b2 ) must be 3-fold higher than the BMP expression rate in the ''hand plate'' (rate r b1 ). Without such a local increase in BMP expression the horseshoe pattern would not split. With time the digit condensations further elongate (Fig. 1D,E,H). In the simulations the digit condensations split while elongating (Fig. 1L). Elongation of patterns becomes more realistic if we solve the model on a uniformly growing domain (Fig. 3E) or if we allow the domain to expand in direction of the highest FGF concentration (Fig. 4). Such deforming growth towards higher FGF concentrations is also observed in cultured limb buds 38,92 . To achieve the characteristic digit shape it is important that the FGF concentrates distally from the BMP-receptor spots (Fig. 4A-C) as also observed in limb buds (Fig. 4D). This is achieved only if BMP-receptor signaling has predominantly a positive impact on FGF expression, i.e. if the BMP concentration is relatively low such that k 1 , R 2 B = k 2 . We therefore predict that at the time of digit formation BMP signaling enhances rather than represses FGF8 expression. We should stress that this is necessary only in order to obtain the correct positioning of the FGF8 expression domain. The BMP-receptor pattern described above can be obtained independently of whether BMP regulates FGF signaling positively or negatively.
In our simulations 5 spots emerge simultaneously in the handplate while in the developing limb digits appear in a sequence (4,2,3/5,1) with digit 1 appearing much later than the other digits. Also the appearance of digit 1 appears to be regulated in a different way from the other digits. Thus digit 1 is still formed in SHH knock-out mutants 10 while digit 1 is lost in FGF8 knock-out mutants 38 . The patterns depend on the specific choice of parameters (Table 1), and four instead of five spots could easily be obtained by altering any but the receptor expression rate (Fig. 5A). Alternatively the rate of FGF expression r F could be lowered 10-fold on the anterior side initially and then increase over time; such asymmetry in expression pattern is indeed observed in the limb bud. Thus while Fgf8 is expressed uniformly in the AER, expression of Fgf4 appears to be biased to the posterior site of the limb bud. FGF8 is expressed already early during limb bud development and remains expressed throughout the patterning process while Fgf4 expression is first detected around E10.5 and ceases around E12. As a result we expect enhanced FGF production in the posterior part of the limb bud at the onset of Sox9 expression. However, in this paper we will continue to analyse homogenous expression patterns.
Regulatory impact of the parameters. The patterns depend on the specific choice of parameters (Table 1) and since the parameters are difficult to determine accurately in experiments we wondered how sensitive our results would be to variations in parameter values. We find that all model parameters affect the pattern (Fig. 5A). Nonetheless the stable parameter ranges are sufficiently wide that molecular noise should not affect the patterning process. Parameters with particular wide stable ranges include d B , the rate of BMP degradation, which can be varied 6-fold without affecting the pattern, as well as the expression rates r b2 , r R and r F for BMP, receptor, and FGF and the threshold for BMP-dependent FGF expression, k 1 .
Since we know the time scale of the process we can relate the nondimensional kinetic rates to their dimensional counterparts. In our simulations splitting of the pattern in the handplate takes about Dt 5 100. The processes of one embryonic day (86400 seconds) therefore corresponds to Dt $ 200 such that t 5 1 corresponds to at least t 5 T 3 t 5 432 seconds in dimensional time. Accordingly the rate of FGF removal would correspond to d F 5 1/T , 2 3 10 23 s 21 in dimensional terms which is similar to the experimentally reported value of 10 23 s 21 87 . Since in our simulations any change in the rate of FGF removal can be counterbalanced by a change in FGF expression (Fig. 5B) a two-fold lower rate in FGF removal could also be obtained if we reduced the FGF expression rate by about 10-fold.
The size of the limb bud changes over time but is about 2 millimeters (mm) at E12.5, the final time of our simulation. Accordingly the length scale is l 5 2 mm/(H 0 1 R 1 1 R 2 )5 2 mm/19.35 5 0.1 mm and the FGF diffusion coefficient in the model corresponds to 34 mm 2 s 21 . We further require that BMP diffuses about 2.7-fold faster than FGF which brings the BMP diffusion constant close to the diffusion constant D , 100 mm 2 s 21 that is typically measured for soluble proteins 87 unless diffusion is impeded by adsorption to the surface 86 . We further require that the BMP receptor diffuses about 100-fold more slowly than the BMP protein. The typical range for diffusional constants of membrane proteins is several orders of magnitude lower than for proteins in solution (i.e. D 5 0.1 2 0.001 mm 2 s 21 76,77 ). A 100-fold reduced diffusion constant may still overestimate the receptor diffusion constant, in particular because the receptor would then diffuse a distance of 3.5 cell lengths (l~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2D R t 1=2 p^3 5 mm) within its receptor half-life time of t 1/2 5 ln(2)/d R 5 960 s even though receptors are mainly restricted to the surface of single cells. We therefore also tested alternative parameterizations in which the receptor diffusion constant is 400-fold lower than the BMP diffusion constant. Importantly we still obtain  analysis. Each parameter as perturbed independently from the reference value given in Table 1 as indicated on the vertical axis and the range was recorded for which pattern were qualitatively preserved on a static domain. The numbers above and below the bar indicate how many digits were either gained (1) or lost (2) as the range was exceeded. (B) The parameter space for which the pattern is preserved is larger (black line) if parameters are changed together as illustrated for the FGF production and degradation rates. The simulations were analyzed at t 5 750. Please note that the FGF degradation rate determines the time scale of the simulations (T 5 1/d F ). Compensation with the FGF production rate thus allow us to simulate the processes on different time scales. qualitative similar, yet smaller and sharper, patterns also with a much lower receptor diffusion constant as long as the FGF and BMP expression are adjusted accordingly. The pattern is also very sensitive to some of the protein expression rates. While there are no reported measurements of protein expression rates and concentrations for the limb bud that these constraints could be related to, there is a large number of mutant phenotypes that can be used for a qualitative comparison. These will be discussed in the following.
FGF. Many different FGFs are expressed during limb development. In the model we consider only those FGFs that are important during digit formation, i.e. FGF4 and FGF8. In mice that are lacking both Fgf4 and Fgf8 function in the forelimb AER limb bud mesenchyme fails to survive and the limb does not develop 93 . FGF4 is expressed only transiently in the limb bud (between E10 and E12 37,43 ) and removal of only FGF4 does not affect digit patterning 39,94 . The effect of Fgf8 removal depends on the developmental stage. Inactivating Fgf8 in early limb ectoderm (Msx2-cre;Fgf8 hindlimbs and in RAR-cre;Fgf8 forelimbs) causes a substantial reduction in limb-bud size, a delay in SHH expression, misregulation of Fgf4 expression, and loss of digit 1 (the thumb) 38 . In contrast, when functional Fgf8 is transiently expressed, as in Msx2-cre;Fgf8 forelimbs, digit 1 is present and digit 2 or 3 is lost 37 . If we remove 50% of Fgf expression in our model we observe less separation between the spots, similar to what would be expected in syndactyly (Fig. 6A). To lose a digit we need to reduce FGF expression to less than 10% of its normal levels. In this case we lose the middle digit as is the case in Msx2-cre;Fgf8 forelimbs without causing syndactyly. One explanation for the quantitative differences could be the dominance of FGF8 expression at later stages of limb bud development such that removal of Fgf8 corresponds to removing 90% of FGF activity.
All of the skeletal defects caused by inactivation of Fgf8 can be rescued when Fgf4 is expressed in place of Fgf8 95 . On the protein level FGF4 can thus functionally replace FGF8 in limb skeletal development. An increase in FGF signaling that occurs when the Fgf4 gain-of-function allele is activated in a wild-type limb bud causes formation of a supernumerary posterior digit (postaxial polydactyly), as well as cutaneous syndactyly between all the digits 95 . Increasing Fgf expression by 50% indeed results in supernumerary digits (Fig. 6B). If the additional expression is restricted to the posterior side (where Fgf4 is predominantly expressed) then the 50% increase results in the posterior placement of one supernumerary digit as well as a merging of spots on the anterior side which would correspond to syndactyly (Fig. 6C).
The model further predicts that the expression and distribution of FGFs remains homogenous initially, but subsequently forms distinct patches as indeed observed in experiments (Fig. S3G-J).
BMP. There are several BMPs expressed in the developing limb, including most importantly BMP2, BMP4, and BMP7 65 . Much as observed for BMP2 in experiments 79 we predict BMP to be first expressed rather homogenously at the boundary of the domain and to be subsequently expressed most strongly outside the digit domains. In our simulation a lower BMP production rate leads to digit loss (Fig. 6E) as indeed observed in the BMP2:BMP4 conditional knock-out 66 . A 30% increased production rate results in polydactyly in our simulation (Fig. 6F) and ectopically  Table 1) for t 5 750. Loss of digits at higher FGF production rates is due to presence of stripes; spots can be recovered if simulations are run longer. (E-H) The impact of changes in the BMP production rate. (E) Digits are lost when the BMP production rate is reduced to 50%. (F) Additional spots emerge as the BMP production rate is increased by 30%. (G) The pattern merge (polysyndactily) as the BMP production rate is increased by 50%. (H) The number of digits at different relative BMP production rates (1 corresponds to the reference rate in Table 1) for t 5 750. Note that at higher production rates pattern merge and digits are no longer observed. The constant but spatially modulated BMP production rates r B1 , r B2 (black and blue lines) and the FGF-dependent BMP production rate r B* (red line) have different effects on the patterning. The spatially modulated BMP r B2 in the stalk has little impact at larger values because most BMP then binds to receptors in the stalk. The spot(s) at the proximal end of the stalk and the two spots appearing at the transition from the stalk to the circular domain were preserved under all conditions and thus not counted.
www.nature.com/scientificreports expressed BMP2 indeed induces duplication of digit 2 and bifurcation of digit 3 84 . Moreover, conditional inactivation of the gene for BMP4 results in polydactyly 67 and other experiments reveal that even though conditional inactivation of a gene for BMP4 downregulates BMP4 expression, the overall amount of BMP signaling increases 68 . As the BMP expression rate is further increased in our simulations ($50% increase) the patterns merge to give rise to polysyndactily (Fig. 6G). Such strong overexpression has not yet been achieved in mutants so that we cannot check this prediction.
BMP receptor. Conditional mutants of BMP receptor type IA lack digits. Since the receptor is an integral part of the proposed Turing mechanism we obtain loss of digits also in the model. Hypomorphic mutants of the BMP receptor have not been described; the model would predict loss of digits (Fig. S3). Mild overexpression of the receptor should lead to a gain of digits while stronger overexpression will lead to a loss of patterning (Fig. S3).
We conclude that the expression rates of BMP and BMP receptor need to be tightly controlled to permit patterning. While FGF is not an essential part of the patterning module FGF enhances BMP production and can thus affect the number of digits.
The impact of domain size and shape. Turing patterns depend on the size of the computational domain and more modes (patterns) emerge as the size of the domain increases 90 . Several mouse mutants have been reported with different shapes and sizes of their limb domain, and larger domains typically correlate with increased digit numbers. Thus GLI3 mutants are polydactylous and the limb domain seems to be widened to 130% while Gremlin mutants are oligodactylous and the limb domain is shrink to about 60% of its original size 10,73 . Previous theoretical models showed that such a correlation can in principle be explained with Turing models 22 , but the domain shapes were not based on real limb buds and no quantitative comparison of domain size and digit number was previously carried out. When we shrink the computational domain uniformly to 60% of its original size we observe a decrease to two BMP-receptor (BR 2 ) patches (Fig. 7A,B), much as is observed in experiments 73 . When we enlarge the computational domain uniformly to 130% of its original size the digit number more than doubles (Fig. 7A). However, in the GLI3 2/2 mutant the expansion is observed only along the anterior-posterior axis 10 . If we increase only the anterior-posterior axis of the limb bud to 130% we observe four additional digits which emerge mainly in the middle of the domain (Fig. 7C,D). A closer analysis of GLI3 2/2 mutants reveals that limb buds are expanded mainly to the anterior side and that extra digits are concentrated in the anterior domain, while the posterior side is about normal. When we analysed a domain that was expanded by 160% to the anterior and constant on the posterior side we observed three additional spots which now emerge in the anterior side (Fig. 7E,F). This corresponds well to the 8 digits observed in the GLI32/2 null mouse 10 .
The impact of growth. The limb bud grows while patterns emerge. Given the importance of domain size and shape we wondered how growth would affect the patterning process. We measured images of limb buds between E9 and E13 and noticed an almost linear growth in proximal-distal direction (Fig. 3A). Limb buds increase in length by about 5-fold within the first four days of development. When we solved the model on a linearly, uniformly growing domain we obtain similar patterns as before (Fig. 3C). However, the number of digits strongly depends on the growth speed (Fig. 3B). For low growth rates we can adjust the number of digits with the help of the other parameters. For increasing growth speeds the speed greatly affects the patterning process and we both lose and gain digits (Fig. 3B). As the growth speed further increases patterns merge and eventually all digits are lost (Fig. 3B). From Figure 3A we can estimate the growth speed as about 7.7 nanometers per second. In dimensionless terms (T 5 432 sec and l 5 0.1 mm) this corresponds to a growth speed of about v g 5 0.03 which is too large to permit patterning (Fig. 3B). This may be the reason why there is a the complex pre-patterning process that we do not consider in this model and which starts already around E9.5 1,2 . The complex regulatory interactions prior to the emergence of digits adjust the expression domains of BMP and FGFs in the limb bud. Such embedding of a Turing mechanism into a sophisticated regulatory network may be important to permit the patterning process to proceed on a rapidly growing domain. Further experiments combined with more detailed modelling also of the other parts of the regulatory network present in the developing limb will, however, be required to either confirm a receptor-ligand based Turing mechanism (embedded in a wider network) or to define an alternative mechanism for digit patterning on the rapidly growing limb bud.

Discussion
Limb development has been studied for decades and much is known about the molecular networks that regulate limb development 1,2 . www.nature.com/scientificreports Based on the available information BMP signaling appears to be the key signaling pathway that regulates the expression of Sox9, the first marker of digit condensations. When modelling BMP signaling we noticed that the regulatory interactions with its receptor constitute a Turing mechanisms. Pattern thus easily emerge. The emerging pattern sensitively depended on the shape and size of the domain, and only when we solved the set of equations on a domain that had been extracted from limb bud images we could reproduce patterns of BMP-receptor activity that resembled the experimentally observed patterns for Sox9. An important test for the suitability of a mathematical model is its consistency of model predictions with a wide range of independent experimental observations. Limb bud development has been studied intensively and a large body of experimental results exists to test the model with. If parameters are chosen appropriately also the phenotype of available mutants can be reproduced. The mechanism is robust to small changes in the parameters as may arise from molecular noise.
Interestingly, while we achieve good correspondence to available experimental data on the static domain, formation of a Turing pattern fails on a domain that grows as fast as measured in the embryo. We propose that the early pre-patterning that involves also SHH/ Ptch/Gli signaling as well as the many other regulatory interactions that we ignore in this simple model are in place to support the patterning mechanism on a rapidly growing domain. We are currently developing a simulation for these processes to better understand the impact of these earlier regulatory processes.
Many models have previously been proposed to explain the emergence of digits, and Turing models in particular have been found previously to be sufficiently flexible to recapitulate the various patterning aspects of digit emergence 30,31,34 . So far, however, it had been difficult to link the components of a potential Turing mechanism to the molecular constituents. We propose that the Turing pattern in the limb arises from a ligand-receptor interaction. For this to work the ligand needs to be a dimer as is the case for BMPs 78 . Moreover, the receptor needs to diffuse, yet at a much lower speed than the ligands as is indeed the case. We notice that the diffusion constant that we are using in this simulation is on the high end of what is plausible; lower receptor diffusion constants can be accommodated but the pattern is then sharper and smaller. Again it is possible and plausible that the additional regulatory interactions in the limb bud help to expand the pattern.
We have previously proposed a mechanism for branch point and branch mode selection in the developing lung 96 . Interestingly, also here receptor-ligand interactions gave rise to a Turing pattern that could recapitulate the different observed branching patterns in the developing lung for physiological parameters and a similar mechanism can be identified also for the developing kidney (Menshykau et al, submitted). While all three systems are based on different signaling proteins (i.e. BMPs in the limb, SHH in the lung, and GDNF in the ureter) they may all exploit the same mechanism to pattern an embryonic field in that these proteins are all multimers that interact with their receptors in a way that induces receptor expression. Interestingly, while BMP can affect its own expression either in a positive or negative feedback via FGFs 52,53 , SHH in the lung clearly engages in a negative feedback via FGF10, while GDNF in the ureter engages in a positive feedback via Wnt signaling 97 . We can obtain the same pattern with a positive and a negative feedback in the limb, lung, and ureter. However, we notice that the possible parameter space is wider in case of a positive feedback. It will be interesting to model further embryonic self-organizing systems to see whether a receptor-ligand interaction may be a general paradigm to enable Turing pattern to emerge in developmental systems.
Further advancements in our understanding of limb development will require the development of three-dimensional models and the inclusion of more signaling factors. The parameterization and validation of such models will require new experimental data that also reveal the three dimensional dynamics. Such information can now be acquired with the help of optical projection tomographs 98 . Further advances in experimental techniques can thus be expected to provide exciting new insights into the regulatory processes of digit formation during limb development.

Methods
Numerical solution of PDEs. The PDEs were solved with finite element methods as implemented in COMSOL Multiphysics 4.1 and 4.2. COMSOL Multiphysics is a wellestablished software package and several studies confirm that COMSOL provides accurate solutions to reaction-diffusion equations both on constant 99 and growing two-dimensional domains [100][101][102] . Both mesh and the time step were refined until further refinement no longer resulted in noticeable improvements as judged by eye. The simulations were optimised for computational efficiency as described in 103 .