延遲性耦合震盪系統

行政院國家科學委員會補助專題研究計畫

（計畫名稱）延遲性耦合震盪系統

Coupled Oscillators with Delays

計畫類別：■個別型計畫 □整合型計畫

計畫編號：

NSC 98-2115-M-009-012-MY3

執行期間：

98 年 08 月 01 日起至民國 101 年 07 月 31 日

執行機構及系所：交通大學應用數學系

計畫主持人：

石至文

共同主持人：

計畫參與人員：廖康伶 蔡澤弘 陳冠維

成果報告類型(依經費核定清單規定繳交)：□精簡報告 ■完整報告

本計畫除繳交成果報告外，另須繳交以下出國心得報告：

□赴國外出差或研習心得報告

□赴大陸地區出差或研習心得報告

■出席國際學術會議心得報告

□國際合作研究計畫國外研究報告

處理方式：

除列管計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年■二年後可公開查詢

中 華 民 國 年 月 日

中文摘要：

本計劃探討具延遲性之耦合震盪系統之動態行為。我們考慮與感興趣的系統包

括一些生物系統與神經細胞模型。目前進行的有一斑馬魚骨節生成的模型、

FitzHugh-Nagumo model on neurons, Morris-Lecar model on neurons及一些artificial

neural models。我們探討的耦合系統行為包括同步化(synchronization)、同步震盪

(synchronous oscillation) 、反相同步化(anti-phase synchronization)等。系統在不

同參數及不同時間延遲下可能收斂到同步常態解、同步周期解、非同步周期解

等。我們已發表三篇論文，另有二篇論文投稿中，一篇論文準備投稿中。

有關 neural networks，去年發表一篇論文 Global synchronization and asymptotic

phases for a ring of identical cells with delayed coupling 於 SIAM J. Mathematical

Analysis。在這工作裡，我們討論系統的同步化及可能的最終同步行為，包括對

稱的單平衡點、多平衡點、同步周期解、非同步周期解。我們分析這些行為與系

統的大小、重要參數、遲滯項的大小之關連。我們也探討由遲滯項的影響而產生

的震盪及同步化的失去。依我們所建立的分析方法，我們回答了文獻中的兩個猜

想。

有關斑馬魚骨節生成之 segmentation clock genes 模型之數學研究已發表兩篇

論文:

K.-L. Liao, C.-W. Shih* and J.-P. Tseng, 2012, “Synchronized oscillations in a

Mathematical model of segmentation in zebrafish”, Nonlinearity 25, 869–904.

K.-L. Liao, C.-W. Shih*, 2012, “A lattice model on somitogenesis of

zebrafish”, Discrete and Continuous Dynamical Systems –Series B, 17 (8).

在前一篇文章裡，我們考慮的是兩個coupled cells, 共八個delay equations 的

系統，主要探討的是在什麼參數條件下，系統能產生周期為30 分鐘之同步周期

解，此為對應 zebrafish 骨節每節生成的時間，及推導oscillation-arrested 的條

件，此對應生成調控之一個階段，進而瞭解參數範圍所對應的動態行為；其中周

期解的存在是用delay Hopf bifurcation theory 做的，計算相當繁複。此外，我們

也做了不少數值模擬，一方面驗證理論，一方面延伸結論。在後一篇文章裡，我

們考慮的是一個N網格模型，在理論基礎之下，設計出參數隨時間及空間的變化，

方程式解的行為對應骨節生成中主要的基因表現動態。我們所建構的方法可以用

來處理其它的 gene regulatory models，並進一步探討生物時鐘、生物韻律的數學

模型。這部份的工作正在進行中，並納入下一期計畫的內容中。

有關 neuron model，我們建構一個 general framework，改良先前的 iteration

scheme，使得我們的方法可以應用在 coupled FitzHugh-Nagumo neurons 及

Morris-Lecar neuron neuron 之 synchronization，under gap junction (linear, diffusive)

coupling 或(nonlinear) synaptic coupling with transmission delay，已完成一篇論文

投稿。

這幾年，我們發展了一個方法(想法)，研究 nonlinear systems 及 coupled

nonlinear systems, with delay or without delay 的 asymptotic behaviors，包含 global

convergence to single equilibrium, existence of multiple equilibrium points，

multistability, convergent dynamics with multiple equilibria (or almost periodic

orbits), synchronization, synchronous oscillations, anti-phase oscillations 等，後二者

配合運用 delay Hopf bifurcation theory。過去，一般而言，研究 global dynamics

多依賴 Lyapunov function technique；其次，某些系統可運用 monotone dynamics

theory。但對複雜的非線性系統或延遲系統(including time-dependent delay), 這兩

種方法頗不可行。我們的方法，稱 sequential upper-lower-dynamics technique，目

前已運用到一些系統，可以解決一些以前的工具沒辦法解決的問題。過去文獻所

使用的數學工具相當有限，其研究 local dynamics 或 bifurcation of spatio-temporal

patterns，多依賴 linearization 及 equivariant bifurcation theory。這些工具很不容易

做到 combined effect of parameters and delays upon the dynamics 結果。Delay 的大

小影響 dynamics，但對非線性延遲系統 with multiple delays，如果系統整個

asymptotic behaviors 不是那麼單純，很少有數學結果明確說多大多小的 delay，

配合什麼樣的參數，是如何改變 dynamics。我們引入了新的方法，也得到新的結

果；我們相信這對非線性動力學中 coupled nonlinear systems with delays 的研究提

供了新的思維；而這些想法也讓我們可以處理其他生物數學的 models，包含一

些 integral-differential equations，multiple- component systems with multiple delays

or distribution delays。

關鍵詞: 收斂動態、震盪凍結、遲滯方程、同步化、同步周期解、多穩定性、

神經網路、分歧性、基因調控、體節生成

英文摘要:

We are interested in the dynamics for systems composed of subsystems with

delayed connection. The subsystems themselves can have time delay of magnitude

different from the connection delay. The systems we study include some gene

regulatory networks such as kinetic models for segmentation clock gene of zebrafish,

and neuronal models.

For neural network models, we published a paper “Global synchronization and

asymptotic phases for a ring of identical cells with delayed coupling” in SIAM J.

Mathematical Analysis. In this work, we consider a neural network which consists

of a ring of identical neurons coupled with their nearest neighbors. Self-feedback

delay and transmission delay are also taken into account in the network. We

investigate synchronization and various synchronous phases for the system. The

possible asymptotical dynamics include convergence to single equilibrium or multiple

equilibria, synchronous oscillations, and asynchronous oscillations in the form of

standing wave. We elucidate the effects from the scale of network, self-decay,

self-feedback strength, coupling strength and delay magnitude upon synchrony,

convergent dynamics and oscillation of the network. The disparity between the

contents of synchrony induced by distinct factors is investigated. Two different types

of multistable dynamics are distinguished. Moreover, oscillation and

desynchronization induced by delays are addressed. We then answer two

conjectures in the literature.

For the kinetic models for segmentation clock gene of zebrafish, we have published

two papers [K.-L. Liao, C.-W. Shih* and J.-P. Tseng, 2012, “Synchronized

oscillations in a Mathematical model of segmentation in zebrafish”, Nonlinearity 25,

869–904], [K.-L. Liao, C.-W. Shih*, 2012, “A lattice model on somitogenesis of

zebrafish”, Discrete and Continuous Dynamical Systems –Series B, 17 (8)].

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 traced out by the segmentation clock genes” which undergo

synchronous oscillation over adjacent cells. In the first work, we analyze the

dynamics for a model on zebrafish segmentation clock-genes which are subject to

direct autorepression by their own products under time delay, and cell-to-cell

interaction through Delta-Notch signaling. For this system of delayed equations, an

ingenious iteration approach is employed to derive the global synchronization and

global convergence to the unique synchronous equilibrium.

On the other hand, by applying the delay Hopf bifurcation theory and the method of

normal form, we derive the criteria for the existence of stable synchronous

oscillations. Our analysis provides the basic range of parameters and delay

magnitudes for stable synchronous, asynchronous oscillation, and oscillation-arrested

dynamics. Based on the derived criteria, further numerical findings on the dynamics

which are linked to the biological phenomena are explored for the considered system.

In the second work, we consider lattice systems which describe the kinetics of the

chief segmentation clock genes in zebrafish under negative feedback regulation with

delay through interaction with the Delta-Notch signaling among neighboring cells.

We first derive the analytical theories for the oscillation-arrested and synchronous

oscillation in an autonomous lattice model. Based on the parameter regimes in the

theories, we design suitable gradients of degradation rates and delays in a

non-autonomous lattice model. Such a lattice system can generate synchronous

oscillations, oscillatory traveling waves, oscillation slowing down, oscillation-arrested,

and high-low expression levels. We further distinguish between different gradient

structures which lead to normal and abnormal segmentations respectively and connect

these structures to the dynamical regimes in the cell-cell model.

Keywords: Convergent dynamics, Oscillation-arrested, Delayed equation,

Synchronization, Synchronous oscillation, Multistability, Neural network,

