Modelling Periodic Oscillations during Somitogenesis

We consider a model of genetic network that has been previously presented by J. Lewis. This model takes the form of delay differential equations with two delays. We give conditions for the local stability of the non-trivial steady state. We investigate the condition underwhich stability is lost and oscillations occur. In particular, we show that when the ratio of the time delays passes a threshold, sustained oscillations occur through a Hopf bifurcation. Through numerical simulations, we further investigate the ways in which various parameters influence the period and the amplitude of the oscillations. In conclusion, we discuss the implications of our results. 1. Introduction. In vertebrate embryos, the most obvious metameric units are the somites. They are aligned symmetrically along both sides of the neural tube and give rise to the repetitive structures including the axial skeleton (vertebrae), the dermis of the back, and the back muscles [3]. The generation and individuation of the somites and their derivatives are fundamental to the evolution and diversification of the vertebrate body plan. As a result, somitogenesis is still a very important yet unresolved problem in developmental biology. Somites are formed in a sequential, temporally ordered process that proceeds from the anterior to the posterior of the embryo. New somites are formed at regular intervals (90 minutes in the chick embryo, 120 minutes in the mouse, and 30 minutes in the zebrafish). To understand the mechanism of the formation of somites, several theoretical models have been proposed. Among these are the clock and wave-front model by Cooke and Zeeman[4], the reaction-diffusion model by Meinhardt[12], and the cell-cycle model by Primmett et al. [15, 16]; see also [22]. Recently, a molecular clock, or somitogenesis clock was identified by Palmeirim et al.[17]. It involves oscillations in both mRNA and protein levels of the basic gene c-hairy1. The expression of c-hairy1(mRNA level, or protein level) sweeps across anteriorly repeatedly, and a new somite is formed every time a wave reaches the border. Moreover, this wave-like propagation of gene expression is caused not by cell

