A General Approach to Synchronization of Coupled Cells

Vol. 12, No. 3, pp. 1354–1393

A General Approach to Synchronization of Coupled Cells∗ Chih-Wen Shih† _{and Jui-Pin Tseng}‡

Abstract. This investigation presents a general framework to establish synchronization of coupled cells and coupled systems. Each individual subsystem is represented by nonlinear differential equations with or without internal or intracellular delay. A general coupling function is employed to depict the communication or interaction between subsystems or cells. Under this framework, the problem of establishing the synchronization for delayed coupled nonlinear systems is transformed to solving a corresponding linear system of algebraic equations. We start by considering a cell-to-cell system under symmetric coupling to present the main idea of the approach. The framework is then extended to theN-cell system under circulant coupling. Delay-dependent, delay-independent, and network-scale–dependent criteria for global synchronization will be established, respectively. The developed scheme can accommodate a wide range of coupled systems. We demonstrate the applications of the present approach to establish synchronization for a gene regulation model, a neuronal model, and some neural networks.

Key words. coupled system, synchronization, delay, gene regulation model, neuronal model, neural network AMS subject classifications. 34K18, 34K20, 92B20, 92C20

DOI. 10.1137/130907720

1. Introduction. Synchronization is a crucial and common phenomenon in various

bio-logical and physical systems. In many regions of the brain, synchronization activity has been observed and implicated as a correlate of behavior and cognition [64]. It is known that synchro-nization encourages the strengthening of mutual connections among neurons. Synchrony and synchronous oscillations are typical activities for gene expressions in cells under interaction. For example, in somitogenesis of vertebrate embryo, the cyclic genes express synchronous os-cillations in neighboring cells at the tail bud of the presomitic mesoderm [25,42,50]. There are many other beautiful examples, including simultaneous flashing of fireflies, crickets chirping in unison, and synchronous activity of pacemaker cells in the heart; see [45,46,60,61]. There is also a large number of studies on synchronization in engineering because of its importance in applications such as synchronized chaos employed in secret communication [14].

Time delays occur in the transcription and translation processes of somitogenesis due to synthesis and trafficking of macromolecules in cells; the lags have been estimated to be around tens of minutes in cell culture [28,42]. For connected neurons, time delay occurs in the prop-agation of action potentials along the axon and the transmission of signals across the synapse ∗_{Received by the editors January 29, 2013; accepted for publication (in revised form) by B. Sandstede May}

21, 2013; published electronically August 6, 2013. The authors were supported in part by the National Center for Theoretical Sciences and National Science Council of Taiwan.

http://www.siam.org/journals/siads/12-3/90772.html

†_{Corresponding author. Department of Applied Mathematics, National Chiao Tung University, National Center}

for Theoretical Sciences, Hsinchu, 300 Taiwan ([email protected]).

‡_{Department of Applied Mathematics, National Pingtung University of Education, Pingtung, 900 Taiwan}

([email protected]).

[7,13,16]. Modeling interacted cells or neurons with delays thus becomes an important con-cern in studying the collective behaviors for coupled-cell systems. On the other hand, as time lag also occurs in transmitting signals among artificial neurons, delay has been incorporated into the neural network modeling [6,9,37,38,51,66]. Indeed, delay can modify the collective dynamics for neural networks; for example, it can induce oscillation or change the stability of the stationary solution [9]. Delay can also induce synchronization [18,39], asynchrony [9], and oscillation death [1]. Crook et al. [13] studied a continuum model of the cortex, with ex-citatory coupling and distance-dependent delays, and found that for small enough delays the synchronous oscillation is stable but for larger delays this oscillation loses stability to a trav-eling wave. Therefore, in addition to synchronization, it is also important to investigate the synchronous phases and their transitions. Developing effective mathematical methodologies and analytic tools elucidating the synchronous activities and collective behaviors with respect to various combinations of parameters, coupling strengths, and delay magnitudes remains an important research task.

It is interesting to explore the dynamical mechanisms underlying the behavior of networks of neurons and biological oscillators. Although the real network architecture can be extremely complicated, rich dynamics arising from the interaction of simple network motifs are believed to provide sources of activities similar to those in real-life systems. Synchronization in coupled dynamical systems has attracted a lot of attention in recent decades. The coupled systems that were most studied for synchronization are various neural network systems and chaotic oscillators. Among this research, some conclude local synchronization which is concerned with the stability of synchronization manifold or solution behavior in a neighborhood of certain synchronous solution, while others obtain global synchronization by showing that all solutions converge to the synchronization manifold or some synchronous solution.

The master stability function, developed by Pecora and Carroll [43, 44], is a well-known approach to studying local synchronization of coupled chaotic systems. This method is based on computing the Lyapunov exponent of the associated variational equation to determine the stability of the synchronization manifold for the coupled systems. However, such an approach leads to a necessary (instead of sufficient) condition for local synchronization; cf. [23].

Methodologies for concluding global synchronization largely involve the notion of Lya-punov functions. For example, Belykh, Belykh, and Hasler employed the “connection graph stability method” combined with the Lyapunov function approach to studying global syn-chronization in small-world networks of chaotic systems [3]. From the viewpoint of feedback control, Nijmeijer and collaborators introduced the notion of passivity and semipassivity and constructed a Lyapunov–Razumikhin function to study global synchronization in coupled sys-tems [47,58,59]. Other works employing the Lyapunov function/functional technique include [7,11,29,33,35,36,48,49,67,68].

Another approach to investigating synchronous oscillation in coupled systems, especially in neural networks, is to apply bifurcation theory to obtain the existence of synchronous periodic solution and use the normal form theory and the center manifold theorem to discuss its stability [6,57]. However, linearization of delayed equations yields another complication, as the linearized system contains exponential functions which make the analysis and computation of characteristic values difficult, especially in the case of multiple delays; see the papers by Campbell and coworkers in [9,10].

Most of the couplings in systems, especially the delayed systems, considered in the lit-erature are linear or linearly diffusive. Such coupling terms are expressed by a summation of connection weights multiplied to the incoming signals from other units, or by coupling strength multiplied to the difference of two corresponding components. Notice that there is a significant difference between diffusive coupling and general nonlinear coupling. For a system comprising identical subsystems under diffusive coupling, its synchronous solution is also a solution for each individual (uncoupled) subsystem, as the coupling parts are annihilated at synchronous states. This is certainly not the case for the general nonlinear coupling scheme. For coupled systems with multiple delays, the analysis is even more complicated. There are nonlinear systems with more intricate coupling and multiple delays, which are distinguished from those systems mentioned above. One such system is the segmentation clock gene model for coupled cells, which describes the gene regulation for vertebrate embryos. In the presomitic mesoderm of zebrafish embryo, neighboring cells interact through delayed, intercellular pos-itive feedback via Delta–Notch signaling [28, 42]. A system modeling the same clock genes was proposed by Uriu, Morishita, and Iwasa in [63], where more complicated nonlinear terms accounting for the Michaelis–Menten-type degradations and general transcription and trans-lation functions with Hill coefficients are considered. Analytic study for such kinetic models is rather difficult, as mentioned in [2]. Other examples include the neuronal models coupled through chemical synapses and neural networks with nonlinear activation functions. We shall discuss these systems in later sections.

As described above, the Lyapunov function technique, used directly or indirectly, is a com-monly and largely adopted approach for tackling synchronization problems, especially global synchronization, in dynamical systems with or without delays. On the one hand, finding a Lyapunov function in highly nonlinear systems, especially with multiple delays, and/or sys-tems with multiple components of different types, seems rather infeasible. On the other hand, for systems with delay, synchronization results concluded from the Lyapunov function ap-proach often reduce to the situation that every solution converges asymptotically to a unique synchronous equilibrium point [9,69]. In addition, typically only delay-independent criteria can be derived under such an approach. Therefore, a mathematical approach to tackling asymptotic behaviors and synchronization for nonlinearly coupled systems, without resort-ing to the Lyapunov function and computations of the characteristic values, is an appealresort-ing advance.

In this investigation, we aim at establishing a general framework based on the idea of “se-quential contracting” to study synchronization for coupled systems. We start by considering the synchronization for a pair of identical subsystems under a general coupling:

(1.1)

 ˙

x(t) = F(xt, t) +G(xt,yt, t),

˙

y(t) = F(yt, t) +G(y_t,x_t, t),

where t≥ t0,x(t), y(t) ∈ Rn, andxt,yt∈ C([−τM, 0];Rn) with τM ≥ 0, are defined by xt(θ) =