Bifurcation, Segmentation clock gene

報告內容:

1. 研究目的、文獻探討

Asymptotic behaviors and collective dynamics for nonlinear systems and coupled

nonlinear systems have been major research subjects in dynamical systems of

mathematics discipline and nonlinear dynamics of physics discipline.

Nevertheless

, mathematical methodologies in concluding global dynamics and

asymptotic behaviors for nonlinear systems and delay equations are quite limited.

The classical approaches are Lyapunov function method and monotone dynamics

theory. However, for complex nonlinear systems, construction of Lyapunov

function or partial order in phase space is often infeasible. In fact, one sees that

Lyapunov functions are expressed by basic functions, whereas the solutions of

differential equations are functions in general sense. The insufficiencies of

Lyapunov method (or Krasovskii theorem) have been revealed in several works; for

example, Slotine from MIT introduced the “Contraction analysis＂[E. Slotine et

al.,1998, On contraction analysis for nonlinear systems, Automatica 34 (6)]。From

there, a series of research works has been reported. However, this approach was

not formulated under mathematical rigorousness: the displacement equation or

variational equation is not precisely defined from the considered equations. If we

extend the interpretation under precise mathematical sense, then that consideration

is too restricted, as using the eigenvalues of certain corresponding linear systems to

obtain the dynamics of the original nonlinear systems leads to local dynamics or

merely partial behaviors. On the other hand, for equations with multiple delays,

analyzing the linearized equations is another complex task, as seen in several works

in the literature (e.g. the ones by S. A. Campbell).

Time delays exist in circuit systems, mechanical systems, physical systems,

biological systems, and traffic flows. Taking delay into account in differential

equations raises the phase space to infinite dimension. While regarded as a

functional differential equation, delay equation admits certain technical difficulty

and its fundamental theory, especially for the case of state-dependent delay, was just

developed recently. The role that time delay plays in mathematical modeling in

application fields has never been overemphasized. For example, time lags,

accounting synthesis and trafficking of macromolecules in cells, are required to

generate oscillation in gene regulation models, as illustrated in [J. Lewis, 2003,

Autoinhibition with transcriptional delay: a simple mechanism for the zebrafish

somitogenesis oscillator, Curr. Biol.].

Nevertheless

, the theory driven by the

application of delay-modeling requires further mathematical ideas and

developments, especially for systems with multiple delays and distribution delay.

Synchronous behavior is ubiquitous in nature, for example, simultaneous

flashing of fireflies, crickets chirping in unison, synchronous activity of pacemaker

cells in the heart, and synchronized bursts of pancreatic beta-cells. There is a large

literature on synchronization in physics and applied mathematics. Concerning

synchronization of coupled neurons, say models of Hodgkin-Huxley, Morris-Lecar,

FitzHugh-Nagumo, Hindmarsh-Rose, there are some mathematical results if the

neurons are under gap junction (linear, diffusive) coupling, but there does not exit

analytical result if the neurons are under synaptic coupling with delayed

transmission (most neurons are actually coupled via chemical synapses), as

commented in [E. Steur et al., 2009, Semi- passivity and synchronization of

diffusively coupled neuronal oscillators, Physica D, 238] .

Mathematical models for gene regulatory networks usually consist of multiple

components if there are several genes involved and their mRNA and proteins are

considered. For example, the basic mathematical cell-cell models on segmentation

clock gene of zebrafish are systems with at least eight components and four delays

which are incurred in the transcription and translation processes. The basic

dynamical phases for such models are synchronous oscillation and

oscillation-arrested. Mathematical analysis had been absent in investigating those

mathematical models. In fact, it was commented in [R. E. Baker, S. Schnell, 2010,

How can mathematics help us explore vertebrate segmentation ? HFSP Journal, 3:1]

that the majority of the models are incredibly difficult to analyze mathematically.

Almost all the existing results are obtained through numerical simulations. Yet the

studies based on numerical computations in the literature contain certain

shortcoming, as remarked in [Uriu, Morishita, Iwasa, 2009, Synchronized

oscillation of the segmentation clock gene in vertebrate development, J. Math. Biol.];

for example, convergence to the numerically observed periodic orbit is slow, since

the dynamical property of such orbit is unknown. Moreover, the basin of attraction

for the oscillation is very small, and hence the oscillation is difficult to observe in

phase space of high dimension.

Mathematical models on biology typically involve many parameters (one or two

dozens). It is a highly nontrivial task to understand the dynamics corresponding to

various parameter values and their various combinations through numerical

computations. Combining analysis and numerical computation is the best

approach in investigating mathematical models. However, developing effective

mathematical approaches to treat complex nonlinear systems with multiple

2. 研究方法:

We have developed a mathematical approach (idea) to investigate the asymptotic

behaviors of nonlinear systems and coupled nonlinear systems, with delay or without

delay, including global convergence to single equilibrium, existence of multiple

equilibrium points, multistability, convergent dynamics with multiple equilibria (or

almost periodic orbits), synchronization, synchronous oscillations, and anti-phase

oscillations. For the latter two on oscillation, delay Hopf bifurcation theory has also

been employed. Our methodology, named “sequential contracting”, has been

applied to a number of systems and resolved some unsolved problems. Designs of

pertinent upper and lower bounds which lead to the targeted asymptotic dynamics is

the first key step in applying sequential contracting. The methodology allows us to

establish global synchronization and multistability under delay-independent or

delay-dependent criteria through different designs of upper and lower dynamics.

For network systems, we can also derive scale-dependent criteria (scale is the

network size) for synchronization. Under our formulation, these nonlinear problems

are reduced to linear problems in a form of Gauss-Seidel iterations. Delays modify

dynamics, but analytic finding depicting how delays and what magnitudes of delays

change the dynamics is rare, especially for systems with nontrivial asymptotic phases

(not just an equilibrium). Our approach provides new thoughts in treating

multiple-component systems with multiple delays or distribution delays, which leads

to developing analytic studies on nonlinear mathematical models on biology.

3. 結果與討論

[K.-L. Liao, C.-W. Shih*, 2012, “A lattice model on somitogenesis of

zebrafish”, Discrete and Continuous Dynamical Systems –Series B, 17 (8)].

DISCRETE AND CONTINUOUS doi:10.3934/dcdsb.2012.17.2789 DYNAMICAL SYSTEMS SERIES B

Volume 17, Number 8, November 2012 pp. 2789–2814

A LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH

Kang-Ling Liao and Chih-Wen Shih* Department of Applied Mathematics

National Chiao Tung University Hsinchu 300, Taiwan

Abstract. Somitogenesis is the process of the development of somites which are segmental structure in vertebrate embryos. This process depends on the expression of segmentation clock genes. In this investigation, we consider lat-tice systems which describe the kinetics of the chief segmentation clock genes in zebrafish under negative feedback regulation with delay through interaction with the Delta-Notch signaling among neighboring cells. We first derive the analytical theories for the oscillation-arrested and synchronous oscillation in an autonomous lattice model. Based on the parameter regimes in the theories, we design suitable gradients of degradation rates and delays in a non-autonomous lattice model. Such a lattice system can generate synchronous oscillations, oscillatory traveling waves, oscillation slowing down, oscillation-arrested, and high-low expression levels. We further distinguish between different gradient structures which lead to normal and abnormal segmentations respectively and connect these structures to the dynamical regimes in the cell-cell model.

1. Introduction. In vertebrate embryos, somites are segmental structure which arises one by one from the presomitic mesoderm (PSM) and lies along the anterior-posterior (AP) axis. Under further development, somites develop into vertebrate, rib, and tail. The whole process of the development of somites, which involves gene regulation in both space and time, is called somitogenesis. During the process, the embryo can be divided into five regions: head, other tissues, determined region (DR), traveling wave region (TWR), and the tail bud (TB), successively.

For zebrafish, somite segmentation depends on the expression of clock genes her 1 and her 7, with neighboring cells interacted through Delta-Notch signaling [15, 16, 21, 23, 25]. It has been observed in biological experiments that gene ex-pression shows synchronous oscillation in the tail bud, and traveling wave pattern arises from the posterior end and moves toward the anterior end of the TWR. The oscillations then slow down, and finally arrest and cells form into somites. In ad-dition, the oscillation period of clock genes in the tail bud of zebrafish is about 30 minutes which match the time interval for somite formation [13,14]. This indicates that the period of the synchronously oscillating gene expression plays an important role in controlling when the somites form. On the other hand, the slowing down of oscillation is controlled by other proteins. Biologists have discovered that the signaling molecule fibroblast growth factor (FGF) regulates differentiation of the PSM cells [10, 17, 32]. Moreover, FGF is only transcribed at the posterior tip of

2000 Mathematics Subject Classification. Primary: 92B25, 37N25, 34K18, 34K13.

Key words and phrases. Somitogenesis, lattice model, oscillation-arrested, synchronous oscil-lation, traveling wave.

* corresponding author.

2790 KANG-LING LIAO AND CHIH-WEN SHIH