1. Introduction.In vertebrate embryos, the most obvious metameric units are the somites.They are aligned symmetrically along both sides of the neural tube and give rise to the repetitive structures including the axial skeleton (vertebrae), the dermis of the back, and the back muscles [3].The generation and individuation of the somites and their derivatives are fundamental to the evolution and diversification of the vertebrate body plan.As a result, somitogenesis is still a very important yet unresolved problem in developmental biology.
Somites are formed in a sequential, temporally ordered process that proceeds from the anterior to the posterior of the embryo.New somites are formed at regular intervals (90 minutes in the chick embryo, 120 minutes in the mouse, and 30 minutes in the zebrafish).
To understand the mechanism of the formation of somites, several theoretical models have been proposed.Among these are the clock and wave-front model by Cooke and Zeeman [4], the reaction-diffusion model by Meinhardt [12], and the cellcycle model by Primmett et al. [15,16]; see also [22].
Recently, a molecular clock, or somitogenesis clock was identified by Palmeirim et al. [17] .It involves oscillations in both mRNA and protein levels of the basic gene c-hairy1.The expression of c-hairy1(mRNA level, or protein level) sweeps across anteriorly repeatedly, and a new somite is formed every time a wave reaches the border.Moreover, this wave-like propagation of gene expression is caused not by cell 2 P FENG AND M NAVARATNA movement but rather by the result of synchronized oscillation of c-hairy1 expression.Subsequent to Palmeirim's study, other cycling genes were identified: Hes1, Hes7, and Hey2 in mouse (Jouve et al. 2000 [8]); Her1 and her7 in zebrafish (Oates and Ho 2002 [14]).Thus, it is proposed that the periodic expression of the cyclic genes is essential for somite segmentation.In fact, it is suggested by Palmeirim et al. [17] that when groups of cells enter a similar phase of the cycle several times, they increase their adhesive properties and aggregate to form a somite.Some synchrony in the cell cycle is observed in the PSM. Figure 1 shows a zebrafish embryo stained by in situ hybridization for gene deltaC mRNA.Oscillation of this gene expression is similar to that of her1 and her7.
However, the cause of the oscillations in gene expression remained to be understood until recently when some mathematical models were proposed to describe such behaviors.In this paper, we study in more detail the mathematical model proposed by Julian Lewis in [9].In Section 2, we discuss the assumptions and the formulation of the model.In Section 3, the steady states and the local stability are considered.In Section 4, we show that Hopf bifurcation may occur when the bifurcation parameter surpasses a threshold.In Section 5, we present our numerical solutions emphasizing the effect of different parameters on the period and amplitude of the periodic solutions.Finally in Section 6, we present a brief conclusion.
Remark 1.The fact that oscillations of a certain gene expression can lead to the formation of somites falls into the scope of one of the most important biological discoveries of the past two decades.Biologists discovered that even though many animals are divergent in forms, they share specific families of genes that regulate the main aspects of their body patterns.For a comprehensive review on this topic, we refer the readers to the wonderful book by Sean Carroll [3].These genes regulate the transcription of cell genes leading to a new state which is defined by the concentration of the gene products.This can be formalized as a concept of genetic network where all genes and their products are members of a network.Methods used to model such genetic regulating networks include the Boolean model that describes the state of genes as either ON or OFF, and the continuous model that uses differential equations to describe the expression of the gene products(concentrations).And there is also the hybrid approach that combines the discrete model with the continuous model.For a comprehensive review on modeling gene networks, we refer the interested readers to [20,21] and references therein.
A simple genetic regulatory network is illustrated in Figure 2. We are interested in the concentration of the proteins and the mRNAs.We shall consider the following simple differential equations, which include the delay effect during the transcription and translation processes.Denoting the number of mRNAs at time t by m and the number of proteins by p, we have the following simple model: where τ p represents the time lag between the initiation of the translation and the appearance of a mature protein and τ m represents the time lag between transcription and the appearance of a mature mRNA molecule; b and c represent the degradation rates of proteins and mRNAs, respectively; a is the translational constant.Also f (p) is a sigmoid function which signifies the switch-like phenomena during the process.Here we shall consider the following type where θ represents the critical protein concentration at which the transcription rate is at half of its basal value and n is the hill coefficient that characterizes the degree of cooperativity.This type of function has often been used in different neural network models.
General mathematical models incorporating delay effect have been studied in [10,11,20].Characteristic of such delay differential equations is that oscillations are generated if the delays surpass a critical value, which is commonly referred to as delay-induced oscillations, see, for example, [2,19].Note that in our present model, we have two delays, and it turns out that when their ratio surpasses a critical value, the oscillations will occur.
Remark 2. The specific gene we are interested here is her1 gene.Overall, τ p could be in the range of one to three minutes and τ m is around ten to twenty minutes; see Table 1. 3. Steady states and local stability.The steady states E * = (p * , m * ) are given by the solutions of am − bp = 0, (3) In fact, p * satisfies bcp n+1 + bcθ n p − kaθ n = 0, which clearly has a unique positive solution.
For τ m = τ p = 0, equation( 1) and equation( 2) read as The characteristic equation of the linearized equation of ( 5)-( 6) around E * is given by (λ + b)(λ + c) − af ′ (p * ) = 0, which clearly has eigenvalues with negative real parts; thus, E * is asymptotically stable.In other words, without delays, the model will not generate sustained oscillations.
Let r = τ m /τ p and we can reduce the system to the following one with dimensionless time t/τ p , which is again denoted by t: To examine the local stability of the steady state E * = (p * , m * ), we linearize equation (7) and equaiton( 8) around E * .To this end, we define or equivalently, p = p * + p * u 1 , m = m * + m * u 2 ; then (7) and (8) become which can be put into the following simpler form where The associated characteristic equation has the following form: Clearly, for each fixed delay τ p , we have a characteristic equation with parameter r.Next we shall study the stability of E * with respect to the ratio parameter r.
The following theorem gives the stability result for the steady state E * .
) is asymptotically stable for all r ≥ 0. When c is small, there exists a critical value r 0 such that the steady state E * is asymptotically stable for r ∈ [0, r 0 ) and unstable for r > r 0 , where and For the proof of Theorem 3.1, we need the following lemma.
Lemma 3.2.[5] Consider the equation where ā, b, c, d and ē are real numbers.Let the hypotheses be as follows: a) If A1, A2 and A3 hold, then all roots of equation( 16) have negative real parts for all r ≥ 0. (b) If A1, A2 and A4 hold, then there exists r 0 > 0 such that when r ∈ [0, r 0 ) all roots of (16) have negative real parts; when r = r 0 , equation( 16) has a pair of purely imaginary roots ±iξ + , and for r > r 0 (16) has at least one root with positive real part.Here r 0 and ξ + are given by The proof of this lemma is quite simple and can be found in [5].Now we are ready to prove Theorem 3.1.
Proof.First we show that for c > ka bθ n−1 n n √ n − 1, E * is asymptotically stable for all r ≥ 0. We shall verify the hypothesis (A1)-(A3) in Lemma 3.2.Comparing characteristic equaiton( 13) with ( 16), we have Thus, E * is asymptotically stable for all r ≥ 0 under the condition c > ka bθ n−1 n n √ n − 1.Now for the change of the stability when c is small.We can verify (A4) as ē2 − d2 := b 2 τ 4 p (c 2 − c 2 1 ) < 0 for c small.From Lemma 3.2, the unique solution of equation( 13) has the following form: and there exists a critical value of the ratio of the time delays so that E * is asymptotically stable for r ∈ [0, r 0 ) and unstable for r > r 0 .This ends the proof of Theorem 3.1.
4. Hopf bifurcation.In this section, we shall study the occurrence of Hopf bifurcation using the ratio of the time delays as the bifurcation parameter.Clearly, time delays do not change the location of the equilibrium in equation ( 1) and equation ( 2), but its stability with the change of a parameter; for example, the ratio of the time delays.If any complex eigenvalue crosses the imaginary axis, then a stable equilibrium loses its stability and oscillations occur locally.We recall the following formulation of the Hopf bifurcation theorem for delay differential equations [7]: Theorem 4.1.We consider the delay differential equation If (a) F is analytic in x and µ in a neighborhood of (0, 0) in R n × R; (b) F (0, µ) = 0 for µ in an open interval containing 0, and x(t) = 0is an isolated stationary solution of ( 17); (c) the characteristic equation of (17) has a pair of complex conjugate eigenvalues λ and λ such that λ(µ) = α(µ) + iw(µ) where w(0) = w 0 > 0, α(0) = 0, α ′ (0) = 0; (d) the remaining eigenvalues of the characteristic equation have strictly negative real part, then the delay differential equation ( 17) has a family of Hopf periodic solutions.More precisely, there is an ǫ 0 > 0 and an analytic function µ(ǫ) = Σ ∞ i=2 µ i ǫ i for 0 < ǫ < ǫ 0 such that for each ǫ ∈ (0, ǫ 0 ) there exists a periodic solution p ǫ (t) occurring for µ(ǫ).If µ(ǫ) is not identically zero, the first nonvanishing coefficient µ i has an even subscript, and there is an ǫ 1 ∈ (0, ǫ 0 ] such that µ(ǫ) is either strictly positive or strictly negative for ǫ ∈ (0, ǫ 1 ).For each L > 2π/w 0 there is a neighborhood V of 0 and an open interval I containing 0 such that for every µ ∈ I the only nonconstant periodic solutions of the delay differential equation (17) with period less than L which lie in V are members of the family p ǫ (t) for values of a satisfying µ(ǫ) = µ, ǫ ∈ (0, ǫ 0 ).For 0 < ǫ < ǫ 0 the period T (ǫ) = 2π[1 + Σ ∞ i=1 τ i ǫ i ]/w 0 of p ǫ (t) is an analytic function.Exactly two of the Floquet exponents of p ǫ (t) approach 0 as ǫ → 0 + .One is 0 for 0 < ǫ < ǫ 0 , and the other is an analytic function β(ǫ) = Σ ∞ i=2 β i ǫ i for 0 < ǫ < ǫ 0 .The periodic solution p ǫ (t) is orbitally asymptotically stable if β(ǫ) < 0 and unstable if β(ǫ) > 0.
Proof.For simplicity, we rewrite the characteristic equation in the following form: where a 1 = (b + c)τ p , e = bcτ 2 p and d = c 1 bτ 2 p .First we recall that at r = r 0 , this equation has a pair of purely imaginary roots ±iξ + .We substitute λ = iξ + and obtain −ξ 2 + + a 1 iξ + + e + de −iξ+(r0+1) = 0, from which we can easily obtain + + e 2 − d 2 = 0.There exists a complex function λ = λ(r) defined in a neighborhood of r 0 such that λ(r 0 ) = ±iξ + .Differentiating on both sides of equation ( 18) with respect to r yields At r = r 0 , λ = iξ + , it follows from equation ( 19) that Note that at r = r 0 , sin(ξ Thus 5. Numerical illustrations.In this section, we would like to see whether her1 could actually generate oscillations by the simple model and how different parameters can affect the period and amplitude of the oscillations.We shall apply the following estimates on the parameters(Table 1).For her7, τ p = 1.7 min and 5.9 < τ m < 20.1.This table is based on [1,6,23,9].If we apply these estimates but keep τ m as a variable, then the steady state is given by (p * , m * ) ≈ (80.8438, 4.1320).Thus we have c 1 ≈ 1.1169 > c, ξ + ≈ 0.5712 and r 0 ≈ 1.8269, and the corresponding critical mRNA delay time is τ m ≈ 5.1153.The period the original system is given by 2π ξ+ • τ p ≈ 31min, which agrees very well with the lab observation that in zebrafish embryo, each somite is formed approximately every 30 minutes. Figure 3 shows that we have sustained oscillations for these parameters.The phase portrait also shows that a limit cycle appears.5.1.Effect of c.In this subsection, we investigate the effect of c (the degradation rate of mRNAs).Our results imply that when c is too large, the steady state E * will always be asymptotically stable.In particular, when we choose the parameters from Table 1 (with τ m ≈ 15min), we arrive at the conclusion that when c > 17.0389, E * is always stable; in other words, the oscillations will be damped and eventually approach E * .On the other hand, the second part of the theorem states that when c < 17.0389, sustained oscillations can be achieved.To clarify the restriction on c, we note from the proof of theorem 3 that only when c < c 1 , there exists a critical value r 0 beyond which we may have sustained oscillations.This condition is equivalent to p * > θ n √ n−1 .From the equation for p * , we have the following restriction on our parameters If we take the parameters in Table 1 but we leave c and Hill coefficient n varying, we obtain the following curve above which the oscillations will be damped and below which the oscillations will be sustained.Figure 5 shows the damping oscillations

Influence of each parameter individually.
In this section, we investigate the effect of each parameter, namely τ m , τ p , b and a, on the period and the amplitude of the periodic solutions.In this part, we computed only the protein level.Our parameters are chosen from Table 1 unless we choose it as a variable.All of our simulations are done through dde23, a delay differential equation solver, which is now part of Matlab official release.We would like to refer the readers to the following website: http://www.runet.edu/thompson/webddes/.
In our simulations, we choose to vary τ m between zeor to twenty minutes, τ p between zero to twenty minutes, the protein degradation rate b between zero to ten molecules/min, the protein synthesis rate between zero to ten molecules/min.The simulations clearly reveal that both τ m and τ p have a strong effect on the total period of the oscillations (Figure 6, upper left and right).Even though one might expect the a reduction of protein synthesis will somewhat affect the period and amplitude.Our simulations clearly show that sustained oscillation is still seen, both the period and the amplitude are almost unchanged (Figure 6, lower left and right).6. Discussion.In this paper we have summarized the mathematical model of J. Lewis for somitogenesis in a zebrafish embryo.We have presented an analysis and numerical simulations of the model equations, which provide a greater understanding of the formation process and allow us to clearly understand which parameters influence the period of the somite-forming process.In particular, we would like to mention the following observations.Our analysis shows that for a certain range of parameters(namely, b, c, a, θ, n and k), Hopf bifurcation occurs when the ration of τ m /τ p crosses a threshold; n other words, we have sustained oscillations.For certain other parameters, we always have damped oscillations.In particular, when the degradation rate of mRNA c is too large, we can only have damped oscillations.The amplitude of such oscillations may still be very large and their decay may be slow.Thus this may be a main factor in explaining a curious observation that in most of the mutants that show a disrupted pattern of somite segmentation, the first few somites appear unaffected, but it fails to appear as the amplitude of oscillation declines to zero.Moreover, our numerical simulation shows that the period of oscillation is determined by the delay τ m and τ p , yet is not so sensitive to the protein synthesis rate and the degradation rate.Note that in our simulations, we have chosen the protein synthesis rate ranging from zero to ten molecules/min while the estimated synthesis rate is 4.5 molecules/min.
It will be interesting to investigate how synchronization can be achieved in multiple cells.All the neighboring cells are communicating through the Notch signaling mechanism.It is our future goal to extend this single cell model to a network of cells.
Finally, we would like to mention that even though the delay differential equation model predicted the period fairly well, it is not the ultimate explanation for somite formation.Other models, such as the clock-and-wavefront model, the reactiondiffusion model and the cell-cycle model, can produce results that are consistent with the observation as well.More specifically, the clock-and-wavefront model can to explain the control of somite number but conflicts with several experiments such as mirosurgical and transplantation experiments.Moreover, it does not explain the formation of anterior and posterior halves of a somite.The reaction-diffusion model agrees with the observation that one full cycle of gene oscillation corresponds to the formation of one somite.However, this model fails to explain the isolation and transplantation experiments.The cell-cycle model, on the other hand, can explain the isolation and transplantation experiments.But it can not explain oscillations of the gene expression and its pattern in the PSM.For more detail, see Schnell and Maini [22] and references therein.

Figure 1 . 2 .
Figure 1.Formation of somites in a zebrafish embryo at the 10-Somites stage: dynamic expression of deltaC mRNA.The oscillations of this gene is similar to that of her1 and her7.(Courtesy of Julian Lewis.)

Theorem 4 . 2 .
When c < The appearance of a limit cycle

Figure 4 .
Figure 4.The dependence of c on the value of Hill coefficient to achieve sustained oscillations.for c = 20 > 17.0389 even when τ m = 40 but a sustained oscillation for c = 15 with τ m = 20.