x(t + θ), yt(θ) = y(t + θ) for θ ∈ [−τ_M, 0], F = (F1, F2, . . . , Fn) is a smooth function which

depicts the intrinsic dynamics of each subsystem, and a smooth functionG = (G₁, G₂, . . . , G_n) expresses the interaction between two coupled subsystems. Herein, time delays in the range [0, τ_M] are considered in the subsystems and coupling terms, (x_t,y_t) denotes the evolution of system (1.1) at time t from (x_t0,y_t0) in C([−τ_M, 0],Rn), and (x(t), y(t)) is the corresponding

solution of system (1.1). The present framework shall cover the ODE case, i.e., when τ_M = 0, and (1.1) reduces to _

˙

x(t) = F(x, t) + G(x, y, t),

˙

y(t) = F(y, t) + G(y, x, t),

where (x, y) lies in R2n or a positively invariant subset of R2n. Let us denote the synchronous set by

(1.2) S := {(x, y) ∈ R2n| x = y}.

We say that a solution of (1.1) is synchronous if it lies inS completely; a solution is

asymp-totically synchronous if its ω-limit set lies in S. The coupled system (1.1) is said to attain

global synchronization if every solution is asymptotically synchronous, i.e., x_i(t)− y_i(t)→ 0, as t → ∞, for all i = 1, . . . , n,

for every solution (x(t), y(t)) = (x₁(t), . . . , x_n(t), y₁(t), . . . , y_n(t)) of system (1.1).

For an arbitrary solution (x(t), y(t)) of system (1.1), where x(t) = (x₁(t), . . . , x_n(t)),

y(t) = (y1(t), . . . , yn(t)), by setting zi(t) = xi(t)− yi(t), we shall consider the following

difference-differential system corresponding to (1.1):

(1.3) z˙_i(t) = F_i(x_t, t) + G_i(x_t,y_t, t)− F_i(y_t, t)− G_i(y_t,x_t, t), i = 1, . . . , n.

System (1.1) attains global synchronization if and only if z_i(t)→ 0, i = 1, . . . , n, as t → ∞, for every (z₁(t), . . . , z_n(t)) satisfying system (1.3), defined from every solution (x(t), y(t)) of (1.1).

In the literature, studying the evolution for the difference of two corresponding compo-nents, such as (1.3), has been a primary target in tackling synchronization problems. The idea of sequential contracting provides a new treatment to analyze such difference-differential sys-tems. This approach unfolds from constructing suitable lower and upper dynamics iteratively for (1.3). Effective designs of lower and upper dynamics can then capture the asymptotic behaviors of the coupled systems (1.1). Under different formulations of lower-upper dynam-ics, delay-dependent criteria and delay-independent criteria for synchronization of (1.1) can be derived, respectively. This approach also leads to a network-scale–dependent criterion for synchronization in network systems. The idea of sequential contracting is quite natural in the following sense. One starts from a preliminary attracting set of S, which usually exists from the dissipative property in coupled systems which admit synchronization. We then formulate a criterion for contraction so that the dynamics converge to S through iteration arguments. Such a formulation imposes mild conditions, as the nonlinear terms in the equations are not overmanipulated by linearization or other treatments.

In section 2, we analyze the asymptotic behavior for a scalar equation associated with the difference-differential equation (1.3). The analysis provides a basis for investigating the synchronization of system (1.1). We present the main theorems for system (1.1) of two cou-pled subsystems in section 3. In subsection3.1, we introduce the main conditions imposed on system (1.1). In subsection 3.2, two synchronization theorems for (1.1), one under a delay-dependent criterion and the other under a delay-indelay-dependent criterion, are established succes-sively through constructing two different lower-upper dynamics for (1.3). We then implement

these theorems to establish the synchronization for coupled FitzHugh–Nagumo neurons under nonlinear coupling with discrete-time delay and distribution delay, respectively, in subsec-tion 4.1. The synchronization for a cell-to-cell kinetic model of segmentation clock genes is demonstrated in subsection 4.2. We extend the framework to treat N -cell (unit) systems under circulant coupling in subsection 5.1. In subsection 5.2, we demonstrate this extension in a K-loop neural network. We compare the present approach with the methodologies for studying synchronization in the literature in section 6.

2. Formulation and component estimate. This section is a preparation for the main

theorems in section 3. Recall the difference-differential equation (1.3), ˙

z_i(t) = F_i(x_t, t) + G_i(x_t,y_t, t)− F_i(y_t, t)− G_i(y_t,x_t, t),

where z_i(t) = x_i(t)−y_i(t), and (x(t), y(t)) with x(t) = (x₁(t), . . . , x_n(t)),y(t) = (y₁(t), . . . , y_n(t)), is an arbitrary solution of system (1.1). The key point of the present approach is to analyze the behavior of z_i(t) through a manipulation of (1.3). To this end, we first consider the following scalar delay-differential equation.

We denote by t₀ the initial time and by τ_M ≥ 0 the upper bound of delay magnitude. Let w(t) be a bounded continuous function for t ≥ t₀, and let h : R × R × R → R and ˜

h :C([−τ_M, 0];R)×C([−τ_M, 0];R)×R → R be continuous functions. Let x_t, y_t∈ C([−τ_M, 0];R)

for t ≥ t0, and set x(t + θ) = xt(θ), y(t + θ) = yt(θ) for θ ∈ [−τM, 0]; we assume that x(t) and y(t) eventually enter and then remain in some closed and bounded interval [ˇq, ˆq]; namely, x(t) and y(t) lie in [ˇq, ˆq] for all t ≥ ˜t0, for some ˜t0 ≥ t0. We suppose that z(t) = x(t)− y(t) satisfies the following scalar equation:

(2.1) z(t) = h(x(t), y(t), t) + ˜˙ h(x_t, y_t, t) + w(t), t≥ t₀.

We shall decompose (1.3), for each i, into an equation of the form (2.1), with the spirit of collecting the instantaneous self-feedback terms in h, delayed self-feedback terms in ˜h, and

cross-coupling terms in w. How (2.1) is connected to (1.3) exactly will be seen in section3.1. We impose the following condition on the argument structure of h and ˜h and the boundedness

of ˜h.

Condition (H0). There exist ˆμ, ˇμ, ˇβ, ˆβ ∈ R, ρh > 0, and 0≤ ¯τ ≤ τM such that for each φ, ψ ∈ {ϕ ∈ C([−τ_M, 0];R) : ϕ(θ) ∈ [ˇq, ˆq], θ ∈ [−¯τ, 0]}, the following properties hold for all t≥ t₀: (H0− i) :  ˇ μ≤ h(φ(0), ψ(0), t)/[φ(0) − ψ(0)] ≤ ˆμ, φ(0) − ψ(0) = 0, h(φ(0), ψ(0), t) = 0, φ(0)− ψ(0) = 0,

(H₀− ii) : |˜h(φ, ψ, t)| ≤ ρh, and there exists a τ = τ (φ, ψ, t)∈ [0, ¯τ] such that

 _ˇ

β≤ ˜h(φ, ψ, t)/[φ(−τ) − ψ(−τ)] ≤ ˆβ, φ(−τ) − ψ(−τ) = 0,

˜

h(φ, ψ, t) = 0, φ(−τ) − ψ(−τ) = 0.

Herein, τ is a function of φ, ψ, and t in (H0-ii); φ and ψ work as xt and yt in (2.1), respectively; in particular, φ(0) (resp., ψ(0)) corresponds to x(t) (resp., y(t)), and φ(θ) (resp.,

ψ(θ)) corresponds to x(t + θ) (resp., y(t + θ)) for θ ∈ [−¯τ, 0]. Thus, in (2.1), condition (H₀)

basically indicates that the dynamics of z(t), or ˙z(t), are controlled by z(t) via some upper and lower factors ˆμ and ˇμ, and by z(t− τ) via some upper and lower factors ˆβ and ˇβ. Moreover,

˜

h(x_t, y_t, t), the delay effect contributed from x_t and y_t on ˙z(t), is bounded. Since condition (H₀) is to be imposed for more than a specific (x(t), y(t)), we describe it by general notation

φ and ψ.

The main result in this section asserts that there exists a bounded and closed interval containing zero to which every solution z(t) of (2.1) converges, under some delay-dependent condition. A variant of this formulation leads to the same conclusion (with a different interval) under a delay-independent condition. We note that the notation t0, ˜t0, [ˇq, ˆq], defined before introducing system (2.1), and ˆμ, ˇμ, ˇβ, ˆβ, ρh, ¯τ , in condition (H₀), will be used throughout this section. For all T ≥ t0, we set