the embryo and progresses from the posterior to the anterior of the PSM with gra-dient. In addition, there exists a threshold such that if the FGF is smaller than this threshold, then cells start to differentiate into somites.

Modeling of the “clock and wavefront” for simitogenesis dates back to Cooke and Zeeman in 1976 [7]; therein, an idea of cyclic motion within the PSM cells along with cell fate determination was advocated. Since then substantial investigations through experiment or modeling have been reported. Based on some experimen-tal evidences [8, 9, 12, 14, 24], Baker et al. [4] proposed a “clock and wavefront” PDE model to investigate the pattern formation mechanism, which depends on the somite factor and the gradient of FGF8 expression. In 2003, Lewis [18] proposed a system of delayed differential equations (DDE) which model the negative feedback regulation and intercellular interaction via Delta-Notch signaling into segmentation clock genes kinetics to generate synchronously oscillatory expression. However, his original model can not produce oscillation slowing down and traveling wave pattern. On the other hand, Cinquin [6] proposed a 13-component multicellular DDE model for zebrafish somitogenesis, that involves heterodimerization of clock proteins, Her1 and Her7, with controlling protein Her13.2 [17,29]. The system was used to model the posterior-anterior slowing of oscillation rate, which leads to formation of clock-wave. Recently, Uriu et al. [30, 31] proposed an ODE model which includes an intermediate process of the gene expression (transport of Her protein from cyto-plasm to nucleus) to replace time delay to generate synchronous oscillations and traveling waves. Campanelli and Gedeon [5] investigated the control mechanisms in the formation of gene expression wave. Multiple transcription binding sites and differential decay rates for the monomers and dimers of the clock protein were con-sidered in generating the waveforms to match experimental observations.

From the observed phenomena, synchronous oscillation, traveling wave pattern, oscillation slowing down, and oscillation-arrested are necessary dynamics for normal segmentation. In addition, FGF or Her13.2 is an important factor to control the oscillation slowing down. Mutations in the Notch cell-cell signaling pathway and the abnormal traveling wave pattern disrupt synchronization, somite formation, and lead to defective somite boundary and irregular segmented body axis [16,18]. Indeed, a number of parts of these phenomena require further analysis to elucidate the whole mechanism for the somitogenesis.

The studies on segmentation clock and the involved internal machinery of the individual cell through mathematical modeling were largely based on numerical simulations, due to the complexity of gene networks. Recently, for Lewis’s two-cell model, analytical theories on synchronous oscillations and oscillation-arrested have been established through delay Hopf bifurcation theory and the sequential-contracting technique respectively, in [19]. The work in [19] sheds a light on the investigation of complex gene networks represented by coupled-cell systems. In this paper, we consider the homodimer model which takes her gene as either her 1 or her 7, and the neighboring cells are interacted through Delta-Notch signaling with delay. We shall extend the theories on synchronous oscillations in the tail bud and oscillation-arrested in the determined region to a lattice model for coupled n cells. Moreover, we show that the oscillation-arrested with high and low expression levels arises for non-identical cells with large degradation rates. Based on the collective behaviors of the multi-cellular model, we shall further modify Lewis’s model to a non-autonomous lattice system to generate traveling waves for the segmentation gene expression in zebrafish. We design a suitable anterior-posterior gradient of

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2791 degradation rates and delays in the lattice model according to the parameter regimes for the corresponding collective dynamics. The cyclic gene expression in this lattice system then exhibits a pattern of traveling wave. We distinguish between different gradient structures which lead to normal and abnormal segmentations respectively and connect these structures to the dynamical regimes in the autonomous model. It is observed that if the differential degradation rates are not posited properly, then the peaks of oscillation will move in incorrect direction or collide and lead to abnormal segmentation, where the size, forming period, and the number of somites are irregular.

The dynamics for the considered spatiotemporal pattern presented in the form of oscillatory traveling wave are rather intricate, as the synchrony among cells is lost, yet certain coherent rhythm remains. Traveling wave formation for the clock gene expression of zebrafish has also been investigated via an ODE model in [30], where the gradient for the reaction parameters was formulated according to numerical studies on the model. The motivation for this investigation is to illustrate that such a traveling wave also exists in delay model. Moreover, in order to elucidate the somitogenesis through mathematical modeling, on the basis of both analytic theories and numerical simulations, we attempt to link the clock-wave patterns to the collective dynamics of the coupled system and learn the underlying mechanism behind the parameter gradients. In addition, it is appealing to map the factors for normal and abnormal segmentations in the considered model to the cell-cell dynamics.

The above-mentioned models of clock gene expression and their interaction with signaling genes are rather difficult to analyze mathematically. As commented in [3], by forgetting about the internal machinery of the segmentation clock in each cell, researchers have turned to employ the notion of phase function to explore the collected dynamics of the coupled oscillators, where an oscillator represents an oscillatory gene in each cell or a group of synchronous cells [22]. Let us mention other related works. A PDE model which describes the movement of cells via cell-cell adhesion mechanism to investigate pattern formation for the aggregation behavior of cell population was reported in [1]. By including an explicit equation for adhesive cell population to the chemical signaling model in [4], the experimentally observed pattern of cellular aggregate in forming somites was reproduced in [2]. On the other hand, periodic traveling wave, as a non-uniform oscillatory spatiotemporal pattern, occurs in cyclic populations and other ecologically relevant scenarios modeled by oscillatory reaction-diffusion equations [26].

The presentation is arranged as follows. In section 2, we extend the analytic theories for oscillation-arrested and synchronous oscillation for segmentation gene to n-cell system and the asynchronous oscillation for the system of two non-identical cells. Then by using the information of these collective behaviors, in section 3, we construct a non-autonomous lattice model with suitable gradients of degradation rates and delays. Through numerical simulations, we illustrate propagation of oscil-latory waves corresponding to normal and abnormal segmentations respectively, and discuss the associated gradient structures. The presentation ends with a conclusion. 2. Autonomous lattice system. In this section, we establish analytical theories for the n-cell kinetic model of clock gene expressions in the cells of zebrafish em-bryo. Oscillation-arrested for the cells at the DR, asynchronous oscillations at the TWR, and synchronous oscillation at the TB shall be discussed. These theories are

2792 KANG-LING LIAO AND CHIH-WEN SHIH

extensions from the investigation of the Lewis’s two-cell model in [19]. The results in this section provide a basis for setting up the gradient structure of the activation parameters and delay magnitudes in a non-autonomous system which exhibits trav-eling wave of gene expression, presented in section 3. We arrange the theories on oscillation-arrested and oscillations in subsections 2.1 and 2.2, respectively. Based on these analytical theories, we summarize the collective behavior of the considered coupled system in subsection 2.3.

We assume that there are n cells located along the AP axis and each cell is identified by index i, i = 1, · · · , n, in succession along the axis. The neighboring cells are interacted through the Delta-Notch signaling. Let x1(i), x2(i), x3(i), x4(i) denote the concentrations of her mRNA, Her protein, delta mRNA, and Delta protein in cell i, respectively. We consider the following kinetic equations:

      

˙x1(i)(t) = gH(x2(i)(t − τ1(i)), ˜x4(i)(t − τ1(i))) − d1(i)x1(i)(t) ˙x2(i)(t) = a2x1(i)(t − τ2(i)) − d2(i)x2(i)(t)

˙x3(i)(t) = gD(x2(i)(t − τ3(i))) − d3(i)x3(i)(t) ˙x4(i)(t) = a4x3(i)(t − τ4(i)) − d4(i)x4(i)(t),

(2.1)

where 1 ≤ i ≤ n, ˜x4(i) = [x4(i−1)+ x4(i+1)]/2, for i = 2, · · · , n − 1 and ˜x4(1) = x4(2), ˜x4(n) = x4(n−1). Herein, a2 (resp., a4) is the protein synthesis rate per mRNA molecule for her (resp., delta) gene; d1(i), d2(i), d3(i), d4(i)are the degradation (decay) rates for her mRNA, Her protein, delta mRNA, and Delta protein for cell i, respectively; τ1(i), τ2(i), τ3(i), τ4(i) are positive numbers which represent the time delays in the processes of her gene transcription, her gene translation, delta gene transcription, and delta gene translation and delivery to cell membrane for cell i, respectively. Functions gH and gD relate the transcription initiation rates to her mRNA and delta mRNA concentrations respectively; they are represented by

gH(u, v) := kH 1 + v P_D0 1 + _Pv D0 + u2 P2 0 , (2.2) gD(u) := kD 1 + u2 P2 0 , (2.3)

for u, v ≥ 0; kH(resp., kD) is the maximal synthesis rate of her (resp., delta) mRNA; P0(resp., PD0) is the critical number of molecules of Her (resp., Delta) protein per cell for inhibition of transcription (resp., activation of Notch). Note that the rate of x1(i), the her mRNA for the cell located at i, 1 < i < n, is promoted by the amount [x4(i−1)+ x4(i+1)]/2 which is the average of the Delta proteins of the neighboring cells at (i−1) and (i+1). For the cell at i = 1 or i = n, there is only one neighboring cell. Both her mRNA and delta mRNA are suppressed by their own Her proteins, as expressed in the definitions of gH and gD.

