Mathematical modeling of chondrogenic pattern formation during limb development: Recent advances in continuous models Journal Pre-proof Mathematical modeling of chondrogenic pattern formation during limb development: Recent advances in continuous models Paramita Chatterjee, Tilmann Glimm, Bogdan Kaźmierczak PII: DOI: Reference: S0025-5564(20)30014-6 https://doi.org/10.1016/j.mbs.2020.108319 MBS 108319 To appear in: Mathematical Biosciences Received date: Revised date: Accepted date: 8 July 2019 17 January 2020 17 January 2020 Please cite this article as: Paramita Chatterjee, Tilmann Glimm, Bogdan Kaźmierczak, Mathematical modeling of chondrogenic pattern formation during limb development: Recent advances in continuous models, Mathematical Biosciences (2020), doi: https://doi.org/10.1016/j.mbs.2020.108319 This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Inc. ฀฀฀฀฀฀฀฀฀ Highlights • We review mathematical and conceptual models of vertebrate chondrogenic pattern formation from the last 50 years. • The last several years have seen the publication of new models based on substantially improved knowledge of the relevant gene regulatory networks, in particular Raspopovic et al.s BSW model (2014), Badugu et al.s model (2012) and Glimm et al.s galectin model (2014) • We include a list of gene products that are relevant to mathematical models of chondrogenic pattern formation in the limb. 1 ฀฀฀฀฀฀฀฀฀ Mathematical modeling of chondrogenic pattern formation during limb development: Recent advances in continuous models Paramita Chatterjee, IPPT PAN, Poland Tilmann Glimm, Western Washington University, USA Bogdan Kaźmierczak, IPPT PAN, Poland Abstract The phenomenon of chondrogenic pattern formation in the vertebrate limb is one of the best studied examples of organogenesis. Many different models, mathematical as well as conceptual, have been proposed for it in the last fifty years or so. In this review, we give a brief overview of the fundamental biological background, then describe in detail several models which aim to describe qualitatively and quantitatively the corresponding biological phenomena. We concentrate on several new models that have been proposed in recent years, taking into account recent experimental progress. The major mathematical tools in these approaches are ordinary and partial differential equations. Moreover, we discuss models with non-local flux terms used to account for cell-cell adhesion forces and a structured population model with diffusion. We also include a detailed list of gene products and potential morphogens which have been identified to play a role in the process of limb formation and its growth. 1 Introduction Organogenesis is one of the most intriguing phenomena in biology. The question ’How does an initially homogeneous and indistinguishable set of cells give rise to subgroups of differentiated cells, tissues and whole organs?’ is an extremely challenging and complicated modeling problem. An important example of organogenesis is vertebrate limb development. In this paper, we give a short overview of the biological background and then discuss mathematical modeling approaches. These models are concerned with the study of limb bud outgrowth and shaping as well as on skeletal pattern formation. While these two processes are dynamically interconnected [111], they are conceptually distinct. Correspondingly, there are two classes of mathematical models of limb development: Those concerned with modeling the growth of the limb buds and those that model the process of pattern formation of pre-skeletal cartilage within the developing buds. We concentrate on continuous models formulated as systems of partial differential equations representing various chemical concentrations as well as cell densities. In the last six years in particular, new classes of reactiondiffusion models [149, 132, 53, 6, 186] have been proposed which are based on new experimental insights into the molecular basis of chondrogenesis and incorporate much more detailed interactions from gene regulatory networks than previous models. While self contained, our survey can also be seen as an update of previous surveys of mathematical models in chondrogenesis (see [197, 196, 52, 189]). The review by Zhang et al., [189] which appeared in 2013, is particularly relevant and covered the history and current state of the art at the time for a mathematical audience. Some of the models that date before 2013 surveyed here have already been covered there. We include them here for their continued relevance and their importance for subsequent models. We also present a list of gene products which have been found to be relevant for limb chondrogenesis, including many which may play important roles as morphogens; see Table 2. The paper is organized as follows: In section 2, we describe in broad strokes the biological background of limb growth and skeletal pattern formation, as well as give a short overview of the complicated history of 2 ฀฀฀฀฀฀฀฀฀ influential conceptual models that have guided thinking and experiments. Section 3 highlights the two main approaches to explain patterning in the limb: reaction-diffusion mechanisms and positional information models. Recent experimental evidence highlights that these two approaches may actually act in tandem, with a reactiondiffusion mechanism laying out periodic patterns that are sculpted through positional information imparted via morphogen gradients as a secondary mechanism. In section 3, the core content of this paper, we indicate several broad approaches modelers have taken and review in depth several mathematical models with special emphasis on a crop of new reaction-diffusion type models. An overview of the models discussed in depth is given in Table 1. 2 Biology of limb bud growth and basic concepts of skeletal pattern formation The process of limb growth is very similar in all tetrapods, but most of the experimental work has been done on chicken and mice embryos. For chicken, the complete developmental process from a fertilize egg to hatching takes about three weeks. Limb buds begin to emerge from the embryonic body at the end of third day and elongate rapidly. On the fourth day the humerus begins to appear in the form of a chondrogenic condensation, that is, a tight aggregate of precartilage cells. The complete pattern of the limb skeleton is laid out in cartilage elements by the seventh day. 2.1 Limb bud outgrowth and shaping The vertebrate limb is an outgrowth from the embryonic body wall. As the outgrowth proceeds, a morphologically distinct ectodermal thickening, consisting of a partially stratified epithelium, forms at the distal tip, known as the apical ectodermal ridge (AER; see Figure 1). The AER is the source of several fibroblast growth factors (FGFs) [189], diffusible morphogens that form a proximal-distal gradient and are instrumental in maintaining growth of the limb bud. (See Table 2.) The AER forms in three steps: induction of precursor cells, migration of precursor cells and compaction of the ridge. Figure 1: Schematic illustration of a limb bud. Positions of the AER (Apical Ectodermal Ridge) and the ZPA (Zone of Polarising Activity) are shown along the limb axes in the developing limb. The image also shows a region which is identified as the ”Progress Zone” in the (now disputed) Progress Zone model. (Modified from [194], page no. 214 ) In the embryo, during the early stages of limb formation, the limb bud is evolving along three axes: the axis from shoulder to hand is known as proximal-distal or PD axis, from thumb to little fingers is known as anterior-posterior or AP axis and the axis from back of hand to palm is called as dorsal-ventral or DV axis as shown in Figure 2. 3 ฀฀฀฀฀฀฀฀฀ Initially, the limb bud is composed of a set of mesodermal cells covered by a relatively thin layer of ectoderm. The limb bud arising from the chick body wall is flat and almost elliptical in cross section with its major and minor axes parallel to the anterior-posterior and dorsal-ventral axes (see Figure 2) respectively [33]. Wing buds start to grow at Hamburger-Hamilton stage 16 1 (HH 16), leg buds at HH17. During limb outgrowth, the proximal-distal length expands quickly [33] and the posterior part of the limb grows faster than the anterior one [70]. Figure 2: Outgrowth and shaping of vertebrate limb. Development of a chicken limb in three axes: Proximal-Distal (humerus in wing, femur in leg), AnteriorPosterior (radius and ulna in wing, tibia and fibula in leg) and Dorsal-Ventral (digits in wing and leg). (Modified from [189]) The AER and underlying mesoderm play an important role in the outgrowth of the limb. If the AER is removed during the growth of the limb, the underlying mesoderm stops dividing and consequently limb outgrowth stops [51, 138, 143]. Also it is observed in tissue cultures that the mesoderm is induced to proliferate if the AER is combined with limb mesoderm [51]. The first to carry out experiments with AER removal was John W. Saunders [142], building on previous ideas and methods by Lillie [80], Peebles [128] and others. In his paper [142], Saunders described the removal of the AER from the tip of the wing bud resulting in truncation of wing depending on the stage at which the AER was removed. He concluded from these results that the AER plays an essential role in growth of a wing and ’the orderly formation of the wing parts’ [169]. The main idea of this experiment was revisited by Rowe and Fallon [137] (for wings and legs). Removal of the AER in the HH18 or early HH19 stage of development results in the absence of some of the digits of the limb. Moreover, the larger the piece cut off from the anterior part of the AER (between the somite levels 15 and 20), the more digits are absent in the final limb. While AER removal leads to immediate cessation of growth, grafting an AER from before HamiltonHamburger stage 29 onto a younger limb results in proper limb development [138]. Likewise, taking a young AER and grafting it onto an older mesoderm does not lead to any limb-elements duplications. This suggests 1 In developmental biology, the Hamburger–Hamilton stages (HH) describe 46 chronological stages in chick development, starting from laying of the fertilized egg and ending with the hatching of the chick. These stages are described, e.g. in [60]. 4 ฀฀฀฀฀฀฀฀฀ that information of the proper sequence of limb development is programmed in the limb’s mesoderm, not in the AER. For many years, the effects of limb morphogens (in particular FGFs) were regarded as mitogenic, i.e., that they promote (isotropic) cell division at the distal tip to drive limb bud outgrowth. Recent works have put the importance of this for limb growth and shaping into question and suggest instead that growth and shape of the limb bud hinges on a chemotactic migratory response of the cells to FGF gradients [79] towards the AER [189, 16, 183]. Cell orientation depends on Wnt signaling, whereas its velocity depends on FGF signaling [55]. Some models (like the model proposed in [18]) suggest that the flow of the limb mesenchyme is mechanically guided by the dorsal and ventral ectoderm, while some indicate that the dorsal ectoderm, excepting the underlying base membrane, is not necessary for normal limb shaping [91, 189]. Another important question concerning limb development is how the identity of limbs (hindlimb and forelimb) is determined. Forelimbs and hindlimbs emerge from the body wall of the embryo, initially composed of undifferentiated mesenchymal cells covered by a layer of ectoderm. The determination of limb identity depends on two transcription factors: Tbx5 and Tbx4 (T-box transcription factor) [50]. Both of these factors accelerate Fgf10 expression [1, 102, 115, 116] during limb bud initiation. The absence of Tbx5 in mice prohibits Fgf10 expression and forelimb skeletal formation [1], while Tbx4 knockout in mouse prevents hindlimb development [109]. 2.2 2.2.1 Skeletal pattern formation in the growing limb Biological background The skeleton is a complex organ which protects a number of internal organs and stores ions, particularly calcium [65]. The skeleton is formed by two different tissues, cartilage and bone, which support each other. Low metabolic rate, avascularity, capability for continued growth and high tensile strength coupled with resilience and elasticity are the distinctive properties of cartilage, whereas constant renovation to meet both mechanical and metabolic demands, vascularity, rigidity and hardness are the properties of bone. Bones differ individually by size and shape, forming discrete arrays of elements. During embryonic development, the limb skeleton develops from the paddle shaped limb bud mesoblast which is covered by epithelium, called ectoderm [130]. In most vertebrates, skeletons progress as a series of primordial cartilage in proximal-distal form [113]. In the case of the chicken, the humerus of the upper wing is established first, then radius-ulna in mid-wing, followed by the wrist bones and the digits lastly. Bones replace cartilage in most vertebrate. Before cartilage forms, mesenchymal cells of the mesoblast form tight aggregates called precartilage mesenchymal condensations. These are the morphological basis of the skeleton [57, 58, 112]. This process is accompanied by the production and secretion of ECM glycoproteins, like fibronectin, which trap the cells at specific positions [45, 46]. The aggregation of cells becomes firm due to cell-cell adhesion which depends on cellsurface molecules such as NCAM [178](neural cell adhesion molecule), N-cadherin [119] or cadherin-11 [74, 85]. In the final step of chondrogenesis, precartilage cells differentiate into cartilage, a process associated with various gene products [120] and cartilage specific ECM [133]. The formation of the limb musculature (a process called myogenesis) happens independently of chondrogenesis [208]. In fact, the precursors of muscles and cartilage derive from two different cell populations: Muscle precursor cells originate in the somites and migrate to the limb buds early in development, whereas precartilage cells originate in the body wall [209, 210]. 2.2.2 Mechanisms of skeletal pattern formation The question of the precise mechanism of pattern formation has a long and complicated history. Several models have been proposed over the years which often varied very much in the level of mechanistic detail they incorporated, ranging from purely conceptual models to very detailed mathematical models. However, two 5 ฀฀฀฀฀฀฀฀฀ types of models have clearly been the most important: Positional information models and models based on the Turing mechanism in reaction-diffusion systems. One positional information model, Wolpert’s Progress Zone model [180], was particularly influential, but starting in the mid-2000s, new experimental results cast serious doubts about its validity [167]. (See discussion below.) Today, the Progress Zone model is widely seen as inadequate to explain all experimental results [163, 160]. Models that are based on self-organizational capacity of the mesoderm have gained in acceptance, in particular reaction-diffusion type models [53, 132, 149, 6]. In section 3, we will describe the mathematical background of some of these new models. 2.2.3 Progress Zone model, Early Specification model, other positional information models Historically, the most influential model has been the Progress Zone model proposed by Lewis Wolpert in [180] (see also [155]), although as mentioned above, it is now regarded as outdated. The Progress Zone (PZ) refers to a region of undifferentiated mesodermal cells approximately 300µm beneath the AER in the developing limb bud (see Figure 1). The PZ model postulates that cells obtain positional information via the amount of time spent within this region. Specifically, the cells ‘measure’ the time they spend in the Progress Zone. When the cells leave this region, their clock stops. Thus the humerus is formed by the cells which are ‘first’ to leave the zone, whereas the digits are formed by the cells that are ‘last’ to leave this zone (see Fig.1 in [167]). This model was motivated by the AER removal experiments by Saunders et. al. [138], and Summerbell et. al. [156, 157] (see Section 2.1). Indeed, in the framework of the model, the removal of the AER stops the change of positional information value in the apical mesenchyme and progress halts in the zone at the tip. Thus cells’ proliferation rate is reduced and they start to differentiate according to the positional value they had just before the AER extirpation. As a result, the outcome limb is shorter and lacks distal elements [155]. The validity of the Progress Zone model was cast in serious doubt by the work of Dudley et. al. in [36], where the effects of extirpation of the AER on the formation of limb elements was revisited. By tracking labeled limb bud cells, it was observed that the removal of the AER resulted in the massive death of underlying distal mesenchymal cells along with the truncation of distal elements. This suggests that distal cells do not take part in the chondrogenetic process after AER removal, and thus that the pattern of truncated elements is produced without reference to events in the most distal zone of the limb bud, i.e. the presumed progress zone [36, 160]. Based on the findings, Dudley et. al. in [36] formulated an alternative ‘Early Specification’ model. This model assumes that the limb pattern is specified at a very early stage of development and the limb segments already have distinct molecular features before the limb bud grows out [36, 167]. Critics of the the Early Specification model contend that it would imply that each of the eight cartilaginous limb bud elements (humerus, radius and ulna, two carpal elements that are initially same size as the radius and ulna but fail to grow and a maximum of three digital elements) would correspond to approximately only four layers of cells [167]. Moreover, recent experiments concerning FGFs (mainly on FGF4 and FGF8) and their effects on pattern formation of cartilage patterns in the limb (see [78, 160]), can not be fully explained by either the Progress Zone or the Early Specification models. A model that incorporates aspects of both these models was recently made by Saiz-Lopez [139], but in essence, it is now widely acknowledged that neither of these two models are fully satisfactory, as acknowledged in a joint paper by Tabin and Wolpert, leading proponents of the two models [163]. It is a common approach to conceptually separate patterning along the three major dimensions of the limb bud, the proximal-distal (PD), the anterior-posterior (AP) and the dorsal-ventral axes (DV). This is likely influenced by the historic popularity of the framework of positional information models, where these three axes could in principle correspond to three different spatial morphogen gradients. Seen from this point of view, the now-superseded Progress Model was proposed to account for proximal-distal patterning. This led researchers to investigate possible independent mechanisms of anterior-posterior patterning. In fact, now-classical experiments [158, 159, 166] show that a region at the posterior end of the growing limb bud, has an important role in organizing AP patterning across the distal limb, known as the zone of polarizing activity (ZPA). The protein Sonic Hedgehog (Shh), produced from ZPA, influences the distinct fates of the limb cells along the anteriorposterior axis. Therefore digit 1 (thumb or big toe) consists of cells with lowest concentrate of Shh in the distal 6 ฀฀฀฀฀฀฀฀฀ limb, whereas the most posterior digits (little finger or small toe) arises from the region close to the ZPA. There is some good evidence that dosage of Shh and duration of exposure to Shh influence digit numbers in the limb [61, 146, 170, 52, 171, 172], although in the absence of Shh expression, digit patterning can still be observed [81]. Thus Shh is dispensable for digit formation. The issue of the effect of Shh on digit identity is arguably even more complex. Indeed, it is still a controversial question whether the different digit morphologies are due to distinct developmental programs among the digits (developmental identities), or are just caused by the spatial positions along the limbs anterior-posterior axis (positional identities); see [177, 154]. In comparison to proximal-distal and anterior-posterior patterning, dorsal-ventral patterning has received much less attention. Experiments done by Chen et. al. in [21] suggest that dorsal-ventral axis is formed by cells derived from both the mesoderm and ectoderm at different stages of development. When parts of limb bud mesoderm were detached and centrifugally compacted limb bud cells were reattached into the ectodermal hull of a three or four day chick wing bud, and then grafted to the flank of the host embryo, it was found that the skeleton and musculature of the distal elements have a dorsal-ventral axis which conforms to that of the ectoderm [65]. Similarly, if intact mesodermal cores were re-merged with rotated ectodermal hulls, the dorsal-ventral axis of the ectoderm was reversed and the skeleton with musculature was also reversed along the dorsal-ventral axis [65, 21]. This indicates that before the appearance of the AER, the ectoderm can specify the dorsal-ventral axis. Other mesoderm and ectoderm recombination experiments revealed that, most probably, the ectoderm acquires dorsal-ventral polarity from the underlying mesoderm prior to limb bud outgrowth at approximately Hamburger-Hamilton stage no 15 [21, 48]. 2.2.3.1 Mathematical Models, in particular reaction-diffusion models The above mentioned models of the formation of the skeletal pattern are conceptual models that synthesize biological, biochemical and physical ideas, but do not use any mathematical specifications of their concepts. In contrast, several mathematical models have been proposed. These models investigate how the interaction of gene expression, cell proliferation, cell movement and adhesion and differentiation can lead to the spontaneous emergence of chondrogenic patterns. The separation of patterning mechanisms by axes (in particular proximaldistal and anterior-posterior) typical of the positional information approach does usually not apply to these models; the same process is responsible for the shape of condensations in all dimensions, although secondary mechanisms may ’refine’ the patterns once they are laid out [114]. Central to most of these models is the Turing mechanism for pattern formation in systems of reacting and diffusing chemicals (see [175]). Among these is Newman and Frisch’s model [110], a reaction-diffusion model by Hentschel et al. [64], the mechanochemical models of Murray and Oster [174, 173] or the boundary model by Meinhardt [94, 96]. In the last six years in particular, new classes of reaction-diffusion models [149, 132, 53, 6, 186] have been proposed which are based on new experimental insights into the molecular basis of chondrogenesis and incorporate much more detailed interactions from gene regulatory networks than previous models. In particular experimental work by Sheth et al. [149] showed that the number of digits is under fine control of the Hox dosage, thus showing that one of the key properties of reaction-diffusion mechanism, a pattern wavelength that is under fine control of the kinetic parameters, was actually realized in the limb. (The findings are summarized in Figure 8.) These findings in particular have greatly increased the acceptance of reaction-diffusion models in the community of developmental biologists [40, 26]. That said, given the very different natures of the conceptual biological models and the more concrete mathematical models, the two are often somewhat hard to compare and contrast. For instance, in the widest sense, Wolpert’s positional information model postulates that a cell’s spatial position determines its fate, but does a priori not stipulate the exact mechanism of how the cell senses its position, or how this information is translated into its behavior. Temporal or spatial gradients of morphogens are central to this mechanism, but there are many different ways in which they may be set up and maintained. (See e.g. the review of Rogers and Schier [136].) In a sense, the Turing mechanism may be regarded as a possible mechanism for the establishment of morphogen gradients, although in a more narrow sense, the positional information concept is that the spatial 7 ฀฀฀฀฀฀฀฀฀ morphogen gradients are set up by regions of specialized cells such as the Zone of Polarizing Activity (ZPA) in the limb; this sense is not compatible with the Turing mechanism, where patterns are set up in an autonomous, self-organizing way that does not require a group of specialized cells [136]. In section 3, we further discuss the relationship between Wolpert’s positional information model and the Turing mechanism, then survey recent mathematical models in the ensuing section. 3 Mathematical ideas behind models of chondrogenic pattern formation Chondrogenesis is one of a plethora of examples of pattern formation in embryogenesis, giving rise to a fundamental questions in developmental biology: What are the fundamental mechanisms by which biological patterns (structures and shapes) form [54]? In the following two subsections, we concentrate on two well known ideas, the reaction-diffusion (RD) mechanism and the Positional Information (PI) approach proposed by Lewis Wolpert. 3.1 Mechanism based on diffusion and interaction between groups of molecules In description of embryological development, reaction-diffusion models are widely used to explain self-regulated pattern formation [76]. In 1952, in the celebrated paper The Chemical Basis of Morphogenesis [175] Alan Turing proposed reaction-diffusion (RD) model addressing the problem of biological patterns formation. Turing showed that a simple system of two equations for interacting morphogens can describe six types of spatial patterns, including traveling waves and oscillations as well as stable periodic patterns, such as stripes or spots, arising from a uniform field of cells [54]. The most important of these findings is the discovery of the mechanism behind the spontaneous formation of spatial patterns in the concentrations of the morphogens. This type of mechanism is known as diffusion driven instability, or Turing instability. Mathematically, the key observation behind the Turing mechanism is that a spatially homogeneous steady state of a dynamical system which is stable to small perturbations in the absence of diffusion, may become unstable to small spatially non-homogeneous perturbations if diffusion is added to the system. This is surprising and unexpected as diffusion usually degrades spatial patterns and leads to uniformity in the long run (see, e.g. [54, 108]). These ideas about pattern formation continue to enjoy popularity among mathematical modelers and developmental biologists. This interest is certainly at least partly due to the fact that reaction and diffusion are two fundamental, ubiquitous processes in development. The first suggestion that the Turing mechanism was responsible for generating the characteristic patterns of the limb skeleton was made by Newman and Frisch in 1979 [110]. (Turing himself never worked on limb development.) Many other different types of patterns in biology, such as animal coat patterns of zebras, leopards, or mollusk shell pigmentation patterns, have been suggested to be RD patterns, arising via various reaction-diffusion systems (see, for example, [95, 108]). However, it is usually difficult to conclusively identify an RD mechanism as the definite patterning mechanism and rule out other possible mechanisms [54, 124]. Indeed, it is often difficult to assign concrete morphogens playing the role of activators or inhibitors in the abstract mathematical model describing the process of pattern formation, although, in several cases, morphogens corresponding to the proposed reaction-diffusion systems have been identified with varying degrees of certainty (see e.g. [4, 40, 89, 100, 105, 132, 151]). In the most popular version of a two-morphogen system, RD mechanisms also requires diffusion coefficients of activating and inhibiting morphogens to be of different sizes, often one or more orders of magnitude. The notion of Turing instability is mainly used to explain the appearance of periodic biological patterns. However, reaction-diffusion mechanisms may be also applied to understand mechanisms of emergence of nonperiodic patterns (see, e.g. [106, 147, 150]). 8 ฀฀฀฀฀฀฀฀฀ 3.2 Positional Information mechanism Wolpert proposed the Positional Information (PI) model in the late 1960s. He was interested in understanding how pattern formation could be directed by existent heterogeneities across the tissue. This is a contrast to pattern formation in reaction-diffusion systems investigated by Turing, where the heterogeneities are not initially given, but rather they arise via a symmetry-breaking process that starts from (in principle arbitrarily small) perturbations within a homogeneous medium. The key idea is that spatial morphogen gradients, i.e. changes in morphogens’ concentrations over space, may result in different cellular behavior, which in turn may lead to the formation of spatial patterns [180, 181, 54]. Compared to RD models, which suggest, for example, that stripes or spots of morphogens directly produce stripes or spots of cell types in the resulting tissue, Wolpert introduced an ‘interpretation’ step, according to which cells can interpret the local concentration of ‘positional molecules’ which in turn determines the fate appropriate for that position [54]. This interpretation step implies an additional freedom that not only allows a smooth, monotonic molecular concentration gradient to give rise to any arbitrary pattern, periodic (like stripes or spots) and non-periodic, but also allows different cell types to interpret the same signal in different ways. Thus the same morphogen gradients may give rise to different patterns depending on the cell types [54]. In the PI framework, the morphogen concentrations effectively work as positional coordinates. In other words, the spatial distribution of the PI molecules may not be isomorphic to the developing or final limb skeletal patterns [189]. A key difference between Turing’s and Wolpert’s ideas, which has certainly had an important impact on their reception, is their intuitive appeals. Turing’s reaction-diffusion mechanism is not intuitive. Self-organizing, morphological patterns arising from ‘nothing’ (to be more precise from arbitrarily small spatial perturbations) are difficult to understand as diffusion usually induces stability. In contrast, Wolpert’s positional information approach is more intuitive and easy to grasp. In Wolpert’s theory, one can regard morphological patterns as a result of an interplay between the spatial concentrations of different morphogens secreted from the corresponding sources, which are usually separated and localized in different positions of the evolving organ. In the case of limb bud growth these sources can be e.g., the AER secreting FGFs, and the ZPA secreting SHH (see Figure1). This does not mean that PI is the ’simpler’ mechanism. In fact, in some sense, the contrary is true. RD patterns may arise spontaneously from the collective behavior of individual agents that obey very simple rules. In contrast, the rules obeyed by the individual agents in a PI mechanism must be by necessity complex, since they incorporate different behaviors according different morphogen concentrations. 3.3 Reaction-diffusion and positional information mechanisms as congruous morphogenetic processes Reaction-diffusion and positional information mechanisms are not mutually exclusive. In fact, gradients of morphogens may locally modulate the parameters of an RD system. Thus a Turing-type bifurcation can explain the emergence of patterns which are fine-tuned through morphogen gradients. These ideas have been around for a long time, e.g. in the pioneering work of Newman and Frisch on reaction-diffusion patterns in the limb [110], but have only relatively recently gained more attention [17, 198, 199, 52, 54, 200]. In [198] and [199] Glimm, Zhang and co-workers investigated a system of activator-inhibitor reaction-diffusion equations where the production of the morphogens was modulated by an external linear gradient, leading to reaction-diffusion equations with explicitly spatially dependent reaction kinetics of the form       ∂ a1 d ∇2 u u + ε · x = F (u, v) + , (1) c a2 ∇2 v v ∂t where x denotes the first component of spatial coordinates (x, y). For the case of a shallow gradient (small εc ) close to the Turing bifurcation, the bifurcating steady state solutions can be investigated via a weakly nonlinear bifurcation analysis in two parameters. These are the parameter εc and a second parameter εd defined via d = d0 − d2 ε2d in the case of a supercritical Turing bifurcation [200, 204], where d0 is the critical 9 ฀฀฀฀฀฀฀฀฀ diffusion coefficient for the appearance of Turing patterns. The resulting patterns are thus perturbations of Turing patterns, influenced by the external gradient. A shallow gradient does not change the wavelength of the patterns, but the amplitude is indeed affected. In fact, as is shown in [199], the amplitude may have a minimum in the interior of the domain. The authors point to an experimental result in the mouse, reported in [205], where a manipulation of the SHH gradient led to the loss of one of the middle digits, digit 2, while all other digits were retained. While this result is hard to explain in a positional information model, it has a simple explanation in the framework of the above toy model: If the digital structures correspond to regions where the activator concentration is above a critical value, then a shallower external gradient can mean that one of the structures in the interior is ”lost” due to the amplitude of the corresponding wave peak falling below the critical value. In a follow-up paper [200], the authors derived explicit conditions on the parameters of the reaction kinetics for the stability of the resulting patterns, showing also that they are stable in the case of a supercritical Turing bifurcation. An interesting observation is that the presence of the gradient relaxes the condition on the relative sizes of the diffusion coefficients of the activator and inhibitor. Patterns can appear if the ratio is closer to one; i.e. the external gradient effectively enlarges the Turing space of the model. 4 Mathematical models related to limb development In the last 50 years, a number of mathematical and computational models have been proposed to investigate the process of vertebrate embryonic limb growth and chondrogenesis. There is enduring interest in the development and refinement of such models due to the ever increasing amount of biological data as well as computational power. To investigate the mechanisms of limb development and pattern formation, mathematical modeling is extremely useful. It allows to rigorously and quantitatively test ideas about mechanisms of growth and pattern formation, potentially revealing which interactions between system components (gene expression, cell motion, cell-cell adhesion, morphogen signaling etc.) are indispensable and which are not. The interplay between experiment and model is a process of reciprocal influence: Models are not only established from the experiments, but often lead to additional experiments verifying new hypotheses. Except for very simple models, it is often difficult to derive quantitative conclusions based on analytical mathematical techniques. Thus simulations, often requiring sophisticated numerical methods, are usually an indispensable tool in the analysis of these models. Mathematical models concerned with limb development can be divided into two groups [189]: (a) Models of the growth and shape of the limb bud. (b) Models of the formation of pre-skeletal patterns An overview of the models we discuss in this paper is given in Table 1. Model Dillon and Othmer Boehm et al. List of models related to Described processes Limb growth described via Navier-Stokes equations coupled with concentration of morphogens. Limb growth model described by Navier-Stokes equation with mesenchyme treated as a viscous incompressible fluid. growth and shape of limb bud References Section Importance/Results [32, 33] 4.1 Modeling the interactions between morphogens and their effects on limb growth. [16] 10 4.1 Directional cell activities (motion and proliferation) are crucial for obtaining realistic limb bud shape. ฀฀฀฀฀฀฀฀฀ List of models related to formation of skeleton patterns via morphogen interaction Dillon et al. Spatiotemporal evolution of Shh [32, 34] 4.2 Non-isomorphic morphogen patgradient terning model. Modeling to validate the experimental findings that ZPA is a source of Shh morphogen. Armstrong Cell-cell adhesion [5] 4.2 Most influential PDE model of et al. cell-cell adhesion. Modeling cell aggregation due to cell-cell adhesion; modeling cell sorting in two interacting populations. Galectin Chondrogenic pattern formation [53] 4.2 Isomorphic morphogen patternmodel by in the chick in vitro via a network ing model. Modeling chicken Glimm et al. of galectins. galectins CG-1A and CG-8 by a system of mixed parabolichyperbolic equations incorporating a non-local flux term describing cell-cell adhesion; based on experimental data on galectin interactions. BMP model Chondrogenic digit pattern for- [6] 4.2 Isomorphic morphogen patternby Badugu mation in the mouse in vivo via ing model. Core mechanism et al. a reaction-diffusion system with is Schnakenberg-like system with BMP, its receptor and FGF. BMP as activator, its receptor as inhibitor. BSW Chondrogenic digit pattern in [149] 4.2 Isomorphic morphogen patternmodel by the mouse in vivo via a reactioning model. Essentially linear reRaspopovic diffusion system with BMP, Sox9 action kinetics; parameters (netet al. and WNT work interactions of BMP, Sox9, WNT) elucidated via matching experimental in-phase and outof-phase expression patterns. Table 1: Overview of models discussed in detail. 4.1 Models of the growth and shape of the limb bud Models that describe the process of vertebrate limb growth and skeletal pattern formation have attracted the interest of researchers since at least the late 1960s. One of the main questions here is what the basic mechanism of limb outgrowth is, and how it is initiated in the early limb [52]. Pioneering work in modeling the growth of the embryonic limb bud was done by Ede and Law in the late 1960s [41, 52]. They proposed a simple model involving cell proliferation and motion to investigate if differences in the rates of cell proliferation at the proximal and distal parts of the limb are necessary to maintain the shape of the limb. Ede and Law concluded from their simulations that while these differences do not affect the early shaping of the limb [52], later on during growth, they are crucial for maintaining the characteristic paddle shape. It was observed by them that during this period of growth, more distal cells are dividing more frequently than the proximal cells with a tendency of cells to move slightly distally. The conclusion that the estimated differences in the rate of proliferation were sufficient to maintain the paddle shape, drawn from Ede and Law’s model, has been widely accepted and used in several subsequent simulations [32, 34, 66, 130]. 11 ฀฀฀฀฀฀฀฀฀ Beginning in the late 2000s, this view was challenged by other authors, namely Murea and Hentschel [107, 189] and Boehm et al. [16]; see the discussion of the latter below. We briefly review here the two arguably most influential PDE models of limb growth, namely the work by Dillon and Othmer [32, 33, 34] and the work by Boehm et al. [16]. (See also Table 1.) In both cases, tissue growth and movement is modeled via Navier-Stokes equations with a moving boundary. Other notable models are those of Murea and Hentschel [107] and Morishita et al. [104]. All these models have been reviewed in [189], and to our best knowledge, there has not been much significant new development since then. • Dillon and Othmer’s model The first mathematical and computational model of limb growth and experimentally observed gene expression patterns using the full apparatus of fluid dynamics modeling was proposed by Dillon and Othmer in [32, 33, 189]. It is a mathematical model of cell fluid flow coupled with elastic boundaries representing the mechanical and biochemical properties of the ectoderm surrounding the limb mesoblast. The fluid motion is described by the Navier-Stokes equations ∇ · u = S(c(x, t)), ρ ∂u 1 + ρ(u · ∇)u = −∇p + µ(∇2 u + ∇S) + ρF. ∂t 3 (2) Here the term S can be interpreted as the local source of growth. It depends on the concentrations c1 and c2 , where c1 corresponds to the morphogen secreted from the AER and c2 corresponds to the morphogen secreted from the ZPA (the pair (c1 , c2 ) is denoted by c above), as well as the location of the tissue within the limb bud, and the age of the limb. x is the position within the limb and t is the age of the limb, ρ is fluid density, p its pressure and µ fluid viscosity. The two morphogens (whose possible identities are FGFs and Shh with sources at the AER and the ZPA, respectively) are governed by a reaction-diffusion-advection system ∂c + ∇ · (uc) = D∇2 c + R(c), ∂t (3) where D is the diffusion matrix for the morphogens. The morphogens are advected with the limb mesoderm, which moves at the local velocity u . The production rates are denoted by R(c) The AER morphogen is only produced in the AER(Ω1 ) and the ZPA morphogen is in the ZPA (Ω2 ). Thus R = (R1 , R2 ) has the form   x ∈ Ωk , rk (c) − κk ck Rk =   −κk ck Otherwise where rk (c) > 0 except at c = 0 and the corresponding Michaelis-Menten kinetics are as follows with rate constants Vk and Kk : c2 c1 r1 (c) = V1 , r2 (c) = V2 . K1 + c 2 K2 + c 1 The source term S has the form: S = s1 c1 + s2 where s1 and s2 are constants. Hence the local growth rate linearly depends on the local concentration of c1 . Even though it was published 20 years ago, Dillon and Othmer’s model is still one of the most thorough model of the interplay of the growth of the limb bud and the spatiotemporal evolution of morphogen concentrations during growth, using current experimental data for estimating parameters. • Model of Boehm et al. In [16], the authors proposed a similar fluid dynamics model like the models in [32, 104], to study numerically the elongation process, using finite element computational method in three dimensions incorporating quantitative 12 ฀฀฀฀฀฀฀฀฀ data on shape changes and proliferation rates. The authors revisited the question of whether a proximal-distal gradient of isotropic cell proliferation rates is crucial for limb outgrowth, as argued first by Ede and Law [41]. Employing a combination of numerical simulations and experimental imaging they came to the conclusion that this is not the case, but rather that directional cell activity (directed cell motion and non-isotropic proliferation) is responsible for limb outgrowth and limb shape. This model, based on Navier-Stokes equations, incorporates the patterns of cell division considering the mesenchyme as a viscous incompressible fluid whose volume increases with s (a distributed material source term, in fact, s represents the proportional volumetric growth per unit time): 1 ∂v + ∇p − ∇ · [∇v] = 0, ∂t Re (4) ∇ · v = s, Here the term v represents velocity and the pressure is represented by the term p. Re is Reynolds number. The authors compared the results of simulations to actual three-dimensional data obtained from imaging of real limb buds. They found that if they used empirical data for the cell proliferation rates, the resulting simulated limbs had a different shape than the real ones: They were more bulbous and less flat than the real ones. Taking the opposite approach as well, the authors also determined the parameters that gave the closest match of the simulated shapes to the real limb shapes. They found that this resulted in source terms that were negative in large parts of the limb bud, corresponding to shrinkage, which is not in accordance with the empirical data. The conclusion from this model is that a gradient of undirected cell proliferation is not sufficient to explain the characteristic paddle shape, and that other, isotropic effects must play a crucial role according to the authors. 4.2 Models of skeletal pattern formation via morphogens We concentrate here on PDE models of skeletal pattern formation based on morphogens. An overview of the models discussed in depth is given in Table 1. In general, one can distinguish two groups of models describing the dynamics of morphogens leading to skeletal pattern formation [189]: ‘isomorphic’ and ‘non-isomorphic’ ones. In isomorphic models, it is assumed that there is an isomorphic (shape preserving) mapping between the spatial distribution of morphogens and the skeletal elements. Thus the cell pattern is a direct ’copy’ of the spatial distribution of the morphogen concentration. This isomorphism can be realized by means of different simple phenomena, e.g. by enhancing the differentiation of mesenchymal cells into cartilage cells (which then differentiate into chondrocytes), attracting (chemotactic) effects, and others. In the non-isomorphic case, there is no a straightforward spatial correspondence between the chemical prepattern and the final location of skeletal elements. In this case, the morphogen dynamics guide patterning mainly by initiating different regulatory pathways [136]. Non-isomorphic models have seen somewhat less interest from mathematical modelers in limb development. This is likely because here the complexity of pattern formation lies almost entirely in the ”program” that the cells execute in response to the morphogen concentration. In isomorphic models, there is an emergent process by which the interplay of many ”simple” agents give rise to complex patterns, which is both more interesting and more intuitive to mathematical modelers. Correspondingly, we review only one non-isomorphic model, that of Dillon et al., which we believe to be the most comprehensive and molecularly detailed model, taking into account both the growth process and the spread of morphogens within the limb bud. (Other notable models of this kind are Hirashima et al. [66] and more recently a model of Shh digit specification by Woolley et al. [207].) The remaining models we discuss are recent models of isomorphic type. The selection of models is by no means complete, and we mostly concentrate on models proposed within the last 15 years or so. Many other influential models are not treated in depth here, e.g. earlier reaction-diffusion models such as to Newman and Frisch’s model [110] and the work by Hentschel et al. [64], the mechanochemical models of Murray and Oster [174, 173] or the boundary model by Meinhardt [94, 96]. 13 ฀฀฀฀฀฀฀฀฀ Among these models, the most molecularly detailed model is that by Hentschel et al. [64], which incorporates a system of activator-inhibitor reaction-diffusion equations together with an explicitly modeled cell density. The activator molecule is identified as TGF-β, and there are three cell types: Two that differ in the expression of two different FGF receptors, and a third one that produces fibronectin. The authors showed that the system was morphostatic, i.e. a prepattern in the concentrations of activator and inhibitor can arise via a Turing bifurcation even if the cell density is kept constant. This model in particular spawned further computational and mathematical investigations [7, 4, 20, 24]; in particular, it was shown that if solved on realistic limb shapes, the system can recapitulate the results of many different experimental manipulations [192]. For a mathematical in-depth treatment of the models mentioned above, as well as others based on discrete models, see the 2013 survey by Zhang et al [189]; for a description for a more general audience, see Glimm et al. [52]. •Model of Dillon et al. Based on their model of the limb growth, Dillon and co-workers incorporated the investigation of the spatiotemporal spread of morphogens, specifically Shh, in a series of papers beginning with [31]. In the framework mentioned in the introduction to this section, we may classify it essentially related to non-isomorphic models, even though the process of condensation of pre-cartilage cells was not explicitly modeled. The relation between Shh receptor Patched (Ptc) and the associated membrane signal transduction factor Smo was discussed by Dillon et al. in [32, 34]. Terms related to the influence of Shh receptor and mediator proteins were coupled (in addition to terms for FGF in 2-D) with Navier-Stokes equations for the growth of the limb bud. Using this model, the authors simulated the effects of ectopic sources of Shh and compares the results with the experiments. The model analyzes the interaction between Shh, Shh transmembrane receptor Patched (Ptc) and Smoothened (Smo) (see [29, 183]), a transmembrane protein mediating Shh signaling through phosphorylation of the Gli family of transcription factors. For example, the interaction between Shh and Ptc can be described, neglecting the other terms, as in [189]: ∂ [Shh] = [diffusion of Shh] − [association of Shh and Ptc] ∂t + [disassociation of Shh-Ptc complex] − [degradation of Shh] + [Shh production], ∂ [Ptc] = −[association of Shh and Ptc] + [disassociation of Shh-Ptc complex] ∂t − [association of Smo and Ptc] + [disassociation and degradation of Smo-Ptc complex] + [Ptc productions by itself and by Smo] − [degradation of Ptc], ∂ [(Shh-Ptc complex)] = [association of Shh and Ptc] ∂t − [disassociation and degradation of Shh-Ptc complex]. 14 ฀฀฀฀฀฀฀฀฀ Figure 3: The computational and experimental results of Patched (Ptc) responses to Sonic hedgehog (Shh) bead implants (upper panels) and ZPA tissue implants (lower panels). (Upper) The rescaled figures of numerical simulations of Ptc concentration 2, 6, and 18 h after bead implants (A, C, and E, respectively), whereas in lower, Ptc concentration 12, 16, and 20 h after tissue implant (A, C, and E, respectively). Experimental results are from Drossopoulou et al. [35] for ptc transcript expression 2, 6, and 16 h post-bead implants and for the ZPA grafts, 4, 8, and 16 h post-implant (B, D, and F, respectively). (Modified from Dillon et al. [34]). It is also assumed in the model that Smo has two forms: active (associated with free Ptc) and inactive (interacting with Shh-Ptc) and that different forms of Shh have the same diffusion constants. The main idea of this work was to compare the effects of implantation of ectopic sources of Shh (Shh beads) and ectopic the ZPA tissue. The excellent agreement between the numerical simulations within the model and experimental findings, seem to prove both the relative validity of the model assumptions as well as the fact that the ZPA is a source of Shh morphogen. A comparison between the model and the experimental results is shown in Fig. 3. • Model of Armstrong et al. In [5], Armstrong et al. designed a nonlinear partial differential equation model to describe the phenomenon of cell-cell adhesion. This model was later analyzed mathematically by Dyson et al. in [38] and [39]. Cell-cell adhesion and cell-extracellular matrix adhesion is modeled via non-local flux term. Cell adhesion is crucial in many biological contexts leading to different kinds of pattern formation. To observe the formation of aggregations of cells or cell clusters for an initially distributed cell population with strong cell-cell adhesion, Armstrong et al. considered first a population of one type of cell with uniform adhesive properties. In one spatial dimension the model reads [39]: ∂u(x, t) ∂u2 (x, t) ∂ = − (uK(u)) 2 ∂t ∂x ∂x where K(u) = α Z (5) R g(u(x + x0 )) ω(x0 ) dx0 . (6) −R Here u(x, t) denotes cell density. Cell-cell adhesion is represented by an adhesion flux term encoded by the term K(u). The magnitude of the effects of adhesion forces on the local cell density are described by the term g(u(x + x0 )). The strength of cell-cell adhesive force is represented by the positive parameter α and ω(x0 ) characterizes the direction and magnitude of the force (changing with x0 ). R is the sensing radius of the cells. As a result, K(u) can essentially be thought of as a weighted sum of adhesive forces exerted by neighboring cells; see Figure 4. Equation (5) is considered for x ∈ (−∞, ∞), with initial condition u(x, 0) = u0 (x). In [5], the term g(u(x + x0 )) is considered to be either linear i.e., g(u(x + x0 )) = u(x + x0 ) or of logistic type: ( u(x + x0 )(1 − u(x + x0 )/M ) if u(x + x0 ) < M g(u(x + x0 )) = 0, otherwise 15 ฀฀฀฀฀฀฀฀฀ Figure 4: Schematic representation of the contributions to the computation of K(u) in (6) due to adhesives forces in the Armstrong et al. model [5]. The case ω(x) = sgn(x) as in (7) is shown. The flux term K(u) at x can be thought of as the net sum of the sketched forces over an interval centered at x whose radius is the sensing radius R. where M represents a maximum density of the population. The form of ω(x0 ) reflects the adhesive nature of cell-cell interaction; in particular, it encodes the attractive influence of the cell density at x0 on a cell at the origin x = 0. The simplest choice for ω(x0 ) is a step function ( −1 − R < x0 < 0 (7) ω(x0 ) = sgn(x0 ) = 1, 0 < x0 < R Other forms for ω(x0 ) can be considered in the model, in particular forms where the magnitude of ω(x0 ) decreases with |x0 |. In the same paper [5], a system describing two populations of different types of cells is set up. The proposed model has the form: ut = uxx − (uKu (u, v))x (8) vt = vxx − (vKv (u, v))x where and Ku (u, v) = Su | Kv (u, v) = Sv | Z Z 1 guu (u(x + x0 ), v(x + x0 ))ωuu (x0 )dx0 + C −1 {z } | Z 1 −1 u−u adhesion 1 gvv (u(x + x0 ), v(x + x0 ))ωvv (x0 )dx0 + C −1 {z } | v−u adhesion Z guv (u(x + x0 ), v(x + x0 ))ωuv (x0 )dx0 {z } (9) u−v adhesion 1 −1 gvu (u(x + x0 ), v(x + x0 ))ωvu (x0 )dx0 . {z } (10) v−v adhesion Here u(x, t) and v(x, t) denote the two populations’ cell densities. The terms Ku (u, v) and Kv (u, v) represent the adhesion of cells. In Ku , the first term describes self-population adhesion of the first type of cells and the other term indicates cross-population adhesion. Similarly, in Kv , the first term represents self-population adhesion of the second type of cells and the other term indicates cross-population adhesion. The self-adhesive strength of the population u and v are represented by the constants Su and Sv , while the cross-adhesive strength between the populations is C. In the same paper [5], the models describing one population and two interacting populations in one spatial dimension were extended to two spatial dimensions. In this case, the one population model takes the form ut = ∇2 u − ∇ · (uK(u)) where K(u) = α Z 1 0 Z (11) 2π g(u(x + r η)) Ω(r) η r dθ dr. 0 16 (12) ฀฀฀฀฀฀฀฀฀ Here x ∈ R2 denotes the position of the cell, x + r η denotes the position of other cells within the sensing disc of radius R scaled to R = 1 and η = η(r, θ) is the unit outward normal to the circle C(x, r). The term η Ω(r) replaces the function ω(x0 ) present in spatially one dimensional model (5)-(6). One can extend the two-population model described above to two spatial dimensions in a straightforward way. In the one population and spatially one dimensional model it is relatively easy to obtain an intuition about the influence of the non-local adhesion terms on solutions. In this case, one can approximate the integrals by local differential terms obtained formally by the expansion of the cell density within the integral [5]. So, substitution of the expansion of u(x, t) u(x + x0 , t) = u(x, t) + x0 ux (x0 , t) + x20 uxx (x0 , t) + · · · 2 into the integral K(u) changes equation (5) for g(u) = u to ut = uxx − Aα[uux ]x − Bα[uuxxx ]x + Φ(x50 ) (13) R1 1 R1 3 x ω(x0 )dx0 are both positive. Let us note that, as ω(x0 ) is odd, where A = −1 x0 ω(x0 )dx0 and B = 6 −1 0 R1 k then −1 x0 ω(x0 )dx0 = 0 for all even integer k. The second order term in equation (13) depends on the first spatial derivative ux as in the models of chemotaxis [4][72]. The above approximation implies cell movement up the gradient (towards higher concentration) of the cell density, therefore cell aggregating may be observed in solutions to equation (5)–(6). On the other hand, the fourth order term has a dampening effect, therefore the non-local term may help in cell aggregations without creating singularities and blows up phenomena. In a similar way, the PDE approximation can be done for two interacting populations. Although the model was constructed initially in one and two dimensions by Armstrong et al. in [5], still higher dimensional studies are important in the context of cancer modeling because of the manner in which cancer cells invade tissue and the possibility of phenomena such as fingering at the invasive front [39]. Therefore the above model was extended in N -dimensions by Dyson et al. in [39] in a straightforward way as follows: ! Z ∂u(x, t) = D∆u(x, t) − ∇ · u(x, t) g(u(x + ξ, t)) ξ ω(|ξ|) dξ + f (u(x, t)) (14) | {z } | {z } ∂t Bρ random motility {z } cell loss and gain | cell adhesion for x ∈ RN , t > 0 and Bρ denoting the N -dimensional ball centered at 0 and of radius ρ with initial condition u(x, 0) = u0 (x), x ∈ RN . Just as in the 1- and 2-dimensional cases in (5) and (8), respectively, cell-adhesion is model via a flux term that represents a weighted average direction of motion taken over an n−dimensional ball. The weighting is given by the cell density. The outcome is an effective flux in a direction of increase of cell density. The function ω(|ξ|) describes how the weighting depends on the distance from the center point of the ball. The simplest case is ω(|ξ|) = 1/|ξ|, but more complicated forms are possible. It is assumed that f (0) = 0 (for biological reason). Also, it is assumed that there is a number P1 > 0 such that f (u) > 0 for u ∈ (0, P1 ), and f (u) < 0 for u > P1 . These assumptions suggest that there is cell gain at lower densities while in higher densities, due to the effects of crowding, cell loss occurs more rapidly in compared to the generation of new cells via division. • Model of Glimm et al. Recently in [53], Glimm et al. proposed a new model related to bone formation based upon the results of an 17 ฀฀฀฀฀฀฀฀฀ experimental paper [15] (see also [84]). The mathematical formulation of this model incorporates a non-local flux term describing cell-cell adhesion forces based on the approach of Armstrong et al. in [5], and has the form of a structured population model with diffusion. In this context, the model proposed by Glimm et al. in [53] differs from models discussed in previous sections, and it has far reaching consequences for its mathematical analysis. This is due to the presence of hyperbolic terms inside the reaction-diffusion equations. (In Eq.(19) ∂ ∂ these are the terms (γ̃(cu1 , cu8 , T1 )R) and (δ̃(cu8 , T8 )R) .) ∂T1 ∂T8 In [15, 84], a crucial role of the galectin family of proteins has been reported to regulate and mediate chondrogenesis during avian limb development. One of these proteins, chicken galectin (CG)-1A acts as a morphogen and a cell adhesion factor, and another, CG-8, is a morphogen and effectively an anti-adhesive factor. CG-1A and CG-8 and their respective counterreceptors are produced by all limb bud mesenchymal cells. The model proposed in [53] not only explains the interactions between CG-1A and CG-8 to form spatial patterns of condensations during cell aggregation and bone formation but also provides the crucial insights of the pattern formation from a physical perspective that the limb skeletal patterning is a morphodynamic process and thus depends on mesenchymal cell motility. In [15] it was found in experiments that CG-1A and CG-8 formed a positive feedback loop, i.e. CG-1A upregulates CG-8 gene expression and CG-8 upregulates CG-1A gene expression. At the same time, CG-1A was found to enhance precartilage condensations and CG-8 suppresses condensations. It was hypothesized that this puzzling behavior could be due to a competition effect: CG-1A and CG-8 may ’compete’ for binding to the same cell-membrane bound counterreceptor. To model this process, Glimm et al. [53] set up a system of equations that involves the cell density, the spread of the diffusible galectins CG-1A and CG-8 and binding and unbinding of the galectins to membranebound counterrecptors. There are two types of counterrecptors: A shared one, to which both CG-1A and CG-8 can bind, and one that is unique to CG-8. The core assumptions are as follows: 1. Mesenchymal cells undergo random motion, biased by cell-cell adhesion. 2. Cell-cell adhesion is mediated by the complex of CG-1A bound to the shared counterreceptor. 3. Binding of CG-1A to the shared counterrecptor also upregulates CG-8 production. 4. Binding of CG-8 to its (unique) counterrecptor upregulates CG-1A production. 5. Binding of CG-8 to the shared counterrecptor has no direct relevant regulatory effect. These assumptions are schematically shown in Figure 5. The fifth point means that CG-8 can effectively ”block”” the shared counterrecptors for CG-1A and thus encodes the competition effect. Figure 5: Schematic illustration of the key players and their basic roles in the galectins model proposed in [53]. (Modified from [53]) 18 ฀฀฀฀฀฀฀฀฀ These assumptions give the following system of integro-differential equations: ∂R = DR ∇2 R − ∇ · (RK(R)) | {z } | {z } ∂t cell diffusion cell-cell adhesion ∂ ∂ ∂ (αR) − 8 (β8 R) − 1 (β1 R) − ∂c1 ∂c8 ∂c8 {z } | binding/unbinding of galectins to counterreceptors ∂ ∂ − [(γ − α − β1 )R] − [(δ − β8 )R] ∂l1 ∂l8 {z } | (15) change in counterreceptors ∂cu1 = D1 ∇2 cu1 | {z } ∂t diffusion +ν̄ | Z c88 RdP {z − | } Z αRdP {z } positive feedback of CG-8 on prod. of CG-1A binding CG-1A to its counterreceptor ∂cu8 = D8 ∇2 cu8 | {z } ∂t diffusion +µ̄ | Z − | c1 RdP {z } Z β1 RdP {z } positive feedback of CG-1A on prod. of CG-8 binding CG-8 to counterreceptor −π̄1 cu | {z 1} (16) degradation −π̄8 cu | {z 8} (17) degradation A term with a bar over it (e.g. µ̄) denotes a constant. Here cu1 = cu1 (t, x) is concentration of freely diffusible CG-1A and cu8 = cu8 (t, x) is concentration of freely diffusible CG-8, whereas R = R(t, x, c88 , c81 , l1 , l8 ) denotes morphogenetic cell density. So the cell density R not only depends on position and time, but also on the concentration of various counterreceptors bound to the cells’ membranes. Thus the counterreceptor concentrations enter the model as structure variables. This approaches is analogous to other structured population models, e.g. models where population density also depends on the age of the individuals. Integration over the space of structure variables is denoted by dP = dc1 dc88 dc81 dℓ1 dℓ8 . As we mentioned above, the cell-cell adhesion term ∇ · (RK(R)) is formulated based on the approach of [5] (see the previous subsection, especially (6)) and is defined as Z Z Z r 8 1 (18) K(R(t, x, c1 , c8 , c8 , l1 , l8 )) = ᾱK c1 c̃1 σ(R(t, x+r, c̃1 , c̃88 , c̃18 , ˜l1 , ˜l8 ))dP̃ dn r. |r| Dρ0 Here ᾱK is a constant which represents the strength of the adhesion and σ(R) has either linear or logistic form. The effective adhesion force on a cell at location x depends on the product of the concentration of bound CG-1A on the cell and the concentration of bound CG-1A at locations x + r, where the distance vector r varies over the n-dimensional (n = 1; 2; 3) ball Dρ0 (x) centered at x. The radius ρ0 is the “sensing” radius, which is a measure of the characteristic distance for adhesion; cells at distance greater than ρ0 do not contribute to the adhesion forces (see figure 4). Thus the integral in the cell adhesion flux represents a weighted average direction of motion, as in the Armstrong model [5]. However, the weighting does not depend on cell density directly, but it is modulated by the concentration of CG-1A on the cells’ surfaces, given by the variable c1 . The term γ −α−β1 models the rate at which the membrane-bound concentration of the shared counterreceptors which are not bound to either galectin changes. The change is due to the expression of new counterreceptors by the cells and degradation (leading to the effective rate γ), the binding and unbinding of the counterreceptor to CG-1A (the rate α) and the binding and unbinding of the counterreceptor to CG-8 (the rate β1 ). Similarly, the term δ − β8 denotes the rate at which the membrane-bound concentration of its own counterreceptors which are not bound to CG-8 galectin changes. (γ, α, β1 and β8 in turn depend on the galectin concentrations. See [53] for the detailed functional relationships.) Assuming “fast galectin binding” to the countereceptors, the following simplified system was obtained in [53] from the full model using two auxiliary variables T1 , T8 : The total concentration of CG-1As counterreceptors is given by T1 = c1 + c18 + l1 19 ฀฀฀฀฀฀฀฀฀ and the total concentration of CG-8s counterreceptors is given by T8 = c88 + l8 . The morphogenetic density (cell density) then becomes a function of T1 and T2 , i.e. R = R(t, x, T1 , T8 ), and the resulting system is as follows: ∂ ∂ ∂R = dR ∇2 R − ∇ · (RK(R)) − (γ̃(cu1 , cu8 , T1 )R) − (δ̃(cu8 , T8 )R), ∂t ∂T1 ∂T8 Z ∞Z ∞ ∂cu1 2 u c88 R dT1 dT8 − cu1 , = ∇ c1 + ν̃ ∂t 0 0 ∂cu8 = ∇2 cu8 + µ̃ ∂t with Z ∞ 0 c88 = c88 (t, x, T8 ) = cu8 T8 , 1 + cu8 c1 = c1 (t, x, T1 ) = cu1 T1 , 1 + f cu8 + cu1 γ̃(cu1 , cu8 , T1 ) =  2cu1 − γ̃2 c1 + c̄˜1 δ̃(cu8 , T8 ) = 1 − δ̃2  Z (19) (20) ∞ 0 c1 R dT1 dT8 − π̃8 cu8 , (21) c1 , cu1 T8 , 1 + cu8 K[R, cu1 , cu8 ](t, x, T1 , T8 ) = Ψ(δ; dist(x, ∂Ω)) Z ∞Z ∞Z s c1 (t, s, T̃1 )σ̃(R(t, s, T̃1 , T̃8 )) ds dT̃1 dT̃8 α̃K c1 (t, x, T1 ) |s| 0 Dr0 (0) 0 Here α̃K is a constant which represents the strength of the adhesion and for some δ > 0 sufficiently small, Ψ(δ; ·) is a smooth cut-off function such that Ψ(δ; y) ≡ 1 for y ≥ 2δ, Ψ(δ; y) ≡ 0 for y ≤ δ. Numerical simulations in [53] show evidence that the system (19)-(20)-(21) can produce spatial patterns in the morphogenetic density R(t; x; T1 ; T8 ) for a wide range of the model parameters. The cell-cell adhesion flux term plays a crucial role in this spatial pattern formation, as can be observed in Figure 6. Moreover, the authors showed that increasing initial CG-1A concentration led to more condensations, in agreement with experiments. Increasing the initial CG-8 concentration led to fewer condensations, but only in certain regions of parameters space. The key to the inhibitory effect of CG-8 is a competition between CG-1A and CG-8 for the shared counterreceptors. It is thus crucial that the percentage of unbound counterrecptors is in some sense low. Simulations show that the critical value lies somewhere between 30% and 10% unbound shared counterrecptors. The model was also used to investigate the role of synchronized oscillations hes1 oscillations in the earliest stages of precartilage condensations [202]. Incorporating oscillating hes1 states in the model, it was shown that asynchronous hes1 states across the domain produced much less regular patterns than synchronously oscillating ones. It was thus argued that these oscillations are a secondary mechanism to ensure regularity of patterns via enforcing simultaneous appearance of the patterns across the developmental field. Furthermore, the model was used in an investigation of the phylogenetic origin of the galectin patterning network. It was found that Gal-8 underwent purifying selection in the sarcopterygians (lobed-fin fishes; the ancestors of all tetrapods). On the level of the model, it was shown that tighter control over the production rates of Gal-8 meant less variability in 20 ฀฀฀฀฀฀฀฀฀ the number of precartilage condensations, thus pointing to a connection with the relatively small variability of the stereotypical limb skeleton within the tetrapods compared to e.g. the fins of ray-finned fishes. RR Figure 6: [53] Distribution of cell density R(t; x; T1 ; T8 )dT1 dT8 at times t = 0 and t = 1 for different values of the cell-cell adhesion constant α̃K . Other values are r0 = 0.04, δ̃2 = 1, γ̃2 = 1, c̄˜1 = 1, f = 0.8, dR = 0.04, π̃8 = 1, ν̃ = 0.8, µ̃ = 2. Initial distributions are represented by dashed lines and distributions at t = 1 are by solid lines. As α̃K is increased, periodic patterns appear as a result of random spatial noise added to the initial distribution. Here periodic boundary conditions are used, so that the positions x = 0 and x = 1 denote the same physical point. (Modified from [53].) • Model of Iber-Badugu In [6], Iber et al. proposed a model for the mechanism of patterning of digits in mouse limb, based on BMPreceptor interaction. BMP signaling along with FGF gradient are important for digits formation [13, 10, 186]. Iber and the coworkers focused on the interactions of BMP (denoted as B in the model), its receptor (denoted as R) and FGF (denoted as F ) under which the digits emerge in the autopod. These interactions are shown in Figure 7. BMP and FGF diffuse throughout the extracellular matrix, whereas the BMP receptors are confined to cell membranes. In the model, the diffusion of the ligand bound receptors, residing mainly within the cell, denoted by C, was considered negligible as they are internalized rapidly [73]. The rate of BMP receptor binding is proportional to R2 B, as BMPs are dimers and can bind two receptors. In the limb bud, BMP2 expression is inhibited by BMP2 signaling. Basing on this fact, the rate of BMP production is assumed to be of the form KB PB . In this way, the BMP and BMP-receptor dynamics model takes the following form: KB + [C] KB ¯ −dB [B] −kon [R]2 [B] + kof f [C] [Ḃ] = D̄B ∆[B] + PB | {z } KB + [C] | {z } | {z } {z } degradation | diffusion complex formation (22) production [Ċ] = kon [R]2 [B] − kof f [C] −dC [C] , {z } | {z } | complex formation degradation where kon and kof f are the binding and dissociation rate constants respectively. D̄B is the diffusion coefficient for BMP molecules. Production of receptor depends on the concentration of C as the signaling of BMP-bound receptors positively regulates receptor production. 21 ฀฀฀฀฀฀฀฀฀ ¯ + pR + pC ([C]) −dR [R] −2kon [R]2 [B] + kof f [C] , [Ṙ] = D̄R ∆[R] {z } | {z } | | {z } | {z } diffusion production degradation (23) complex formation where pR and pC are constants. The receptor ligand complex assumes its quasi steady-state almost instantaneously as the dynamics of receptors ligands complex are much faster than the dynamics of BMP, hence the concentration of bound receptors are proportional to R2 B, i.e., [C] ∼ kon kon [R]2 [B] = KC [R]2 [B]; KC = kof f + dC kof f + dC It was shown in [6] that the system 22–23 was sufficient to produce pattern and it can be reduced to classical Schnakenberg type if pC = 2dC and dB = 0. Expression of BMP is induced by the FGF signaling, so the model was extended by Badugu et al. in [6] by introducing the production rate PB as a function of FGF concentration F , PB (F ) = pb + p∗B [F ]n KB n K + [C] , [F ]n + KBF B where KBF and n are the Hill constant and Hill coefficient repectively. It has found that BMP-bound receptors signaling may stimulate as well as inhibit FGF-dependent processes [37, 125]. To account for this, FGF activity was chosen as [C]n KFn 2 PF ([C]) = pF , [C]n + KFn 1 [C]n + KFn 2 where KF 1 ≪ KF 2 are the Hill constants for the activation and inhibition impacts of BMP signaling. So the dynamics of FGF is as follows, KFn 2 [C]n − dF [F ] , [Ḟ ] = D̄F ∆[F ] + pF n n n | {z } | {z } [C] + KF 1 [C] + KFn 2 {z } degradation | diffusion (24) production where D̄j (j = B, R, F ) are the diffusion coefficients with D̄R ≪ D̄B , D̄F . The shape of the domain was extracted from limb bud images at E12.5 and hence the system of equations were solved on a growing domain (see Figure. 7). Zero-flux boundary conditions for B and R were assumed, while FGF production was implemented as a flux boundary condition, ~n · ∇F = ρF (R2 B)n κn2 , n 2 + κ1 (R B)n + κn2 (R2 B)n where ~n is the unit normal vector. An additional remark should be made here concerning a specific role of Sox9 gene which is not explicitly taken into account in this model. Sox9 serves an important role in digit pattern formation, but according to [6], only as a marker of endochondral differentiation. To be more precise, BMP-2 signalling stimulates Sox9 expression and this enhances Noggin expression, which has a negative impact on BMP signalling by changing BMP into an inactive complexes [6, 193]. One criticism of the model by Badugu et al., raised in Raspopovic et al. [132], is that it includes diffusion of BMP receptors R throughout the ECM, although in fact the receptors are membrane-bound and do not diffuse, but are rather advected with the cells. It is unclear whether pattern formation is possible without this diffusion assumption. However, Kurics et. al [83] showed for a related case that pattern formation is possible in receptor-ligand reaction-diffusion systems if the receptor is confined to diffuse only on certain subdomains, interpreted as cell, and that this even enlarges the Turing space, i.e. the region of parameter space for which Turing pattern are possible. 22 ฀฀฀฀฀฀฀฀฀ Figure 7: (A) Schematics of interactions in the model of Badugu et al. [6]. The interactions of BMP, its receptor and FGF form a loop. In fact, BMP receptor complexs are formed due to the reversible binding of BMP and its receptor, which induce the production of receptors and enhances the FGF activity, while FGF induces BMP expression. (B) Therefore, BMP has a positive impact on receptor as well as on FGF, both via BMP-receptor complexs (excluded from the figure). It should be mentioned here that receptors act as auto-activatory when they are bounded by BMP, whereas BMP enhances self-decay by receptor binding, and hence they are autoinhibitory, and a mutual enhancment is observed between BMP and FGF. (C) The domain of computation, based on the shape of a limb bud, at E11.5. The radial axes of the elliptic bud are denoted as Ri , (i = 1, 2, 3, 4). Height and width of the stalk are represented by H0 and W0 , respectively. In the stalk, the height of the domain is denoted by H1 . The expression of BMP is upgraded at the height of the domain in the stalk. (Modified from [6]) • BSW Model (Raspopovic et al.) In Turing-type reaction-diffusion systems, the wavelength of the patterns produced typically shows a strong dependence on the parameters. So changes in the parameters will lead to changes in the wavelength. This observation is one of the principal objections to the applicability of the Turing mechanism to bone pattern formation in the vertebrate limb. After all, the number of elements is very stable and e.g. derivation from the number of digits (in the form of extra digits (polydactyly) or missing or fused digits (syndactyly)) relatively rare. Proponents of reaction-diffusion mechanisms have pointed out that the prevalence of congenital limb defects is high relative to other defects, and higher than e.g. either Down syndrome or cleft palate, with prevalence reported as 22.7 per 10,000 birth in a study in Thailand [69]; in a large study in Hungary, this number was reported as 1 in 1816, or about 5.5 per 10,000 births ([42] as cited in [52]). Still, these incidences are quite low on an absolute scale and do not put in question the basic argument of the robustness of the limb patterning network. A more convincing reply was presented in a remarkable study by Sheth et al. (2012) [149]. The authors generated a series of mouse mutants which lacked alleles for three genes that have beeen shown to be important in digit formation, namely the distal Hox genes Hoxd11-13, Hoxa13, and Gli3, the major mediator of Sonic Hedgehog sigmaling in limb development. Sheth et al present a total of 15 mutant types: 5 combinations of different Hox gene deletions; namely Hoxa13+/+ ; Hoxd11-13+/+ , Hoxa13+/- ; Hoxd11-13+/- , Hoxa13+/+ ; Hoxd11-13-/- , Hoxa13+/- ; Hoxd11-13-/- , and Hoxa13-/- ; Hoxd11-13-/- . Each of these types was combined with either the normal GLi3 dose, Gli3+/+ , or the heterozygous dose Gli3+/XtJ , or the null dose Gli3XtJ/XtJ . With progressive removal of Hox and Gli3, phenotypes show more and more digits, from the control number of 5 to 13 for the Hoxd11-13+/- ; Hoxa13-/- ;Gli3XtJ/XtJ mutant. Here the number of digits depends on the Hox dose, with finely graduated steps (see Figure 8). 23 ฀฀฀฀฀฀฀฀฀ Sheth et al. created a simple linear reaction-diffusion model. The parameters of the reaction kinetics involving a generic activator and a generic inhibitor were kept constant except for the activator-dependent production rate of the inhibitor, which was assumed to be under the joint control of FGFs and Hox genes. The model showed that indeed, a reduction of Hox dose led to a decrease of the wavelength of the Turing pattern; this wavelength could be tuned through control of the Hox dose. These experimental results were then incorporated into a much more detailed model by Raspopovic et al. (2014) [132], the so called BSW model. This model takes into account Sox9, the earliest skeletal marker in the mouse, bone morphogenetic proteins (BMPs) and WNT. It was found that all three show spatially periodic expression patterns. Sox9 was exactly out of phase with BMP and WNT, i.e. the peaks of concentration of Sox9 coincided with the concentration troughs of BMP and WNT, and vice versa. Two of the regulatory interactions between the three components are known: WNT signaling inhibits Sox9 and BMP upregulates Sox9. The other relationships in a Turing reaction-diffusion network were chosen in such a way that the linearized solutions show the same phase pattern as the experiments, with BMP and WNT being in-phase and Sox9 being exactly out-of-phase relative to them. In the linear reaction kinetics, a third order term was added to prevent blow up of concentrations. This yields the so-called BSW model: ∂s = αs + k2 b − k3 w − (s − s′ )3 ∂t ∂b = α b − k 4 s − k 5 b + d b ∇2 b ∂t ∂w = α w − k 7 s − k 9 w + d w ∇2 w ∂t Here s(x, t), b(x, t) and w(x, t) are the concentrations of Sox9,BMP, and WNT, respectively. The parameters αs , αw , αb , k2 , k3 , k4 , k5 , k7 , k9 are positive constants, and db , dw are the diffusion coefficients of BMP and WNT, respectively. Note that Sox9 does not diffuse. The condition that BMP and WNT be in-phase and Sox9 out-of-phase means that the system is qualitatively different from activator-inhibitor systems, where the activator and inhibitor are necessarily in-phase relative to each other, i.e. the peaks of the activator concentration coincide with the peaks of the inhibitor concentration. (This is a simple consequence of the sign pattern of the linearization matrix, see [206].) Out-of-phase concentrations are a feature of so-called substrate-depletion systems, which are two component Turing reaction-diffusion systems with the opposite sign patterns on the off-diagonal of the linearization matrix relative to the more familiar activator-inhibitor system. This is because the self-activating molecule inhibits the other one, and the self-inhibiting one activates the other. While the three-component BSW system does not fall directly into this classification of two-component systems, it does share similarities with the substrate-depletion model. Raspopovic et al. solved the BSW equations on realistic, growing two-dimensional limb shapes. To get the correct proximal-distal orientation of the stripes of the Turing patterns, they had to include FGF and Hoxd13 gradients, which they postulated to modulate the inhibitory effect of Sox9 on BMP and WNT. With these additions, the system could robustly recapitulate the main features of experimental Sox9 expression. 24 ฀฀฀฀฀฀฀฀฀ Figure 8: In experiments reported by Sheth et al. [149] for a set of different mouse mutants, the number of digital elements finely depends on the Hox dose. The images show Sox9 expression, a marker for precartilage condensations, in the mouse at embryonic days E12.5 (or E13.5 if indicated). (Modified from [149].) 5 Discussion and Outlook Recent years have seen the proposal of several new mathematical models of pattern formation in the vertebrate limb, with the the Iber-Badugu model [6], the Glimm-Bhat-Newman galectin model [53] and the BSW model [132]. These models rely on a much more in-depth understanding of relevant gene regulatory networks than previous models. One of the limitations of current modeling approaches is the relative lack of data for species other than the mouse and the chicken. Thus it is presently not clear to which extent the proposed models transcend the specific details of the model organism to reflect generic mechanisms that apply to all tetrapods. For instance, Sox9 expression patterns are quite distinct for the mouse, the chicken and turtles, respectively [103], and Sox9-null limb mesenchyne still exhibits precartilage condensations [8]. In turn, galectin-1 null mice have been shown to have normally developed limb skeletons [49], so that questions about the generality of both the galectin model and the BSW model remain. All mathematical models presented in this survey have taken into account only a small number of components, in contrast with the hundreds of molecules that have been shown to play a role in limb chondrogenesis. The process itself is characterized by a great robustness and redundancy of many components. Newman et al. argue in [114] that these models are not to be understood as necessarily competing explanations, but rather represent 25 ฀฀฀฀฀฀฀฀฀ different modules of a multisystem complex which each are capable of generating patterns in a self-organizing way, but whose interplay yields the redundancy which is the source of the extraordinary robustness of the overall patterning process. For instance, they speculate that the BSW network evolved from a differentiation-inducing module that served as a ‘readout’ of an existing prepattern to a self-organizing patterning system in its own right. Arguably, despite the progress in recent years with new mathematical models, the conceptual understanding of the mechanisms of chondrogenic pattern formation in the limb is lagging far behind experimental investigations, which have generated huge amounts of data through increasingly sophisticated visualization and experimentation techniques. Besides expanding the investigation of component mechanisms to more species than the mouse and the chicken, the analysis of how these different components may interact, reinforce each other and yield robustness is a crucial future task. Addressing these problems will certainly encompass sophisticated integrated computational multi-scale models carefully vetted against data. However, the mathematical investigation of simpler, more analytically tractable models will remain relevant, for instance in addressing the question of how two independent self-organizing systems acting in concert may enhance the robustness of the overall patterning process. Acknowledgements: P.C. and B.K. were supported by the National Science Centre (Poland) grant 2016/21/B/ST1/03071. The authors would like to thank an anonymous reviewer for comments that have substantially improved the quality of the paper. Declaration of interest All authors declare that they have no competing interest. References [1] Agarwal, P., Wylie, J.N., Galceran, J., Arkhitko, O., Li, C., Deng, C., Grosschedl, R., Bruneau, B.G., Tbx5 is essential for forelimb bud initiation following patterning of the limb field in the mouse embryo. Development, 130(3), (2003), 623-633. [2] Akiyama, H., Chaboissier, M.-C., Martin, J.F., Schedl, A., de Crombrugghe, B., The transcription factor Sox9 has essential roles in successive steps of the chondrocyte differentiation pathway and is required for expression of Sox5 and Sox6. Genes Dev 16,(2002), 28132828. [3] Alber, M., Hentschel, H.G.E., Kazmierczak, B., Newman, S. A., Existence of solutions to a new model of biological pattern formation. Journal of mathematical analysis and applications, 308(1), (2004), 175-194. [4] Alber, M., Glimm, T., Hentschel, H.G.E., Kazmierczak, B., Newman, S.A., Stability of n-dimensional patterns in a generalized Turing system: implications for biological pattern formation. Nonlinearity, 18(1), (2004), 125. [5] Armstrong, N.J., Painter, K.J., Sherratt, J.A., A continuum approach to modelling cell-cell adhesion. Journal of Theoretical Biology, 243(1), (2006), 98-113. [6] Badugu, A., Kraemer, C., Germann, P., Menshykau, D., Iber, D., Digit patterning during limb development as a result of the BMP-receptor interaction. Scientific reports, 2, (2012), 991. [7] Alber, M., Glimm, T., Hentschel, H.G.E., Kazmierczak, B., Zhang, Y.-T., Zhu, J., and Newman, S.A. (2008). The morphostatic limit for a model of skeletal pattern formation in the vertebrate limb. Bull. Math. Biol. 70, 460483. [8] Barna, M., Niswander, L., Visualization of Cartilage Formation: Insight into Cellular Properties of Skeletal Progenitors and Chondrodysplasia Syndromes. Developmental Cell, 12(6), (2007), 931-941. 26 ฀฀฀฀฀฀฀฀฀ [9] Bastida, M. F., Sheth, R., Ros, M. A. , A BMP-Shh negative-feedback loop restricts Shh expression during limb development. Development, 136(22), (2009), 3779-3789. [10] Beenken, A., Mohammadi, M. The FGF family: biology, pathophysiology and therapy. Nature reviews Drug discovery, 8(3), (2009), 235-253. [11] Benson, D. L., Sherratt, J. A., Maini, P. K., Diffusion driven instability in an inhomogeneous domain. Bulletin of mathematical biology, 55(2), (1992), 365-384. [12] Bénazet, J.D., Bischofberger, M., Tiecke, E., Gonalves, A., Martin, J.F., Zuniga, A., Naef, F., Zeller, R., A self-regulatory system of interlinked signaling feedback loops controls mouse limb patterning. Science, 323(5917), (2009), 1050-1053. [13] Bénazet, J.D., Zeller, R., Vertebrate limb development: moving from classical morphogen gradients to an integrated 4-dimensional patterning system. Cold Spring Harbor perspectives in biology, 1(4), (2009), a001339. [14] Bénazet, J.D., Pignatti, E., Nugent, A., Unal, E., Laurent, F., Zeller, R., Smad4 is required to induce digit ray primordia and to initiate the aggregation and differentiation of chondrogenic progenitors in mouse limb buds. Development, 139(22), (2012), 4250-4260. [15] Bhat, R., Lerea, K.M., Peng, H., Kaltner, H., Gabius, H.J., Newman, S.A., A regulatory network of two galectins mediates the earliest steps of avian limb skeletal morphogenesis. BMC developmental biology, 11(1), (2011), 1. [16] Boehm, B., Westerberg, H., Lesnicar-Pucko, G., Raja, S., Rautschka, M., Cotterell, J., Swoger, J. and Sharpe, J., The role of spatially controlled cell proliferation in limb bud morphogenesis. PLoS biology, 8(7), (2010), 1512. [17] Borckmans, P., Dewel, G., De Wit, A., Walgraef, D., Turing bifurcations and pattern selection. In Chemical waves and patterns, Springer Netherlands. (1995), 323-363. [18] Borkhvardt, V. G., The growth and form development of the limb buds in vertebrate animals. Ontogenez, 31(3),(1999): 192-200. [19] Boulet, A.M., Moon, A.M., Arenkiel, B.R. and Capecchi, M.R., The roles of Fgf4 and Fgf8 in limb bud initiation and outgrowth. Developmental biology, 273(2), (2004), 361-372. [20] R. Chaturvedi, C. Huang, B. Kazmierczak, T. Schneider, J.A. Izaguirre, T. Glimm, H.G.E. Hentschel, J.A. Glazier, S.A. Newman, M.S. Alber, On multiscale approaches to three-dimensional modelling of morphogenesis. J. R. Soc. Interface 2, (2005), 237. [21] Chen, H., Johnson, R.L., Dorsoventral patterning of the vertebrate limb: a process governed by multiple events. Cell and tissue research, 296(1), (1999), 67-73. [22] Chen, Y., Knezevic, V., Ervin, V., Hutson, R., Ward, Y., Mackem, S., Direct interaction with Hoxd proteins reverses Gli3-repressor function to promote digit formation downstream of Shh. Development, 131(10), (2004), 2339-2347. [23] Church, V. L., Francis-West, P., Wnt signalling during limb development. International Journal of Developmental Biology, 46(7), (2002), 927-936. [24] Cickovski, T., Aras, K., Alber, M.S., Izaguirre, J.A., Swat, M., Glazier, J.A., Merks, R.M.H., Glimm, T., Hentschel, H.G.E., and Newman, S.A. (2007). From Genes to Organisms Via the Cell A Problem-Solving Environment for Multicellular Development. Comput Sci Eng 9, 5060. 27 ฀฀฀฀฀฀฀฀฀ [25] Clase, K.L., Mitchell, P.J., Ward, P.J., Dorman, C.M., Johnson, S.E., Hannon, K., FGF5 stimulates expansion of connective tissue fibroblasts and inhibits skeletal muscle development in the limb. Developmental Dynamics, 219(3), (2000), 368-380. [26] Cooper, K.L. (2015). Self-organization in the limb: a Turing mechanism for digit development. Current Opinion in Genetics & Development 32, 9297. [27] Crossley, P.H., Minowada, G., MacArthur, C.A., Martin, G.R., Roles for FGF8 in the induction, initiation, and maintenance of chick limb development. Cell, 84(1), (1996), 127-136. [28] De Lise, A.M., Fischer, L., Tuan, R.S., Cellular interactions and signaling in cartilage development. Osteoarthritis and cartilage, 8(5), (2000), 309-334. [29] N. Denef, D. Neubuser, L. Perez, S.M. Cohen, Hedgehog induces opposite changes in turnover and subcellular localization of patched and smoothened, Cell, 102, (2000), 521. [30] Dillon, R., Maini, P. K., Othmer, H. G., Pattern formation in generalized Turing systems. Journal of Mathematical Biology, 32(4), (1994), 345-393. [31] Dillon, R. H., A mathematical model of vertebrate limb development with modulated reaction and diffusion. University of Utah, (1993). [32] Dillon, R., Othmer, H. G., A mathematical model for outgrowth and spatial patterning of the vertebrate limb bud. Journal of theoretical biology, 197(3), (1999), 295-330. [33] Dillon, R. H., ”Mathematical modeling of vertebrate limb development”. Mathematical Models for Biological Pattern Formation. Springer New York, (2001), 39-57. [34] Dillon, R., Gadgil, C., Othmer, H. G., Short-and long-range effects of Sonic hedgehog in limb development. Proceedings of the National Academy of Sciences, 100(18), (2003), 10152-10157. [35] Drossopoulou, G., Lewis, K. E., Sanz-Ezquerro, J. J., Nikbakht, N., McMahon, A. P., Hofmann, C., Tickle, C., A model for anteroposterior patterning of the vertebrate limb based on sequential long-and short-range Shh signalling and Bmp signalling. Development, 127(7), (2000), 1337-1348. [36] Dudley, A.T., Ros, M.A., Tabin, C.J., A re-examination of proximodistal patterning during vertebrate limb development. Nature, 418(6897), (2002), 539-544. [37] Duprez, D.M., Kostakopoulou, K., Francis-West, P.H., Tickle, C., Brickell, P.M., Activation of Fgf-4 and HoxD gene expression by BMP-2 expressing cells in the developing chick limb. Development, 122(6), (1996), 1821-1828. [38] Dyson, J., Gourley, S.A., Villella-Bressan, R. and Webb, G.F., Existence and asymptotic properties of solutions of a nonlocal evolution equation modeling cell-cell adhesion. SIAM Journal on Mathematical Analysis, 42(4), pp.1784-1804, 2010. [39] Dyson, J., Gourley, S.A. and Webb, G.F., A non-local evolution equation model of cellcell adhesion in higher dimensional space. Journal of biological dynamics, 7(sup1), pp.68-87, 2013. [40] Economou, A.D., Green, J.B., Modelling from the experimental developmental biologists viewpoint. In Seminars in cell & developmental biology, Academic Press, 35, (2014), 58-65. [41] Ede, D.A., Law, J.T., Computer simulation of vertebrate limb morphogenesis. Nature, 221, (1969), 244-248. 28 ฀฀฀฀฀฀฀฀฀ [42] Evans, J.A., Vitez, M. and Czeizel, A., Congenital abnormalities associated with limb deficiency defects: a population study based on cases from the Hungarian Congenital Malformation Registry (19751984). American journal of medical genetics, 49(1), (1994), 52-66. [43] Fallon, J.F., Lopez, A., Ros, M.A., Savage, M.P., Olwin, B.B. and Simandl, B.K., FGF-2: apical ectodermal ridge growth signal for chick limb development. Science, 264(5155), (1994), 104-107. [44] Farshi, P., Ohlig, S., Pickhinke, U., Höing, S., Jochmann, K., Lawrence, R., Dreier, R., Dierker, T. and Grobe, K., Dual roles of the Cardin-Weintraub motif in multimeric Sonic hedgehog. Journal of Biological Chemistry, 286(26), (2011), 23608-23619. [45] Frenz, D. A., Akiyama, S. K., Paulsen, D. F.,Newman, S. A., Latex beads as probes of cell surfaceextracellular matrix interactions during chondrogenesis: evidence for a role for amino-terminal heparinbinding domain of fibronectin. Developmental biology, 136(1), (1989), 87-96. [46] Frenz, D. A., Jaikaria, N. S., Newman, S. A., The mechanism of precartilage mesenchymal condensation: a major role for interaction of the cell surface with the amino-terminal heparin-binding domain of fibronectin. Developmental biology, 136(1), (1989), 97-103. [47] Ganan, Y., Macias, D., Duterque-Coquillaud, M., Ros, M. A., Hurle, J. M., Role of TGF beta s and BMPs as signals controlling the position of the digits and the areas of interdigital cell death in the developing chick limb autopod. Development, 122(8), (1996), 2349-2357. [48] Geduspan, J.S., Maccabe, J.A., Transfer of dorsoventral information from mesoderm to ectoderm at the onset of limb development. The Anatomical Record, 224(1), (1989),79-87. [49] Georgiadis, V. , Stewart, H. J. , Pollard, H. J. , Tavsanoglu, Y., Prasad, R. , Horwood, J. , Deltour, L. , Goldring, K. , Poirier, F., Lawrence-Watt, D. J., Lack of galectin-1 results in defects in myoblast fusion and muscle regeneration, Dev. Dyn.,236(4), (2007), 1014–1024. [50] Gibson-Brown, J.J., Agulnik, S.I., Chapman, D.L., Alexiou, M., Garvey, N., Lee, S.M., Papaioannou, V.E., Evidence of a role for T- genes in the evolution of limb morphogenesis and the specification of forelimb/hindlimb identity. Mechanisms of development, 56(1), (1996), 93-101. [51] Gilbert, S. F., Developmental biology: The anatomical tradition. (2000). [52] Glimm, T., Headon, D., Kiskowski, M.A., Computational and mathematical models of chondrogenesis in vertebrate limbs. Birth Defects Research Part C: Embryo Today: Reviews, 96(2), (2012), 176-192. [53] Glimm, T., Bhat, R. Newman, S.A., Modeling the morphodynamic galectin patterning network of the developing avian limb skeleton. Journal of theoretical biology, 346, (2014), 86-108. [54] Green, J.B., Sharpe, J., Positional information and reaction-diffusion: two big ideas in developmental biology combine. Development, 142(7), (2015), 1203-1211. [55] Gros, J., Hu, J. K. H., Vinegoni, C., Feruglio, P. F., Weissleder, R., Tabin, C. J., WNT5A/JNK and FGF/MAPK pathways regulate the cellular events shaping the vertebrate limb bud. Current Biology, 20(22), (2010), 1993-2002. [56] Hajihosseini, M.K., Heath, J.K., Expression patterns of fibroblast growth factors-18 and-20 in mouse embryos is suggestive of novel roles in calvarial and limb development. Mechanisms of development, 113(1), (2002), 79-83. [57] Hall, B. K., Miyake, T. S., Divide, accumulate, differentiate: cell condensation in skeletal development revisited. International Journal of Developmental Biology, 39, (1995), 881-94. 29 ฀฀฀฀฀฀฀฀฀ [58] Hall, B. K., Miyake, T., All for one and one for all: condensations and the initiation of skeletal development. Bioessays, 22(2), (2000), 138-147. [59] Hamburger, V., Morphogenetic and axial selfdifferentiation of transplanted limb primordia of 2day chick embryos. Journal of Experimental Zoology Part A: Ecological Genetics and Physiology, 77(3), (1938), 379-399. [60] Hamburger, V., Hamilton, H.L., A series of normal stages in the development of the chick embryo. Journal of morphology, 88(1), (1951), 49-92. [61] Harfe, B.D., Scherz, P.J., Nissim, S., Tian, H., McMahon, A.P., Tabin, C.J., Evidence for an expansionbased temporal Shh gradient in specifying vertebrate digit identities. Cell, 118(4), (2004), 517-528. [62] Hartmann, C., Tabin, C. J., Dual roles of Wnt signaling during chondrogenesis in the chicken limb. Development, 127(14), (2000), 3141-3159. [63] Heine, U.L., Munoz, E.F., Flanders, K.C., Ellingsworth, L.R., Lam, H.Y., Thompson, N.L., Roberts, A.B., Sporn, M.B., Role of transforming growth factor-beta in the development of the mouse embryo. The Journal of Cell Biology, 105(6), (1987), 2861-2876. [64] H.G.E. Hentschel, T. Glimm, J.A. Glazier, S.A. Newman, Dynamical mechanisms for skeletal pattern formation in the vertebrate limb. Proc. R. Soc. B 271, (2004), 1713 -1722. [65] Herrero, M. A., López, J. M., Bone formation: Biological aspects and modelling problems. Journal of Theoretical Medicine, 6(1), (2005), 41-55. [66] Hirashima, T., Iwasa, Y., Morishita, Y., Distance between AER and ZPA is defined by feed-forward loop and is stabilized by their feedback loop in vertebrate limb bud. Bulletin of mathematical biology, 70(2), (2008), 438-459. [67] Hopyan, S., Sharpe, J., Yang, Y., Budding behaviors: Growth of the limb as a model of morphogenesis. Developmental Dynamics, 240(5), (2011), 1054-1062. [68] Hung, I.H., Yu, K., Lavine, K.J., Ornitz, D.M., FGF9 regulates early hypertrophic chondrocyte differentiation and skeletal vascularization in the developing stylopod. Developmental biology, 307(2), (2007), 300-313. [69] Jaruratanasirikul, S., Tangtrakulwanich, B., Rachatawiriyakul, P., Sriplung, H., Limpitikul, W., Dissaneevate, P., Khunnarakpong, N. and Tantichantakarun, P., Prevalence of congenital limb defects: Data from birth defects registries in three provinces in Southern Thailand. Congenital anomalies, 56(5), 203-208, 2016. [70] Javois, L.C., Pattern specification in the developing chick limb. Pattern formation, (1984), 557-579. [71] Kawakami, Y., Tsuda, M., Takahashi, S., Taniguchi, N., Esteban, C. R., Zemmyo, M., ... & Asahara, H. Transcriptional coactivator PGC-1α regulates chondrogenesis via association with Sox9. Proceedings of the National Academy of Sciences of the United States of America, 102(7), (2005), 2414-2419. [72] Keller, E.F. and Segel, L.A., Model for chemotaxis. Journal of theoretical biology, 30(2), (1971), 225-234. [73] Kicheva, A., Pantazis, P., Bollenbach, T., Kalaidzidis, Y., Bittig, T., Jülicher, F., González-Gaitán, M., Kinetics of morphogen gradient formation. Science, 315(5811), (2007), 521-525. [74] Kimura, Y., Matsunami, H., Inoue, T., Shimamura, K., Uchida, N., Ueno, T., Miyazaki, T., Takeichi, M., Cadherin-11 expressed in association with mesenchymal morphogenesis in the head, somite, and limb bud of early mouse embryos. Developmental biology, 169(1), (1995), 347-358. 30 ฀฀฀฀฀฀฀฀฀ [75] Kook, S.H., Jeon, Y.M., Lim, S.S., Jang, M.J., Cho, E.S., Lee, S.Y., Choi, K.C., Kim, J.G. and Lee, J.C., Fibroblast growth factor-4 enhances proliferation of mouse embryonic stem cells via activation of c-Jun signaling. PloS one, 8(8), 2013, e71641. [76] Kondo, S., Miura, T., Reaction-diffusion model as a framework for understanding biological pattern formation. Science, 329(5999), (2010), 1616-1620. [77] Leonard, C. M., Fuld, H. M., Frenz, D. A., Downie, S. A., Massague, J., Newman, S. A. , Role of transforming growth factor-β in chondrogenic pattern formation in the embryonic limb: stimulation of mesenchymal condensation and fibronectin gene expression by exogenenous TGF-β and evidence for endogenous TGFβ-like activity. Developmental biology, 145(1), (1991), 99-109. [78] Lewandoski, M., Sun, X., Martin, G.R., Fgf8 signalling from the AER is essential for normal limb development. Nature genetics, 26(4), (2000), 460-463. [79] Li, S., Muneoka, K., Cell migration and chick limb development: chemotactic action of FGF-4 and the AER. Developmental biology, 211(2), (1999), 335-347. [80] Lillie, F. R., Experimental studies on the development of the organs in the embryo of the fowl (Gallus domesticus). The Biological Bulletin, 5(2), (1903), 92-124. [81] Litingtung, Y., Dahn, R. D., Li, Y., Fallon, J. F., & Chiang, C. Shh and Gli3 are dispensable for limb skeleton formation but regulate digit number and identity. Nature, 418(6901), (2002), 979-983. [82] Liu, Z., Xu, J., Colvin, J.S., Ornitz, D.M., Coordination of chondrogenesis and osteogenesis by fibroblast growth factor 18. Genes & development, 16(7),(2002), 859-869. [83] Kurics, T., Menshykau, D., Iber, D., Feedback, receptor clustering, and receptor restriction to single cells yield large Turing spaces for ligand-receptor-based Turing models. Phys. Rev. E 90, (2014) 022716. https://doi.org/10.1103/PhysRevE.90.022716 [84] Lorda-Diez, C.I., Montero, J.A., Diaz-Mendoza, M.J., Garcia-Porrero, J.A., Hurle, J.M., Defining the earliest transcriptional steps of chondrogenic progenitor specification during the formation of the digits in the embryonic limb. PLoS One, 6(9), (2011), e24546. [85] Luo, Y., Kostetskii, I. and Radice, G.L., N-cadherin is not essential for limb mesenchymal chondrogenesis. Developmental dynamics, 232(2), (2005), 336-344. [86] Macias, D., Ganan, Y., Sampath, T. K., Piedra, M. E., Ros, M. A., Hurle, J. M., Role of BMP-2 and OP-1 (BMP-7) in programmed cell death and skeletogenesis during chick limb development. Development, 124(6), (1997), 1109-1117. [87] Maden, M., Intercalary regeneration in the amphibian limb and the rule of distal transformation. Development, 56(1), (1980), 201-209. [88] Maden, M., Mustafa, K., The structure of 180 supernumerary limbs and a hypothesis of their formation. Developmental biology, 93(1), (1982), 257-265. [89] Marcon, L., Sharpe, J., Turing patterns in development: what about the horse part?. Current opinion in genetics & development, 22(6), (2012), 578-584. [90] Mariani, F.V., Ahn, C.P., Martin, G.R., Genetic evidence that FGFs have an instructive role in limb proximaldistal patterning. Nature, 453(7193), (2008), 401-405. [91] Martin, P. Lewis, J., Normal development of the skeleton in chick limb buds devoid of dorsal ectoderm. Developmental biology, 118(1), (1986), 233-246. 31 ฀฀฀฀฀฀฀฀฀ [92] Martin, G. R., The roles of FGFs in the early development of vertebrate limbs. Genes & Development, 12(11), (1998), 1571-1586. [93] Martin, G., Making a vertebrate limb: new players enter from the wings. Bioessays, 23(10), (2001), 865-868. [94] Meinhardt, H., Models for the ontogenetic development of higher organisms. In Reviews of Physiology, Biochemistry and Pharmacology, Volume 80, (1978), 47-104, Springer Berlin Heidelberg. [95] Meinhardt, H., Models of biological pattern formation. Vol. 6. London: Academic Press, (1982). [96] Meinhardt, H., A boundary model for pattern formation in vertebrate limbs. Journal of embryology and experimental morphology, 76(1), (1983), 115-137. [97] Meinhardt, H., Gierer, A., Pattern formation by local self-activation and lateral inhibition. BioEssays, 22(8), (2000), 753-760. [98] Meinhardt, H., Models of biological pattern formation: from elementary steps to the organization of embryonic axes. Current topics in developmental biology, 81, (2008), 1-63. [99] Mic, F. A., Sirbu, I. O., Duester, G., Retinoic acid synthesis controlled by Raldh2 is required early for limb bud initiation and then later as a proximodistal signal during apical ectodermal ridge formation. Journal of Biological Chemistry, 279(25), (2004), 26698-26706. [100] Michon, F., Forest, L., Collomb, E., Demongeot, J., Dhouailly, D., BMP2 and BMP7 play antagonistic roles in feather induction. Development, 135(16), (2008), 2797-2805. [101] Min, H., Danilenko, D.M., Scully, S.A., Bolon, B., Ring, B.D., Tarpley, J.E., DeRose, M., Simonet, W.S., Fgf-10 is required for both limb and lung development and exhibits striking functional similarity to Drosophila branchless. Genes & development, 12(20), (1998), 3156-3161. [102] Minguillon, C., Del Buono, J., Logan, M.P., Tbx5 and Tbx4 are not sufficient to determine limb-specific morphologies but have common roles in initiating limb outgrowth. Developmental cell, 8(1), (2005), 75-84. [103] Montero, J. A., Lorda-Diez, C. I. , Francisco-Morcillo, J., Chimal-Monroy, J. , Garcia-Porrero, J. A. , Hurle, J. M. , Sox9 Expression in Amniotes: Species-Specific Differences in the Formation of Digits. Front Cell Dev Biol,5 (2017), 23 pages. [104] Morishita, Y., Iwasa, Y., Growth based morphogenesis of vertebrate limb bud. Bulletin of mathematical biology, 70(7), (2008), 1957-1978. [105] Mou, C., Jackson, B., Schneider, P., Overbeek, P.A., Headon, D.J., Generation of the primary hair follicle pattern. Proceedings of the National Academy of Sciences, 103(24), (2006), 9075-9080. [106] Müller, P., Rogers, K.W., Jordan, B.M., Lee, J.S., Robson, D., Ramanathan, S., Schier, A.F., Differential diffusivity of Nodal and Lefty underlies a reaction-diffusion patterning system. Science, 336(6082), (2012), 721-724. [107] Murea, C., Hentschel, H., A finite element method for growth in biological development. Mathematical Biosciences and Engineering, 4(2), (2007), 339. [108] Murray, J.D., Mathematical Biology II: Spatial Models and Biomedical Applications. (third ed.) Springer, New York (2003). [109] Naiche, L.A., Papaioannou, V.E., Loss of Tbx4 blocks hindlimb development and affects vascularization and fusion of the allantois. Development, 130(12), (2003), 2681-2693. 32 ฀฀฀฀฀฀฀฀฀ [110] Newman, S. A. & Frisch, H. L., Dynamics of skeletal pattern formation in developing chick limb. Science, 205(4407), (1979), 662-668. [111] Newman, S. A., Müller, G. B., Origination and innovation in the vertebrate limb skeleton: an epigenetic perspective. Journal of Experimental Zoology Part B: Molecular and Developmental Evolution, 304(6), (2005), 593-609. [112] Newman, S. A., Tomasek, J. J., Morphogenesis of Connective Tissues. Extracellular matrix, 2, (1996), 335. [113] Newman, S.A., Christley, S., Glimm, T., Hentschel, H.G.E., Kazmierczak, B., Zhang, Y.T., Zhu, J., Alber, M., Multiscale models for vertebrate limb development. Current topics in developmental biology, 81, (2008), 311-340. [114] Newman, S. A. , Glimm, T., Bhat, R., The vertebrate limb: An evolving complex of self-organizing systems, Prog. Biophys. Mol. Biol., 137 (2018) 12–24. [115] Ng, J.K., Kawakami, Y., Büscher, D., Raya, Á., Itoh, T., Koth, C.M., Esteban, C.R., Rodrguez-Len, J., Garrity, D.M., Fishman, M.C., Belmonte, J.C.I., The limb identity gene Tbx5 promotes limb initiation by interacting with Wnt2b and Fgf10. Development, 129(22), (2002), 5161-5170. [116] Nishimoto, S., Wilde, S.M., Wood, S., Logan, M.P., RA acts in a coherent feed-forward mechanism with Tbx5 to control limb bud induction and initiation. Cell reports, 12(5), (2015), 879-891. [117] Niswander, L., Tickle, C., Vogel, A., Booth, I., Martin, G.R., FGF-4 replaces the apical ectodermal ridge and directs outgrowth and patterning of the limb. Cell, 75(3), (1993), 579-587. [118] Niswander, L., Pattern formation: old models out on a limb. Nature Reviews Genetics, 4(2), (2003), 133-143. [119] Oberlender, S. A., Tuan, R. S., Expression and functional involvement of N-cadherin in embryonic limb chondrogenesis. Development, 120(1), (1994), 177-187. [120] Okazaki, K., Sandell, L. J. Extracellular matrix gene regulation. Clinical orthopaedics and related research, 427, (2004), S123-S128. [121] Ohbayashi, N., Shibayama, M., Kurotaki, Y., Imanishi, M., Fujimori, T., Itoh, N., Takada, S., FGF18 is required for normal cell proliferation and differentiation during osteogenesis and chondrogenesis. Genes & Development, 16(7), (2002), 870879. [122] Ohuchi, H., Nakagawa, T., Yamamoto, A., Araga, A., Ohata, T., Ishimaru, Y., Yoshioka, H., Kuwana, T., Nohno, T., Yamasaki, M., Itoh, N., The mesenchymal factor, FGF10, initiates and maintains the outgrowth of the chick limb bud through interaction with FGF8, an apical ectodermal factor. Development, 124(11), (1997), 2235-2244. [123] Ohuchi, H., Kimura, S., Watamoto, M., Itoh, N., Involvement of fibroblast growth factor (FGF) 18-FGF8 signaling in specification of left-right asymmetry and brain and limb development of the chick embryo. Mechanisms of development, 95(1), (2000), 55-66. [124] Oster, G.F., Lateral inhibition models of developmental processes. Mathematical Biosciences, 90(1-2), (1988), 265-286. [125] Pajni-Underwood, S., Wilson, C.P., Elder, C., Mishina, Y., Lewandoski, M., BMP signals control limb bud interdigital programmed cell death by regulating FGF signaling. Development, 134(12), (2007), 2359-2368. 33 ฀฀฀฀฀฀฀฀฀ [126] Pan, Y., Liu, Z., Shen, J., Kopan, R., Notch1 and 2 cooperate in limb ectoderm to receive an early Jagged2 signal regulating interdigital apoptosis. Developmental biology, 286(2), (2005), 472-482. [127] Pearse, R.V., Scherz, P.J., Campbell, J.K., Tabin, C.J., A cellular lineage analysis of the chick limb bud. Developmental biology, 310(2), (2007), 388-400. [128] Peebles, F., On the interchange of the limbs of the chick by transplantation. The Biological Bulletin, 20(1), (1910), 14-18. [129] Petit, F., Sears, K.E., Ahituv, N., Limb development: a paradigm of gene regulation. Nature Reviews Genetics. 2017. [130] Poplawski, N. J., Swat, M., Gens, J. S., Glazier, J. A., Adhesion between cells, diffusion of growth factors, and elasticity of the AER produce the paddle shape of the chick limb. Physica A: Statistical Mechanics and its Applications, 373, (2007), 521-532. [131] Probst, S., Kraemer, C., Demougin, P., Sheth, R., Martin, G.R., Shiratori, H., Hamada, H., Iber, D., Zeller, R., Zuniga, A., SHH propagates distal limb bud development by enhancing CYP26B1-mediated retinoic acid clearance via AER-FGF signalling. Development, 138(10), (2011), 1913-1923. [132] Raspopovic, J., Marcon, L., Russo, L., Sharpe, J., Digit patterning is controlled by a Bmp-Sox9-Wnt Turing network modulated by morphogen gradients. Science, 345(6196), (2014), 566-570. [133] Ratcliffe, A., Mow, V. C., Articular Cartilage. Extracellular matrix, 1, (1996), 234. [134] Riddle, R.D., Johnson, R.L., Laufer, E., Tabin, C., Sonic hedgehog mediates the polarizing activity of the ZPA. Cell, 75(7), (1993), 1401-1416. [135] Riley, B.B., Savage, M.P., Simandl, B.K., Olwin, B.B. and Fallon, J.F., Retroviral expression of FGF-2 (bFGF) affects patterning in chick limb bud. Development, 118(1), (1993), 95-104. [136] Rogers, K.W., and Schier, A.F. (2011). Morphogen Gradients: From Generation to Interpretation. Annual Review of Cell and Developmental Biology 27, 377407. [137] Rowe, D.A., Fallon, J.F., The effect of removing posterior apical ectodermal ridge of the chick wing and leg on pattern formation. Development, 65(Supplement), (1981), 309-325. [138] Rubin, L.,Saunders, J. W. J., Ectodermal-mesodermal interactions in the growth of limb buds in the chick embryo: constancy and temporal limits of the ectodermal induction. Developmental Biology 28, (1972),94112. [139] Saiz-Lopez, P., Chinnaiya, K., Campa, V.M., Delgado, I., Ros, M.A., Towers, M., An intrinsic timer specifies distal structures of the vertebrate limb. Nature communications, 6, (2015). [140] Saiz-Lopez, P., Chinnaiya, K., Towers, M., Ros, M.A., Intrinsic properties of limb bud cells can be differentially reset. Development, 144(3), (2017), 479-486. [141] Sato, K., Koizumi, Y., Takahashi, M., Kuroiwa, A., Tamura, K., Specification of cell fate along the proximal-distal axis in the developing chick limb bud. Development, 134(7), (2007), 1397-1406. [142] Saunders, J. W., The proximodistal sequence of origin of the parts of the chick wing and the role of the ectoderm. Journal of Experimental Zoology Part A: Ecological Genetics and Physiology, 108(3), 363-403. [143] Saunders, J. W., Cairns, J. M., Gasseling, M. T., The role of the apical ridge of ectoderm in the differentiation of the morphological structure and inductive specificity of limb parts in the chick. Journal of morphology, 101(1), (1957), 57-87. 34 ฀฀฀฀฀฀฀฀฀ [144] Saunders, J.W., Is the Progress Zone Model a Victim of Progress? https://doi.org/10.1016/S0092-8674(02)00936-4 Cell 110, 541543. (2002) [145] Scherz, P.J., Harfe, B.D., McMahon, A.P., Tabin, C.J., The limb bud Shh-Fgf feedback loop is terminated by expansion of former ZPA cells. Science, 305(5682), (2004), 396-399. [146] Scherz, P. J., McGlinn, E., Nissim, S., Tabin, C. J., Extended exposure to Sonic hedgehog is required for patterning the posterior digits of the vertebrate limb. Developmental biology, 308(2), (2007), 343-354. [147] Schier, A.F., Nodal morphogens. Cold Spring Harbor perspectives in biology, 1(5), (2009), a003459. [148] Sekine, K., Ohuchi, H., Fujiwara, M., Yamasaki, M., Yoshizawa, T., Sato, T., Yagishita, N., Matsui, D., Koga, Y., Itoh, N., Kato, S., Fgf10 is essential for limb and lung formation. Nature genetics, 21(1), (1999), 138-141. [149] Sheth, R., Marcon, L., Bastida, M.F., Junco, M., Quintana, L., Dahn, R., Kmita, M., Sharpe, J. , Ros, M.A., Hox genes regulate digit patterning by controlling the wavelength of a Turing-type mechanism. Science, 338(6113), (2012), 1476-1480. [150] Shiratori, H., Hamada, H., The left-right axis in the mouse: from origin to morphology. Development, 133(11), (2006), 2095-2104. [151] Sick, S., Reinker, S., Timmer, J., Schlake, T., WNT and DKK determine hair follicle spacing through a reaction-diffusion mechanism. Science, 314(5804), (2006), 1447-1450. [152] St-Jacques, B., Hammerschmidt, M., McMahon, A.P., 1Indian hedgehog signaling regulates proliferation and differentiation of chondrocytes and is essential for bone formation. Genes & development, 13(16), (1999), 2072-2086. [153] Stricker, S., & Mundlos, S. Mechanisms of digit formation: human malformation syndromes tell the story. Developmental Dynamics, 240(5), (2011), 990-1004. [154] Stewart, T.A., Liang, C., Cotney, J.L., Noonan, J.P., Sanger, T.J., and Wagner, G.P. (2019). Evidence against tetrapod-wide digit identities and for a limited frame shift in bird wings. Nat Commun 10, 3244. [155] Summerbell, D., Lewis, J. H., Wolpert, L., Positional information in chick limb morphogenesis. Nature, 244, (1973), 492-496. [156] Summerbell, D., A quantitative analysis of the effect of excision of the AER from the chick limb-bud. Journal of embryology and experimental morphology. 32(3), (1974),651-660. [157] Summerbell, D., The zone of polarizing activity: evidence for a role in normal chick limb morphogenesis. Journal of embryology and experimental morphology. 50, (1979), 217-233. [158] Summerbell, D., The effect of local application of retinoic acid to the anterior margin of the developing chick limb. Journal of embryology and experimental morphology, 78(1), (1983), 269-289. [159] Summerbell, D., Harvey, F., Vitamin A and the control of pattern in developing limbs. Progress in clinical and biological research, 110, (1983), 109. [160] Sun, X., Mariani, F. V., Martin, G. R., Functions of FGF signalling from the apical ectodermal ridge in limb development. Nature, 418(6897), (2002), 501-508. [161] Sun, M., Ma, F., Zeng, X., Liu, Q., Zhao, X.L., Wu, F.X., Wu, G.P., Zhang, Z.F., Gu, B., Zhao, Y.F., Tian, S.H., Triphalangeal thumbpolysyndactyly syndrome and syndactyly type IV are caused by genomic duplications involving the long range, limb-specific SHH enhancer. Journal of medical genetics, 45(9), (2008), 589-595. 35 ฀฀฀฀฀฀฀฀฀ [162] Szeto, D.P., Ryan, A.K., O’Connell, S.M., Rosenfeld, M.G., P-OTX: a PIT-1-interacting homeodomain factor expressed during anterior pituitary gland development. Proceedings of the National Academy of Sciences, 93(15),(1996), 7706-7710. [163] Tabin, C.,Wolpert, L., Rethinking the proximodistal axis of the vertebrate limb in the molecular era. Genes development, 21(12), (2007),1433-1442. [164] Thaller, C., Eichele, G., Retinoid Signaling in Vertebrate Limb Development. Annals of the New York Academy of Sciences, 785(1), (1996), 1-11. [165] Tickle, C., Summerbell, D., Wolpert, L., Positional signalling and specification of digits in chick limb morphogenesis. Nature, 254(5497), (1975), 199-202. [166] Tickle, C., Alberts, B., Wolpert, L., Lee, J., Local application of retinoic acid to the limb bond mimics the action of the polarizing region. Nature, 296, (1982), 564-566. [167] Tickle, C., Wolpert, L., The progress zonealive or dead?. Nature cell biology, 4(9), (2002), E216-E217. [168] Tickle, C., How the embryo makes a limb: determination, polarity and identity. Journal of Anatomy, 227(4), (2015), 418-430. [169] Tickle, C., An historical perspective on the pioneering experiments of John Saunders. Developmental Biology. (2017). [170] Towers, M., Mahood, R., Yin, Y., Tickle, C., Integration of growth and specification in chick wing digitpatterning. Nature, 452(7189), (2008), 882-886. [171] Towers, M., Tickle, C., Growing models of vertebrate limb development. Development, 136(2), (2009), 179-190. [172] Towers, M., Tickle, C., Generation of pattern and form in the developing limb. International Journal of Developmental Biology, 53(5), (2009), 805-812. [173] Oster, G.F., Murray, J.D., Maini, P.K., A model for chondrogenic condensations in the developing limb: the role of extracellular matrix and cell tractions. J Embryol Exp Morphol 89, (1985), 93112. [174] Oster, G.F., Murray, J.D., Harris, A.K., Mechanical aspects of mesenchymal morphogenesis . Development 78, (1983), 83125. [175] Turing, A.M., The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 237(641), (1952), 37-72. [176] Varjosalo, M., Taipale, J., Hedgehog: functions and mechanisms. Genes & development, 22(18), (2008), 2454-2472. [177] Wada, N., Spatiotemporal changes in cell adhesiveness during vertebrate limb morphogenesis. Developmental Dynamics, 240(5), (2011), 969-978. [178] Widelitz, R. B., Jiang, T. X., Murray, B. A., Chuong, C. M., Adhesion molecules in skeletogenesis: II. Neural cell adhesion molecules mediate precartilaginous mesenchymal condensations and enhance chondrogenesis. Journal of cellular physiology, 156(2), (1993), 399-411. [179] Willier, B. H., Rawles, M. E., The control of feather color pattern by melanophores grafted from one embryo to another of a different breed of fowl. Physiological Zoology, 13(2), (1940), 177-201. 36 ฀฀฀฀฀฀฀฀฀ [180] Wolpert, L., Positional information and the spatial pattern of cellular differentiation. Journal of theoretical biology, 25(1), (1969), 1-47. [181] Wolpert, L., Positional information and pattern formation. Current topics in developmental biology, 6, (1971), 183-224. [182] Wolpert, L., Tickle, C., Sampford, M., Lewis, J.H., . The effect of cell killing by X-irradiation on pattern formation in the chick limb. Development, 50(1), (1979), 175-198. [183] Wyngaarden, L.A., Vogeli, K.M., Ciruna, B.G., Wells, M., Hadjantonakis, A.K., Hopyan, S., Oriented cell motility and division underlie early limb bud morphogenesis. Development, 137(15), (2010), 2551-2558. [184] Yang, Y., Wnts and wing: Wnt signaling in vertebrate limb development and musculoskeletal morphogenesis. Birth Defects Research Part C: Embryo Today: Reviews, 69(4), (2003), 305-317. [185] Yonei-Tamura, S., Endo, T., Yajima, H., Ohuchi, H., Ide, H., Tamura, K., FGF7 and FGF10 directly induce the apical ectodermal ridge in chick embryos. Developmental biology, 211(1), (1999), 133-143. [186] Zeller, R., The temporal dynamics of vertebrate limb development, teratogenesis and evolution. Current opinion in genetics & development, 20(4), (2010), 384-390. [187] Zeng, W., Thomas, G.L., Glazier, J.A., Non-Turing stripes and spots: a novel mechanism for biological cell clustering. Physica A: Statistical Mechanics and its Applications, 341, (2004), 482-494. [188] Zenjari, C., Boilly-Marer, Y., Desbiens, X., Oudghir, M., Hondermarck, H., Boilly, B., Experimental evidence for FGF-1 control of blastema cell proliferation during limb regeneration of the amphibian Pleurodeles waltl. International Journal of Developmental Biology, 40(5), (1996), 965-971. [189] Zhang, Y. T., Alber, M. S., Newman, S. A. , Mathematical modeling of vertebrate limb development. Mathematical biosciences, 243(1), (2013), 1-17. [190] Zhu, J., Zhang, Y. T., Newman, S. A., Alber, M., Application of discontinuous Galerkin methods for reaction-diffusion systems in developmental biology. Journal of Scientific Computing, 40(1-3), (2009), 391418. [191] Zhu, J., Zhang, Y. T., Newman, S. A., Alber, M. S., A finite element model based on discontinuous Galerkin methods on moving grids for vertebrate limb pattern formation. Mathematical Modelling of Natural Phenomena, 4(04), (2009), 131-148. [192] Zhu, J., Zhang, Y. T., Alber, M. S., Newman, S. A., Bare bones pattern formation: a core regulatory network in varying geometries reproduces major features of vertebrate limb development and evolution. PLoS One, 5(5), (2010), e10892. [193] Zimmerman, L. B., De Jess-Escobar, J. M., & Harland, R. M. The Spemann organizer signal noggin binds and inactivates bone morphogenetic protein 4. Cell, 86(4), (1996), 599-606. [194] Bamshad, Michael, et al. Reconstructing the history of human limb development: lessons from birth defects. Pediatric research 45.3 (1999): 291. [195] Scott F. Gilbert, Developmental Biology, 9th Edition. [196] Newman, S., Bhat, R., Activator-inhibitor dynamics of vertebrate limb pattern formation. Birth Defects Res C Embryo Today. (2007) Dec;81(4):305-19 [197] Maini PK1, Solursh M. Cellular mechanisms of pattern formation in the developing limb. Int Rev Cytol. (1991);129:91-133. 37 ฀฀฀฀฀฀฀฀฀ [198] Glimm, T., Zhang, J., Shen, Y.-Q., Interaction of Turing patterns with an external linear morphogen gradient., Nonlinearity (2009) 22, 2541-2560. [199] Glimm, T., Zhang, J., Shen, Y.-Q., Newman, S.A., Reaction-diffusion systems and external morphogen gradients: the two-dimensional case, with an application to skeletal pattern formation, Bull. Math. Biol. (2012) 74, 666687. [200] Glimm, T., Zhang, J., Shen, Y.-Q., Stability of Turing-Type Patterns in a ReactionDiffusion System with an External Gradient, Int. J. Bifurcation Chaos (2017) 27, 1750003. [201] Bhat, R., Chakraborty, M., Glimm, T., Stewart, T.A., and Newman, S.A.Deep phylogenomics of a tandem-repeat galectin regulating appendicular skeletal pattern formation. BMC Evolutionary Biology 16, 162 (2016). [202] Bhat, R., Glimm, T., Linde-Medina, M., Cui, C., and Newman, S.A. Synchronization of Hes1 oscillations coordinates and refines condensation formation and patterning of the avian limb skeleton. Mech. Dev. 156,(2019), 4154. [203] Rutkowski, T.P., Kohn, A., Sharma, D., Ren, Y., Mirando, A.J., and Hilton, M.J. HES factors regulate specific aspects of chondrogenesis and chondrocyte hypertrophy during cartilage development. J. Cell. Sci. 129, (2016), 21452155. [204] Ermentrout, B. Stripes or spots? Nonlinear effects in bifurcation of reaction-diffusion equations on the square.(1991). Proc. Roy. Soc. London Ser. A 434, 413417. [205] Harfe, B.D., Scherz, P.J., Nissim, S., Tian, H., McMahon, A.P., and Tabin, C.J. Evidence for an expansionbased temporal Shh gradient in specifying vertebrate digit identities. Cell 118, 517528. [206] Edelstein-Keshet, Leah. (2005). Mathematical Models in Biology (Society for Industrial and Applied Mathematics). [207] Woolley, T.E., Baker, R.E., Tickle, C., Maini, P.K., and Towers, M. (2014). Mathematical modelling of digit specification by a sonic hedgehog gradient. Developmental Dynamics 243, 290298. [208] Swalla, B.J., and Solursh, M. (1986). The independence of myogenesis and chondrogenesis in micromass cultures of chick wing buds. Developmental Biology 116, 3138. [209] Murphy, M., and Kardon, G. (2011). Origin of vertebrate limb muscle: the role of progenitor and myoblast populations. Curr. Top. Dev. Biol. 96, 132. [210] Duprez, D. (2002). Signals regulating muscle formation in the limb during embryonic development. Int. J. Dev. Biol. 46, 915925. 38 ฀฀฀฀฀฀฀฀฀ List of various gene products relevant to models in limb chondrogenesis Family Sub-family Role References FGFs (Fi- FGF-1 Involved in embryonic development, [188] broblast known as cell growth, tissue repair Growth FGF acidic Modifies endothelial cell migration and Factors ) proliferation. Controls blastema cell proliferation at the time of limb regeneration of the amphibians. FGF-2, Controls the patterning along Proximo- [43, 92, 135] [12, 55, known as Distal and Anterior-Posterior axes. 79, 117, 130] FGF basic Overexpression of FGF-2 up-regulates proliferation of mesenchymal cells, leads to a duplications along the Anterior-Posterior axis. FGF-4, also Induces limb-bud initiation, growth and [19, 75, 118, 160] known as patterning. FGF-K or Promotes stem cell proliferation. K-FGF FGF-5 Enhances expansion of connective tis- [25] sue fibroblasts. Suppresses skeletal muscle development in the limb. FGF-7, also Induces the formation of an AER in [185] known as dorsal median ectoderm. KGF FGF-8 Stimulates the activities of AER. [19, 27, 78, 160] Participates in the initiation of Shh expression in the mesoderm. Maintains mesoderm outgrowth and Shh expression in the established limb bud. FGF-9, also Involved in formation of proximal skele- [68, 160, 171] known as tal element in the developing limb. HBGF-9 and Regulates early stages of chondrogenGAF esis and promotes skeletal vascularization and osteogenesis. FGF-10 Regulates fgf10 gene expression in the [93, 101, 122, 148] lateral plate mesoderm and may be involved in the determination process of the limb territories. Acts as an endogenous initiator for limb formation. Involved in communication between limb mesenchyme and AER. 39 ฀฀฀฀฀฀฀฀฀ FGF-18 Hedgehog Family Indian hedgehog (Ihh) Shh (Sonic hedgehog) Notch family Notch-1 HES family HES1 TGFβ (transforming growth factor beta) Galectins Galectin-1 Galectin-8 RA (retinoic acid, a metabolite of vitamin A) Plays a negative up-regulating role in skeletal development and bone homeostasis. Acts for specification of L-R asymmetry on limb development. Lack of FGF18 in mice results in expanded zones of proliferating and hypertrophic chondrocytes and increased chondrocyte proliferation, differentiation, and Indian hedgehog signaling. Involved in the growth of the endochondral skeleton, but is not directly involved in limb development. Regulates vertebrate organogenesis, especially the growth of digits of limbs. Shh is secreted from the ZPA. Regulates interactions between physically adjacent cells. Contributes in the regulation of mesenchymal apoptosis during digit formation. Involved in limb mesenchymal development, especially has an impact on autopod from the dorsal and/or ventral ectoderm. Notch downstream target, suppresses chondrogenesis, oscillates during precartilage condensation, may coordinate condensation formation. Regulates chondrocyte formation, proliferation and differentiation during limb development. [82, 123, 121] Mediates cell-cell adhesion, enhances precartilage condensations, upregulates expression of Gal-8; Chicken-Gal-1A and -Gal-8 are the earliest marker of sites of pre-cartilage condensations Suppresses precartilage condensations, upregulates Gal-1 Stimulates the growth of the posterior end of the limb. [15, 53] 40 [152, 176] [176] [160, 9, 22, 28, 44, 61, 134, 145, 146, 161] [126] [202, 203] [47, 63, 77, 86] [131] ฀฀฀฀฀฀฀฀฀ Wnt/betacatenin TBox TBX4 TBX5 Hox genes HoxA HoxB BMP(Bone morphogenetic protein) SOX Family & Acts during: 1. limb initiation, 2. limb patterning, 3. late limb morphogenesis, 4. myoblast differentiation in the limb, 5. long bone development. Induces the growth of hindlimb. Accelerates the expression of FGF10 and growth of forelimb Control the body plan along the craniocaudal axis of an embryo. At a specific position, Hox genes are sequentially activated in a rostrocaudal pattern and this is crucial for the induction of limb growth. Upregulates FGF expression. SOX9 [23, 62, 184] [129] [129] [129] [6, 2, 132] Early marker of precartilage condensa- [8, 71, 132] tions, necessary for formation of cartilage, but not for the formation of precartilage condensations. Table 2: Gene products relevant to models in limb chondrogenesis. 41