|w|max_{(T ) := sup{|w(t)| : t ≥ T }, |w|}max_{(∞) := lim}

T →∞|w|

max_{(T ).}

To capture the dynamics of system (2.1), we define the following functions: ˆ h(ξ) =  (ˆμ + ˆβ)ξ + 3ρh+|w|max(t₀) for ξ ≥ 0, (ˇμ + ˇβ)ξ + 3ρh+|w|max(t₀) for ξ < 0, ˇ h(ξ) =−ˆh(−ξ).

Obviously, ˆμ + ˆβ ≥ ˇμ + ˇβ. If ˆμ + ˆβ < 0, then ˆh(ξ) ≥ ˇh(ξ) for all ξ ∈ R; moreover, ˆh and ˇh are

piecewise linear, are decreasing, and have unique zeros at ˆAh and ˇAh, respectively; see Figure 1. Notably,

(2.2) ˆh( ˇAh) =−ˇh( ˆAh) = (ˆμ + ˇμ + ˆβ + ˇβ)(3ρh+|w|max(t₀))/(ˆμ + ˆβ) > 0.

The following lemma asserts that functions ˆh(·) − ρh and ˇh(·) + ρh provide preliminary upper and lower bounds for the dynamics of (2.1).

Lemma 2.1. Assume that condition (H₀) holds and ˆμ + ˆβ < 0. If z(t) satisfies (2.1), then (2.3) ˇh(z(t)) + ρh ≤ ˙z(t) ≤ ˆh(z(t)) − ρh for all t≥ ˜t0+ ¯τ .

Consequently, there exists a T_x,y ≥ ˜t0+ 2¯τ such that z(t) ∈ [ ˇAh, ˆAh], and | ˙z(t)| < ˆh( ˇAh), for

all t≥ T_x,y− ¯τ.

Proof. Let us verify (2.3). Recall that x(t), y(t)∈ [ˇq, ˆq] for all t ≥ ˜t0. For all t≥ ˜t0+ ¯τ , z(t) = x(t)− y(t) satisfies

˙

z(t) = h(x(t), y(t), t) + ˜h(˜x_t, ˜y_t, t) + w(t) + ˜h(x_t, y_t, t)− ˜h(˜x_t, ˜y_t, t),

where ˜x_t(θ) := x(t), ˜y_t(θ) := y(t) for all θ∈ [−τ_M, 0] are constant in θ. Note that x_t, y_t, ˜x_t, ˜y_t∈ {φ ∈ C([−τM, 0];R) : φ(θ) ∈ [ˇq, ˆq], θ ∈ [−¯τ, 0]}. If x(t) ≥ y(t), then h(x(t), y(t), t) ≤ ˆμ · z(t),

˜

h(˜x_t, ˜y_t, t) ≤ ˆβ · z(t), and ˜h(x_t, y_t, t)− ˜h(˜x_t, ˜y_t, t) ≤ 2ρh by condition (H₀). Consequently, ˙z(t)≤ (ˆμ + ˆβ) · z(t) + 2ρh+|w|max(t0) =: ˆh(z(t))− ρh. On the other hand, if x(t) < y(t), then h(x(t), y(t), t) ≤ ˇμ · z(t), and ˜h(˜x_t, ˜y_t, t) ≤ ˇβ · z(t). Therefore, ˙z(t) ≤ ˆh(z(t)) − ρh. Similar arguments lead to ˇh(z(t)) + ρh ≤ ˙z(t). From (2.3), we obtain ˙z(t)≤ −ρh if z(t) ≥ ˆAh, and

Figure 1. Configurations of functions ˆh, ˇh, ˆh(0)(·, T ), and ˇh(0)(·, T ).

˙z(t) ≥ ρh if z(t) ≤ ˆAh, for t ≥ ˜t₀+ ¯τ . Subsequently, there exists a T_x,y ≥ ˜t₀+ 2¯τ such that z(t)∈ [ ˇAh, ˆAh], which in turn yields | ˙z(t)| < ˆh( ˇAh) for all t≥ T_x,y− ¯τ by (2.3).

Herein, the notation T_x,y indicates its dependence on x(t) and y(t). Now, let us consider the following condition for (2.1).

Condition (A1). ˆμ + ˆβ < 0 and ¯β ¯τ < 3ρh(ˆμ + ˆβ)/[(ˆμ + ˇμ + ˆβ + ˇβ)(3ρh+|w|max_(t

0))], where ¯

β := max{| ˇβ|, | ˆβ|}.

The latter inequality in condition (A1) requires that if the delayed effect ˜h in (2.1) exists, i.e., ¯β = 0, the allowable maximal magnitude of time lag ¯τ should be small enough. From

(2.2), condition (A1) yields ¯β ¯τ ˆh( ˇAh) < 3ρh, and there exists an ε₀> 0 such that

(2.4) β ¯¯τ ˆh( ˇAh) + ε₀ < 3ρh.

For each T ≥ t0, we further introduce the following functions: ˆ h(0)(ξ, T ) =  (ˆμ + ˆβ)ξ + ¯β ¯τ ˆh( ˇAh) +|w|max(T ) + ε0 for ξ ≥ 0, (ˇμ + ˇβ)ξ + ¯β ¯τ ˆh( ˇAh) +|w|max(T ) + ε0 for ξ < 0, ˇ h(0)(ξ, T ) =−ˆh(0)(−ξ, T ). Notably, condition (A1) implies (2.4), and

ˇ

h(ξ) < ˇh(0)(ξ, T ) < ˆh(0)(ξ, T ) < ˆh(ξ) for all ξ∈ R. (2.5)

Let ˇm(0)(T ) (resp., ˆm(0)(T )) be the unique solution of ˇh(0)(·, T ) = 0 (resp., ˆh(0)(·, T ) = 0) lying in interval [ ˇAh, ˆAh]; see Figure1. Notably, ˆm(0)(T ) =− ˇm(0)(T ) > 0, and [− ˆm(0)(T ), ˆm(0)(T )]⊂ [ ˇAh, ˆAh], by (2.5). Recall T_x,y introduced in Lemma 2.1. The following lemma reveals that ˇ

h(0)(·, T ) + ε0 and ˆh(0)(·, T ) − ε0 provide lower and upper bounds finer than ˇh(·) + ρh and ˆ

h(·) − ρh, respectively, for the dynamics of system (2.1), as time gets larger.

Lemma 2.2. Assume that conditions (H₀) and (A1) hold. If z(t) satisfies (2.1), then for

each T ≥ T_x,y, we have

(2.6) hˇ(0)(z(t), T ) + ε₀ ≤ ˙z(t) ≤ ˆh(0)(z(t), T )− ε₀ for all t≥ T. Consequently, z(t) eventually enters and stays afterward in [− ˆm(0)(T ), ˆm(0)_{(T )].}

Proof. By condition (H0), there exist some μ(t) with ˇμ≤ μ(t) ≤ ˆμ, β(t) with ˇβ ≤ β(t) ≤ ˆβ, and τ (t) := τ (x_t, y_t, t)≤ ¯τ, such that the terms h(x(t), y(t), t) and ˜h(x_t, y_t, t) in (2.1) become

h(x(t), y(t), t) = μ(t)[x(t)− y(t)] = μ(t)z(t),

˜

h(x_t, y_t, t) = β(t)[x(t− τ(t)) − y(t − τ(t))] = β(t)z(t − τ(t)).

Thus, (2.1) can be rewritten as follows:

(2.7) ˙z(t) = μ(t)z(t) + β(t)z(t− τ(t)) + w(t). For t≥ T ≥ T_x,y, applying the mean value theorem to (2.7) yields

˙z(t) = μ(t)z(t) + β(t)[z(t)− τ(t) ˙z(s)] + w(t),

where s≥ t − ¯τ ≥ T_x,y− ¯τ; hence | ˙z(s)| < ˆh( ˇAh) by Lemma 2.1. Consequently, if z(t) ≥ 0, then ˙z(t)≤ (ˆμ + ˆβ)z(t) + ¯β¯τˆh( ˇAh) +|w|max(T ) =: ˆh(0)(z(t), T )− ε0; if z(t) < 0, then ˙z(t)≤ (ˇμ + ˇβ)z(t) + ¯β ¯τ ˆh( ˇAh) +|w|max(T ) =: ˆh(0)(z(t), T )− ε₀. Hence, the right-hand inequality of (2.6) is verified. The left-hand one can be treated similarly. Since ˆm(0)(T ) and ˇm(0)_{(T ) are} the unique zeros of ˆh(0)(·, T ) and ˇh(0)(·, T ), respectively, and ˆm(0)(T ) = − ˇm(0)(T ) > 0, we conclude that z(t) eventually enters and stays afterward in [− ˆm(0)(T ), ˆm(0)(T )], as depicted in Figure 1.