We call system (2.1) the one for “n identical cells” if τj(i) = τj(1), for all i = 2, 3, · · · , n and j = 1, 2, 3, 4, and

dj(i)= dj(1), for all i = 2, 3, · · · , n, and j = 1, 2, 3, 4; (2.4) i.e., (2.1) becomes       

˙x1(i)(t) = gH(x2(i)(t − τ1), ˜x4(i)(t − τ1)) − d1x1(i)(t) ˙x2(i)(t) = a2x1(i)(t − τ2) − d2x2(i)(t)

˙x3(i)(t) = gD(x2(i)(t − τ3)) − d3x3(i)(t) ˙x4(i)(t) = a4x3(i)(t − τ4) − d4x4(i)(t),

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2793 for 1 ≤ i ≤ n, where dj = dj(1) and τj = τj(1), j = 1, 2, 3, 4. Herein, we call (2.1) the system for “n non-identical cells”, if (2.4) does not hold. We interpret that cells having different degradation rates is a consequence of gradient FGF.

We consider the evolution Ψ(t, φ) of (2.1) from initial condition φ = (φ1, · · · , φ4n) ∈ C([−τM, 0], R4n

+) at initial time t0= 0, where

τM := max {τ1(i), τ2(i), τ3(i), τ4(i), 1 ≤ i ≤ n}, R4n_{+ := {(x1, · · · , x4n)|xj ≥ 0, j = 1, · · · , 4n}.}

Let X(t; φ) be the solution of (2.1) defined by X(t + θ; φ) = Ψ(t, φ)(θ), θ ∈ [−τM, 0], for t > 0. We further denote X(t) = (x(1)(t), · · · , x(n)(t)) = X(t; φ) = (x(1)(t; φ), · · · , x(n)(t; φ)), where x(i) = (x1(i), x2(i), x3(i), x4(i)) for each i, if φ is not specified.

The basic dynamical properties for the system of two identical cells, i.e., n = 2 in (2.5), have been established in [19]. We present herein the extension of those prop-erties to the general coupled n-cell system (2.1). The following two propositions, derived from the idea of sequential component-estimates, show that every solution X_{(t; φ) of (2.1) exists and remains nonnegative for t ∈ [0, ∞), for any initial} con-dition φ ∈ C([−τM, 0], R4n

+), given any positive delays τj(i) and degradation rates dj(i), and other parameters. Their proofs resemble the ones for n = 2 in [19] and are omitted.

Proposition 2.1. C([−τM, 0], R4n

+) is positively invariant under the flow generated by system (2.1).

Proposition 2.2. There exists a compact set Q = Πn

i=1(Q1(i)× Q2(i)× Q3(i)× Q4(i)) ⊂ R4n

+, such that X(t; φ) converges to Q as t → ∞, for arbitrary φ ∈ C([−τM, 0], R4n

+), where Q1(i), Q2(i), Q3(i), and Q4(i) are defined by Q1(i)= [ˇq1(i), ˆq1(i)] := [0, kH

d1(i)], Q2(i) = [ˇq2(i), ˆq2(i)] := [0, a2kH d1(i)d2(i)], Q3(i)= [ˇq3(i), ˆq3(i)] := [cˇ3(i)

d3(i), kD

d3(i)], Q4(i) = [ˇq4(i), ˆq4(i)] := [ a4ˇc3(i) d3(i)d4(i), a4kD d3(i)d4(i)], and ˇc3(i):= kDP2 0/[P02+ ( a2k H d1(i)d2(i)) 2_{], for 1 ≤ i ≤ n.}

2.1. Oscillation-arrested. In this subsection, we study the kinetics of gene ex-pressions for the n cells located in the determined region. For the system of n non-identical cells (2.1), we shall show that their gene expressions tend to a non-uniform steady state corresponding to high and low expression levels as all degradation rates are large enough.

For the system of n identical cells, it is straightforward to verify that there exists a synchronous equilibrium. For system (2.1) of non-identical cells, we utilize the implicit function theorem to deduce the existence of an asynchronous equilibrium. We then employ the sequential-contracting technique and the Gauss-Seidel itera-tion to derive the criterion for the global convergence to this equilibrium. Such a dynamical scenario is linked to the state of oscillation-arrested for the cells at the DR.

First, analogous to the two-cell case in [19], there exists a positive synchronous equilibrium for the system (2.5) of n identical cells.

Proposition 2.3. There exists a positive synchronous equilibrium point ¯X_{= (¯x, · · · ,} ¯

x) for system (2.5), where ¯x= (¯x1, ¯x2, ¯x3, ¯x4), ¯ x1= d2x2¯ a2 , ¯x3= kDP2 0 d3(P2 0 + ¯x22) , ¯x4= a4kDP 2 0 d3d4(P2 0 + ¯x22) ,

2794 KANG-LING LIAO AND CHIH-WEN SHIH and ¯x2 is a positive solution to the equation:

P (ξ) := c5ξ5_{+ c3ξ}3_{− c2ξ}2_{+ c1ξ = a2kHP}4 0(d3d4PD0+ a4kD), with c1 := d1d2P4 0(d3d4PD0 + a4kD) > 0, c2 := a2d3d4kHP 2 0PD0 > 0, c3 := 2d1d2d3d4P2 0PD0 > 0, and c5:= d1d2d3d4PD0> 0.

Proposition2.3follows from solving the stationary equations of (2.5). With the existence of the synchronous equilibrium ¯X_{, we further define}

¯

γ1:= −∂gH

∂u (u, v)|u=¯x2,v=¯x4> 0, ¯γ2:= ∂gH ∂v (u, v)|u=¯x2,v=¯x4 > 0, ¯ γ3:= −dgD du (u)|u=¯x2 > 0, (2.6) M1:=     −d1 −¯γ1 0 0 a2 −d2 0 0 0 −¯γ3 −d3 0 0 0 a4 −d4     , M2:=     0 0 0 ¯γ2 0 0 0 0 0 0 0 0 0 0 0 0     , M3:=     0 0 0 ¯γ2/2 0 0 0 0 0 0 0 0 0 0 0 0     , (2.7) M :=             M1 M2 0 · · · 0 M3 M1 M3 . .. ... 0 M3 M1 M3 . .. ... .. . . .. ... ... ... 0 .. . . .. M3 M1 M3 0 · · · 0 M2 M1             4n×4n . (2.8)

We derive the existence of equilibrium for system of non-identical cells (2.1) as follows. Let int(R4n

+) denote the interior of R4n+ .

Proposition 2.4. For given degradation rate ¯d := (d1, d2, d3, d4, · · · , d1, d2, d3, d4) ∈ int(R4n

+) for system (2.5) and other fixed parameters and delays in (2.5), if

det(M ) 6= 0, (2.9)

then system (2.1) admits a positive equilibrium point X∗_{= (x}∗ (1), x ∗ (2), · · · , x ∗ (n)) with x∗_(i) = (x∗ 1(i), x ∗ 2(i), x ∗ 3(i), x ∗

4(i)), i = 1, 2, · · · , n, provided that the degradation rate (d1(1), d2(1), d3(1), d4(1), · · · , d1(n), d2(n), d3(n), d4(n)) is close to ¯d.

Proof. We define a function F : R8n_{→ R}4n _{by F (X, d) = (F1(X, d), F2(X, d), · · · ,} F4n(X, d)) with

F4(i−1)+1(X, d) := gH(x2(i), ˜x4(i)) − d1(i)x1(i), F4(i−1)+2(X, d) := a2x1(i)− d2(i)x2(i), F4(i−1)+3(X, d) := gD(x2(i)) − d3(i)x3(i),

F4i(X, d) := a4x3(i)− d4(i)x4(i),

for X = (x(1), x(2), · · · , x(n)) and d = (d(1), d(2), · · · , d(n)) with x(i)= (x1(i), x2(i), x3(i), x4(i)) and d(i) = (d1(i), d2(i), d3(i), d4(i)), for i = 1, 2, · · · , n. For any

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2795 ¯

d = (d1, d2, d3, d4, · · · , d1, d2, d3, d4) ∈ int(R4n

+), there exists ¯X ∈ int(R4n+) such that F ( ¯X_{, ¯d) = 0, thanks to Proposition2.3. According to (2.9), we have}

J = |∂Fk

∂x` |(X,d)=(¯X,¯d)6= 0, for k, ` = 1, 2, · · · , 4n;

therefore, by the implicit function theorem, if d is close to ¯d, then there exists a positive equilibrium X∗ _{= (x}∗

(1), x ∗

(2), · · · , x ∗

(n)) close to ¯Xsuch that F (X ∗_{, d) =} 0.

Remark 2.1. (i) If system (2.1) is designated for n non-identical cells, then the positive equilibrium X∗ _{in Proposition2.4is asynchronous in general, see Example} 2.1. (ii) For n = 2, the condition (2.9) can be replaced by d1d2d3d4+ a2¯γ1d3d4− a2a4¯γ2¯γ36= 0.

