Hybrid reaction-diffusion and clock-and-wavefront model for the arrest of oscillations in the somitogenesis segmentation clock
Jesús Pantoja-Hernández, Víctor F. Breña-Medina, Moisés Santillán
RREACTION-DIFFUSION MODEL FOR THE ARREST OF OSCILLATIONS IN THESOMITOGENESIS SEGMENTATION CLOCK
JES ´US PANTOJA-HERN ´ANDEZ , V´ICTOR F. BRE ˜NA-MEDINA , AND MOIS´ES SANTILL ´AN Abstract.
The clock and wavefront model is one of the most accepted models for explaining the embryonicprocess of somitogenesis. According to this model, somitogenesis is based upon the interaction between agenetic oscillator, known as segmentation clock, and a moving wavefront, which provides the positionalinformation indicating where each pair of somites is formed. Recently, Cotterell et al. (2015) reported aconceptually different mathematical model for somitogenesis. The authors called it a progressive oscillatoryreaction-diffusion (PORD) model. In this model, somitogenesis is driven by short-range interactions andthe posterior movement of the front is a local, emergent phenomenon, which is not controlled by globalpositional information. With the PORD model, it was possible to explain some experimental observationsthat are incompatible with the clock and wavefront model. However the PORD model has the disadvantageof being quite sensitive to fluctuations. In this work, we propose a modified version of the PORD modelin order to overcome this and others inconveniences. By means of numerical simulations and a numericalstability analysis, we demonstrate that the modified PORD model achieves the robustness characteristic ofsomitogenesis, when the effect of the wavefront is included. Introduction
Somitogenesis, the process by which somites are formed, is an essential developmental stage in manyspecies. Somites are bilaterally paired blocks of mesoderm cells that form along the anterior-posterior axisof the developing embryo in segmented animals (Maroto et al., 2012). In vertebrates, somites give rise toskeletal muscle, cartilage, tendons, endothelial cells, and dermis. Somites form with a strikingly regularperiodicity, that is preserved among embryos of a single species. From this, and other reasons, scientistshave been attracted to somitogenesis for decades. One of the earliest conceptual attempts to explain thisregularity is the so-called clock and wavefront model, originally proposed by Cooke and Zeeman (1976).According to this model, somitogenesis occurs due to the interaction between: (i) autonomous oscillationsof a network of genes and gene products, which causes presomitic-mesoderm cells to oscillate between apermissive and a non-permissive state, in a consistently (clock-like) timed fashion; and (ii) a wavefront (alsoknown as determination front) of signaling that slowly progresses in an anterior-to-posterior direction. Asthe wavefront comes in contact with cells in the permissive state, they undergo a mesenchymal-epithelialtransition, forming a somite boundary, and resetting the process for the next somite.The clock and wavefront model gained relevance when the expression of several genes under the Notch,Wnt, and FGF pathways was discovered to oscillate cyclically, with the same period as that of somiteformation (Palmeirim et al., 1997; Pourqui´e, 2001; Gibb et al., 2010; Pourqui´e, 2011). This, together with theexistence of substances, like Wnt3a and FGF8, whose concentrations vary along the presomitic mesoderm incharacteristic patterns that travel at constant velocity in an anterior to posterior direction, seemed to confirmthe general assessments of the clock and wavefront model (Dubrulle et al., 2001; Gibb et al., 2010; Pourqui´e,2011). As a matter of fact, this model has been very successful because it agrees with numerous experimentalobservations, and because much progress has been done in elucidating the segmentation-clock clockwork, andthe way it interacts with the Wnt3a and FGF8 gradients. Despite its success, some experimental observations
Date : September 8, 2020. a r X i v : . [ q - b i o . CB ] S e p J. Pantoja-Hern´andez, V. Bre˜na–Medina and M. Santill´an AR β Figure 1.
Schematic representation of the gene regulatory network introduced by (Cotterellet al., 2015) and studied in the present work. This network consists of two genes: anactivator A and a repressor R . Solid lines ending in arrowheads (hammerheads) denotepositive (negative) regulation. β represents external regulation of gene A via a wavefrontsubstance. Dashed lines correspond to passive diffusion into the extracellular medium. Thediffusion process depicted by the blue arrow was not originally included in the Cotterellet al. model, but is included in the modified model version studied in this work.have been reported that are incompatible with the clock and wavefront paradigm. For instance, somites canform in the absence of Wnt or FGF gradients, albeit in a disorderly fashion (Naiche et al., 2011; Dias et al.,2014).Recently, Cotterell et al. (2015) introduced a novel mathematical model for somitogenesis. Contrarilyto the clock and wavefront model, in the model by Cotterell et al.—which the authors call a progressiveoscillatory reaction-diffusion (PORD) system—somitogenesis is driven by short-range interactions. Hence,the posterior movement of the front is a local, emergent phenomenon, that is not controlled by globalpositional information. As shown by Cotterell et al., their PORD model is compatible with several importantfeatures of somitogenesis (like size regulation), that former reaction-diffusion models were unable to explain.Furthermore, the PORD model makes predictions regarding FGF-inhibition and tissue-cutting experimentsthat are more consistent with experimental observations than those of clock and wavefront models.At the core of the Cotterell et al. model there is a gene regulatory motif schematically depicted in Fig. 1.In this network, a hypothetical gene codes for an activator protein (A), which enhances its own expression, aswell as that of another gene that codes for a hypothetical repressor protein (R). In turn, the repressor proteindown-regulates the gene coding for the activator, and is able to diffuse into the extracellular medium andaffect neighboring cells. Cotterell et al. did not identify the network genes, but showed that in case it existed,the gene network in Fig. 1 can explain most of the experimental observations regarding somitogenesis. Inparticular, with the correct parameter values, such gene network can generate sustained oscillations, andwhen repressor diffusion is included, it gives rise to a stationary pattern of gene expression like that observedin somitogenesis, even in the absence of a FGF wavefront. eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 3 Even though the hypothetical genes in the Cotterell et al. model have not been identified, we find thismodel interesting for the following reasons: (i) as far as we have observed, it accounts for some importantfeatures of somitogenesis that traditional clock and wavefront models fail to explain; (ii) a gene networkwith a similar architecture as that of Fig. 1—albeit with a different dynamic behavior—has been invokedto explain oscillation arrest in somitogenesis (Santill´an and Mackey, 2008; Zavala and Santill´an, 2012); andfinally, (iii) the architecture of the gene network in Fig. 1 is an ubiquitous motif in the intricate transcription-factor regulatory network of the genes under the Wnt, Notch and FGF pathways (Gibb et al., 2010; Zavalaand Santill´an, 2012).In opposition to its multiple virtues, the Cotterell et al. model is quite sensitive to fluctuations. Forinstance, when random perturbations of initial conditions or intrinsic noise are present, a disordered patternof gene expression arises. This behavior resembles what happens with PSM cells that are not under theinfluence of a wavefront, but contrasts with the observed robustness of somite formation under normalconditions. We hypothesize from the above discussion that, although a wavefront is not strictly necessary forthe Cotterell et al. model to generate a gene expression pattern consistent with somitogenesis, the interactionwith a wavefront is essential to explain the observed robustness of this phenomenon. The present work isaimed at testing this hypothesis, and discussing the corresponding biological implications.The manuscript is organized as follows. In section 2 we develop an extension of the Cotterell et al. model,which is aimed at avoiding some of the former model inconveniences. In section 3 we define the parameterset where the dynamics of the non-diffusion system give place to a limit cycle family. In section 4 we explainthe numerical methods used as well as the initial and boundary conditions. In section 5 we show resultsof the simulations for each one of the selected conditions. Finally, we discuss the relevance of the resultsobtained as well as the limitations of the model in section 6.2.
Mathematical Model
In this section, we introduce a modified version of the Cotterell et al. model. As the functions taken intoaccount in the original model are unrealistic from a biochemical point of view, we propose key modifications,which consist in replacing the gene-expression regulatory functions and introducing a diffusion transportprocess for both activator and repressor protein concentrations.Consider the gene network schematically represented in Fig. 1. Assume that the half life of mRNAmolecules corresponding to genes A and R is much shorter than that of the corresponding proteins. Then,a quasi-stationary approximation can be made for the equations governing mRNA dynamics, which yieldsthe following reaction-diffusion system for the concentration of proteins A and R : ∂A∂t = α A P A ( A, R ) − µ A A + D A ∇ A , (1a) ∂R∂t = α R P R ( A ) − µ R R + D R ∇ R , (1b)where α A and α R are the maximum possible rate of activator and repressor production; P A and P B are theprobabilities that the promoters of genes A and R are active; degradation rate constants of each proteinconcentration are denoted by µ A and µ R , respectively; and D A and D R are the corresponding diffusioncoefficients.To take into account the roles of the activator and the repressor, P A ( A, R ) must be a monotonic increasing(decreasing) function of A ( R ), while P R ( A ) ought to be a monotonic increasing function of A . Indeed, J. Pantoja-Hern´andez, V. Bre˜na–Medina and M. Santill´an
Cotterell et al. (2015) proposed the following functions: P A ( A, R ) = Φ (cid:18) l A − l R + β l A − l R + β (cid:19) , (2a) P R ( A ) = l A l A , (2b)where l , l , and l define the strengths of regulatory interactions between A and R , and β is the backgroundregulatory input of A . To prevent negative values, Cotterell et al. introduced the function Φ( x ) = x.H ( x ),where H ( x ) is the standard Heaviside function ( H ( x ) = 1 for x ≥ H ( x ) = 0 for x < P A ( A, R ) defined in (2a) fulfills the requirement of being a monotonic decreasingfunction of R and a monotonic increasing function of A , it shows features that are biologically challenging: • There is neither biological nor biochemical motivation for the introduction of function Φ. • Function Φ( x ) is non smooth at the origin. This feature is unusual in biologically inspired mathe-matical models, and may cause unexpected complications while studying the dynamical system. • Since the term l R appears with negative sign in the denominator of the argument of function Φ, inthe right hand side of Eq. (2a), function P A ( A, R ) may be divergent for certain values of A and R .In order to address the above discussed issues, we propose a slightly novel approach by substituting thepromoter-activation probabilities in (2) with functions that are consistent with the assumptions that theactivator and the repressor compete for the same binding site in the promoter region of gene A , and thatboth the activator and the repressor interact with their corresponding binding sites in the promoter regionsof genes A and R in a cooperative fashion (Santill´an, 2008). In so doing, we have P A ( A, R ) = β + ( A/K ) n A/K ) n + ( R/K ) n , (3a) P R ( A ) = ( A/K ) n A/K ) n , (3b)where K denotes the half saturation constant for the binding reaction between the activator and thepromoter of gene A ; the Hill coefficient that accounts for a cooperative interaction between the activatorand gene A promoter is given by n ; as in the original model, β is the background regulatory input of A ;the half saturation constant and Hill coefficient of the interaction between the repressor and the promoterof gene A are represented by K and n , respectively; and K and n , equivalently, are the half saturationconstant and Hill coefficient for the interaction between the activator and the promoter of gene R .Notice that (1) along with (3) constitute a reaction-diffusion system for the gene expression networkdepicted in Fig. 1, which accounts for spatio-temporal interactions between these two protein concentrations.Upon re-scaling position and time as x (cid:48) = x/L , y (cid:48) = y/L , z (cid:48) = z/L , and t (cid:48) = tµ A and substituting thefollowing dimensionless variables and parameters a = Aµ A α A , d a = D A µ A L , k = K µ A α , k = K µ R α ,r = Bµ r α r , d r = D R µ A L , k = K µ R α R , µ = µ R µ A , eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 5 where L is a characteristic length on system (1)-(3), we obtain the reaction-diffusion system in a dimensionlessform ∂a∂t (cid:48) = P a ( a, r ) − a + d a ∇ (cid:48) a, (4a) ∂r∂t (cid:48) = µ (cid:0) P r ( a ) − r + d r ∇ (cid:48) r (cid:1) , (4b)where ∇ (cid:48) is the Laplacian with respect to ( x (cid:48) , y (cid:48) , z (cid:48) ), and P a ( a, r ) = β + ( a/k ) n a/k ) n + ( r/k ) n , (4c) P r ( a ) = ( a/k ) n a/k ) n . (4d)To ease notation, we suppress symbol ( · ) (cid:48) from this point onward along the paper.3. Parameter estimation
Since the genes here modeled are hypothetical, it is impossible to estimate the model parameter valuesfrom experimental data. Instead, we performed a bifurcation analysis of the system with no diffusion( d a = d r = 0), employing a continuation method implemented in XPPAUT (Ermentrout, 1987). The resultsof this analysis are presented in Fig. A.1 of Appendix A. From this analysis, we determined the parameterintervals (see Table 1) for which the system shows sustained oscillations in the absence of diffusion, whichis crucial to set the periodicity of somitogenesis. In our simulations, we consider parameter values in themiddle of the intervals reported in Table 1. Following (Cotterell et al., 2015), we fixed n = µ = 1, andassumed that n = n = n . After this, the only free parameters in the model are the diffusion coefficients. k ∈ [0 . , . k ∈ [0 . , . k ∈ [1 . , . n ∈ [2 . , . β ∈ (0 , . Table 1.
Parameter intervals for which the system with no diffusion shows sustained stableoscillatory behavior. 4.
Numerical methods
Under the supposition that the presomitic mesoderm (PSM) can be regarded as one-dimensional, weconsidered a single spacial dimension, x , with boundaries at x = 0 and x = 1. These boundaries set anobservation window in the PSM, where x = 0 corresponds to the posterior extreme. To numerically solvesystem (4), we implemented a standard finite-difference three-point stencil and Euler’s algorithm in Julia ;homogeneous Neumann boundary conditions were included: ∂a∂x (cid:12)(cid:12)(cid:12)(cid:12) (0 ,t ) = ∂a∂x (cid:12)(cid:12)(cid:12)(cid:12) (1 ,t ) = ∂r∂x (cid:12)(cid:12)(cid:12)(cid:12) (0 ,t ) = ∂r∂x (cid:12)(cid:12)(cid:12)(cid:12) (1 ,t ) = 0 , (5)with, unless otherwise stated, the following initial conditions: a ( x,
0) = (cid:26) .
05 for x = 1 , ≤ x < , r ( x,
0) = 0 , for 0 ≤ x ≤ . (6) J. Pantoja-Hern´andez, V. Bre˜na–Medina and M. Santill´an x t A x t B x t C Figure 2.
Spatio-temporal evolution of a from system (4), with boundary condition asin (5), for different initial conditions: (A) r ( x,
0) = 0 for all x ∈ [0 , a ( x,
0) = 0 for0 ≤ x <
1, and a (1 , t ) = 0 .
05; (B) r ( x,
0) = a ( x,
0) = 0 for all x ∈ [0 , r ( x,
0) = 0and a ( x,
0) = 0 for all x ∈ [0 ,
1] and a ( x (cid:48) ,
0) = 0 .
05 at x (cid:48) = 0 . , . , . k = 0 . k = 0 . k = 2, n = 3, β = 0 . d a = 0, and d r = 2 . × − .That is, the system is assumed to be initially homogeneous, except for a perturbation at the anterior extremeof the observation window.We also performed stochastic simulations in which additive white noise was added to the system. Forthese simulations we substituted equations (4a) and (4b) by ∂adt = P a ( a, r ) − a + d a ∇ a + c v a dW, (7a) ∂rdt = µ ( P r ( a ) + r − d r ∇ r ) + c v r dW, (7b)where a and r respectively denote the stationary values of a and r , c v is the coefficient of variation ofthe added white noise, and dW is a normally-distributed white noise term with mean zero and variance 1.To solve this system of stochastic partial differential equations we employed the Euler-Maruyama method,implemented in Julia . 5.
Results
We started by reproducing the results in (Cotterell et al., 2015). To this end, we set d a = 0, d r = 2 . × − , β = 0 .
5, and numerically solved the model equations as described in Section 4. The simulation results areshown in Fig. 2A. Observe that an oscillatory behavior gradually gives rise to a steady pattern consisting ofalternated high and low gene-expression regions. Notice as well that the pattern formation dynamics startat the initial perturbation position, and propagate with constant speed. The initial perturbation is essentialfor the appearance of the pattern. If the system is initially homogeneous, the oscillatory behavior continuesindefinitely—see Fig. 2B. On the other hand, when more than one initial perturbations are present, each oneof them originates a pattern-formation wave, and when two such waves collide, they cancel out—see Fig. 2C.These results suggest that the transition from an oscillatory behavior to a steady pattern of gene expressionis due to a diffusion driven instability interacting with a limit cycle. To verify this, we investigated thespatial stability of the dynamical system in (4) and (5) in Appendix A. We were able to confirm that thelimit cycle is unstable with the current parameter values set; see Fig. A.3, panel (c), at κ = 0.The results described in the previous paragraph are consistent with the reported experimental observationsthat somites can form in the absence of a wavefront. As a matter of fact, they explain why somites formalmost simultaneously and irregularly; in other words, any initial perturbation in the mesoderm tissue rapidlyoriginates a somite boundary, and the emerging pattern formation wave almost immediately collides withneighboring waves. However, even though the wavefront is not strictly necessary for the formation of somites, eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 7 x t A x t B Figure 3.
Spatio-temporal evolution of a from system (4), with boundary condition asin (5), for: (A) spatially stable with β = 0 .
5; (B) spatially unstable with β = 0 .
01. In bothcases, the initial conditions for a were selected from random uniform distributions in theinterval [0 , . r ( x,
0) = 0 for all x ∈ [0 , k = 0 . k = 0 . k = 2, n = 3, d a = 5 × − and d r = 2 . × − .there are several reports that confirm its importance (Sawada et al., 2001; Naiche et al., 2011) . In particular,the regularity of somite formation has been demonstrated to be extremely robust to many different kinds ofperturbations on both the mesoderm tissue and the differentiation wavefront, which is not explained by thecurrent PORD model.We hypothesize that somitogenesis robustness can be achieved via the interaction between the gene net-work and a wavefront. This wavefront is originated by substances, like FGF8 (in what follows we refer toFGF8 only, but the discussion applies to other substances as well, like Wnt3a), which are produced in theembryo tail bud and diffuse to the rest of the PSM. In consequence, FGF8 concentration decreases in theposterior to anterior direction. On the other hand, as the embryo grows, the tail bud recedes leaving PSMcells behind. Hence, the spatial FGF8 distribution moves in the anterior to posterior direction, like a wave-front, as time passes. In other words, we expect that a stable limit cycle emerges for large enough values of β (that accounts for the PSM interaction with the wavefront), which maintains an oscillatory gene expressiondespite perturbations, thus preventing somite formation close to the tail bud. Furthermore, if the limit cycleturns unstable below a given β threshold, any perturbation would lead to the formation of a somite at thePSM position where the threshold is reached. This dynamic behavior may explain why somites are formedin an irregular fashion in the absence of FGF8, and how a FGF8 wavefront robustly drives somitogenesis.In Appendix, A we analyze the stability of the PDE system when depending of β and d a parameters,and found that for a fixed value for d a = 5 × − the system undergoes a bifurcation from a stable to anunstable limit cycle as the value of β decreases. To illustrate these findings, we present in Fig. 3 the resultsof two simulations: one in which the limit cycle is stable and another in which it is unstable. In Fig. 4 wepresent the results of two more runs where, instead of considering a constant value of β , we assume that itis given by β ( x, t ) = β K n K n + ( x − vt ) n . (8)Notice that this expression corresponds to a external regulation of gene A, which is a temporarily andspatially dependent profile that decays in a sigmoidal fashion in the posterior to anterior direction, andtravels with speed v > K = 0 . n = 16 (these values were also employed for the simulation inFig. 5), and considered two different speed values v = 0 .
02, and v = 0 .
04. Observe that, in both simulations,the pattern formation waves have the same periodicity, but the regions corresponding somites are largerfor the faster wave. This result agrees with the experimental observation that somites are larger when thevelocity of FGF8 profile is increased (Sawada et al., 2001).
J. Pantoja-Hern´andez, V. Bre˜na–Medina and M. Santill´an x t A x t B x t C x t D Figure 4.
Spatio-temporal evolution of a from system (4) for β -wave speed as in (8), withboundary condition as in (5). (A) β -wave evolution and (B) somite-pattern formation forspeed v = 0 .
02; (C) β -wave evolution and (D) somite-pattern formation for speed v = 0 . a consists of random uniform distributions in the interval [0 , . r ( x,
0) = 0 for all x ∈ [0 , k = 0 . k = 0 . k = 2, n = 3, d a = 5 × − and d r = 2 . × − . x t A x t B x t C Figure 5.
Spatio-temporal evolution of A from system (7) for β -wave speed as in (8), withboundary condition as in (5): (A) β -wave evolution for speed v = 8 for somite-patternformation with noise intensity (B) c v = 0 .
05 and (C) c v = 0 .
1. Initial conditions for A consists of random uniform distributions in the interval [0 , . R ( x,
0) = 0 for all x ∈ [0 , k = 0 . k = 0 . k = 2, n = 3, d a = 5 × − and d r = 2 . × − .Finally, we performed simulations in which additive white noise was added to both variables ( a and r ),and confirmed the robustness of the system behavior to this kind of variability. As can be seen in Fig. 5,results of a typical simulation show that, when the noise coefficient of variation is c v = 0 .
05, somite formationproceeds in a very precise way, despite the relatively large fluctuations in the gene expression level. On thecontrary, when c v = 0 .
1, although somites continue sequentially emerging, the periodicity of somitogenesisis lost, as well as the regularity of somite sizes.To summarize, the results presented in the paragraphs above confirm the hypothesis that a PORD model(modified to include diffusion of both the repressor and the activator) undergoes a bifurcation from a stable eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 9 to an unstable limit cycle, as a result of its interaction with the wavefront, and that this bifurcation isenough to explain the observed robustness of somitogenesis. To the best of our understanding, the modifiedPORD model has a couple of quintessential features: it explains why somitogenesis may occur in the absenceof a wavefront and its inherent somites irregular prompting, and how the process achieves robustness as aconsequence of the wavefront. 6.
Concluding Remarks
We have studied a modified version of the somitogenesis PORD model originally introduced by Cotterellet al. (2015). The most important of the introduced modifications is the assumption that the activatorprotein can also diffuse, although with a much smaller diffusion coefficient than the repressor. With thismodification, the model undergoes a bifurcation, from a stable to an unstable limit cycle, as the value ofthe parameter accounting for the background regulatory input of the activator decreases. From a biologicalperspective, the bifurcation just described allows the model to explain why somites can form in the absenceof a wavefront (which traditional front and wavefront models failed to explain), reassesses the role of thewavefront as a conductor for somitogenesis, and makes the model behavior robust to random fluctuations;notice that the latter is one of the weak points of the original PORD model.In the clock and wavefront models, there is a consensus that the oscillatory behavior is originated by a genenetwork with time-delayed negative feedback regulation (Monk, 2003; Lewis, 2003). This claim is supportedby multiple experimental reports, which have elucidated some of the underlying regulatory mechanisms(Schr¨oter et al., 2012). In contrast, the present and the Cotterell et al. (2015) models not only rely on adifferent gene network architecture to generate oscillations, but the genes in the network are hypothetical.From this perspective, the weight of experimental evidence seems to favor clock and wavefront models.Nonetheless, upon taking this into consideration, we believe that we have provided convincing evidence thatreaction-diffusion and positional information (wavefront) mechanisms could work together in somitogenesis.Further investigating this possibility may allow a better understanding of such a fascinating phenomenon.
Acknowledgements
JP-H thanks CONACYT-M´exico for granting him a doctoral scholarship. VFBM thanks the financialsupport by Asociaci´on Mexicana de Cultura AC. MS acknowledges the financial support of CONACYT-M´exico, grant INFRA-302610.
AppendixAppendix A. Somite pattern-formation dynamical features
System (4), along with homogenous Neumann boundary conditions, gathers the essential ingredients ofsomite pattern-formation dynamics. Particularly, the kinetic terms consist of Hill functions, whose co-efficients corresponding to a non-negative cooperative interaction. Thus, we assume that n = 1 and n = n = n ≥ µ = 1. Initially, we set d a = d r = 0. In so doing, we get the purelykinetic system: dadt = f ( a, r ) , drdt = g ( a, r ) , (A.1a)where the field components are given by f ( a, r ) = β + ( a/k ) n a/k ) n + ( r/k ) n − a , g ( a, r ) = ak + a − r . (A.1b)This system steady states satisfy the relation g ( a ) = β , where (a) (b) (c) (d) (e) Figure A.1.
Bifurcation diagrams of system (A.1) for slowly varying: (a) background reg-ulatory input of activator β , effective half-saturation parameters (b) k , (c) k and (d) k ,and (e) Hill coefficient n . Blue (orange) solid lines correspond to stable (unstable) branchesof steady states, and green (red) lines are stable (unstable) limit-cycle branches. Squares(filled circle) indicate supercritical (subcritical) Hopf bifurcations, and triangle is for indi-cating fold bifurcations. Other parameter set values are k = 0 . k = 0 . k = 2 . n = 3 .
0, and β = 0 .
5, respectively. eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 11 g ( a ) := a − a n /k n + a n +1 /k n + a (cid:18) a/k k + a (cid:19) n . (A.2)From (A.2), notice that g ( a ) satisfies that: (i) g (0) = 0, and (ii) g ( a ) → + ∞ , when a → + ∞ . In consequence,there exists a ∗ > g ( a ∗ ) = β >
0, which further implies that w ∗ = ( a ∗ , r ∗ ) is a steady-state ofsystem (A.1), with a ∗ > r ∗ >
0, since r ∗ = a ∗ k + a ∗ . In other words, the existence of at least one steady-state in the first quadrant is guaranteed. We canstraightforwardly prove, from the Poincare-Bendixon theorem, that a limit-cycle family emerges as a resultof this steady-state undergoing a Hopf bifurcation (HB). In order to disclose this implication, we performed anumerical continuation by using the algorithms implemented in XPPAUT (Ermentrout, 1987). The resultingbifurcation diagrams obtained by slowly varying each parameter of system (A.1) are depicted in Fig. A.1.Notice that the system undergoes HBs for all parameters values in Table 1. In Fig. A.1(a), a single HBfor parameter β takes place, where periodic orbits vanish at the bifurcation point HB. This suggests thatthe lower the β -input value, the larger the amplitude and the longer the period of stable orbits for thehomogeneous system; and that no periodic orbits occur for β (cid:29)
1, nonetheless. In contrast, as is shown inFigs. A.1(b)-(d), two supercritical HB points occur for parameters k , k and k , which determine a finiteinterval for each parameter wherein a family of periodic orbits exist. Interestingly, a bi-stability interval forparameter n is delimited by n -values where a subcritical HB and fold bifurcation (LP) points happen. Thatis, a branch consisting of unstable limit cycles emanates from the subcritical HB, which stabilizes at theLP points. Hence, stable steady-states coexist with stable limit-cycles for this interval, where the unstableperiodic branch plays a critical role for initial states. This result, as is depicted in Fig. A.1(e), indicatesthat an on-and-off gene switch is present, which provides hysteretical features triggered by key values of thecharacteristic Hill coefficient. We may argue from this that the cooperativity in the gene regulatory networkfavors system robustness to parameter variations.On the other hand, as has been discussed above, the external regulation of gene A is an essential ingredientof the somite pattern formation dynamics. Recall that this element is captured by β . In addition, theactivator diffusion also plays a crucial role in the somite pattern-formation dynamics. We study now theinterplay between β and d a , which triggers the somite formation mechanism that we have proposed. To doso, we include diffusion terms in (A.1) to get the reaction-diffusion system ∂a∂t = f ( a, r ) + d a ∇ a , ∂r∂t = g ( a, r ) + d r ∇ r , (A.3)where the kinetic terms are as in (A.1b).Now, upon defining w = ( a, r ) T , we obtain that system (A.3) can be set up in vector notation as w t = F ( w ) + D ∇ w , where F ( w ) = ( f ( a, r ) , g ( a, r )) T and D = diag ( d a , d r ). In so doing, for an isolated root of F ( w ) = 0 given by w ∗ = ( a ∗ , r ∗ ), system (A.3) has a local solution of the form w ( x, t ) = ∞ (cid:88) m =0 γ m e − λ ( κ ) t w m ( x ) , (A.4)where w m ( x ) satisfies the Helmholtz equation ∇ w m + κ w m = 0, in which the so-called wave mode is denotedby κ . The Fourier coefficients γ m are determined by the initial conditions, and λ determines whether (A.4)converges, and hence is bounded, as t → + ∞ . These three parameters not only shape solution (A.4), butalso are intrinsic to the geometry and boundary conditions of the system into consideration. Moreover, theyalso depend on the wave number m ∈ N ; see (Murray, 1989) for further details. β d a pattern classi fi cation I Turing PatternII Turing-Hopf PatternIII Hopf PatternIV No Pattern
I II IIIIV IIV
Figure A.2.
Two-parameter space for d a and β . Each pattern region is plotted in colorsaccordingly to the table in the right-hand side, and transition lines correspond to each regionboundary. Other parameter set values are d r = 2 . × − , k = 5 × − , k = 10 − , k =2 . , n = 3 . w ∗ to get, in vectorial form, w t = Jw + D ∇ w , (A.5)where J is the Jacobian matrix at w ∗ . Now, as our interest lies on the dynamics in one spatial dimension,we have that w m ( x ) = cos( κx ), where κ = mπ/L satisfies the Helmholtz equation for homogenous Neumannboundary conditions as in (5). Thus, (A.5) is satisfied by (A.4), when the dispersion relation λ = λ ( κ ) isgiven by | J − Dκ − λI | = 0 , I ∈ R × , (A.6)which can be seen by substituting (A.4) into (A.5). As a result, it relates the temporal growth rate and thespatial wave mode κ , which parametrises the finite spatial domain. Note that (A.6) leads to λ + b ( κ ) λ + c ( κ ) = 0 , (A.7a)where b ( κ ) = ( d a + d r ) κ − ( f a + g r ) , (A.7b) c ( κ ) = d a d r κ − ( d r f a + d a g r ) κ + ( f a g r − f r g a ) . (A.7c)As we are interested in the linear stability of the steady state w ∗ , characterized by parameters β and d a , wenotice that the parameter space for spatial instability of Turing type is given by conditions: f a + g r < , f a g v − f r g a > , (A.8a) d r f a + d a g r > , d r f a + d a g r > (cid:112) d a d r ( f a g r − f r g a ) . (A.8b)These conditions provide the ingredients that give place to non homogeneous spatially extended stationarystates. We are however interested in a mechanism that triggers sustained oscillations of the gene networkthat give place to a stationary pattern in the long term. Such a process is dynamically provided by theTuring-Hopf bifurcation (THB); see, for instance, (Castillo et al., 2016; Liu et al., 2007). From (A.7a), thisbifurcation is prompted by obtaining purely imaginary eigenvalues by slowly varying d a or β . In so doing,we notice that two key conditions must be met: (i) f a + g r = 0 at the THB point, and (ii) dλ ( κ ; p ∗ ) /dp (cid:54) = 0, eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 13 k − − − λ ( κ ) IV β =2.8, d a =5e-05 (a) k − − − − λ ( κ ) I β =2.8, d a =5e-06 (b) k − − − − λ ( κ ) II β =0.3, d a =5e-05 (c) k − − − − λ ( κ ) III β =1.0, d a =1e-04 (d) Figure A.3.
Samples of dispersion relations for each region depicted in Fig. A.2. The realparts of the eigenvalues are in solid lines, and the imaginary parts in dashed lines. Panel (a)corresponds to region IV, where the real part is negative for all κ ; in panel (b) the real parthas two roots, which gives place to a Turing type stationary pattern; in panel (c), Turingand Hopf bifurcations occur for different wave modes: the real part line is positive for wavemodes corresponding to non-zero imaginary part, and a Turing instability occurs as in panel(b); in panel (d), a typical Hopf bifurcation feature is exhibited as the real part is positiveonly for non-zero portions of the imaginary part. The values of d a and β are respectivelyshown in each panel, and other parameters values are k = 5 × − , k = 10 − , k = 2 . , n =3 . p ∗ = d ∗ a or p ∗ = β ∗ represent the THB parameter value. Thus,in order to obtain the parameter regions where a Turing bifurcation, HB and THB occur in system (A.3), wesolve (A.7) for λ ( κ ). In Fig. A.2, the parameter space on scope is portrayed, where four different stabilityfeatures for the selected range values of parameters β and d a are identified in four regions, labeled accordingly:Turing pattern, Turing-Hopf pattern, Hopf pattern, and no pattern. Notice that, at the transition curvebetween region II and III, parameters β and d a follow an inverse relation; in other words, the larger parameter β is, the lower diffusivity d a is needed for the bifurcation to take place. In addition, notice that for a fixedvalue of d a = 5 × − , slowly varying β from 0 up to 1 drives the system through two bifurcations, whichin turn originate three different mechanisms. That is, no pattern is formed for 0 ≤ β < . × − , which is x t (a) I t a , r (b) I ar x t (c) II t a , r (d) II ar Figure A.4.
Time-step simulations of (A.3) with homogeneous Neumann boundary con-ditions and parameter set values as in region I and II plotted in Fig. A.2. An azimuthalview of the spatio-temporal solution shows the pattern formation dynamics for a Turingtype (panel a) and a Turing-Hopf type (panel c). Temporal evolution for each case in theleft-hand side column of the activator and repressor at x = 0 . x = 0 . . × − < β < .
5, to get in the Hopf patternregion, which is held by 0 . < β ≤ . eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 15 x t (a) x t (b) Figure A.5.
Simulations with parameter β as function of x and t as indicated in the text.For (a), the initial conditions are perturbations of the steady state, while in (b) the initialconditions are the final state of (a).by having a negative real part of the dispersion relation for all wave modes, and hence the homogeneoussteady state suffers no symmetry breaking. In other words, the combination of real and imaginary parts of λ ( κ ) in (A.4), and provided by (A.6), determines each pattern type as well as transitions between regions.For a detailed discussion about the former approach, see (Liu et al., 2007).Regions I and II are particularly relevant as a stationary pattern is formed, although a key mark lies ontransitory dynamical behavior for each scenario. To illustrate this distinguished mark, we perform time-stepruns for a setting in both regions, by having a perturbed steady state as an initial condition; see top panelsin Fig. A.4. As can be seen in panel (a), the system is initially in a homogeneous steady state with a slightperturbation. As time goes by, a heterogeneous pattern arises. This is a consequence of the unstable wavemodes as is shown in Fig. A.3, panel (b), which corresponds to the Turing pattern, region I, in Fig. A.2. Inaddition, in panel (b), a time evolution is observed for the activator and repressor states at x = 0 . λ ( κ ). Such oscillatory dynamics goes on until unstable wave modes allowa stationary pattern to arise. Notice that this mark is clearly observed in panel (d); see bottom panels inFig. A.4. In other words, even though both dynamical configurations give place to stationary patterns, thecrucial oscillating feature previous to finally get a fixed pattern is added by having a Turing-Hopf mechanismin play.In Fig A.5, we show additional results, where a spatio-temporal dependent parameter β is taken intoconsideration. There, we take β as in (8) where v = 0 .
02 and v = − .
02 for the run in panels (a) and(b), respectively. The initial conditions for panel (b) corresponds to the final profile of panel (a). Noticethat once the pattern is completely formed, even though β varies in a inverse direction, the pattern is notdestroyed. This is typical trait of a hysteretical process. In other words, this result indicates that theproposed mechanism is robust since, once the wave front depicted by β in (8) prompts the formation ofsomites, this process cannot be undo. References
J. A. Castillo, F. S´anchez-Gardu˜no, and P. Padilla. A Turing–Hopf Bifurcation Scenario for Pattern Forma-tion on Growing Domains.
Bulletin of Mathematical Biology , 78(7):1410–1449, 2016. ISSN 15229602. doi:10.1007/s11538-016-0189-6.J. Cooke and E. Zeeman. A clock and wavefront model for control of the number of repeated structuresduring animal morphogenesis.
Journal of Theoretical Biology , 58(2):455–476, jan 1976. doi: 10.1016/s0022-5193(76)80131-2. URL https://doi.org/10.1016/s0022-5193(76)80131-2 .J. Cotterell, A. Robert-Moreno, and J. Sharpe. A local, self-organizing reaction-diffusion model can explainsomite patterning in embryos.
Cell Systems , 1(4):257–269, oct 2015. doi: 10.1016/j.cels.2015.10.002. URL https://doi.org/10.1016/j.cels.2015.10.002 .A. S. Dias, I. de Almeida, J. M. Belmonte, J. A. Glazier, and C. D. Stern. Somites without a clock.
Science ,343(6172):791–795, 2014. doi: 10.1126/science.1247575.J. Dubrulle, M. J. McGrew, and O. Pourqui´e. FGF signaling controls somite boundary position and regulatessegmentation clock control of spatiotemporal hox gene activation.
Cell , 106(2):219–232, 2001. doi: 10.1016/s0092-8674(01)00437-8.B. Ermentrout.
Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Re-searchers and Students (Software, Environments and Tools) . Society for Industrial and Applied Mathe-matics, 1987. ISBN 0898715067.S. Gibb, M. Maroto, and J. K. Dale. The segmentation clock mechanism moves up a notch.
Trends in CellBiology , 20(10):593–600, oct 2010. doi: 10.1016/j.tcb.2010.07.001. URL https://doi.org/10.1016/j.tcb.2010.07.001 .J. Lewis. Autoinhibition with transcriptional delay.
Current Biology , 13(16):1398–1408, 2003. doi: 10.1016/s0960-9822(03)00534-7.R.-T. Liu, S.-S. Liaw, and P. K. Maini. Oscillatory Turing Patterns in a Simple Reaction-Diffusion System.
J. Kor. Phys. Soc. , 50(1):234–238, 2007.M. Maroto, R. A. Bone, and J. K. Dale. Somitogenesis.
Development , 139(14):2453–2456, jun 2012. doi:10.1242/dev.069310. URL https://doi.org/10.1242/dev.069310 .N. A. Monk. Oscillatory expression of hes1, p53, and NF- κ b driven by transcriptional time delays. CurrentBiology , 13(16):1409–1413, 2003. doi: 10.1016/s0960-9822(03)00494-9.J. D. Murray.
Mathematical Biology . Oxford, second edition, 1989. ISBN 9783540572046. doi: 10.1007/978-3-662-08542-4.L. A. Naiche, N. Holder, and M. Lewandoski. FGF4 and FGF8 comprise the wavefront activity that controlssomitogenesis.
Proceedings of the National Academy of Sciences , 108(10):4018–4023, 2011. doi: 10.1073/pnas.1007417108.I. Palmeirim, D. Henrique, D. Ish-Horowicz, and O. Pourqui´e. Avian hairy gene expression identifies amolecular clock linked to vertebrate segmentation and somitogenesis.
Cell , 91(5):639–648, 1997. doi:10.1016/s0092-8674(00)80451-1.O. Pourqui´e. Vertebrate somitogenesis.
Annual Review of Cell and Developmental Biology , 17(1):311–350,2001. doi: 10.1146/annurev.cellbio.17.1.311.O. Pourqui´e. Vertebrate segmentation: From cyclic gene networks to scoliosis.
Cell , 145(5):650–663, may2011. doi: 10.1016/j.cell.2011.05.011. URL https://doi.org/10.1016/j.cell.2011.05.011 .M. Santill´an. On the use of the hill functions in mathematical models of gene regulatory networks.
Mathematical Modelling of Natural Phenomena , 3(2):85–97, 2008. doi: 10.1051/mmnp:2008056. URL https://doi.org/10.1051/mmnp:2008056 .M. Santill´an and M. C. Mackey. A proposed mechanism for the interaction of the segmentation clock and thedetermination front in somitogenesis.
PLoS ONE , 3(2):e1561, feb 2008. doi: 10.1371/journal.pone.0001561.URL https://doi.org/10.1371/journal.pone.0001561 . eaction-diffusion model for the arrest of oscillations in the somitogenesis segmentation clock 17 A. Sawada, M. Shinya, Y. J. Jiang, A. Kawakami, A. Kuroiwa, and H. Takeda. Fgf/MAPK signalling isa crucial positional cue in somite boundary formation.
Development , 128(23):4873–4880, 2001. ISSN09501991.C. Schr¨oter, S. Ares, L. G. Morelli, A. Isakova, K. Hens, D. Soroldoni, M. Gajewski, F. J¨ulicher, S. J. Maerkl,B. Deplancke, and A. C. Oates. Topology and dynamics of the zebrafish segmentation clock core circuit.
PLoS Biology , 10(7):11, 2012. doi: 10.1371/journal.pbio.1001364.E. Zavala and M. Santill´an. An analysis of overall network architecture reveals an infinite-period bifurcationunderlying oscillation arrest in the segmentation clock.
Mathematical Modelling of Natural Phenomena , 7(6):95–106, 2012. doi: 10.1051/mmnp/20127605. URL https://doi.org/10.1051/mmnp/20127605 . Centro de Investigaci´on y de Estudios Avanzados del IPN, Unidad Monterrey, V´ıa del Conocimiento 201,Parque PIIT, 66628 Apodaca NL, M´exico.2