Lemmas 2.1and 2.2 demonstrate the formulation of lower and upper bounds for the dy-namics of (2.1) in succession. In the same spirit, we shall formulate finer lower and upper bounds iteratively to capture the asymptotic dynamics of (2.1). Now, let {ε_k}∞_k=1 be a de-creasing sequence with ε1 < ε0, and let εk → 0 as k → ∞. For k ∈ N and T ≥ t0, we define ˆ h(k)(ξ, T ) :=  (ˆμ + ˆβ)ξ + ¯β ¯τ ˆh(k−1)( ˇm(k−1)(T ), T ) +|w|max(T ) + ε_k, ξ ≥ 0, (ˇμ + ˇβ)ξ + ¯β ¯τ ˆh(k−1)( ˇm(k−1)(T ), T ) +|w|max(T ) + ε_k, ξ < 0, ˇ h(k)(ξ, T ) :=−ˆh(k)(−ξ, T ),

where ˇm(k)(T ) ≤ 0 is the unique zero of ˇh(k)(·, T ), and ˆm(k)(T ) = − ˇm(k)(T ) ≥ 0 is the unique zero of ˆh(k)(·, T ). By arguments similar to Lemma 2.3 in [53], it can be shown that under condition (A1), for any fixed T ≥ t₀, {ˆh(k)(·, T )|_{[ ˇAh, ˆAh]}}_k≥1 are uniformly bounded and equicontinuous; in addition, ˆh(k)(·, T ) is decreasing with respect to k. There exists a continuous function ˆh(∞)(·, T ) defined on [ ˇAh, ˆAh] such that

(2.8) hˆ(k)(·, T ) ↓ ˆh(∞)(·, T ) uniformly on [ ˇAh, ˆAh], as k→ ∞,

by the Ascoli–Azela theorem. Since ˆh(k)(·, T ) is decreasing with respect to k, there exists an

m(T ) ∈ R such that ˆm(k)(T )→ m(T ) decreasingly as k → ∞, where ˆm(k)(T ) is the unique

zero of ˆh(k)(·, T ). With (2.8), ˆm(k)(T )→ m(T ), and the continuity of ˆh(k) and ˆh(∞), we can derive that ˆh(∞)(·, T ) is a vertical shift of ˆh(k)(·, T ) and satisfies

(2.9) ˆh(∞)(ξ, T ) = 

(ˆμ + ˆβ)ξ + ¯β ¯τ ˆh(∞)(−m(T ), T ) + |w|max(T ), ξ≥ 0,

(ˇμ + ˇβ)ξ + ¯β ¯τ ˆh(∞)(−m(T ), T ) + |w|max(T ), ξ < 0,

where m(T ) is the unique zero to ˆh(∞)(·, T ). Moreover, from the configuration of ˆh(∞)(·, T ), it follows that

(2.10) 0≤ m(T ) = |w|max(T )/{−ˆμ − ˆβ + ¯β¯τ(ˇμ + ˆμ + ˇβ + ˆβ)}.

The detailed computation for (2.10) is arranged in AppendixA. By (2.8) and that|w|max(T ) decreases with respect to T , we conclude that ˆm(k)(T )→ m(T ) decreasingly, as k → ∞. In addition, there exists an ¯m≥ 0, such that m(T ) → ¯m decreasingly, as T → ∞. Thus,

∩k≥0,T ≥t0[− ˆm(k)(T ), ˆm(k)(T )] =∩T ≥t0[−m(T ), m(T )] = [− ¯m, ¯m].

It can be argued by induction that for arbitrarily fixed T ≥ T_x,y, and n ∈ N, there exists an increasing sequence{T_k}n_k=0 with T_k+1≥ T_k+ ¯τ for k = 0, 1, . . . , n− 1, and T₀ ≥ T + ¯τ, such

that  _ˇ

h(k)(z(t), T ) + ε_k≤ ˙z(t) ≤ ˆh(k)(z(t), T )− ε_k for t≥ T_k+ ¯τ , k = 0, 1, . . . , n− 1, z(t)∈ [− ˆm(k)(T ), ˆm(k)(T )] for t≥ T_k+1, k = 0, 1, . . . , n− 1.

This then leads to the fact that z(t) which satisfies (2.1) eventually enters and remains in [− ˆm(k)(T ), ˆm(k)(T )] for each T ≥ T_x,y and k ∈ N and hence converges to interval [− ¯m, ¯m] as t→ ∞. Based on these arguments, we conclude the following proposition.

Proposition 2.3. Assume that conditions (H₀) and (A1) hold. If z(t) satisfies (2.1), then

z(t) converges to some interval [− ¯m, ¯m] as t→ ∞. Moreover, 0≤ ¯m≤ |w|

max_(∞)

−ˆμ − ˆβ + ¯β¯τ(ˇμ + ˆμ + ˇβ + ˆβ).

The conclusion in this proposition is ¯τ -dependent. By recomposing the upper and lower

functions (see AppendixB) and using arguments similar to those for Proposition2.3, we can derive the following ¯τ -independent conclusion.

Proposition 2.4. If z(t) satisfies (2.1), then z(t) converges to interval [− ˜m, ˜m], as t→ ∞, under condition (H₀) and

condition (A2) : 0≤ ¯β < −ˆμ/[1 + |w|max(t0)/ρh].

Moreover,

0≤ ˜m≤ |w|

max_(∞) −ˆμ − ¯β .

Remark 2.1. When introducing (2.1), x(t) and y(t) are assumed to enter and remain in [ˇq, ˆq] eventually. Such an assumption can be weakened to that x(t) and y(t) converge to [ˇq, ˆq]

for the delay-independent result in Proposition 2.4.

3. Synchronization for coupled system (1.1). We shall derive a delay-dependent criterion and a delay-independent criterion for the global synchronization of system (1.1), based on Propositions 2.3and2.4, respectively. Let (x_t,y_t) be the solution evolved from an arbitrarily fixed initial condition, and let (x₁(t), . . . , x_n(t), y₁(t), . . . , y_n(t)) be the corresponding entire solution for system (1.1). Recall from (1.3) that z_i(t) := x_i(t)− y_i(t) satisfies the following difference-differential system:

˙

z_i(t) = F_i(x_t, t) + G_i(x_t,y_t, t)− F_i(y_t, t)− G_i(y_t,x_t, t)

=: H_i(x_t,y_t, t), i = 1, . . . , n.

(3.1)

Our aim is to show that z_i(t) → 0 for all i = 1, . . . , n, as t → ∞, to conclude the global synchronization for system (1.1). In subsection 3.1, we introduce two basic assumptions on system (1.1). In subsection 3.2, we establish the synchronization theorems for (1.1).

3.1. Dissipative and argument conditions. We make two basic assumptions on system

(1.1). The first is associated with the dissipative property of system (1.1), and the second is related to the argument structure for H_i in (3.1).

Assumption (D). All solutions of system (1.1) eventually enter and then remain in some compact set Q × Q, where Q := [ˇq1, ˆq1]× · · · × [ˇqn, ˆqn]⊂ Rn.

Notably, under assumption (D), all solutions of system (1.1) exist on [t₀,∞). To introduce

the second assumption, we decompose function H_i in (3.1) as

(3.2) H_i(Φ, Ψ, t) = h_i(φ_i(0), ψ_i(0), t) + ˜h_i(φ_i, ψ_i, t) + w_i(Φ, Ψ, t),

where Φ = (φ₁, . . . , φ_n), Ψ = (ψ₁, . . . , ψ_n) ∈ C([−τ_M, 0];Rn). Herein, h_i (resp., ˜h_i) refers to the instantaneous (resp., delayed) part of H_i contributed from φ_i and ψ_i, and w_i collects all cross-coupling terms. Such a decomposition for H_i is always achievable since a trivial case is h_i = ˜h_i ≡ 0 and w_i ≡ H_i. A nontrivial decomposition for the coupled FitzHugh–Nagumo system (4.1) is illustrated in section 4.1. The following second assumption is associated with the argument structure of h_i, ˜h_i and the boundedness of ˜h_i and w_i.

Assumption (H). For i = 1, . . . , n, there exist ˇμ_i, ˆμ_i, ˆβ_i, ˇβ_i ∈ R, ρh_i, ρw_i ≥ 0, ¯μ_ij, ¯β_ij ≥ 0, and