While system (2.1) admits an equilibrium X∗_{, we can discuss the global} conver-gence to X∗_{. Let us translate system (2.1) from X}∗ _{to the origin, and retain the} same notation:           

˙x1(i)(t) = −d1(i)x1(i)(t) − d1(i)x∗ 1(i) +gH(i)(x2(i)(t − τ1(i)) + x∗

2(i), ˜x4(i)(t − τ1(i)) + ˜x ∗ 4(i)) ˙x2(i)(t) = −d2(i)x2(i)(t) + a2x1(i)(t − τ2(i))

˙x3(i)(t) = −d3(i)x3(i)(t) − d3(i)x∗

3(i)+ gD(x2(i)(t − τ3(i)) + x ∗ 2(i)) ˙x4(i)(t) = −d4(i)x4(i)(t) + a4x3(i)(t − τ4(i)),

(2.10)

where 1 ≤ i ≤ n and ˜x∗

4(i) = [x∗4(i−1)+ x∗4(i+1)]/2, for i = 2, · · · , n − 1 and ˜x∗4(1)= x∗

4(2), ˜x ∗ 4(n) = x

∗

4(n−1). Note that after the translation, X(t) eventually converges to Q − X∗_{= Π}n i=1(Q∗1(i)× Q ∗ 2(i)× Q ∗ 3(i)× Q ∗ 4(i)), as t → ∞, where Q ∗ j(i):= [ˇqj(i)− x∗ j(i), ˆqj(i)− x ∗

j(i)], for j = 1, 2, 3, 4 and i = 1, 2, · · · , n.

We shall show that every solution of system (2.10) converges to the origin, as t → ∞. By the mean value theorem, we put each component of (2.10) in the form ˙xj(i)(t) = −dj(i)xj(i)(t) + wj(i)(t), j = 1, 2, 3, 4, 1 ≤ i ≤ n, (2.11) where w1(1)(t) := ∂gH ∂u (u(1), v(1)) · x2(1)(t − τ1(1)) + ∂gH ∂v (u(1), v(1)) · x4(2)(t − τ1(1)), w1(i)(t) := ∂gH

∂u (u(i), v(i)) · x2(i)(t − τ1(i)) +∂gH

∂v (u(i), v(i)) · [x4(i−1)(t − τ1(i)) + x4(i+1)(t − τ1(i))]/2,

for i = 2, · · · , n − 1, w1(n)(t) := ∂gH ∂u (u(n), v(n)) · x2(n)(t − τ1(n)) + ∂gH ∂v (u(n), v(n)) · x4(n−1)(t − τ1(n)), and for 1 ≤ i ≤ n,

w2(i)(t) := a2x1(i)(t − τ2(i)), w3(i)(t) := dgD

du (η(i)) · x2(i)(t − τ3(i)), w4(i)(t) := a4x3(i)(t − τ4(i)),

with u(i):= u(i)(t − τ1(i)) (resp., η(i) := η(i)(t − τ3(i))) lying between x2(i)(t − τ1(i)) + x∗

2(i)and x∗2(i)(resp., x2(i)(t−τ3(i))+x∗2(i)and x∗2(i)) and v(i) := v(i)(t−τ1(i)) between x4(2)(t − τ1(1)) + x∗

2796 KANG-LING LIAO AND CHIH-WEN SHIH x∗

4(i+1))/2 and (x∗4(i−1)+ x∗4(i+1))/2, and x4(n−1)(t − τ1(n)) + x∗4(n−1) and x∗4(n−1), for i = 1, i = 2, 3, · · · , n − 1, and i = n, respectively.

It can be shown that there exist 4n intervals Ij(i) := [−δj(i), δj(i)] to which the (4i + j − 4)th component of X(t) converges respectively, for j = 1, 2, 3, 4 and 1 ≤ i ≤ n. Moreover, the lengths of the intervals can be estimated as

0 ≤ δj(i)≤ |wj(i)|max(∞)/dj(i), (2.12) where |wj(i)|max_{(t) := sup{|wj(i)(s)| : s ≥ t}, and}

|wj(i)|max_{(∞) := limt→∞|wj(i)|}max_{(t), for j = 1, 2, 3, 4, 1 ≤ i ≤ n; cf. [27,28]. Let} us define

ρ1(1):= max u∈Q2(1), v∈Q4(2)

{|∂gH(u, v)

∂u |}, ρ2(1):=u∈Q2(1)max, v∈Q4(2)

{|∂gH(u, v) ∂v |}, ρ1(i):= max

u∈Q2(i), v∈ ˜Q4(i)

{|∂gH(u, v)

∂u |}, ρ2(i):=u∈Q2(i)max, v∈ ˜Q4(i)

{|∂gH(u, v) ∂v |}, for 1 < i < n, ρ1(n):= max u∈Q2(n), v∈Q4(n−1) {|∂gH(u, v)

∂u |}, ρ2(n):=u∈Q2(n)max, v∈Q4(n−1)

{|∂gH(u, v) ∂v |}, ρ3(i):= max u∈Q2(i) {|dgD(u) du |}, for 1 ≤ i ≤ n, (2.13)

where ˜Q4(i):= [(ˇq4(i−1)+ ˇq4(i+1))/2, (ˆq4(i−1)+ ˆq4(i+1))/2], for 1 < i < n. Now, we present the main result of this subsection.

Theorem 2.5. Assume that an equilibrium X∗ _{of system (2.1) exists. Then every} solution of (2.1) converges to the equilibrium X∗ _provided

d1(i)> ρ1(i)+ ρ2(i), d2(i) > a2, d3(i)> ρ3(i), d4(i)> a4, (2.14) for all i = 1, 2, · · · , n.

Proof. By utilizing the estimate in (2.12) to each component in (2.11) successively, it can be shown that there exists a sequence of estimates for δj(i), namely, δj(i)≤ δ(k)_j(i), j = 1, 2, 3, 4, i = 1, · · · , n, where for each k ∈ N,

0 ≤ δ_1(i)(k) :=        (ρ1(1)δ₂₍₁₎(k−1)+ ρ2(1)δ₄₍₂₎(k−1))/d1(1), if i = 1,

(ρ1(i)δ_2(i)(k−1)+ ρ2(i)(δ(k)_4(i−1)+ δ_4(i+1)(k−1))/2)/d1(i), if 1 < i < n, (ρ1(n)δ(k−1)_2(n) + ρ2(n)δ_4(n−1)(k) )/d1(n), if i = n,

0 ≤ δ_2(i)(k) := (a2/d2(i))δ(k)_1(i), 0 ≤ δ_3(i)(k) := (ρ3(i)/d3(i))δ(k)_2(i), 0 ≤ δ_4(i)(k) := (a4/d4(i))δ(k)_3(i),

(2.15) and δ(0)_2(i) := max{|ˇq2(i)− x∗

2(i)|, |ˆq2(i)− x ∗ 2(i)|}, δ

(0)

4(i) := max{|ˇq4(i)− x ∗

4(i)|, |ˆq4(i)− x∗

4(i)|}, with ρ1(i), ρ2(i), and ρ3(i) defined in (2.13). Moreover, the (4i + j − 4)th component for the solution X(t) of system (2.10) converges to I_j(i)(k) := [−δ(k)_j(i), δ(k)_j(i)], as t → ∞, for all k ∈ N. The detailed arguments are similar to those in [19]. From the iterative estimates (2.15), we observe that the sequence {δ(k)_j(i)|j = 1, 2, 3, 4, 1 ≤

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2797 i ≤ n, k ∈ N} is in the form of the Gauss-Seidel iteration for solving the linear system

(ML + E)y = 0, where M := [mpq]1≤p,q≤4n with

m4(p−1)+2,4(p−1)+1= −a2, m4p,4(p−1)+3= −a4, for 1 ≤ p ≤ n,

m4(p−1)+1,4(p−1)+2= −ρ1(p), m4(p−1)+3,4(p−1)+2 = −ρ3(p), for 1 ≤ p ≤ n, m4p+1,4p= m4p+1,4(p+2)= −ρ2(p+1)/2, for 1 ≤ p ≤ n − 2, m1,8= −ρ2(1), m4(n−1)+1,4(n−1)= −ρ2(n), mpq= 0, otherwise, L := I4n×4n, E := diag{d1(1), d2(1), d3(1), d4(1), · · · , d1(n), d2(n), d3(n), d4(n)}.

Since ML + E is strictly diagonal-dominant under condition (2.14) (cf. [27,33]), we obtain δ_j(i)(k) → 0, as k → ∞, for each 1 ≤ j ≤ 4 and 1 ≤ i ≤ n.

