This content has been downloaded from IOPscience. Please scroll down to see the full text.
Download details:
IP Address: 140.113.38.11
This content was downloaded on 28/04/2014 at 21:22
Please note that terms and conditions apply.
Synchronized oscillations in a mathematical model of segmentation in zebrafish
View the table of contents for this issue, or go to the journal homepage for more 2012 Nonlinearity 25 869
(http://iopscience.iop.org/0951-7715/25/4/869)
Nonlinearity 25 (2012) 869–904 doi:10.1088/0951-7715/25/4/869
Synchronized oscillations in a mathematical model of
segmentation in zebrafish
Kang-Ling Liao1, Chih-Wen Shih1,3and Jui-Pin Tseng2
1Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan 2Department of Applied Mathematics, National Pingtung University of Education, Pingtung 900,
Taiwan
E-mail:[email protected]
Received 30 July 2010, in final form 29 August 2011 Published 28 February 2012
Online atstacks.iop.org/Non/25/869 Recommended by J A Sherratt
Abstract
Somitogenesis is a process for the development of somites which are transient, segmental structures that lie along the anterior–posterior axis of vertebrate embryos. The pattern of somites is governed by the segmentation clock and its timing is controlled by the clock genes which undergo synchronous oscillation over adjacent cells in the posterior presomitic mesoderm (PSM). In this paper, we analyze a mathematical model which depicts the kinetics of the zebrafish segmentation clock genes subject to direct autorepression by their own products under time delay, and cell-to-cell interaction through Delta–Notch signalling. Our goal is to elucidate how synchronous oscillations are generated for the cells in the posterior PSM, and how oscillations are arrested for the cells in the anterior PSM. For this system of delayed equations, an iteration technique is employed to derive the global convergence to the synchronous equilibrium, which corresponds to the oscillation-arrested. By applying the delay Hopf bifurcation theory and the center manifold theorem, we derive the criteria for the existence of stable synchronous oscillations for the cells at the tail bud of the PSM. Our analysis provides the basic parameter ranges and delay magnitudes for stable synchronous, asynchronous oscillation and oscillation-arrested. We exhibit how synchronous oscillations are affected by the degradation rates and delays. Extended from the analytic theory, further numerical findings linked to the segmentation process are presented.
Mathematics Subject Classification: 37N25, 34K18, 92D15 (Some figures may appear in colour only in the online journal)
3 Author to whom any correspondence should be addressed.
1. Introduction
Somites are transient, segmental structures that lie along the anterior–posterior axis of vertebrate embryos. They arise from the presomitic mesoderm (PSM) repeatedly and develop into the vertebrae, rib and tail, under further differentiation. This whole process, termed somitogenesis, is controlled by the segmentation clock which comprises a multicellular genetic network of oscillators acting in the PSM to drive periodic expression of the cyclic genes [11,18,22,26,28]. Somitogenesis involves tight gene regulation in both space and time. The cyclic genes express synchronous oscillations in the tail bud of the PSM. Underlying the morphogenetic rhythm of somitogenesis, repeated waves of oscillatory gene expression propagate from posterior to anterior; the oscillations then slow down and finally arrest and cells are formed into somites. Key cycling genes for the segmentation clock have been identified, namely, Hes1 and Hes7 in mouse, hairy1 and hairy2 in chick and her1 and her7 in zebrafish [6,18,27]. In particular, the oscillation period of the cyclic genes in the tail bud of zebrafish is about 30 min, which matches the time taken to generate a somite [17].
Theoretical modelling of somite segmentation starts from a ‘clock and wavefront’ model proposed by Cooke and Zeeman [5]. The interface between the gene expression oscillation and the oscillation-arrested is called ‘wavefront’ which determines where the somites form. It has been discovered in experiments that the signalling molecule fibroblast growth factor (FGF) regulates the differentiation of PSM cells [9,19,38]. FGF8 is only transcribed at the posterior tip of the embryo and spreads anteriorly with gradient. In addition, there exists a threshold for FGF8: if the FGF8 is smaller than the threshold, then the oscillation starts to slow down and cells will then form into somites. Based on the experimental observations of FGF8 [7,8,13,17], Baker et al [1] performed a mathematical study on the pattern formation mechanism for somites under moving gradient of FGF8 expression along the anterior–posterior axis of vertebrate embryos. The segmentation clock genes were not taken into account in that model.
While somite oscillators in chick and mouse are considered more complicated than fish, zebrafish has been a standard model organism for investigating somitogenesis. For zebrafish,
her1 and her7 are the genes that satisfy the conditions to be central components of the
somitogenesis oscillator. Modelling the kinetics of segmentation clock genes for zebrafish was proposed by Lewis who advocated that direct autorepression of her1 and/or her7 by their own products provides a mechanism for the intracellular oscillator. Accordingly, negative feedback with time delay in gene regulation was modelled to generate the oscillatory gene expression [21]. In zebrafish, the synchronization is realized by intercellular interaction via Delta–Notch signalling. These oscillatory genes, her1, her7, and the Notch ligand DeltaC all belong to the Notch signalling pathway.
Notch signalling pathway has been shown to play a critical role in somitogenesis in a number of experiments. The hypotheses of the Notch signalling pathway include oscillation-generator, boundary-formation, synchronization and polarity-formation [22,25]. For zebrafish, synchronization and polarity-formation are the only functions of Notch signalling pathway in the tail bud and the somites, respectively. This was justified in [22,25] via the mutants caused by blocking the Notch signalling pathway through adding the chemical inhibitor DAPT or overactivating the Notch signalling pathway by driving the expression of NICD (activated form of Notch). In mathematical modelling, it is assumed that NICD generated is proportional to the amount of Delta protein in the neighbouring cells. Hence, one still takes into account the expressions of her1, her7 and delta genes to investigate the function of Notch signalling pathway in zebrafish.
In this investigation, we aim at establishing analytical theories for the collective behaviour in the cell–cell kinetic model for the somitogenesis clock genes in zebrafish. We shall consider the homodimer system and take her gene as either her1 or her7; the two neighbouring cells interact through Delta–Notch signalling with delay. Let x1, x2, x3, x4 (respectively,
y1, y2, y3, y4) represent the concentrations of her mRNA, Her protein, delta mRNA and Delta protein of the first cell (respectively, the second cell), respectively. The kinetics of these genes evolve according to the following equations:
˙ x1(t )= gH(x2(t− τ1), y4(t− τ1))− d1x1(t ) ˙ x2(t )= a2x1(t− τ2)− d2x2(t ) ˙ x3(t )= gD(x2(t− τ3))− d3x3(t ) ˙ x4(t )= a4x3(t− τ4)− d4x4(t ) ˙ y1(t )= gH(y2(t− τ1), x4(t− τ1))− d1y1(t ) ˙ y2(t )= a2y1(t− τ2)− d2y2(t ) ˙ y3(t )= gD(y2(t− τ3))− d3y3(t ) ˙ y4(t )= a4y3(t− τ4)− d4y4(t ). (1.1)
Herein, the first equation describes the evolution of her mRNA, with its first and second terms depicting the transcription and the degradation (with rate d1), respectively. The second equation depicts the kinetics of Her protein, where the first and second terms express the translation with protein synthesis rate a2per mRNA molecule, and the degradation (with rate d2), respectively. Similar interpretations apply to delta mRNA in the third equation and Delta protein on the cell membrane in the fourth equation of (1.1). τ1, τ2, τ3, τ4are positive numbers which represent the time delays, incurred in the following processes:
τ1= τmh: her gene transcription,
τ2= τph: her gene translation,
τ3= τmd : delta gene transcription,
τ4= τpd : delta gene translation and delivery to cell membrane.
The functions gH and gD relate the transcription initiation rates to her and delta protein concentrations, respectively, gH(u, v)= kH 1 + bPv D0 + ν 1 + bPv D0 + ν + u2 P2 0 , for u, v 0, (1.2) gD(u)= kD 1 +Pu22 0 , for u 0, (1.3)
where kH(respectively, kD) is the maximal synthesis rate of her (respectively, delta) mRNA;
P0(respectively, PD0)is the critical number of molecules of Her (respectively, Delta) protein
per cell. In addition, b = 1 corresponds to the normal functioning of Notch signalling pathway; b = 0 represents the blockade of Notch signalling pathway (by adding DAPT);
ν > 0 means adding NICD to over-express the Notch signalling pathway. Since we only consider normal segmentation, we take b= 1 and ν = 0 throughout this paper. Note that the Delta protein of the neighbouring cell stimulates the expression of her mRNA, as indicated in function gH.
System (1.1) describes the kinetics of the gene expressions of two cells in contact; the locations of these two cells are regarded as fixed. The first four and last four equations are
coupled in a symmetric manner through the first and the fifth equations. We are interested in seeing the collective behaviour of the coupled-cell system (1.1) corresponding to various combinations of parameters and delay magnitudes. In particular, we shall investigate how synchronous oscillation and oscillation-arrested are achieved in this system. By synchronous
periodic solution, we mean a periodic solution of (1.1) with xi(t )= yi(t ), i= 1, 2, 3, 4. If a synchronous periodic orbit is stable, then we say that the system admits (stable) synchronous
oscillation. We will also address asynchronous or anti-phase oscillation which is related to
mutation of segmentation. The parameter and delay magnitudes which yield synchronous oscillations and clock-wave formation are important in understanding the models and the defects of somitogenesis. As pointed out in [34], it is not obvious how synchronized oscillations are achieved via Delta–Notch signalling, as Delta proteins stimulate the expression of the her mRNA in neighbouring cells via Notch receptors, but the active expression of her in these cells suppresses the delta gene within them. Hence, how synchronization can be achieved, while sustaining the oscillation, requires a careful and detailed analysis in the study of somitogenesis models.
Synchronized oscillations for the her1/her7 heterodimer oscillator model and the homodimer oscillator based on her1 or her7 (system (1.1)) with a different transcription function gH have been investigated numerically in [21]. It was indicated therein that the behaviour for her1/her7 heterodimer oscillator model is in general similar to that of a homodimer oscillator based on her1 or her7 alone. On the other hand, function gHin (1.2) was adopted in [25] to describe the mutant caused by the absence or over-expression of Notch signalling. This modified form of gHhas the advantage that it can be reduced to the decoupled (single-cell) model when b = 0. However, if the cells are coupled (i.e. b = 0), then the dynamics for the system with either gHare quite similar. Since our goal is to see how the Notch signalling promotes the synchronization, we only consider normal somite formation, and hence, the Notch signalling should be present. While the main focus is on system (1.1) in this work, our approach can be extended to treat the heterodimer model and systems with general transcription functions gHand gD.
Mathematical analysis truly confirms the dynamics of the considered systems, and hence provides a solid ground for comparison among various models and further extensions. The stability of nontrivial steady state and periodic orbit through the Hopf bifurcation theory for Lewis’s single-cell model has been studied in [10]. Bifurcation analysis for the oscillatory gene expressions of other single-cell models was reported in [36,37]. Analytic results on how synchronous oscillations between adjacent cells are achieved have been lacking, to the best of our knowledge. How delays affect the collective behaviour and oscillation in biological clocks has been an interesting issue to be tackled [16]. Yet, nonlinear systems with several components and multiple delays make mathematical analysis a difficult task [3,31]. An effective way to analyze the characteristic equation which contains exponential functions, is to allocate the purely imaginary eigenvalues and compute the exchange of stability for the equilibrium as the first step. In applying the delay Hopf bifurcation theory, the computation in determining the direction and stability of bifurcated periodic orbit through the center manifold theory and the normal form theory is rather involved [36]. On the other hand, while stable synchronous periodic orbit is local dynamics, ruling out possible oscillation as an indication of oscillation-arrested requires a global result. Our analytic approach in concluding the regimes of the synchronous oscillation and oscillation-arrested is especially illuminating for biological systems with many parameters and multiple delays. Our investigation sheds light on certain important issues; for example, that large coupling strength benefits the synchronous oscillation, asserted in [28], actually also depends on appropriate magnitudes of degradation rates and delays.
The binding and dissociation of a gene protein to and from its site on DNA are stochastic processes. The deterministic system (1.1) is, in fact, modelled under the assumption that the random flickering behaviour of her1 and her7 can be replaced by their instantaneous time average. It was observed in [21] that the noisy system oscillates in the same way as the deterministic system, if the protein synthesis is at the full normal rate; therein, further numerical observation on the oscillations for the noisy system when the protein synthesis is attenuated was also reported. The effects from noise on segmentation have also been studied numerically in [2,24,28,33,35].
There were some extensions of Lewis’s model. Her13.2 has been discovered to interact with clock proteins and control the rate of oscillation in zebrafish [19,32]. Cinquin [4] proposed a multicellular delay system for zebrafish somitogenesis that involves heterodimerization of clock proteins Her1 and Her7, with protein Her13.2, to model the posterior-to-anterior slowing of oscillation rate, which leads to formation of clock-wave. Campanelli and Gedeon [2] generalized Lewis’s model to a system with two different genetic control mechanisms to initiate the gene-expression wave, one on the number of clock protein transcription binding sites, the other on differential decay rates for clock protein monomers and dimers.
In the somitogenesis, time delays which have been estimated in the range of tens of minutes in cell culture are due to synthesis and trafficking of macromolecules in cells. Modelling with delays raises some mathematical technicality in differential equations. Replacing the time delay by taking into account certain intermediate process in the cell, namely, the translocation of Her protein from cytoplasm to nucleus, an ODE model has been studied in [33,34]. Therein, Michaelis–Menten type reaction for the degradation and more general transcription function with Hill coefficients were employed. The investigation was mainly based on numerical computations. In fact, in that ODE system, as some of those over ten parameters varies, the equilibrium value, the linearization and its eigenvalues all change, which make bifurcation analysis a difficult task.
Another modelling of vertebrate segmentation is to depict the phase dynamics of coupled oscillators, where each of the oscillators represents a cell or a group of synchronous cells in the PSM. It was shown that disruption of Delta–Notch intercellular coupling increases the period of zebrafish somitogenesis and the embryonic segment length. An instability resulting from decreased coupling delay time was predicted by the theory of phase oscillator and justified by the experiments [16,24]. A phase equation depicting the peak gene expression in the travelling wave for a lattice system of ODEs was discussed in [33]. On the one hand, the parameters in the kinetic models are more difficult to measure than those for the phase oscillators in vivo [12,16,34]. On the other hand, those phase models are built under certain assumptions; how they are linked to the solution behaviour of the kinetic equations which depict the interaction of cyclic genes, such as (1.1), remains to be investigated. This present work hopes to provide a basis for this linkage. Our analytic theory not only delineates qualitatively how the degradation rates, the other focal parameters and delays affect synchronous oscillations, but also provides the parameter regimes for oscillation-arrested, synchronous and asynchronous oscillations.
For mouse, the oscillating network of signalling genes underlying the segmentation clock is more complex. It was reported through experiments by Dequ´eant and coworkers [6] that the inhibitors of the Notch, FGF and Wnt signalling pathways are involved in the segmentation clock. In addition, numerical simulation on an ODE model for Notch, FGF and Wnt signalling pathways which contain the inhibitors discovered in [6], was performed to investigate the oscillation and the interactions among these three pathways in [13]. On the other hand, delay models for gene regulation of Axin2, Hes1 and Lfng, which belong to the Notch and Wnt signalling pathway, were employed to investigate the segmentation clock in mouse and chick in [29].
The paper is organized as follows. In section2, we study the positive invariance of the considered domain and the dissipative property for system (1.1), for the validity of the model equations. We then show the existence of a synchronous equilibrium and apply an iteration technique developed in [30,31] to derive the global convergence to the synchronous equilibrium. We regard this regime of parameters as the one for oscillation-arrested which corresponds biologically to the formed somites. Although there are other proteins involved in slowing down the oscillation, we interpret that these effects are reflected in altering the parameters, such as degradation rates, of the system. In section3, we employ the delay Hopf bifurcation theorem to prove the existence of synchronous periodic solutions and use the center manifold theorem and the normal form method to analyze the stability of the bifurcated periodic solutions. In section4, we summarize the collective behaviour of the coupled system (1.1), and give some numerical simulations to demonstrate our analytic results. In section5, the paper ends with a conclusion.
2. Global convergence to steady state
In this section, we first examine the basic properties for system (1.1) to ensure that it is suitable to model gene regulation. We then show that (1.1) has a synchronous equilibrium ¯X. A criterion for global convergence to ¯X is then derived.
For the coupled system (1.1) with four time delays τ1, τ2, τ3, τ4, in this paper, we consider the evolution (t, φ) of (1.1) from initial condition φ = (φ1, . . . , φ8)∈ C([−τM,0],R8+)at initial time t0= 0, where
τM := max {τ1, τ2, τ3, τ4}, R8+:= {(x1, . . . , x4, y1, . . . , y4)| xi, yi 0, i = 1, . . . , 4}. Let X(t; φ) denote the solution induced from (1.1), which is defined by X(t + θ; φ) =
(t, φ)(θ ), θ ∈ [−τM,0], for t > 0. We further denoteX(t) = (x(t), y(t)) = X(t; φ) =
(x(t; φ), y(t; φ)), where x = (x1, x2, x3, x4), y = (y1, y2, y3, y4), when φ is not specified. To ensure that (1.1) is proper in modelling the gene regulation, one first needs to ensure that the solutionX(t; φ) is always nonnegative, for any φ ∈ C([−τM,0],R8+).
Proposition 2.1. C([−τM,0],R8+) is positively invariant for system (1.1).
Proof. We shall show that xi(t ) 0, yi(t ) 0, i = 1, 2, 3, 4, for solution X(t) evolved from arbitrary φ ∈ C([−τM,0],R8+). Let I be the maximal interval of existence for solutionX(t). First,
˙
x3(t )= −d3x3(t )+ gD(x2(t− τ3)) >−d3x3(t ), (2.1) for all t ∈ I, since gD(u)= kD/(1 + u2/P02) >0, for any u. It follows that x3(t ) 0, for all
t ∈ I, through comparison arguments. Similarly, we can show that y3(t ) 0, for all t ∈ I, by the symmetry of system (1.1). With x3(t ) 0 for all t ∈ I, we use similar argument to derive
x4(t ), y4(t ), x1(t ), y1(t ), x2(t ), y2(t ) 0, for all t ∈ I, successively. We next derive the existence of global attracting set for system (1.1).
Proposition 2.2. There exists a closed and bounded setQ= 4
i=1Qi× 4i=1Qi ⊂ R8+, such
thatX(t; φ) converges to Q for arbitrary φ ∈ C([−τM,0],R8+), where Qiare defined by
Q1= [ ˇq1,ˆq1] := [0, kH/d1], Q2= [ ˇq2,ˆq2] := [0, a2kH/(d1d2)], Q3= [ ˇq3,ˆq3] := [ˇc3/d3, kD/d3], Q4= [ ˇq4,ˆq4] := [a4ˇc3/(d3d4), a4kD/(d3d4)] and ˇc3:= kDP02 P2 0 + ( a2kH d1d2) 2. (2.2)
Proof. The assertion follows from estimating dxi/dt, dyi/dt of system (1.1), and hence
xi(t ), yi(t ), successively. From the definition of functions gHand gD, we have 0 < gH(u, v) kH, 0 < gD(u) kD, for all u, v 0.
Let X(t; φ) = X(t) = (x1(t ), . . . , x4(t ), y1(t ), . . . , y4(t )) be the solution evolved from
φ∈ C([−τM,0],R8+). Note that, from proposition2.1, xi(t ), yi(t ) 0, for all i = 1, 2, 3, 4,
t ∈ I. Hence,
−d1x1(t ) <x˙1(t ) −d1x1(t )+ kH, for all t ∈ I.
Consequently, x1(t ) exists on [0,∞) and converges to [0, kH/d1] =: ˜Q1, as t → ∞. Subsequently, x2(t )also exists on [0,∞), with respect to the second component of system (1.1); in addition, for any > 0, there exists some t
1 >0 such that
−d2x2(t ) ˙x2(t ) −d2x2(t )+ a2kH/d1+ , for t t1;
thus x2(t ) converges to [0, a2kH/(d1d2)+ /d2] for all > 0, and hence converges to [0, a2kH/(d1d2)]=: ˜Q2, as t→ ∞. By similar arguments, we can justify that x3(t )and x4(t ) exist on [0,∞) and converge to ˜Q3:= [0, kD/d3] and ˜Q4 := [0, a4kD/(d3d4)], respectively. Moreover, we can also prove that yi(t )exists on [0,∞) and converges to set ˜Qi, i= 1, 2, 3, 4, by the symmetry of system (1.1).
So far, we have constructed a region, 4i=1Q˜i × 4i=1Q˜i, to which every solutionX(t) converges. We can continue to construct an attracting region smaller than 4i=1Q˜i× 4i=1Q˜i, for system (1.1). Note that the convergence of x2(t )to ˜Q2= [0, a2kH/(d1d2)] leads to that for any > 0, there exists some t2> t1such that gD(x2(t )) > kDP02/[P02+(a2kH/(d1d2))2]− =:
ˇc3− , for all t t2− τM >0. It follows from the third component of system (1.1) that −d3x3(t )+ ˇc3− < ˙x3(t ) −d3x3(t )+ kD, for all t t2.
Subsequently, x3(t )converges to [ˇc3/d3− /d3, kD/d3] for all > 0, and hence converges to [ˇc3/d3, kD/d3]=: Q3 = [ ˇq3,ˆq3], as t → ∞. By repeating the similar process, we can show that x4(t )converges to interval Q4 = [ ˇq4,ˆq4] := [a4ˇc3/(d3d4), a4kD/(d3d4)], as t → ∞. Setting Q1 = [ ˇq1,ˆq1] := [0, kH/d1] and Q2 = [ ˇq2,ˆq2] := [0, (a2kH)/(d1d2)], we also conclude that xi(t )converges to set Qi, for i = 1, 2. Symmetrically, yi(t )converges to set
Qi, i= 1, 2, 3, 4.
Remark 2.1. (i) By continuing the estimates in the proof of proposition2.2iteratively, we can confine the asymptotical behaviour of system (1.1) to smaller attracting set. For simplicity of presentation, we merely demonstrate the idea via two iteration steps. (ii) The assertion of proposition2.2indicates that every solution of (1.1) exists on [0,∞). (iii) Proposition2.1can be extended to conclude that the Delta protein is always positive, if its initial value is positive. Moreover, according to proposition2.2, the Delta protein becomes larger than a positive value after certain time. This indicates that the Notch signalling is always present, if b= 0.
Next, we shall show that system (1.1) always admits one positive synchronous equilibrium. Herein, an equilibrium (x1∗, . . . , x4∗, y1∗, . . . , y4∗)of (1.1) is synchronous if xi∗ = y∗i, for
i= 1, 2, 3, 4.
Proposition 2.3. There exists a synchronous equilibrium point
for system (1.1), where ¯x1= d2¯x2 a2 , ¯x3= kDP 2 0 d3(P02+¯x22) , ¯x4 = a4kDP 2 0 d3d4(P02+ ¯x22) , and ¯x2is a positive solution to the equation:
p(ξ ):= c5ξ5+ c3ξ3− c2ξ2+ c1ξ = a2kHP04(d3d4PD0+ a4kD), (2.3) and c5 = d1d2d3d4PD0 > 0, c3 = 2d1d2d3d4P 2 0PD0 > 0, c2 = a2d3d4kHP 2 0PD0 > 0, c1= d1d2P04(d3d4PD0+ a4kD) >0.
Proof. Finding the synchronous equilibrium for (1.1) amounts to solving
gH(x2, x4)− d1x1 = 0,
a2x1− d2x2= 0,
gD(x2)− d3x3 = 0,
a4x3− d4x4= 0.
(2.4)
Then (¯x1,¯x2,¯x3,¯x4)satisfies (2.4) if and only if ¯x1= d2¯x2 a2 , ¯x3= kDP02 d3(P02+¯x22) , ¯x4 = a4kDP02 d3d4(P02+ ¯x22) ,
and¯x2satisfies (2.3). Since p(ξ )→ ∞ as ξ → ∞ and p(0) = 0 < a2kHP04(d3d4PD0+ a4kD),
there exists a solution ξ0>0 satisfying p(ξ0)= a2kHP04(d3d4PD0+ a4kD).
Note that every component of ¯X is positive. In the rest of this section, we shall discuss global convergence and synchronization. A criterion will be derived to yield both global synchronization and global convergence to the synchronous equilibrium point ¯X for system (1.1). To this end, we let ˜x(t) = x(t) − ¯x, ˜y(t) = y(t) − ¯x, and still denote ˜x, ˜y by x, y respectively. System (1.1) becomes
˙ x1(t )= −d1x1(t )+ gH(x2(t− τ1)+ ¯x2, y4(t− τ1)+ ¯x4)− d1¯x1, ˙ x2(t )= −d2x2(t )+ a2x1(t− τ2), ˙ x3(t )= −d3x3(t )+ gD(x2(t− τ3)+ ¯x2)− d3¯x3, ˙ x4(t )= −d4x4(t )+ a4x3(t− τ4), ˙ y1(t )= −d1y1(t )+ gH(y2(t− τ1)+ ¯x2, x4(t− τ1)+ ¯x4)− d1¯x1, ˙ y2(t )= −d2y2(t )+ a2y1(t− τ2), ˙ y3(t )= −d3y3(t )+ gD(y2(t− τ3)+ ¯x2)− d3¯x3, ˙ y4(t )= −d4y4(t )+ a4y3(t− τ4). (2.5)
We shall show that every solution converges to the origin for system (2.5). Consider an arbitrary solution X(t) = (x1(t ), . . . , x4(t ), y1(t ), . . . , y4(t )) of (2.5); thenX(t) exists on [0,∞) according to remark2.1(ii). As xi(t ), yi(t )satisfy (2.5), by the mean value theorem, we have ˙ xi(t )= −dixi(t )+ wi(t ), ˙ yi(t )= −diyi(t )+ wi+4(t ), (2.6)
where i= 1, 2, 3, 4, and w1(t ):= ∂gH ∂u (u1(t− τ1), v1(t− τ1))x2(t− τ1)+ ∂gH ∂v (u1(t− τ1), v1(t− τ1))y4(t− τ1), w2(t ):= a2x1(t− τ2), w3(t ):= dgD du (u3(t− τ3))x2(t− τ3), w4(t ):= a4x3(t− τ4), w5(t ):= ∂gH
∂u (˜u1(t− τ1),˜v1(t− τ1))y2(t− τ1)+ ∂gH ∂v (˜u1(t− τ1),˜v1(t− τ1))x4(t− τ1), w6(t ):= a2y1(t− τ2), w7(t ):= dgD du (˜u3(t− τ3))y2(t− τ3), w8(t ):= a4y3(t− τ4),
where u1(t− τ1)(respectively, v1(t − τ1), u3(t − τ3), ˜u1(t− τ1), ˜v1(t− τ1), ˜u3(t − τ3)) is between x2(t− τ1)+ ¯x2and ¯x2(respectively, y4(t− τ1)+ ¯x4 and ¯x4, x2(t− τ3)+ ¯x2and ¯x2, y2(t− τ1)+ ¯x2 and ¯x2, x4(t − τ1)+ ¯x4 and ¯x4, y2(t− τ3)+ ¯x2 and ¯x2). Note that the solutionX(t) of system (2.5) eventually converges toQ− ¯X = 4
i=1Q∗i × 4i=1Q∗i, where
Q∗i := [ ˇqi− ¯xi,ˆqi− ¯xi], for i= 1, 2, 3, 4. Every component of (2.6) is of the form ˙
x(t )= −cx(t) + w(t), (2.7)
where c > 0 is a constant and w(t) is a continuous scalar function. For any t 0, we denote |w|max(t ):= sup{|w(s)| : s t}, |w|max(∞) := lim
t→∞|w| max(t ). It is not difficult to derive the following property; cf [30].
Lemma 2.4. Every solution of (2.7) converges to an interval [−˜δ, ˜δ], where
0 ˜δ |w|max(∞)/c.
According to lemma2.4, there exist eight intervals Ii := [−δi, δi], i= 1, · · · , 8, to which the ith component ofX(t) converges respectively. Moreover,
0 δi |wi|max(∞)/di, for i= 1, . . . , 4, 0 δi |wi|max(∞)/di−4, for i= 5, . . . , 8.
Next, we construct a sequence of nonnegative numbers{δ(k)i }∞k=1through an iteration process to derive sharper estimates for δi, i= 1, . . . , 8. We introduce
ρ1:= max ∂gH(u, v) ∂u : u ∈ Q2, v∈ Q4 , ρ2:= max∂gH(u, v) ∂v : u ∈ Q2, v∈ Q4 , (2.8) ρ3:= max dgD(u) du : u ∈ Q2 .
Proposition 2.5. For each i = 1, . . . , 8, there exists a sequence of nonnegative numbers
{δ(k)
i }∞k=1with δ (k)
i δisuch that for each k, the ith component for the solutionX(t) of system
(2.5) converges to Ii(k):= [−δ (k) i , δ (k) i ], as t→ ∞, and δ (k) i satisfies 0 δ1(k)= δ5(k):= (ρ1δ2(k−1)+ ρ2δ4(k−1))/d1, 0 δ2(k)= δ6(k):= (a2/d2)δ (k) 1 , 0 δ3(k)= δ7(k):= (ρ3/d3)δ (k) 2 , 0 δ4(k)= δ8(k):= (a4/d4)δ3(k), (2.9)
where δ2(0):= max{| ˇq2− ¯x2|, | ˆq2− ¯x2|} and δ4(0):= max{| ˇq4− ¯x4|, | ˆq4− ¯x4|}, k 1, and ρi
are defined in (2.8).
Proof. From lemma 2.4, x1(t ) and y1(t )converge to intervals [−δ1, δ1] and [−δ5, δ5], as
t → ∞, respectively, where
0 δi |wi|max(∞)/d1, i= 1, 5.
Note that both x2(t )and y2(t )eventually converge to Q∗2 = [ ˇq2 − ¯x2,ˆq2− ¯x2], and both
x4(t )and y4(t )eventually converge to Q∗4 = [ ˇq4 − ¯x4,ˆq4 − ¯x4]; therefore u1(t− τ1)and ˜u1(t− τ1)eventually converge to Q2and v1(t− τ1)and˜v1(t− τ1)eventually converge to Q4. Accordingly, for i= 1, 5,
|wi|max(∞) ρ1max{| ˇq2− ¯x2|, | ˆq2− ¯x2|} + ρ2max{| ˇq4− ¯x4|, | ˆq4− ¯x4|} =: ρ1δ(0)2 + ρ2δ(0)4 .
We may say that x1(t ) and y1(t ) converge to the closed and bounded intervals I1(1) := [−δ(1)
1 , δ (1)
1 ] ⊃ I1 and I5(1) := [−δ5(1), δ(1)5 ] ⊃ I5, respectively, where δ1(1) = δ (1)
5 :=
(ρ1δ2(0)+ ρ2δ4(0))/d1.It then follows that for i = 2, 6,
δi |wi|max(∞)/d2 (a2/d2)δ1(1),
since|w2(t )| = |a2x1(t − τ2)|, |w6(t )| = |a2y1(t − τ2)|, and thus |wi|max(∞) a2δ(1)1 ,
i = 2, 6. We then conclude that x2(t )and y2(t )converge to closed and bounded intervals
I2(1) := [−δ2(1), δ2(1)]⊃ I2and I6(1) := [−δ(1)6 , δ6(1)] ⊃ I6, respectively, where δ2(1) = δ(1)6 :=
(a2/d2)δ1(1). By continuing this process, we obtain (2.9).
In the following theorem, with the help of proposition2.5, we derive the condition for δi(k) to converge to zero as k→ ∞, which yields δi = 0, for i = 1, . . . , 8.
Theorem 2.6. System (1.1) achieves global convergence to the synchronous equilibrium point ¯X; that is, xi(t ), yi(t ) → ¯xi, as t → ∞, i = 1, 2, 3, 4, for solution
(x1(t ), . . . , x4(t ), y1(t ), . . . , y4(t )) of system (1.1), evolved from any initial condition in C([−τM,0],R8+), if
ρ1a2d3d4+ ρ2ρ3a2a4 < d1d2d3d4=: β4. (2.10)
Proof. To justify the assertion, it suffices to show that δ(k)i converges to zero as k→ ∞ for all
i= 1, · · · , 8. From the iterative estimate in (2.9), we obtain
δ1(k)= (ρ1δ (k−1) 2 + ρ2δ (k−1) 4 )/d1 = ρ1a2 d1d2 δ1(k−1)+ ρ2a4ρ3a2 d1d4d3d2 δ1(k−1) = ρ1a2d3d4+ ρ2ρ3a2a4 d1d2d3d4 δ1(k−1) =: Rδ(k−1) 1 .
Thus δ1(k)→ 0, as k → ∞, since 0 < R < 1, by (2.10). Subsequently, δ(k)i → 0, as k → ∞,
Remark 2.2. (i) In fact, we can compute the upper bounds for the derivatives of gHand gD: 0∂gH(u, v) ∂u 2kHP02PD0ˆq2(PD0+ ˆq4) P04(PD0+ ˇq4)2 =: ˆρ1, for u∈ Q2, v∈ Q4, 0∂gH(u, v) ∂v kHP02PD0ˆq 2 2 P04(PD0+ ˇq4)2 =: ˆρ2, for u∈ Q2, v∈ Q4, (2.11) 0dgD(u) du 2a2kDkH P02d1d2 =: ˆρ3, for u∈ Q2.
Note that ˆρi is computable and ˆρi ρi, i = 1, 2, 3. Therefore, the convergent result in theorem2.6also holds under the following condition
ˆρ1a2d3d4+ ˆρ2ˆρ3a2a4 < β4, (2.12) which is easier to verify from computation view point. (ii) Roughly speaking, condition (2.10) in theorem 2.6 favors large decay rates and small synthesis rates. (iii) Recalling remark 2.1(i), we can establish smaller attracting region Q if further iterative estimate is performed. Subsequently, the value ρi (respectively, ˆρi) can be lowered to relax condition (2.10) (respectively, (2.12)). (iv) Under condition (2.10) or (2.12), ¯X obtained in proposition2.3 is the unique positive equilibrium for system (1.1), according to theorem2.6.
Theorem2.6asserts that
xi(t )→ ¯xi and yi(t )→ ¯xi, as t→ ∞, i = 1, 2, 3, 4,
for solution X(t) = (x1(t ), . . . , x4(t ), y1(t ), . . . , y4(t ))of system (1.1), evolved from any initial condition inC([−τM,0],R8+). This corresponds to non-oscillatory or oscillation-arrested phase for system (1.1), which is associated with the state of formed somites. This scenario also corresponds to the global synchronization of system (1.1), as
xi(t )− yi(t )→ ¯xi− ¯xi = 0, as t → ∞, i = 1, 2, 3, 4.
In fact, by analyzing the difference systemdtd(xi(t )− yi(t )), and employing arguments similar to those in proposition2.5, we can establish the global synchronization for system (1.1) directly, also under condition (2.10).
Corollary 2.7. Under condition (2.10) or (2.12), system (1.1) attains global synchronization: xi(t )− yi(t )→ 0, as t → ∞, i = 1, 2, 3, 4, for solution (x1(t ), . . . , x4(t ), y1(t ), . . . , y4(t ))
of system (1.1), evolved from any initial condition inC([−τM,0],R8+).
It is usually a challenging task to investigate global dynamics for nonlinear systems and delay coupled systems. Without knowing and using Lyapunov function, the above iteration formulation works well in deriving global synchronization and global convergence to the synchronous equilibrium.
3. Synchronous oscillations
Although corollary2.7concludes the global synchronization, it does reduce to the situation that every solution tends to the synchronous steady state. Thus, there does not exist any oscillation for system (1.1), under condition (2.10) or (2.12). In this section, we turn to target the local dynamics and employ bifurcation theory to investigate the synchronous oscillations.
We shall apply the delay Hopf bifurcation theory to prove the existence of nontrivial synchronous periodic solutions for system (1.1), in section 3.1. We then use the center manifold theorem and the normal form method to compute the stability of the bifurcated
periodic solutions in section3.2. With positive invariance of the synchronous manifold, the result then leads to synchronous oscillations for (1.1). To depict the loss of synchrony and the phase transition as parameters or delay magnitudes vary, the asynchronous oscillation will also be discussed.
3.1. Bifurcation of periodic orbits
In this section, we shall investigate the periodic solutions bifurcated from the synchronous equilibrium ¯X of system (1.1) via Hopf bifurcation theorem. By observing the structure of the characteristic equation, while holding s= τ3+ τ4fixed, we take r= τ1+ τ2as the bifurcation parameter. We will work on system (2.5) which is a translation of (1.1) from ¯X to the origin. The linearization of system (2.5) at the origin is given by
˙ x1(t )= −d1x1(t )+ 11x2(t− τ1)+ 12y4(t− τ1), ˙ x2(t )= −d2x2(t )+ a2x1(t− τ2), ˙ x3(t )= −d3x3(t )+ 31x2(t− τ3), ˙ x4(t )= −d4x4(t )+ a4x3(t− τ4), ˙ y1(t )= −d1y1(t )+ 51y2(t− τ1)+ 52x4(t− τ1), ˙ y2(t )= −d2y2(t )+ a2y1(t− τ2), ˙ y3(t )= −d3y3(t )+ 71y2(t− τ3), ˙ y4(t )= −d4y4(t )+ a4y3(t− τ4), (3.1) where 11= 51:= ∂gH ∂u (¯x2,¯x4), 12= 52:= ∂gH ∂v (¯x2,¯x4), 31 = 71:= dgD du(¯x2). (3.2) For convenience, we set
γ1:= − ∂gH ∂u (u, v)|u= ¯x2,v= ¯x4= 2kHP02PD0¯x2(PD0+ ¯x4) [PD0(P 2 0 + ¯x22)+ P02¯x4]2 , γ2:= ∂gH ∂v (u, v)|u= ¯x2,v= ¯x4= kHP02PD0¯x 2 2 [PD0(P 2 0 + ¯x22)+ P02¯x4]2 , γ3:= − dgD du (u)|u= ¯x2 = 2kDP02¯x2 (P02+ ¯x22)2.
Indeed, γ1 = − 11 = − 51, γ2 = 12 = 52, γ3 = − 31 = − 71 and γ1, γ2, γ3 > 0. The characteristic equation for (3.1) is
(λ, τ1, τ2, τ3, τ4)= 0, (3.3) where (λ, τ1, τ2, τ3, τ4)is given by det λ+ d1 γ1e−τ1λ 0 0 0 0 0 −γ2e−τ1λ −a2e−τ2λ λ+ d2 0 0 0 0 0 0 0 γ3e−τ3λ λ+ d3 0 0 0 0 0 0 0 −a4e−τ4λ λ+ d4 0 0 0 0 0 0 0 −γ2e−τ1λ λ+ d1 γ1e−τ1λ 0 0 0 0 0 0 −a2e−τ2λ λ+ d2 0 0 0 0 0 0 0 γ3e−τ3λ λ+ d1 0 0 0 0 0 0 0 −a4e−τ4λ λ+ d4 .
By letting r := τ1+ τ2and s := τ3+ τ4, the characteristic equation (3.3) can be factored as +(λ, τ1, τ2, τ3, τ4)· −(λ, τ1, τ2, τ3, τ4)= +(λ, r, s)· −(λ, r, s)= 0, (3.4) where ±(λ, τ1, τ2, τ3, τ4):= λ4+ β1λ3+ β2λ2+ β3λ+ β4 + a2γ1[λ2+ (d3+ d4)λ+ d3d4]e−(τ1+τ2)λ± a2a4γ2γ3e−(τ1+τ2+τ3+τ4)λ, ±(λ, r, s):= λ4+ β1λ3+ β2λ2+ β3λ+ β4 + a2γ1[λ2+ (d3+ d4)λ+ d3d4]e−rλ± a2a4γ2γ3e−(r+s)λ, (3.5) β1:= d1+ d2+ d3+ d4, β2:= d2d4+ d3d4+ d2d3+ d1d4+ d1d2+ d1d3, β3:= d2d3d4+ d1d2d4+ d1d3d4+ d1d2d3,
and β4has been introduced in (2.10). Note that each βiis positive.
We now investigate, for fixed s 0, the existence of r-induced periodic solutions bifurcated from the origin of (2.5). We shall start from finding a purely imaginary root iw of the characteristic equation (3.3) and its corresponding bifurcation value r = r(w) so that
(iw, r(w), s)= 0. To this end, we substitute λ = iw, with w > 0, into ±(λ, r, s)= 0 and
collect the real and imaginary parts:
w4− β2w2+ β4+ a2γ1(d3d4− w2)cos (rw) + a2γ1w(d3+ d4)sin (rw) ±a2a4γ2γ3cos ((r + s)w)= 0,
−β1w3+ β3w+ a2γ1w(d3+ d4)cos (rw) + a2γ1(w2− d3d4)sin (rw)
∓a2a4γ2γ3sin ((r + s)w)= 0. (3.6)
By the properties of trigonometric functions, (3.6) can be written as
w4− β2w2+ β4+ [a2γ1(d3d4− w2)± a2a4γ2γ3cos (sw)] cos (rw) +[a2γ1w(d3+ d4)∓ a2a4γ2γ3sin (sw)] sin (rw)= 0, −β1w3+ β3w+ [a2γ1w(d3+ d4)∓ a2a4γ2γ3sin (sw)] cos (rw) − [a2γ1(d3d4− w2)± a2a4γ2γ3cos (sw)] sin (rw)= 0, or, equivalently, L±(w)· sin (φ±+ rw)= −w4+ β2w2− β4, L±(w)· cos (φ±+ rw)= β1w3− β3w, (3.7) where L±(w) := [a2γ1(d3d4 − w2) ± a 2a4γ2γ3cos (sw)]2 + [a2γ1w(d3 + d4) ∓ a2a4γ2γ3sin (sw)]2>0, if w is a solution of (3.7); φ ±∈ [0, 2π) and satisfies sin (φ±)= [a2γ1(d3d4− w2)± a2a4γ2γ3cos (sw)]/ √ L±(w), cos (φ±)= [a2γ1w(d3+ d4)∓ a2a4γ2γ3sin (sw)]/ √ L±(w).
Summing up the square of equations (3.7) gives
Q±(w)= (a2d3d4γ1)2=: , (3.8) where Q±(w)= w8+ (β12− 2β2)w6+ (−a22γ12+ β22− 2β1β3+ 2β4)w4 + [−a22(d32+ d42)γ12+ β32− 2β2β4± 2a22a4γ1γ2γ3cos (sw)]w2 ± [2a2 2a4(d3+ d4)γ1γ2γ3sin (sw)]w + [β42− a22a42γ22γ32∓ 2a22a4d3d4γ1γ2γ3cos (sw)].
Figure 1. Graphs of functions z= Q−(w), z= Q+(w), z= , with the parameter data in [25].
We summarize the properties of functions Q±(w)as follows; cf figure1. (P1) : Q±(0)= β2
4− (a2a4γ2γ3)2∓ 2a22a4d3d4γ1γ2γ3. Obviously,
Q+(0) < Q−(0).
(P2) : Q+(w)and Q−(w)are dominated by the leading term w8, since| sin u|, | cos u| 1. Consequently, both Q+(w)and Q−(w)are strictly increasing eventually, and increase to ∞, as w → ∞.
According to properties (P1) and (P2), we derive the following lemma.
Lemma 3.1.
(i) For any fixed s 0, there exists at least a positive solution to Q+(w)= if Q+(0) < ;
namely,
0 < β4 < a2a4γ2γ3+ a2d3d4γ1. (3.9)
(ii) For any fixed s 0, each of Q−(w)= and Q+(w)= has at least a positive solution
if Q−(0) < ; namely,
0 < β4 <|a2a4γ2γ3− a2d3d4γ1|. (3.10)
There may exist multiple solutions to Q+(w)= or Q−(w)= . In the following, we denote by w+(respectively, w−) a positive solution to Q+(·) = (respectively, Q−(·) = ). Now, we find the value of r such that iw±is a purely imaginary root of ±(·, r, s) = 0. We
divide the first equation of (3.7) by the second and obtain tan (φ±+ rw)= S(w)/C(w),
S(w):= −w4+ β2w2− β4,
C(w):= β1w3− β3w.
Let σ = + or −. For a fixed wσ, there exists a sequence{rσ(k)(wσ)}k∈Z:
rσ(k)(wσ):= 1 wσ tan−1 S(wσ) C(wσ) − φσ+ 2kπ , if C(wσ) >0 , 1 wσ tan−1 S(wσ) C(wσ) − φσ+ (2k− 1)π , if C(wσ) <0 , 1 wσ 3π 2 − φσ+ 2kπ , if C(wσ)= 0, S(wσ) <0, 1 wσ π 2 − φσ+ 2kπ , if C(wσ)= 0, S(wσ) >0, (3.11)
such that σ(iwσ, rσ(k)(wσ), s) = 0. To simplify the notation, we denote rσ(k) := rσ(k)(wσ) as the bifurcation value, and we shall consider the case r(k)
σ > 0. Accordingly, a positive solution w+ (respectively, w−) of Q+(·) = (respectively, Q−(·) = ) corresponds to a pair of purely imaginary roots±iw+(respectively,±iw−) of +(·, r+(k), s)= 0 (respectively,
−(·, r−(k), s)= 0).
Next, to apply Hopf bifurcation theory, we further impose the conditions of simple root and transversality. For σ = +, −, we consider
Condition (C1)σ: Qσ(wσ)= 0, and all other positive solutions to Q+(·) = and Q−(·) = are not integer multiples of wσ.
Condition (C2)σ: [Rσ(wσ, rσ(k))]2+ [Iσ(wσ, rσ(k))]2= 0, W1,σ2 (wσ, rσ(k))+ W2,σ2 (wσ, rσ(k))= 0, and Q1,σ(wσ, rσ(k))W1,σ(wσ, rσ(k))+ Q2,σ(wσ, rσ(k))W2,σ(wσ, rσ(k))= 0, where
R±(w, r):= −3w2β1+ β3+ a2γ1(d3+ d4− d3d4r+ w2r)cos (wr)
−a2γ1w(−2 + d3r+ d4r)sin (wr)∓ a2a4γ2γ3(r+ s) cos ((r + s)w),
I±(w, r):= −4w3+ 2wβ2− a2γ1w(−2 + d3r+ d4r)cos (wr)
−a2γ1(d3+ d4− d3d4r+ w2r)sin (wr)± a2a4γ2γ3(r+ s) sin ((r + s)w),
Q1,±(w, r):= −a2w[(d3+ d4)γ1wcos (wr) + γ1(−d3d4+ w2)sin (wr) ∓a4γ2γ3sin ((r + s)w)],
Q2,±(w, r):= a2w[γ1(d3d4− w2)cos (wr)± a4γ2γ3cos ((r + s)w) +(d3+ d4)γ1wsin (wr)],
W1,±(w, r):= −3w2β1+ β3+ a2γ1(d3+ d4− d3d4r+ w2r)cos (wr)
∓a2a4γ2γ3(r+ s) cos ((r + s)w)− a2γ1w(−2 + d3r+ d4r)sin (wr),
W2,±(w, r):= −[4w3− 2wβ2+ a2γ1w(−2 + d3r+ d4r)cos (wr)
+a2γ1(d3+ d4− d3d4r+ w2r)sin (wr)∓ a2a4γ2γ3(r+ s) sin ((r + s)w)]. These terms are obtained from the proof of the following theorem.
Theorem 3.2. For a fixed s 0, assume that there exists a positive solution w+(respectively,
w−) to Q+(w) = (respectively, Q−(w) = ) satisfying conditions (C1)+ and (C2)+
(respectively, (C1)− and (C2)−) for some r+(k) >0 (respectively, r (k)
− >0) and k∈ Z. Then
Hopf bifurcation occurs at r= r+(k)(respectively, r = r (k)
− ), and a periodic orbit is bifurcated
from the zero solution of system (2.5). Moreover, the periodic orbit bifurcated at r+(k) >0 is
synchronous, while the one bifurcated at r−(k)>0 is asynchronous.
Proof. Under the assumption, it suffices to verify the transversality for the assertion of Hopf
bifurcation. We only prove the first case. We compute that
∂
∂λ+(λ, r, s)|λ=iw+,r=r+(k) = {4λ
3+ 3β
1λ2+ 2β2λ+ β3+ a2γ1(2λ + d3+ d4)e−rλ
−ra2γ1[λ2+ (d3+ d4)λ+ d3d4]e−rλ− (r + s)a2a4γ2γ3e−(r+s)λ}|λ=iw+,r=r+(k)
= R+(w+, r+(k))+ iI+(w+, r+(k)),
where R+, I+ are defined as above. Thus, ∂λ∂+(λ, r, s)|λ=iw+,r=r+(k) = 0, due to condition
(C2)+. Hence, there exist a δ (k)
+ >0 and a smooth function λ (k) + : (r (k) + − δ (k) + , r (k) + + δ (k) + )→ C such that +(λ (k) + (r (k) + ), r (k) + , s) = +(λ (k) + (r (k) + )) = 0 and λ (k) + (r (k) + ) = iw+, by the implicit function theorem. Differentiating +(λ
(k) + (r))= 0 with respect to r at r = r (k) + , we obtain (λ(k)+ )(r+(k))= Q1,+(w+, r (k) + )+ iQ2,+(w+, r+(k)) W1,+(w+, r+(k))+ iW2,+(w+, r+(k)) ,
where Q1,+, Q2,+, W1,+, W2,+are defined as above. Thus, Re((λ(k)+ )(r (k)
+ ))= 0, by condition
(C2)+. Therefore, Hopf bifurcation occurs at the origin of system (2.5) at r = r+(k), according to [15]. Next, note that
S := {(z1, z2, z3, z4, z1, z2, z3, z4): zi ∈ C([−τM,0],R+), i= 1, 2, 3, 4} (3.12) is positively invariant under (1.1) and (2.5). In addition, the characteristic equation for the linearized system of (2.5) restricted to the synchronous manifoldS at the origin, is exactly
+(λ, r, s)= 0. Therefore, by the uniqueness in Hopf bifurcation theorem and the positive invariance ofS, the periodic orbit generated by the bifurcation at r= r+(k)is synchronous. On the other hand, the periodic orbit generated by the bifurcation at r = r−(k) is asynchronous;
cf [14,20,39] for details.
Remark 3.1. (i) In lemma3.1, (3.9) and (3.10) are sufficient conditions for Q+(w)= and
Q±(w)= to have a solution, respectively. These conditions provide a basic parameter range
for ±(·, r, s) = 0 to admit purely imaginary roots. Therefore, theorem3.2can be used with condition (3.9) or (3.10). For parameters not satisfying (3.9) or (3.10), we still can compute all positive solutions to Q±(w)= numerically, if they exist. (ii) It is also possible to use s
as the bifurcation parameter, as observed from the structure of the characteristic equation for (3.1). The computation will be similar to using r. (iii) We have seen that β4, the product of all decay rates di, i = 1, 2, 3, 4, has appeared in the conditions for theorem2.6and lemma3.1. The parameters corresponding to the stable synchronous oscillations lie outside the range for global convergence to the equilibrium in theorem2.6, namely,
β4 ρ1a2d3d4+ ρ2ρ3a2a4.
Indeed, β4serves as an organizer for the dynamics of system (1.1). We shall pursue this further in the next section.
3.2. Stability of periodic solution
Under the situation that there exist solutions w+to Q+(·) = and/or w−to Q−(·) = , the following value rcwhich stands for the first bifurcation value of r can be defined:
rc:= min{r : r = r+(k)(w+) >0 or r= r (k)
− (w−) >0, k∈ Z,
Q+(w+)= , Q−(w−)= }. (3.13)
Moreover, we set wcas the positive solution to Q+(·) = with rc= r (k)
+ (wc)or to Q−(·) = with rc= r
(k)
− (wc), for some k ∈ Z. We note that rcand wcare well defined under (3.9) or (3.10).
In this subsection, we shall analyze the stability of the periodic solution induced by
r= τ1+ τ2for system (2.5) at r near rc, by the normal form theory. To this end, we will first discuss the stability for the equilibrium of the system as r < rc. It will be shown that in some parameter range, the origin of system (2.5) is asymptotically stable for r < rcand any τ3 0 and τ4 0.
We consider (3.4) with r= 0; namely,
+(λ,0, s)· −(λ,0, s)= 0. (3.14)
We aim at deriving the conditions for parameters under which all roots of (3.14) have negative real parts, for all s 0. We first consider (3.14) with s= 0 as follows:
+(λ,0, 0)· −(λ,0, 0)= 0, (3.15)
±(λ,0, 0)= λ4+ α1λ3+ α2λ2+ α3λ+ α4±,
where α1:= β1, α2:= β2+a2γ1, α3:= β3+a2γ1(d3+d4)and α4±:= β4+a2d3d4γ1±a2a4γ2γ3. For later use, we denote
D1 := det(α1)= α1>0, D2 := det α1 α3 1 α2 = d2 4(d1+ d2+ d3)+ d4(d1+ d2+ d3)2+ (d1+ d2)[d32+ d1d3+ d2(d1+ d3)+ a2γ1] > 0, D±3 := det α11 αα32 α0±4 0 α1 α3 = α1α2α3− α32− α21α4±, D±4 := det α1 α3 0 0 1 α2 α±4 0 0 α1 α3 0 0 1 α2 α±4 = α4±· D±3.
Note that D1and D2are both positive. Therefore, by the well-known Routh–Hurwitz criterion, all roots of (3.15) have negative real parts if and only if
D3+>0, D3−>0, D+4 >0 and D4−>0. (3.16)
Now, we give the following result on the stability of the origin for (2.5) with r= s = 0.
Lemma 3.3. All roots of (3.15) have negative real parts, and hence the origin of system (2.5) with r = s = 0 is asymptotically stable if and only if
α1α2α3− (α3)2
(α1)2
− a2a4γ2γ3− a2d3d4γ1> β4> a2a4γ2γ3− a2d3d4γ1. (3.17)
Proof. Note that α+
4 > α−4 and D+3 < D−3. Thus α4− >0 and D3+ >0 if and only if (3.16) holds. It can be verified that (3.17) holds if and only if α4− >0 and D+
3 >0. The assertion
thus follows.
Following lemma3.3, we shall show that the origin is an asymptotically stable equilibrium of (2.5), for any s 0 and r = 0 under some additional condition. We substitute λ = iν with
ν 0 into ±(λ,0, s)= 0 and collect the real and imaginary parts to obtain ∓a2a4γ2γ3cos (νs)= ν4+ a2γ1(d3d4− ν2)− β2ν2+ β4, ±a2a4γ2γ3sin (νs)= ν[a2(d3+ d4)γ1− ν2β1+ β3].
(3.18) Summing up the square of equations (3.18) gives
where
E(ν)= ν8+ (−2a2γ1+ β12− 2β2)ν6
+ [a22γ12+ β22+ 2a2γ1(d3d4− (d3+ d4)β1+ β2)− 2β1β3+ 2β4]ν4 + [a22(d32+ d42)γ12+ β32− 2β2β4− 2a2γ1(d3d4β2− (d3+ d4)β3+ β4)]ν2 + (a2d3d4γ1+ β4)2.
A direct computation gives E(ν)= 8ν7+ E5ν5+ E3ν3+ E1ν,where
E5= 6(d12+ d 2 2+ d 2 3 + d 2 4)− 12a2γ1, E3= 4[d12d32+ d22(d12+ d32)+ 2a2d1d2γ1+ a22γ12+ d42(d12+ d22+ d32)]− 8a2(d32+ d42)γ1, E1= 2[d32(d1d2+ a2γ1)2+ d42(d22(d12+ d32)+ 2a2d1d2γ1+ a22γ12+ d12d32)]− 4d32d42a2γ1.
Lemma 3.4. There does not exist any real solution to (3.19) under conditions:
E1>0, E3>0, E5 >0, (3.20)
β4> a2a4γ2γ3− a2d3d4γ1. (3.21)
Proof. Obviously, E(ν) is strictly increasing on [0,∞), under condition (3.20). Moreover, it is straightforward to verify that E(0) = (a2d3d4γ1+ β4)2 > (a2a4γ2γ3)2 if (3.21) holds. Consequently, E(ν) > (a2a4γ2γ3)2, for all ν 0, and hence for all ν ∈ R, as E is an even
function.
Remark 3.2. By the definition of Ei and that γ1 is bounded, we find that, roughly speaking, inequalities (3.20) favor larger d1and d2and smaller a2, d3 and d4. Note that (3.17) implies (3.21).
By lemmas3.3and3.4, for arbitrarily fixed s 0 and r = 0, all roots of (3.14) have negative real parts under conditions (3.17) and (3.20). Moreover, under such assumptions, there does not exist any solution λ with nonnegative real part to ±(λ, r, s)= 0 for r < rc and any s 0. Accordingly, we conclude the following result.
Proposition 3.5. Assume that (3.9) or (3.10) holds. For any fixed s 0, the origin is an asymptotically stable equilibrium of system (2.5) for r < rc, under (3.17) and (3.20).
Remark 3.3. (i) In lemma3.1, condition (3.9) or (3.10) is used to guarantee the existence of purely imaginary roots iw±for ±(·, r, s) = 0. Herein, condition (3.10) is compatible with (3.21) as a2a4γ2γ3− a2d3d4γ1<0. If a2a4γ2γ3− a2d3d4γ1>0, we still can use (3.9) to find the solution of Q+(w)= , when β4satisfies (3.21). (ii) If both Q+(w)= and Q−(w)= do not have any solution, then for any fixed s 0 and r 0, the origin is an asymptotically stable equilibrium of system (2.5), under (3.17) and (3.20).
Setting σ = + (respectively, σ = −) if Q+(wc)= (respectively, Q−(wc) = ), we introduce the following conditions:
Condition (B1): Qσ(wc) = 0, and all other positive solutions to Q+(·) = and
Q−(·) = satisfy w = mwcfor any integer m;
Condition (B2): [Rσ(wc, rc)]2+ [Iσ(wc, rc)]2 = 0, W1,σ2 (wc, rc)+ W2,σ2 (wc, rc)= 0, and
![Figure 1. Graphs of functions z = Q − (w) , z = Q + (w) , z = , with the parameter data in [ 25].](https://thumb-ap.123doks.com/thumbv2/9libinfo/7934216.157300/15.892.154.638.123.357/figure-graphs-functions-z-q-q-parameter-data.webp)
![Figure 5. Dynamics of x 4 (t ), y 4 (t ) , for (a) t ∈ [0, 600], (b) t ∈ [2600, 3000], with](https://thumb-ap.123doks.com/thumbv2/9libinfo/7934216.157300/28.892.153.700.128.349/figure-dynamics-x-t-y-t-t-b.webp)
![Figure 7. Evolutions of x4 (t ), y 4 (t ) for t ∈ [2600, 3000], with (d ∗ , s) , τ 3 as in figure 6, and](https://thumb-ap.123doks.com/thumbv2/9libinfo/7934216.157300/29.892.152.697.132.346/figure-evolutions-x-t-y-t-τ-figure.webp)
![Figure 8. Evolutions of x4 (t ), y 4 (t ) for t ∈ [2600, 3000], in example 4.2; (a) asynchronous](https://thumb-ap.123doks.com/thumbv2/9libinfo/7934216.157300/30.892.145.703.131.339/figure-evolutions-x-t-y-t-example-asynchronous.webp)