0≤ ¯τ_i, ¯τ_ij ≤ τ_M, j= i, such that for each Φ, Ψ ∈ C_Q:={˜Φ = (φ1, . . . , φn)∈ C([−τM, 0];Rn) : φ_i(θ) ∈ [ˇq_i, ˆq_i], for all θ ∈ [−¯τ_i, 0], i = 1, . . . , n}, the following three properties hold for all t≥ t₀: (H− i) :  ˇ μ_i ≤ h_i(φ_i(0), ψ_i(0), t)/[φ_i(0)− ψ_i(0)]≤ ˆμ_i, φ_i(0)− ψ_i(0)= 0, h_i(φ_i(0), ψ_i(0), t) = 0, φ_i(0)− ψ_i(0) = 0, (H− ii) : |˜h_i(φ_i, ψ_i, t)| ≤ ρh_i and there exists τ_i = τ_i(φ_i, ψ_i, t)∈ [0, ¯τ_i], such that

 _ˇ

β_i ≤ ˜h_i(φ_i, ψ_i, t)/[φ_i(−τ_i)− ψ_i(−τ_i)]≤ ˆβ_i, φ_i(−τ_i)− ψ_i(−τ_i)= 0, ˜

h_i(φ_i, ψ_i, t) = 0, φ_i(−τ_i)− ψ_i(−τ_i) = 0, (H− iii) : |w_i(Φ, Ψ, t)| ≤ ρw_i and there exists τ_ij = τ_ij(Φ, Ψ, t)∈ [0, ¯τ_ij], j = i, such that

|wi(Φ, Ψ, t)| ≤ Σj=i{¯μij|φj(0)− ψj(0)| + ¯βij|φj(−τij)− ψj(−τij)|}.

In practical application, assumption (H) can be realized by suitable manipulation on (3.2) through some estimates, as all solutions of system (1.1) eventually stay in compact setQ × Q

inR2n, under assumption (D). ThusC_Qdepends on the compact setQ in assumption (D), and (H-i) and (H-ii) are multiple-component versions of conditions (H₀-i) and (H₀-ii) in section2, respectively. For later use, we set

(3.3) L¯_ij := ¯μ_ij+ ¯β_ij.

Let us explain the connection of assumption (H) to system (3.1). In assumption (H), Φ (resp., Ψ) plays the role of x_t(resp., y_t), φ_i(0) (resp., ψ_i(0)) of x_i(t) (resp., y_i(t)), and φ_i(θ) (resp.,

ψ_i(θ)) of x_i(t + θ) (resp., y_i(t + θ)) in system (3.1). Assumption (H) basically asserts that ˙z_i(t) in system (3.1) is dominated by z_i(t) via some lower and upper factors ˇμ_i and ˆμ_i (see (H-i)), by z_i(t− τ_i) via some lower and upper factors ˇβ_i and ˆβ_i (see (H-ii)), and by |z_j(s)|, j = i, when s = t and s = t− τ_ij (t− τ_ij is certain uniformly bounded delayed time) via some upper factors (see (H-iii)). Actually, assumption (H) strongly relies on assumption (D) under which ˜

h_i and w_i are bounded on set C_Q, and hence such argument conditions on h_i, ˜h_i, and w_i can be verified by applying the mean value theorem basically. On the other hand, formulating proper forms of h_i, ˜h_i, and w_i for a considered system is also important for assumption (H) to be met.

3.2. Synchronization for system (1.1). In this section, we shall establish the global

synchronization of system (1.1) under assumptions (D) and (H). By (3.2), we can rewrite the difference-differential system (3.1) as follows:

(3.4) z˙_i(t) = h_i(x_i(t), y_i(t), t) + ˜h_i((x_t)_i, (y_t)_i, t) + w_i(t),

where we regard w_i(x_t,y_t, t) as a function of t, i.e., w_i(t) := w_i(x_t,y_t, t), as (x_t,y_t) is the solution evolved from a fixed initial condition, mentioned in section 3.1. Under assumptions (D) and (H), each ith component in (3.4) is in the form of (2.1) and satisfies condition (H₀) with ˇμ = ˇμ_i, ˆμ = ˆμ_i, ρh = ρh_i, ˇβ = ˇβ_i, ˆβ = ˆβ_i, ¯τ = ¯τ_i. Now, let us introduce the multi-component versions of conditions (A1) and (A2).

Condition (S1). ˆμ_i+ ˆβ_i < 0 and ¯β_iτ¯_i < τ_iS for all i = 1, . . . , n, where ¯

β_i:= max{| ˇβ_i|, | ˆβ_i|}, τ_iS:= 3ρ h

i(ˆμi+ ˆβ_i)

(ˆμ_i+ ˇμ_i+ ˆβ_i+ ˇβ_i)(3ρh_i + ρw_i ). Condition (S2). ¯β_i<−ˆμ_i/(1 + ρw_i /ρh_i) for all i = 1, . . . , n.

Note that condition (S1) is delay-dependent, while condition (S2) is delay-independent. Assume that condition (S1) holds; by Proposition 2.3, for each i = 1, . . . , n, there exists an interval I_i:= [− ¯m_i, ¯m_i], to which z_i(t) converges, as t→ ∞; moreover,

(3.5) 0≤ ¯m_i≤ |w_i|max(∞)/η_i,

where

(3.6) η_i:=−ˆμ_i− ˆβ_i+ ¯β_iτ¯_i(ˇμ_i+ ˆμ_i+ ˇβ_i+ ˆβ_i).

The following proposition shows that ¯m_i can be further estimated iteratively.

Proposition 3.1. Assume that condition (S1) holds. Then for each i = 1, . . . , n, there exists a sequence {m(k)_i }∞_k=1 which satisfies

(3.7) m¯_i ≤ m(k)_i = ⎡ ⎣i−1 j=1 ¯ L_ijm(k)_j + n  j=i+1 ¯ L_ijm(k−1)_j ⎤ ⎦  η_i

for k≥ 1, where m(0)_i := ρw_i /η_i and ¯L_ij is defined in (3.3).

Proof. We prove the proposition by induction and sketch the main process. Under

as-sumption (H), |w_i|max_{(∞) ≤ ρ}w

i ; consequently, ¯mi ≤ m(0)i for all i = 1, . . . , n. Assume that m(k)_i in (3.7) have been defined, and hence z_i(t) converges to [−m(k)_i , m(k)_i ], as t → ∞, for

k = 1, . . . , k0− 1, i = 1, . . . , n, and k = k0, i = 1, . . . , − 1 < n. By condition (H-iii), |w(t)| = |w(xt,yt, t)| ≤ Σj={¯μj|zj(t)| + ¯βj|zj(t− τj(xt,yt, t))|}; then, |w_|max(∞) ≤ (Σ−1_j=1L¯_jm_j(k0)+ Σn_j=+1L¯_jm(k_j 0−1)); hence

0≤ ¯m_≤ |w_|max(∞)/η_≤ [Σ−1_j=1L¯_jm_j(k0)+ Σ_j=+1n L¯_jm_j(k0−1)]/η_ =: m(k_ 0).

This completes the proof.

We observe that {m(k)_i | i = 1, 2, . . . , n} in Proposition 3.1 is exactly the Gauss–Seidel iteration for solving the linear system

(3.8) Mv = 0,

where

(3.9) M := D_M− L_M− U_M= [m_ij]_1≤i,j≤n, m_ii= η_i, m_ij =−¯L_ij, for i= j,

and D_M,−L_M, and−U_Mrepresent the diagonal, strictly lower-triangular, and strictly upper-triangular parts of M, respectively; ¯L_ij and η_i are defined in (3.3) and (3.6), respectively. For each i = 1, 2, . . . , n, z_i(t) = x_i(t)− y_i(t) converges to [− ¯m_i, ¯m_i] as t → ∞, and ¯m_i ≤ m(k)_i

for all k. Thereby, the problem of synchronization for system (1.1) reduces to solving the linear problem (3.8). Restated, system (1.1) achieves global synchronization if m(k)_i → 0, as

k→ ∞, for all i = 1, 2, . . . , n. One sufficient condition for the convergence of the Gauss–Seidel

iteration for (3.8) is the strict diagonal-dominance of M, which is straightforward to verify. However, for some systems, such as Example 4.1, such a condition is too strong a criterion for synchronization. On the other hand, it is well known that the necessary and sufficient condition for the convergence of the Gauss–Seidel iteration for (3.8) is that the absolute magnitudes of all eigenvalues of the iteration matrix (D_M− L_M)−1U_M are less than unity; see [27]. Based on such a condition, we obtain the main results in this investigation. The assertion shall be derived by computing the eigenvalues for certain corresponding matrices. Other criteria for the convergence of the Gauss–Seidel method [22,27] may provide conditions which are easier to verify, without computing the eigenvalues of these matrices. We assume that system (1.1) satisfies assumptions (D) and (H).