Remark 2.2. (i) Theorem2.5shall be adopted to indicate that if the cells located in the DR have large degradation rates, then the gene expressions of her and delta in these cells eventually converge to some steady state X∗_{. For non-identical cells} (with different degradation rates) in this region, X∗ _{is non-uniform among different} cells. We shall present an example below to show high and low expressions in the components of X∗_{, which correspond to the anterior and posterior of a formed} somite. We interpret this scenario as that when the peak and foot of the oscillatory wave arrive at the anterior end of the TWR, there generate distinct reactions which vary the degradation rates and form high and low expression levels, respectively. (ii) If we consider system (2.5) of n identical cells, then there exist an attracting set Q = (Π4

i=1Qi)n⊂ R4n+, with

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)], where ˇc3:= kDP2 0/[P02+(a2k H d1d2)

2_{], and a positive synchronous equilibrium ¯X, thanks} to Proposition2.3. By applying arguments analogous to the ones in Theorem2.5, it can be shown that every solution of system (2.5) converges to the synchronous equilibrium ¯X_{, provided} ˆ ρ1a2d3d4+ ˆρ2ρ3a2a4ˆ < d1d2d3d4, (2.16) where ˆρ1 := 2kHP_D0qˆ2(PD0+ˆq4) P2 0(PD0+ˇq4)2 , ˆρ2 := kHP_D0qˆ2₂ P2 0(PD0+ˇq4)2, and ˆρ3 := 2a2kDkH P2 0d1d2 . Note that condition (2.16) is weaker than condition (2.14), and the assertion reveals that the gene expressions eventually tend to a uniform steady state ¯X_.

The following example shows the existence of an asynchronous equilibrium in the system of two non-identical cells.

Example 2.1. We consider the system of two non-identical cells, i.e., system (2.1) with n = 2; we set

2798 KANG-LING LIAO AND CHIH-WEN SHIH First, if these two cells have close degradation rates

d1(1)= 2.3, d2(1)= 4.1, d3(1)= 2.3, d4(1) = 4.1, d1(2)= 2.31, d2(2)= 4.2, d3(2) = 2.29, d4(2)= 4.2, then the system admits an asynchronous equilibrium

X∗_{≈ (12.82, 14.07, 12.77, 14.01, 12.82, 13.74, 12.89, 13.81),} and solutions tend to this equilibrium X∗ _{if the initial value is close to X}∗_.

Next, if we set the degradation rates of these two cells apart from each other, and satisfy condition (2.14):

d1(1)= 2.3, d2(1) = 5, d3(1) = 2.3, d4(1)= 5, d1(2)= 4, d2(2) = 4.6, d3(2)= 4, d4(2) = 4.6, then the system admits a globally stable asynchronous equilibrium

X∗_{≈ (13.2, 11.88, 13.18, 11.87, 7.96, 7.78, 7.95, 7.78).}

Note that the gene expression of cell 1 (resp., cell 2) tends to a higher level (resp., lower level) of expression which corresponds to the anterior (resp., posterior) of a somite, as these two cells are located at the DR.

2.2. Oscillations. In zebrafish embryo, clock genes for the cells at the TB show synchronous oscillations with period about 30 minutes. Such a scenario can be depicted by considering coupled system (2.5) of n identical cells. Employing the delay Hopf bifurcation theory to establish synchronous periodic solutions has been completed for the two-cell case in [19]. Note that the synchronous manifold