Theorem 3.2. Assume that condition (S1) holds. Then system (1.1) achieves global

syn-chronization if the Gauss–Seidel iteration for linear system (3.8) converges to zero, the unique

solution of (3.8), or, equivalently,

max

1≤i≤n{|λi| : λi : eigenvalue of (DM− LM)

−1_U

M} < 1.

By Proposition 2.4and arguments similar to those for Proposition 3.1 and Theorem3.2, we can derive the delay-independent criterion for the synchronization of system (1.1).

Theorem 3.3. Assume that condition (S2) holds. Then system (1.1) achieves global

syn-chronization if the Gauss–Seidel iteration for linear system ˜

Mv = 0, (3.10)

˜

M := [ ˜m_ij]1≤i,j≤n, ˜mii=−ˆμi− ¯βi, ˜mij =−¯Lij, f or i= j,

converges to zero, the unique solution of (3.10), or, equivalently, max

1≤i≤n{|λi| : λi: eigenvalue of (DM˜ − LM˜)

−1_U

˜ M} < 1,

where D_M˜, −L_M˜, and −U_M˜ are the diagonal, strictly lower-triangular, and strictly upper-triangular parts of ˜M, respectively.

Remark 3.1. (i) For the delay-independent result in Theorem3.3, assumption (D) can be relaxed to that all solutions of system (1.1) converge toQ×Q, as t → ∞; see Remark2.1. (ii) The contents of matricesM and ˜M actually reflect the structure of the coupling configuration. (iii) Let us translate the notation and theory to the ODE case. Consider z_i(t) = x_i(t)− y_i(t), which satisfies

˙

z_i(t) = F_i(x, t) + G_i(x, y, t) − F_i(y, t) − G_i(y, x, t), = h_i(x_i, y_i, t) + w_i(x, y, t),

where ˜h_i = 0 in (3.2). In assumption (H), C_Q is replaced by Q, (H-i) becomes ˇμ_i ≤

h_i(x_i, y_i, t)/[x_i − y_i] ≤ ˆμ_i if x_i − y_i = 0 and h_i(x_i, y_i, t) = 0 if x_i − y_i = 0, (H-ii) is not needed, and (H-iii) is adjusted to |w_i(x, y, t)| ≤ ρw_i and |w_i(x, y, t)| ≤ Σ_j=iμ¯_ij|x_j− y_j|. The matrix M in (3.9) becomes identical to matrix ˜M in (3.10), and Theorem 3.2 reduces to Theorem 3.3.

The idea of sequential contracting was applied to study the global synchronization and asymptotic phases in a basic neural network with nearest-neighbor coupling in [53]. Therein, each unit of the coupled system is a scalar equation. In this paper, we have extended this idea to coupled systems in the form (1.1), where each unit itself is a system of differential equations and may contain intrinsic delays. Herein, we have established a general framework to accommodate a variety of nonlinearly coupled systems for studying synchronization. Un-der this framework, the problem of establishing synchronization for systems unUn-der delayed and nonlinear coupling was transformed to solving a corresponding linear system of algebraic equations. In the process, we have improved the formulation and analysis so that the con-vergence of the corresponding Gauss–Seidel iteration is determined by the optimal condition (both sufficient and necessary). This has enhanced the applicability of the synchronization theory, as shown in the following sections.

4. Implementation of approach. We shall apply the theory developed in section 3 to establish synchronization for two models. We shall examine assumptions (D) and (H), and condition (S1) or (S2), to apply Theorem 3.2 or 3.3. We illustrate the applications with the classical FitzHugh–Nagumo neuronal model and a representative gene regulation model on the segmentation clock genes in zebrafish in subsections4.1 and 4.2, respectively.

4.1. Coupled FitzHugh–Nagumo neurons. The FitzHugh–Nagumo model was first

sug-gested by FitzHugh in 1961 [19], and its equivalent circuit was created by Nagumo, Arimoto, and Yoshizawa in 1962 [40] to describe a prototype of excitable systems. FitzHugh–Nagumo equations, while modified from the van der Pol equation, capture the essence of the cubic nullcline nature of the voltage-component in the simplified Hodgkin–Huxley equations; see [17].

Let us consider the excitable FitzHugh–Nagumo system coupled with time delay [6],

(4.1) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ˙ x₁(t) =−x3₁(t) + (a + 1)x2₁(t)− ax₁(t)− x₂(t) + cf (y₁(t− τ)), ˙ x2(t) = bx1(t)− γx2(t), ˙ y1(t) =−y13(t) + (a + 1)y12(t)− ay1(t)− y2(t) + cf (x1(t− τ)), ˙ y₂(t) = by₁(t)− γy₂(t),

where a, b, γ > 0, and c > 0 is the coupling strength; the sigmoidal coupling function f lies in the following class:

(4.2) {f ∈ C1: f (0) = 0, δ := f (0) > f (ξ) > 0,|f(ξ)| < ρ for ξ = 0}.

In system (4.1), the individual dynamics are governed by the FitzHugh–Nagumo neuron [6, 19,40]:

(4.3)



˙u(t) =−u3(t) + (a + 1)u2(t)− au(t) − v(t), ˙v(t) = bu(t)− γv(t).

In referring to the notation in (1.1),F = (F1, F2) is now

F₁(Φ, t) =−φ3₁(0) + (a + 1)φ2₁(0)− aφ₁(0)− φ₂(0), (4.4)

F₂(Φ, t) = bφ₁(0)− γφ₂(0); (4.5)

the two subsystems are connected via a sigmoidal coupling, a simplification of synaptic cou-pling, with time delay, i.e., G = (G₁, G₂), and

G1(Φ, Ψ, t) = cf (ψ1(−τ)), (4.6)

G₂(Φ, Ψ, t) = 0, (4.7)

where Φ = (φ₁, φ₂), Ψ = (ψ₁, ψ₂) ∈ C([−τ, 0]; R2). In (4.6), the fixed time delay τ is of discrete-time type. In reality, time delay is likely varying each time an action potential is propagated from neurons, and incorporating a distribution of delays to represent the time lags in some range of values with some associated probability distribution is an alternative formulation [7]. In this situation, the term f (ψ₁(t− τ)) in (4.6) can be modified to

(4.8)

_τ 0

f (ψ1(−σ))K(σ)dσ,

where function K is the kernel of the distribution representing the probability density func-tion of time delay. Equafunc-tion (4.3), a single FitzHugh–Nagumo neuron, can exhibit excitable behavior, in the sense that a small perturbation away from its quiescent state can result in a large excursion of its potential before returning to the quiescent state [24]. It was indicated in [6] that the paradigmatic example of the FitzHugh–Nagumo system in the form of (4.3) does not admit periodic solutions for any parameters a, b, γ and exhibits excitable behaviors clearly for certain parameter ranges, for instance,

(4.9) b > γ2, a b, a  γ;

in particular, the only attractor is in the form of a stable equilibrium at the origin if

(4.10) 4b/γ > (a− 1)2.

The investigation of the behavior of the coupled FitzHugh–Nagumo system, which takes into account time delays in signal transmission, has been a subject of considerable interest. The previous works [6,4,5,30] focus on the stability of the trivial equilibrium and delay-induced or coupling-induced bifurcation, which gives rise to synchronous or asynchronous oscillation. The stable synchronous periodic solution for system (4.1) was investigated in [6]. Through nu-merical simulation, it was shown that the system exhibits global convergence to this periodic solution. However, analytical evidence for this global dynamics has been lacking. In [70], via the method of the Lyapunov functional, synchronization conditions for the system consisting of three FitzHugh–Nagumo neurons with delayed coupling and smooth sigmoidal amplification functions were derived. However, the arguments strongly relied on additional consideration of the instantaneous self-feedback term in the coupling and hence provided a delay-independent criterion. Indeed, the existing analytical tools for studying global dynamics and synchro-nization for neuronal models with nonlinear and delayed coupling are rather limited. Herein, we shall derive a delay-dependent criterion and establish the global synchronization for sys-tem (4.1). Our approach can also establish delay-independent and delay-dependent global synchronization for the model considered in [70].

A nontrivial decomposition in the form of (3.2) for the coupled FitzHugh–Nagumo system (4.1) is formulated as follows: From (4.4)–(4.7),

H₁(Φ, Ψ, t) =−φ3₁(0) + (a + 1)φ2₁(0)− aφ₁(0)− [−ψ₁3(0) + (a + 1)ψ2₁(0)− aψ₁(0)] + c[f (ψ₁(−τ)) − f(φ₁(−τ))] − [φ₂(0)− ψ₂(0)]; consequently, we set h₁(φ₁(0), ψ₁(0), t) = p(φ₁(0))− p(ψ₁(0)), (4.11) ˜ h₁(φ₁, ψ₁, t) = c[f (ψ₁(−τ)) − f(φ₁(−τ))], (4.12) w₁(Φ, Ψ, t) =−[φ₂(0)− ψ₂(0)], (4.13)

where p(ξ) :=−ξ3+ (a + 1)ξ2− aξ. On the other hand, from

(4.14) H₂(Φ, Ψ, t) =−γ[φ₂(0)− ψ₂(0)] + b[φ₁(0)− ψ₁(0)],

we have h2(φ2(0), ψ2(0), t) =−γ[φ2(0)− ψ2(0)], (4.15) ˜ h₂(φ₂, ψ₂, t)≡ 0, (4.16) w₂(Φ, Ψ, t) = b[φ₁(0)− ψ₁(0)]. (4.17)

Notably, by the mean value theorem, h1 in (4.11) and ˜h1 in (4.12) can be written in the following form: h1(φ1(0), ψ1(0), t) = [−3s12+ 2(a + 1)s1− a][φ1(0)− ψ1(0)], (4.18) ˜ h1(φ1, ψ1, t) =−cf (s2)[φ1(−τ) − ψ1(−τ)] (4.19)

for some s₁ between φ₁(0) and ψ₁(0), and s₂ between φ₁(−τ) and ψ₁(−τ). In observing (4.11)–(4.19), we see that h_i (resp., ˜h_i) can be transformed into a multiple of φ_i(0)− ψ_i(0) (resp., φ_i(−τ) − ψ_i(−τ)), and the ratio can be further estimated. Notice that the terms

−3s2

1+ 2(a + 1)s1− a in (4.18) and −cf (s2) in (4.19) are bounded when s1, s2 are restricted to some compact set in R. On the other hand, roughly speaking, w_i can be transformed into a linear combination of φ_j(·) − ψ_j(·), j = i.

Now, we show that system (4.1) satisfies assumptions (D) and (H). We define, for k∈ N, (4.20) P(k)(ξ) :=−ξ4+ (a + 1)ξ3− aξ2+|cρ(k−1)ξ|, where ρ(0) := ρ, and ρ(k):= max{|f(ξ)| : ξ ∈ [−γ2_{+ b¯q}(k)_/γ,_γ2_{+ b¯q}(k)_/γ]}, (4.21) ¯ q(k):= max{|ξ| : P(k)(ξ) = 0}. (4.22)

Herein, the parameters a, b, c, γ, ρ and function f were introduced in (4.1) and (4.2). Lemma 4.1. All solutions of system (4.1) eventually enter and then remain in ˜Q(k)× ˜Q(k), for each k ∈ N, where

˜

Q(k):= [−γ2_{+ b¯q}(k)_/γ,_γ2_{+ b¯q}(k)_{/γ]× [−b¯q}(k)_{/γ, b¯q}(k)_/γ].

The proof of Lemma4.1is arranged in AppendixC. Actually, ¯q(k) are well defined for all

k∈ N and are strictly decreasing with respect to k. Subsequently, for larger k, ˜Q(k)provides a smaller attracting region for the dynamics of system (4.1) and hence relaxes the conditions for our synchronization formulation. Throughout this subsection, we consider that system (4.1) satisfies assumption (D) with Q = ˜Q(k)_{=: Q}∗ _{for some fixed k. In some cases (see Example} 4.1), one does need larger k to meet the synchronization criterion. Accordingly, the evolutions for each subsystem in (4.1) will eventually enter and remain in the set:

C∗

Q ={Φ = (φ1, φ2)∈ C([−τ, 0]; R2) : φi(θ)∈ [−q_i∗, q_i∗], i = 1, 2, θ∈ [−τ, 0]}, where

q₁∗:=γ2_{+ b¯q}(k)_{/γ, q}∗

2 := b¯q(k)/γ,

and ¯q(k)is defined in (4.22). Below, let us show that system (4.1) actually satisfies assumption (H). For all t≥ t₀, and Φ = (φ₁, φ₂), Ψ = (ψ₁, ψ₂)∈ C_Q∗,

−q_i∗≤ φi(θ), ψ_i(θ)≤ q_i∗ for θ∈ [−τ, 0], i = 1, 2. (4.23)

Accordingly, by the definitions of ˜h_i and w_i in (4.12), (4.16) and (4.13), (4.17) for system (4.1), we obtain

|˜h1(φ1, ψ1, t)| ≤ 2cMf, |˜h2(φ2, ψ2, t)| = 0, |w1(Φ, Ψ, t)| ≤ 2q2∗, |w2(Φ, Ψ, t)| ≤ 2bq∗1, where

(4.24) M_f := max{|f(ξ)| : ξ ∈ [−q₁∗, q₁∗]}.

This yields the boundedness of ˜h_i and w_i in assumption (H). The argument conditions for functions h_i, ˜h_i, and w_i formulated in (4.13) and (4.15)–(4.19) can then be confirmed. Let us examine these conditions for h₁ and ˜h₁, as the other cases are simpler. Note that the terms

s_i, i = 1, 2, in h₁ and ˜h₁ defined in (4.18) and (4.19) both satisfy−q₁∗ ≤ s_i≤ q∗₁ due to (4.23). It follows from a direct computation that

λ≤ −3s2₁+ 2(a + 1)s₁− a ≤ (a2− a + 1)/3, d_f ≤ f (s2)≤ δ, where λ :=−3(q₁∗)2− 2(a + 1)q₁∗− a, (4.25) d_f := min{f (ξ) : ξ∈ [−q∗₁, q∗₁]} > 0. (4.26)

Consequently, h1 in (4.18) and ˜h1 in (4.19) satisfy, respectively,

λ≤ h₁(φ₁(0), ψ₁(0), t)/[φ₁(0)− ψ₁(0)]≤ (a2− a + 1)/3 if φ₁(0)− ψ₁(0)= 0,

−cδ ≤ ˜h1(φ1, ψ1, t)/[φ1(−τ) − ψ1(−τ)] ≤ −cdf if φ1(−τ) − ψ1(−τ) = 0. From these arguments, we conclude the following lemma.

Lemma 4.2. System (4.1) satisfies assumption (H) with ˇμ₁ = λ, ˆμ₁ = (a2 − a + 1)/3, ˇ

μ2 = ˆμ2 = −γ, ˇβ1 = −cδ, ˆβ1 = −cdf, ¯τ1 = τ , ρ1h = 2cMf, ˜h2 ≡ 0, ρw1 = 2q∗2, ρw2 = 2bq1∗,

¯

μ₁₂ = 1, ¯μ₂₁ = b, and ¯β₁₂= ¯β₂₁= 0, where M_f, λ, and d_f are defined in (4.24), (4.25), and (4.26), respectively.

Assumptions (D) and (H) are thus satisfied for system (4.1) by Lemmas 4.1 and 4.2. By applying Theorem 3.2, we derive a delay-dependent criterion for synchronization of system (4.1).

Theorem 4.3. System (4.1), the two FitzHugh–Nagumo neurons under delayed sigmoidal

coupling, achieves global synchronization if

(4.27) c > [b/γ + (a2− a + 1)/3]/d_f ≥ 0 and τ < min{τ₁F, τ₂F},

where τ₁F := 3Mf[(a 2_{− a + 1)/3 − cd} f] δ[(a2_{− a + 1)/3 + λ − c(d} f+ δ)](3cMf + q∗2) , τ₂F := b/γ− cdf + (a 2_{− a + 1)/3} cδ[(a2_{− a + 1)/3 + λ − c(d} f + δ)],

and M_f, λ, and d_f are defined in (4.24), (4.25), and (4.26), respectively.

Proof. Notice that (4.27) implies c > (a2− a + 1)/(3d_f) and τ < τ₁F, which in turn lead to meeting condition (S1). Moreover, the corresponding matrices in (3.9) are

M =  η₁ −¯L₁₂ −¯L21 η2  , D_M=  η₁ 0 0 η2  , L_M=  0 0 ¯ L21 0  , U_M =  0 L¯₁₂ 0 0  ,