S := {(x(1), · · · , x(n)) : x(i) = x(1)∈ C([−τM, 0], R4+), for all i = 2, · · · , n} is positively invariant under the flow generated by system (2.5). Therefore, we can consider the dynamics of (2.5) restricted to S:

       ˙x1(t) = gH(x2(t − τ1), x4(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). (2.17)

The computation and analysis for the two-cell case in [19] can then be adapted to the n-cell system (2.5).

On the other hand, oscillatory wave of gene expression traveling from the poste-rior end to the anteposte-rior end in the TWR has been observed in experiments. This spatiotemporal pattern indicates that gene expressions for the cells at the TWR remain oscillatory. These oscillations are out of synchrony, but sustain a suitable spatiotemporal pace, so that the peaks of oscillations can move anteriorly without collision.

In this subsection, we consider the cells located in the TWR with different degra-dation rates, and show that such asynchronous oscillation can be generated. We also illustrate that cell with larger degradation rate has smaller amplitude of os-cillation. For simplicity, we consider the coupled system of two cells with different

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2799 degradation rates                            ˙x1(1)(t) = gH(x2(1)(t − τ1), x4(2)(t − τ1)) − d1(1)x1(1)(t) ˙x2(1)(t) = a2x1(1)(t − τ2) − d2(1)x2(1)(t) ˙x3(1)(t) = gD(x2(1)(t − τ3)) − d3(1)x3(1)(t) ˙x4(1)(t) = a4x3(1)(t − τ4) − d4(1)x4(1)(t) ˙x1(2)(t) = gH(x2(2)(t − τ1), x4(1)(t − τ1)) − d1(2)x1(2)(t) ˙x2(2)(t) = a2x1(2)(t − τ2) − d2(2)x2(2)(t) ˙x3(2)(t) = gD(x2(2)(t − τ3)) − d3(2)x3(2)(t) ˙x4(2)(t) = a4x3(2)(t − τ4) − d4(2)x4(2)(t), (2.18)

where dj(1)6= dj(2), for some j ∈ {1, 2, 3, 4}. Note that the synchronous manifold S is no longer positively invariant under system (2.18). Our goal is to employ delay Hopf bifurcation theory to explore the existence of asynchronous periodic solution for system (2.18).

Following the discussions in section 2, we assume that system (2.18) has an asynchronous equilibrium X∗_{= (x}∗

(1), x∗(2)). We then translate system (2.18) from X∗ _{to the origin; the linearized system at the origin is}

                       ˙x1(1)(t) = −d1(1)x1(1)(t) + `11x2(1)(t − τ1) + `12x4(2)(t − τ1) ˙x2(1)(t) = −d2(1)x2(1)(t) + a2x1(1)(t − τ2) ˙x3(1)(t) = −d3(1)x3(1)(t) + `31x2(1)(t − τ3) ˙x4(1)(t) = −d4(1)x4(1)(t) + a4x3(1)(t − τ4) ˙x1(2)(t) = −d1(2)x1(2)(t) + `51x2(2)(t − τ1) + `52x4(1)(t − τ1) ˙x2(2)(t) = −d2(2)x2(2)(t) + a2x1(2)(t − τ2) ˙x3(2)(t) = −d3(2)x3(2)(t) + `71x2(2)(t − τ3) ˙x4(2)(t) = −d4(2)x4(2)(t) + a4x3(2)(t − τ4), (2.19) where `11:= ∂gH ∂u (x ∗ 2(1), x ∗ 4(2)), `12:= ∂gH ∂v (x ∗ 2(1), x ∗ 4(2)), `31:= dgD du (x ∗ 2(1)), `51:= ∂gH ∂u (x ∗ 2(2), x∗4(1)), `52:= ∂gH ∂v (x ∗ 2(2), x∗4(1)), `71:= dgD du (x ∗ 2(2)). (2.20) For convenience, we set γ1 = −`11, γ2 = `12, γ3 = −`31, γ4 = −`51, γ5 = `52, γ6= −`71 and γi> 0, i = 1, 2, · · · , 6.

The characteristic equation for (2.19) is given by 0 = 4(λ, τ1, τ2, τ3, τ4) := det  A(1) B(1) B(2) A(2)  , (2.21) where A(1):=   λ + d1(1) γ1e−τ1 λ 0 0 −a2e−τ2 λ λ + d2(1) 0 0 0 γ3e−τ3 λ λ + d3(1) 0 0 0 −a4e−τ4 λ λ + d4(1)  , B(1) := 0 0 0 −γ2e−τ1 λ 0 0 0 0 0 0 0 0 0 0 0 0 ! , A(2):=   λ + d1(2) γ4e−τ1 λ 0 0 −a2e−τ2 λ λ + d2(2) 0 0 0 γ6e−τ3 λ λ + d3(2) 0 0 0 −a4e−τ4 λ λ + d4(2)  , B(2) := 0 0 0 −γ5e−τ1 λ 0 0 0 0 0 0 0 0 0 0 0 0 ! .

2800 KANG-LING LIAO AND CHIH-WEN SHIH We compute to obtain 4(λ, τ1, τ2, τ3, τ4) = λ8+ β7λ7+ β6λ6+ β5λ5+ β4λ4+ β3λ3+ β2λ2+ β1λ + β0 +e−(τ1+τ2)λ_(α16λ6_{+ α15λ}5_{+ α14λ}4_{+ α13λ}3_{+ α12λ}2_{+ α11λ + α10)} +e−2(τ1+τ2)λ_(α24λ4_{+ α23λ}3_{+ α22λ}2_{+ α21λ + α20) − e}−2(τ1+τ2+τ3+τ4)λ_α30, where αij and βi are positive constants. We set r := τ1+ τ2, s := τ3+ τ4, and write 4(λ, τ1, τ2, τ3, τ4) = 4(λ, r, s). To employ the bifurcation theory, we take r as a bifurcation parameter, while holding s fixed. Finding purely imaginary roots iw to the characteristic equation (2.21) can be proceeded by solving

eλr· 4(λ, r, s) = 0. (2.22)

The real and imaginary parts of equation (2.22) read as 

Re(w; r) := A0(w) + A1s(w) cos (rw) + A2s(w) sin (rw) = 0

Im(w; r) := B0(w) + B1s(w) cos (rw) + B2s(w) sin (rw) = 0, (2.23) where A0(w) := −α16w6+ α14w4− α12w2+ α10, A1s(w) := w8_{− β6w}6_{+ (α24+ β4)w}4_{− (α22+ β2)w}2_{+ (α20+ β0)} −α30cos (2sw), A2s(w) := β7w7_{− β5w}5_{+ (−α23+ β3)w}3_{+ (α21− β1)w + α30sin (2sw),} B0(w) := α15w5− α13w3+ α11w, B1s(w) := −β7w7+ β5w5− (α23+ β3)w3+ (α21+ β1)w + α30sin (2sw), B2s(w) := w8_{− β6w}6_{+ (−α24+ β4)w}4_{+ (α22− β2)w}2_{+ (−α20+ β0)} +α30cos (2sw). For any fixed s ≥ 0, we define

P0(w) := A2s(w)B1s(w) − A1s(w)B2s(w), P1(w) := A0(w)B2s(w) − A2s(w)B0(w), P2(w) := A1s(w)B0(w) − A0(w)B1s(w), and set Q1(w) := P12(w) + P22(w), Q2(w) := P02(w), for w ∈ R, Q(w) := Q2(w) − Q1(w), for w ∈ R.

We derive the criterion for the existence of purely imaginary eigenvalue iw∗and the corresponding bifurcation values {r(k)∗ (w∗)}k∈Zin the following lemma.

Lemma 2.6. If

(−α20+ α30+ β0)2(α20− α30+ β0)2< α210(−α20+ α30+ β0)2, (2.24) then Q(·) = 0 has at least one positive root w∗. Moreover, if

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2801

then there exists a corresponding sequence {r(k)∗ (w∗)}k∈Z defined by

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, (2.26) such that w∗ and r∗(k)(w∗) satisfy 4(iw∗, r

(k)

∗ (w∗), s) = 0, for k ∈ Z, where C(w) := P1(w)/P0(w), S(w) := P2(w)/P0(w), for w ∈ U.

Proof. Since the leading term of Q(w) is w32and Q(0) = (−α20+α30+β0)2(−α210+ (α20−α30+β0)2), through arguments similar to those in [19], we derive the existence of solution w∗ to Q(·) = 0 under condition (2.24). Next, according to (2.25), we can further define

C(w) := P1(w)/P0(w), S(w) := P2(w)/P0(w), for w ∈ U,

such that C2_{(w∗) + S}2_{(w∗) = 1, for all w∗ ∈ U , due to Q(w∗) = 0. Moreover,} taking cos (rw) = C(w∗) and sin (rw) = S(w∗) with w = w∗ _{and r defined by} (2.26) into (2.23) yields Re(w∗; r∗(k)(w∗)) = 0 and Im(w∗; r

(k)

∗ (w∗)) = 0. Hence, 4(iw∗, r∗(k)(w∗), s) = 0, for k ∈ Z.

Finally, to apply Hopf bifurcation theory, we need the following simple-root con-dition and a computable concon-dition which yields transversality. For simplicity, we denote r(k)∗ := r(k)∗ (w∗).

Condition (C1): Q0(w∗) 6= 0, and all other positive solutions to Q(·) = 0 are not integer multiples of w∗. Condition (C2): ∂ ∂λ(e λr_{· 4(λ, r, s))|} λ=iw∗,r=r (k) ∗ 6= 0 and Re(λ 0_(r(k) ∗ )) 6= 0. We thus obtain the following oscillation theorem for two non-identical cells. Theorem 2.7. For a fixed s ≥ 0, assume that there exists a positive solution w∗ to Q(·) = 0 satisfying (2.25) and conditions (C1) and (C2) for some r(k)∗ > 0 and k ∈ Z. Then Hopf bifurcation occurs at r = r∗(k)and an asynchronous periodic orbit is bifurcated from X∗ _{in system (2.18).}

In what follows, we provide an example to show that if system (2.18) of two non-identical cells satisfies the conditions in Theorem 2.7, then the amplitudes of the bifurcated asynchronous periodic solution for the two cells are distinct. Example 2.2. For system (2.18), we set

a2= a4= 4.5, kH= kD= 33, P0= 40, PD0 = 400,

τ2= 2.8, τ3= 40, τ4= 20, (2.27)

and choose

d1(1)= d2(1)= d3(1)= d4(1)= 1, d1(2)= d2(2) = d3(2)= d4(2)= 1.2. (2.28) Then the system satisfies the conditions in Theorem 2.7 and admits a bifurcated asynchronous periodic solution when τ1 is near rc− τ2≈ 4.13373 − 2.8 = 1.33373.

2802 KANG-LING LIAO AND CHIH-WEN SHIH

In Fig. 1, the periods of the oscillation for cell 1 and cell 2 are the same, but the amplitude of the oscillation for cell 1 is larger than the one of cell 2. We interpret that admitting smaller amplitude for cell 2 corresponds to the degradation rates closer to the values where the oscillations arrest.

(a) (b)

Figure 1. Time series of x2(1)(t) and x2(2)(t) in system (2.18) with parameters (2.27) and degradation rates (2.28), from initial value φ = (10, 180, 7, 130, 8, 170, 6, 138), for (a) t ∈ [0, 600], (b) t ∈ [2600, 3000]. The period of the solution is close to 12 minutes, but the amplitudes of x2(1)(t) and x2(2)(t) are about 9.34 and 2.57, respectively.

Remark 2.3. (i) Theorem 2.7 implicates that two coupled cells with different degradation rates exhibit asynchronous oscillations under certain conditions. Hence, the considered model (2.1) has the potential to generate moving oscillatory wave in a lattice system, if the cells located in the TWR possess pertinent gradient of the degradation rates along the AP axis. (ii) If the solutions w∗ to Q(·) = 0 satisfy (2.25), then the stability of the bifurcated asynchronous periodic solution can be analyzed by the center manifold theorem and the normal form method, as performed in [19] for the system of two identical cells.

2.3. Collective behaviors. So far, we have studied global convergence of dynam-ics for the n-cell lattice model (2.1), synchronous oscillations over identical cells, and asynchronous oscillation for non-identical cells. In this subsection, we sum-marize these dynamics and connect them to the corresponding biological scenarios, depicted in Fig. 2.

(D1): Identical cells located at the TB exhibit synchronous oscillations. This dy-namics can be examined via Hopf bifurcation theory for the system (2.5) of n identical cells. The parameter and delay conditions for such dynamics are as the ones for the system of two identical cells.

(D2): Cells located at the TWR with different degradation rates admit asynchro-nous oscillations. Moreover, if we consider two coupled non-identical cells, then this asynchronous oscillation can be examined via Theorem2.7.

(D3): The oscillatory gene expressions arrest and tend to high and low expression levels for non-identical cells with large degradation rates, located at the DR. This dynamics is obtained in Theorem2.5for general coupled n-cell system (2.1).

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2803

1 2

" "

i



1

_i

i



1

" "

N



1

N

Determined region (D3) Traveling wave region (D2) Tail bud (D1)

AP axis Anterior Posterior coupled coupled Cells coupled

Figure 2. N cells aligned in the AP axis, each cell is coupled with its nearest neighbors. The PSM which contains the traveling wave region and the tail bud moves posteriorly as the embryo grows.

Based on these theoretical results, in section 3, we construct a non-autonomous lattice system with gradients of degradation rates and delays to exhibit the biological scenarios in the TB, TWR, and DR.

3. Traveling wave in non-autonomous system. It has been observed in bio-logical experiments that normal segmentation requires a pertinent traveling wave of gene expressions propagating from the posterior toward the anterior of the TWR, with increasing period, decreasing amplitude, wave length, and wave speed. Herein, the wave speed is referred to the speed of peaks in the oscillations.

In this section, we shall explore the traveling wave patterns in the N -cell model equations. Based on the theoretical results and the understanding on the dynamics of the two-cell system in [19] and the n-cell system (2.1) established in section 2, we shall design suitable gradients of degradation rates and delay magnitudes to produce traveling wave patterns. This will lead to a non-autonomous lattice system where some degradation rates and delay magnitudes vary with respect to space and time. More precisely, we take into account the dynamics regimes in (D1) to gen-erate synchronous oscillations in the TB, in (D3) to arrest the oscillations in the DR, and combine the dynamics in (D2) and the cell-cell system to devise the gra-dients of degradation rates and delays. Under suitable design on the gragra-dients, this lattice system can generate various traveling wave patterns, but only some of them represent normal segmentation. As one of the main issues in current investigations of somitogenesis is to explore the cause of mutant, we shall pursue the factors that result in inappropriate propagation of waves, including wrong direction of waves and collision of oscillations and peaks.

We assume that there are N cells located along the AP axis, which are divided into three parts: DR, TWR, and TB, successively, and the PSM contains the TB and TWR. The gene expressions for cells at the TB show synchronous oscillations with period about 30 minutes; the oscillations then slow down with increasing period in the TWR, where traveling wave of gene expression proceeds from the posterior to the anterior. Finally, the oscillations arrest and cells develop into somites in the DR.

According to the experiments, during the growth of the embryo, the cells located at the anterior end of the TWR are removed from the TWR and become mature somites. Meanwhile, the cells located behind the posterior end of the TWR enter the TWR to compensate the loss of anterior cells. Moreover, the growth rate of the embryo is close to a constant. Therefore, this scenario indicates that the overall

2804 KANG-LING LIAO AND CHIH-WEN SHIH

length of the TWR maintains a constant value and the PSM moves posteriorly with a constant speed, along with the growth of embryo. According to these observations, we assume that the length of the TWR is constant, denoted by L, and the PSM moves posteriorly with a constant speed v := 1/m which matches the growth rate of the embryo, where L is an integer and m is a positive real number. When t = 0, the DR lies in 0 ≤ i ≤ K − 1, the TWR is in K ≤ i ≤ K + L − 1, and the TB is in i ≥ K + L, where K is an integer. Furthermore, for each t with m(˜k − 1) < t ≤ m˜k, the DR lies in i ≤ K + ˜k − 2, and the TWR is in K + (˜k − 1) ≤ i ≤ (K + L − 1) + (˜k − 1) with the anterior and posterior ends of the TWR at K + (˜k − 1) and at (K + L − 1) + (˜k − 1) respectively; the TB is in i ≥ K + L + (˜k − 1). Fig. 3 displays the regions TB, TWR, and DR along the AP axis, for ˜k ∈ N. In addition, the dynamics in regions TB, TWR, and DR correspond to (D1), (D2), and (D3) regimes, depicted in section 2.3, respectively.

"

#

#

"

"

#

#

"

m

t

i

TWR

DR

TWR

TWR

TWR

DR

DR

DR

Figure 3. The region anterior (resp., posterior) to the TWR is the DR (resp., the TB). The blue blocks represent the TWR which moves posteriorly one grid per m minutes as the embryo grows with speed v = 1/m.

To generate oscillatory traveling waves, we propose the following non-autonomous lattice system modified from (2.1):

      

˙x1(i)(t) = gH(x2(i)(t − τ1(i)(t)), ˜x4(i)(t − τ1(i)(t))) − d(i)(t)x1(i)(t) ˙x2(i)(t) = a2x1(i)(t − τ2) − d(i)(t)x2(i)(t)

˙x3(i)(t) = gD(x2(i)(t − τ3)) − d(i)(t)x3(i)(t) ˙x4(i)(t) = a4x3(i)(t − τ4) − d(i)(t)x4(i)(t),

(3.1)

for 1 ≤ i ≤ N , where gH and gDare defined as (2.2) and (2.3) respectively, and the parameters and delays are adopted from Lewis [11,18,23]:

a2= a4= 4.5, kD= kH= 33, P0= 40, PD0 = 400, τ2= 2.8, τ4= 20, (3.2) τ1= 3.8 ± 1.0 + Tinit−her, τ3= 8.4 ± 1.2 + Tinit−delta, (3.3) where Tinit−her (resp., Tinit−delta) is the time for bound inhibitory protein and for the mRNA polymerase to produce sufficient transcription for her1 (resp., delta)

LATTICE MODEL ON SOMITOGENESIS OF ZEBRAFISH 2805 gene. Note that the degradation rates

d1= 0.252525, d2= 0.23, d3= 0.273224, d4= 0.23,

adopted in a two-cell system in [23] for the wild type, are close to each other. Therefore, for simplicity, we set all degradation rates of a single cell identical at any time (but different cells can have different degradation rates) in system (3.1). This will allow us to focus on the effect from the gradient of degradation rates with respect to space (i).

From the numerical simulation on system (3.1), if we set the same constant degra-dation rates for all cells (d(i) = d for all i) and let τ1(i) decrease (resp., increase) as i decreases, then the system generates oscillatory wave from the anterior (resp., posterior) to the posterior (resp., anterior) of the TWR. Hence, the delay magni-tudes play an important role in the direction of wave propagation. On the other hand, the molecular gradients, such as FGF and Wnt signaling, control the period of oscillation for individual cell and cause oscillation-arrested. Hence, similar to the formulation in [30], we set the degradation rates to vary as a reaction to these molecular gradients, during the process. By combining these two factors, we shall design the non-autonomous system (3.1) with degradation and delay gradients to generate various oscillatory wave patterns.

Let us first classify the degradation rates. We denote d∗ = dj = dj(i), j = 1, 2, 3, 4, i = 1, 2, in the system of two identical cells. According to our analysis (Corollary 3.7 and Fig. 3 of [19]), with the parameters in (3.2), if we fix s = τ3+ τ4, then the ranges for d∗ to have stable synchronous periodic solution and to have stable asynchronous periodic solution appear in successive subintervals, one next to the other. We denote these subintervals as IS, IAS, and ISS, where IS ∪ ISS (resp., IAS) contains the degradation rates d∗with which the system admits a stable synchronous (resp., asynchronous) periodic solution when dj(i) = dj = d∗, j = 1, 2, 3, 4, i = 1, 2, and r near the minimal bifurcation value rc, via Hopf bifurcation. Moreover, from Hopf bifurcation theory, only synchronous periodic solutions appear when d∗ ∈ ISS, and either synchronous or asynchronous oscillation can occur when d∗ ∈ IS. In addition, the interval ISS is usually located on the right of IS ∪ IAS. For example, when s = 60 (τ3= 40, τ4= 20),

IS∪ IAS = [0.15, 1.3865] and ISS = [sa, sb] := [1.38652, 1.38824]. (3.4) Next, we employ the function ps := ps(d∗) to delineate the periods of these oscil-lations, for a given s ≥ 0. In Fig. 4, for s = 60, we depict ps(d∗) for d∗ ∈ Sd := {0.15, 0.16, · · · , 1.37, 1.38, 1.3865}. It is interesting to see that the period ps(d∗) decreases (resp., increases) with respect to d∗∈ {0.15, 0.16, · · · , 0.94, 0.95}, (resp., {0.96, 0.97, · · · , 1.37, 1.38, 1.3865}). Based on this variation of period in Fig. 4, for convenience, we set

Ipd:= [0.15, 0.95] and Ipi:= [0.96, 1.3865]. (3.5) We specify two degradation rates, dc and ds, for the setting of degradation gradient in (3.1). Let dc be a quantity satisfying condition (2.16), so that there is a globally asymptotically stable synchronous equilibrium ¯X _{for system (2.5) with dj = dc,} j = 1, 2, 3, 4, for any fixed r = τ1+ τ2≥ 0 and s = τ3+ τ4≥ 0. On the other hand, with dj= ds, j = 1, 2, 3, 4, there exists a stable periodic solution bifurcated at rcvia Hopf bifurcation theory. Furthermore, we can find the her transcription delay, τs, such that τs+ τ2is the minimal bifurcation value rcfor the system with parameters

2806 KANG-LING LIAO AND CHIH-WEN SHIH *

d

)