where η1 := −(a2− a + 1)/3 + cdf + cδτ [(a2 − a + 1)/3 + λ − c(df + δ)], η2 := γ ¯L12 = 1, ¯

L₂₁= b. A direct computation reveals that the corresponding matrix (D_M− L_M)−1U_M =



0 L¯₁₂/η₁

0 L¯₁₂L¯₂₁/(η₁η₂) 

admits the eigenvalues 0 and ¯L₁₂L¯₂₁/(η₁η₂) = b/(η₁η₂). We obtain b/(η₁η₂) < 1 under (4.27) using c > [b/γ + (a2 − a + 1)/3]/d_f and τ < τ₂F. The assertion thus follows from Theorem 3.2.

The present approach can also be applied to the case of distribution delay in (4.8). Now ˜ h1 is modified to (4.28) h˜1(φ1, ψ1, t) = c _τM 0 [f (ψ₁(−σ)) − f(φ₁(−σ))]K(σ)dσ,

where functionK is the kernel of the distribution. One of the commonly used distributions is the uniform distribution: for some τmin_{> 0, > 0,}

(4.29) K(σ) := ⎧ ⎨ ⎩ 0 if 0≤ σ < τmin, 1/ if τmin≤ σ ≤ τmin+ , 0 if τmin+ < σ≤ τ_M. By the definition of K, we obtain

˜ h₁(φ₁, ψ₁, t) = c  _τmin_+ τmin [f (ψ1(−σ)) − f(φ1(−σ))]dσ  =−cf(ς)[φ1(−s) − ψ1(−s)]

for some s∈ (τmin, τmin+ ), and ς between φ₁(−s) and ψ₁(−s). Thus, −cδ ≤ −cf(ς)≤ −cd_f. Basically, the arguments are valid for kernel K with compact support. Accordingly, we can verify that the assertions in Lemmas 4.1and 4.2hold, but with ¯τ1 = τmin+ instead. In the following, (4.1) denotes system (4.1) with the coupling f (φ₁(−τ)) replaced by the distribution

delay (4.8) with kernel K in (4.29). Indeed, Lemma 4.1 also holds for system (4.1) . By arguments similar to those in Theorem 4.3, we conclude the following theorem.

Theorem 4.4. System (4.1) achieves global synchronization, provided

c > [b/γ + (a2− a + 1)/3]/d_f ≥ 0 and τmin+ < min{τ₁F, τ₂F},

where d_f, τ₁F, and τ₂F are defined as in Theorem 4.3.

It was derived in [6] that coupled system (4.1) with no time delay (i.e., τ = 0) undergoes a supercritical Hopf bifurcation at c = c₀ := a + γ, which gives rise to a stable synchronous periodic solution. Numerical simulation shows that the oscillation remains a global attractor of the system for a large range of c greater than c0. Here, we note that f in system (4.1) is merely required to satisfy f (0) = 0 and f (0) > 0 for the local dynamics derived from the bifurcation analysis in [6], but we require f to be bounded for the consideration of global dynamics herein. The following example demonstrates that, under our synchronization framework, system (4.1) with parameters a, b, γ satisfying (4.9), (4.10), and τ = 0 admits global synchronization with stable synchronous oscillation as c is larger than and near c₀. This gives an analytical support to the numerical finding in [6]. By considering τ as a bifurcation parameter, and with fixed parameters a, b, γ, c, this example illustrates that the stable oscillation sustains as τ is small enough so that no further bifurcation occurs. The system loses synchrony as τ is larger than a certain critical bifurcation value and yields to a stable antiphase periodic solution.

Example 4.1. Consider (4.1) with a = 0.5, b = 0.00126, γ = 0.02, and f (ξ) = 5 tanh(0.2ξ). Choose Q∗ = ˜Q(2). This system with τ = 0 undergoes a supercritical Hopf bifurcation at

c = c₀ = a + γ = 0.52, which gives rise to a stable synchronous periodic solution. Let us consider the system with fixed c = 0.52001 which is slightly larger than c₀. In this situation, by the bifurcation analysis with respect to τ in [6] (cf. Figure 3 in [6]), the system undergoes a subcritical Hopf bifurcation at the first critical value τ_1,−0 ≈ 0.0008, and a supercritical Hopf bi-furcation at the second critical value τ_2,+0 ≈ 19.505, where a stable antiphase periodic solution emerges. On the other hand, a direct computation gives [b/γ + (a2− a + 1)/3]/d_f ≈ 0.5049,

τ₁F ≈ 0.0027, τ₂F ≈ 0.0004; the system satisfies (4.27) in Theorem 4.3 and hence achieves global synchronization if τ < min{τ₁F, τ₂F} ≈ 0.0004. Here, we note that the synchronization

criterion in Theorem 4.3 does not hold if we choose Q∗ = ˜Q(1) instead. Figure 2 illustrates that the system with τ = 0 and 0.0002, respectively, which are smaller than min{τ₁F, τ₂F, τ_1,−0 }

near 0.0004, admits a stable synchronous oscillation. Figure 3 demonstrates that the system with τ = 20, which is slightly larger than τ_2,+0 , admits a stable antiphase oscillation.

Remark 4.1. The present framework can also accommodate coupled FitzHugh–Nagumo

systems under diffusive coupling, such as system (4.1) with the coupling terms cf (y1(t− τ))

and cf (x₁(t− τ)) replaced by c[y₁(t− τ) − x₁(t)] and c[x₁(t− τ) − y₁(t)]. To apply the present synchronization theories, one needs to examine assumption (D), say, via the approach in [41,58]. Assumption (H) can then be verified subsequently.

4.2. Cell-to-cell kinetic model. In this subsection, we consider a cell-to-cell model on the

kinetics of the segmentation clock genes in zebrafish, proposed by Uriu, Morishita, and Iwasa

9000 9200 9400 9600 9800 10000 −8 −6 −4 −2 0 2 4 6 8x 10 −4 Delay τ=0.0002 time t solution x 1 ,x2 ,y1 ,y2 x 1(t) y₁(t) x₂(t) y 2(t) 9000 9200 9400 9600 9800 10000 −8 −6 −4 −2 0 2 4 6 8x 10 −4 _{Delay τ=0} time t solution x 1 ,x2 ,y1 ,y2 x₁(t) y 1(t) x 2(t) y₂(t) 9000 9200 9400 9600 9800 10000 −8 −6 −4 −2 0 2 4 6 8x 10 −4 _{Delay τ=0.0002} time t solution x 1 ,x2 ,y1 ,y2 x 1(t) y 1(t) x₂(t) y₂(t)

Figure 2.Time series of the solution of system (4.1), evolved from (0.001+t, −0.002·t, 0.001·sin t, −0.001·t)

at initial time t0 = 0, tends to a synchronous (in-phase) oscillation. Here, a = 0.5, b = 0.00126, γ = 0.02,

c = 0.52001, f(ξ) = 5 tanh(0.2ξ), and τ = 0, τ = 0.0002, respectively.

9500 9600 9700 9800 9900 10000 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1x 10 −3 _{Delay τ=20} time t solution x 1 ,y1 x 1(t) y₁(t) 9500 9600 9700 9800 9900 10000 −1 0 1x 10 −4 _{Delay τ=20} time t solution x 2 ,y2 x₂(t) y 2(t) 9500 9600 9700 9800 9900 10000 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1x 10 −3 _{Delay τ=20} time t solution x 1 ,y1 x₁(t) y 1(t) 9500 9600 9700 9800 9900 10000 −1 0 1x 10 −4 _{Delay τ=20} time t solution x 2 ,y2 x 2(t) y 2(t)

Figure 3. Time series of the solution of system (4.1), evolved from (−0.0001, −0.0002, 0.0001, 0.0002) at

initial time t0 = 0, tends to an antiphase oscillation, with a = 0.5, b = 0.00126, γ = 0.02, c = 0.52001,

f(ξ) = 5 tanh(0.2ξ), and τ = 20. in [62,63]: (4.30) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˙ x₁(t) = g_H(x₃(t), y₄(t))− f₁(x₁(t)), ˙ x2(t) = ν3x1(t)− f2(x2(t)), ˙ x₃(t) = ν₅x₂(t)− f₃(x₃(t)), ˙ x₄(t) = g_D(x₃(t))− f₄(x₄(t)), ˙ y1(t) = gH(y3(t), x4(t))− f1(y1(t)), ˙ y₂(t) = ν₃y₁(t)− f₂(y₂(t)), ˙ y₃(t) = ν₅y₂(t)− f₃(y₃(t)), ˙ y4(t) = gD(y3(t))− f4(y4(t)).