(

d

p

Figure 4. The periods of oscillations, ps(d∗), corresponding to a sequence of d∗ ∈ Sd, in the system of two identical cells with degradation rates dj(i) = dj = d∗, j = 1, 2, 3, 4, i = 1, 2, when r = τ1+ τ2 is close to the minimal bifurcation value rc.

(3.2) and degradation rate ds. Accordingly, we choose ds∈ IS∪ ISS to generate os-cillations, and take dc > dsas larger degradation rates to create oscillation-arrested. Next, we choose some value d`∈ (ds, dc) so that the system with parameters (3.2) and degradation rate d` has a stable periodic solution bifurcated at the minimal bifurcation value rc = τ`+ τ2.

We set the degradation rates equal to dcand ds, for the cells at the DR and TB respectively, so that oscillation-arrested in the DR and synchronous oscillation with period ps(ds) in the TB take place, according to regimes (D3) and (D1) respectively. Moreover, we let the degradation rates vary from ds at the posterior end to d` at the anterior end of the TWR. In addition, we divide the TWR into three parts R3, R2, and R1, from the anterior to the posterior, so that the degradation rates are continuous and decreasing linearly, delay τ1(i)is continuous and linear in each region, and the absolute values of the slope of degradation rate and τ1(i) are largest (resp., smallest) in R3(resp., R1). The arrangement is according to the dependence of the gradient slopes in d(i) and τ1(i) on the rate of period variation of the oscillations. Based on this setting, we design d(i)(t) and τ1(i)(t) for system (3.1) in the following five regions. We denote ˜d = (d`− ds)/3 and ˜τ = (τ`− τs)/3, and set d(i)(t) and τ1(i)(t) for cell i at time t in the respective region as follows:

(i) 1 ≤ i ≤ K − 1, i.e., the determined region initially,

d(i)(t) = dc, for t ≥ 0, τ1(i)(t) = τ`, for t ≥ 0;