From approximate synchronization to identical synchronization in coupled systems

We establish a framework to investigate approximate synchronization of coupled systems under general coupling schemes. The units comprising the coupled systems may be nonidentical and the coupling functions are nonlinear with delays. Both delay-dependent and delay-independent criteria for approximate synchronization are derived, based on an approach termed sequential contracting. It is explored and elucidated that the synchronization error, the distance between the asymptotic state and the synchronous set, decreases with decreasing difference between subsystems, difference between the row sums of connection matrix, and difference of coupling time delays between different units. This error vanishes when these factors decay to zero, and approximate synchronization becomes identical synchronization for the coupled system comprising identical subsystems and connection matrix with identical row sums, and with identical coupling delays. The application of the present theory to nonlinearly coupled heterogeneous FitzHugh-Nagumo neurons is illustrated. We extend the analysis to study approximate synchronization and asymptotic synchronization for coupled Lorenz systems and show that for some coupling schemes, the synchronization error decreases as the coupling strength increases, whereas in another case, the error remains at a substantial level for large coupling strength.


1.
Introduction. Synchronization has appeared as a common coherent behavior in many biological and physical systems [29,47,50,62]. Complex networks of coupled nodes, in which each node is a dynamical system, have been widely exploited to model many complex behaviors in sciences, engineering, and society, and have attracted much attention in recent decades, see [4,5,61] and the references therein. In addition, there is a large number of studies on synchronization of coupled neural networks because of its potential engineering applications such as secret communication [12,57], pattern recognition [21], and parallel image processing [28].
The collective behavior of a coupled system is determined by not only the dynamics of individual subsystems, but also the coupling structure between subsystems. A network of coupled systems can exhibit a variety of interesting behaviors that are qualitatively different from their isolated behaviors. In addition, coupling time delays or communication delays can modify the collective behaviors. Indeed, time delays have been incorporated into neural network modeling and artificial systems [8,9,11,14,33,40,65].
To investigate synchronization phenomena of coupled systems and complex networks, various models and scenarios have been studied, including complete or identical synchronization [30,33,45,46,54,58,63,64,65], phase synchronization [43,44], lag synchronization [43], partial synchronization [56], generalized synchronization [26,34], and almost synchronization [15]. Among these works, identical synchronization of homogeneous systems consisting of identical units is the simplest form of synchronization, and has been studied extensively. However, due to various possible deviations, natural imperfection, and wider consideration such as real-world complex networks, heterogeneous systems which consist of nonidentical units are more practical, see for instance [37,50,60] and the references therein. For heterogeneous systems, it is natural to relax the notion of identical synchronization to more general sense. The description of synchronous behaviors for heterogeneous systems is more complicated, and the formulations and treatments for identical synchronization in homogeneous systems may no longer be valid. Recently, there were some papers focusing on global dynamics of synchronization for heterogeneous networks. Several notions centering around nonidentical synchronization were considered therein, according to the respective approach employed. They include the so-called approximate synchronization, quasisynchronization, bounded synchronization, and practical synchronization. In [18], by introducing a virtual target trajectory, some conditions of quasi-synchronization were derived, and synchronization error was estimated for heterogeneous networks under linear diffusive coupling. In [66], bounded synchronization of some dynamical networks was established by introducing the average dynamics of all nodes. A criterion of global synchronization subject to bounded maximum state deviation between nodes was then derived. Bounded synchronization, in a sense different from the one in [66], for asymmetrical coupled networks was investigated in [59]. In [41,42], an approach was developed to analyze approximate synchronization of a coupled system comprising two nonidentical systems of nonlinear differential equations, where a uniform estimate of the synchronization error with respect to the parameters of the systems was derived. Chaotic approximate synchronization between nonidentical master-slave systems was addressed in [20,24,52]. Applications of control schemes to studying approximate synchronization in heterogeneous networks were shown in [19,25,48]. The studies on practical synchronization of coupled nonidentical systems with linearly diffusive couplings can be found in [36,38]. Frequency synchronization and phase synchronization of phase models, such as the Kuramoto model, have been investigated in [16,22,53].
To summarize, the equations considered in the studies of approximate synchronization originated from motivations in physics, networks, or engineering, largely carry linear, and/or diffusive (or akin to diffusive) coupling functions. This can be seen from the above-mentioned papers [18,24,41,52,59,66] for the ones without delays, and [19,20] for the ones with delays. On the other hand, depending on the methodologies employed, some of those investigations conclude local dynamics, whereas others established global dynamics. In particular, the works which analyze the stability of synchronous manifold or near-synchronous manifold via linearization and the ones involving master stability function conclude local dynamics, see also [1,2,49,51]. On the other hand, Lyapunov function theory and ideas developed from this theory are typical tools for concluding global synchronization. In addition, the criteria for synchronization in some of these papers, such as [59,66], involve certain differential matrix inequalities, and a treatment is needed to examine the criteria.
Another direction of investigations on synchronization comes from biological perspectives, where systems of equations that are used to model the biological processes are highly nonlinear. For example, Hill-type functions are used to describe the repressions of transcriptions in gene regulation. Such functions appear in the coupling terms to describe interactions among cells or proteins in those models [27,55]. In addition, transcriptional and translational time delays, or transmission delays are often taken into account in these models. This makes complicated the coupled systems modeling these biological processes, and new mathematical treatments are required, see [3,10].
From the literature survey, it is tempting to develop a framework which can accommodate the basic content of approximate synchronization for nonlinearly coupled systems and heterogeneous complex networks originated form various motivations. And this framework should be able to cover general forms of coupling with delays. It is also important to know what factors contribute to synchronization error and to identify whether the cause of error is resulted from the differences between individual units or coupling structure. The goal of this work is to establish such a general framework.
In this paper, we shall adopt the term approximate synchronization to describe the common nonidentical synchronization of heterogeneous coupled systems and network systems with a bound on the synchronization error. Let us consider N subsystems which may be nonidentical: where the dynamics of the ith subsystem is governed by F i = (F i,1 , . . . , F i,K ) which is a continuously differentiable function of x t i ∈ C([−τ M , 0]; R K ), and x t i (θ) = x i (t + θ) for θ ∈ [−τ M , 0]. Herein, we denote by C([−τ M , 0]; R n ) the set of continuous functions from [−τ M , 0] to R n , for an n ∈ N, and consider time delays in the range [0, τ M ], τ M ≥ 0. A general coupled network system comprising these N subsystems can be described aṡ where x i (t) = (x i,1 (t), . . . , x i,K (t)) ∈ R K . Herein, coupling functions G ij are smooth functions, each ω ij (t) is a bounded smooth function of t, which represents the time-dependent connection weight from node j to node i, and W (t) := [ω ij (t)] N ×N is the connection matrix. We denote by (x t 1 , . . . , x t N ) the evolution of system (2) at time t from the initial condition (x t0 1 , . . . , x t0 N ) in C([−τ M , 0], R KN ) at t = t 0 , and (x 1 (t), . . . , x N (t)) is the corresponding solution of system (2). System (2) reduces to ordinary differential equations when τ M = 0.
While our approach can treat system (2), the formulation is very cumbersome. So in this paper, we illustrate our approach by considering system (2) , ω ij (t) ≡ ω ij , and function G ij = (G ij,1 , . . . , G ij,K ) of the following form: for k ∈ K := {1, . . . , K}, where G k : R → R are continuously differentiable and nondecreasing functions, and τ ij ∈ [0, τ M ]. In this situation, the component form (the kth component in the ith subsystem) of system (2) becomeṡ System in the form of (4) is still very general, and is commonly considered in the studies of synchronization problems in the literature. For example, the systems that were studied in [18,30,33,36,41,54,58,59,63,65,66] are in the form of (4), or similar to (4). Usually, network system (2) or (4) is called homogeneous if it comprises identical subsystems, i.e., F i ≡ F j , for all i, j ∈ N , and called heterogeneous otherwise. It is important to note that for homogeneous systems, the synchronous set S : and where κ i is the ith row sum of the connection matrix: Herein, δ τ measures the difference of coupling delays between different components, and δ κ measures the difference between row sums of connection matrix. Notably, δ τ = 0 if system (4) has identical coupling delays or does not admit coupling delays; δ κ = 0 if the connection matrix has identical row sums. For homogeneous systems without coupling delays or with identical coupling delay, there have been a large amount of investigations on global identical synchronization, especially in the diffusive coupling case: κ i = 0 for all i ∈ N , see [30,33,45,58,63,65] and the references therein. Herein, system (4) is said to attain global identical synchronization if for every solution (x 1 (t), . . . , x N (t)) of system (4). There were also extensive studies on local synchronization, where the stability of synchronous set S is discussed, see [23,33,39,58]. Basically, the invariance of synchronous set under the flow generated by system (4) is a prerequisite for identical synchronization of (4). However, for a heterogeneous system (4) or a homogeneous system (4) without condition (5), or without (6), the synchronous set is not invariant in general, and thus the realization of identical synchronization is unpractical. We note that for coupled systems with mismatched parameters, some studies on the stability of near-synchronous states have been reported in [1,2,49,51]. Another interesting point is that the asymptotic collective dynamics for (4) is distinct from the ones for individual subsystems when isolated, if some κ i is nonzero. This leads to the fact that there can be more varieties for the collective behaviors of non-diffusively coupled systems than the ones of diffusively coupled systems.
Let δF ij := F i − F j , and where · is certain norm on R K . Then δ F := max{ δF ij : i, j ∈ N } measures the maximal difference between subsystems. In this paper, the coupled system (4) is said to attain global (ε-)approximate synchronization if there exists an ε ≥ 0 such that lim sup for every solution (x 1 (t), . . . , x N (t)) of system (4) (i.e., in the sense of global dynamics). The expression in (10) indicates an estimation of the synchronization error, and ε herein stands for a bound of the error. Moreover, ε = 0 corresponds to identical synchronization. Certainly, this definition with (10) will be meaningful in the sense of approximate synchronization if ε decreases with decreasing inhomogeneity of the system, i.e., ε decreases as δ F decreases, and ε → 0 as δ F → 0. We shall show that this is indeed the case if the row sums of connection matrix and coupling delays are identical, i.e., (5) and (6) hold. When the row sums of connection matrix or coupling delays are not identical, network connection may also contribute to the error ε. This is quite natural, but has not been elucidated mathematically, likely because the existing approaches are insufficient for such a task. It will be revealed that ε decreases as each of δ F , δ κ , and δ τ decreases, and ε → 0 as these factors decay to zero. Our results are much more thorough than the related studies in the literature. For instance, in [20,24,25], a bound of synchronization error similar to (10) was formulated for master-slave systems with parameter mismatches in subsystems. The estimated bound was shown to depend on the magnitude of parameter mismatches. The present approach can also establish approximate synchronization if (10) is replaced by lim sup t→∞ x i (t) − x j (t) ≤ ε for all i, j ∈ N , cf. [59], or lim sup t→∞ x i (t) − s(t) ≤ ε for all i ∈ N , where s(t) := 1 N j∈N x j (t), cf. [66]. If there is a coupling strength c considered in the coupling term, say if the summation term in (4) is modified to one may also discuss asymptotic synchronization. System (4) is said to attain asymptotic synchronization if it attains ε-approximate synchronization and ε → 0 as c → ∞. With the invariant manifold theory, such asymptotic synchronization has been discussed by Hale [17]. However, it is not true that approximate synchronization always yields asymptotic synchronization. With the present framework of approximate synchronization, we are able to explore this interesting issue.
When system (4) is homogeneous, i.e., F i = F j for all i, j ∈ N , and (5) and (6) hold, the formulation in this investigation reduces to the one for identical synchronization. Recently, an iteration argument named "sequential contracting" has been applied to study identical synchronization to conclude global dynamics in [45,46,54]. More precisely, the framework for identical synchronization of coupled systems in a general form was established in [45], via the sequential contracting approach. When considering a coupled system with connection in additive form similar to (4), the condition of circulant coupling in [45] was lifted in [54]. This was completed by performing a suitable transformation on the coupling terms in the difference-differential equations associated with the coupled system. In this paper, to take into account distinct subsystems, nonidentical row sums of connection matrix, and variant coupling delays between different components, a new formulation on the associated difference-differential equations is developed. In particular, we are able to derive an efficacious decomposition on the difference-differential equations, which allows us to extract and collect the above-mentioned inhomogeneous factors from the decomposition. Combining with such a formulation, we then perform a more general transformation on the coupling terms in the difference-differential equations.
Actually, in the somehow limited publications on approximate synchronization of heterogeneous system (4), the couplings are commonly linear and diffusive (i.e., G k is linear for all k and κ i = 0 for all i) and without coupling delays (i.e., τ M = 0), cf. [18,36,41,59,66]. To the best of our knowledge, there has not been an analytical result on approximate synchronization for a heterogeneous system (4), with nonlinear and delayed coupling, in the sense of (10) or similar to (10). We remark that basically, a criterion is needed to attain synchronization, even for a coupled system comprising identical subsystems. Thus, deriving a criterion for approximate synchronization is also the key objective in our framework, in addition to indicating how synchronization error depends on the difference between subsystems and on the coupling structure.
Dissipativeness which implies the existence of a globally attracting set, is a basic property often imposed on the systems considered for synchronization and other asymptotic behaviors [17,30,45,54]. It was utilized directly or indirectly in various methodologies for synchronization theory [13,18,19,20,24,48,52]. In those papers, the synchronization criteria involve the ranges or bounds of solutions of the systems, and the dissipative property are actually utilized, although this might not have been spelled out. In fact, dissipative property for coupled systems is a highly nontrivial issue, and was a chief concern in [17,41,42]. The idea of sequential contracting [45,54] is to start with a preliminary estimation on the attracting region, and then make use of the dissipative property inside the coupled systems via an iterative argument. We therefore make the following assumption: Assumption (D): All solutions of system (4) eventually enter, and then remain in some compact set Q N := Q × . . . × Q ⊂ R N K , where Under assumption (D), every solution (x 1 (t), . . . , x N (t)) of system (4) exists on The remainder of this paper is organized as follows. In Section 2, we present the sequential-contracting framework to derive both delay-dependent and delayindependent criteria for approximate synchronization of system (4). In Section 3, with the criteria derived in Section 2, we study the approximate synchronization of nonlinearly and nondiffusively coupled FitzHugh-Nagumo neurons. We discuss approximate synchronization and asymptotic synchronization for coupled Lorenz equations in Section 4. This paper then ends with a conclusion. Some technical parts are arranged in the Appendices.
2. Approximate synchronization via sequential contracting. Our approach to establishing the approximate synchronization of system (4) relies on a new setting of attracting set, an efficacious decomposition of the associated differencedifferential system, and the sequential contracting technique.
Recall that the evolution (x t 1 , . . . , x t N ) of (4) is an element in C([−τ M , 0]; R N K ), and (x 1 (t), . . . , x N (t)) is the corresponding solution of (4), where x i (t) = (x i,1 (t), . . . , for θ ∈ [−τ M , 0]. Let us denote by (x t i,k ) ′ the derivative of function x t i,k . It follows from (11) that and thus (x i,k ) ′ andẋ i,k both represent the derivative of function x i,k . Notably, τ M ≥ 0 is the upper bound of delay magnitude, and τ M = 0 reduces to the ODE case. In the following expressions, we take (Φ 1 , . . . , , and it indicates the evolution of (4) when (Φ 1 , . . . , respectively, when they represent the components of a solution of (4). In addition, we take Ξ = (ξ 1 , . . . , ξ K ) as a general element in R K .
Via Q and [q k ,q k ] defined in assumption (D) and the connection matrix W = [ω ij ] N ×N , we define the following quantities successively: We further define the following subset of C([−τ M , 0]; R N K ): where Remark 1. Assumption (D) was also imposed to study the identical synchronization of homogeneous coupled systems in [45,54]. Therein, a notion similar to C Q was adopted, and can be described by Compared to the set in (20), C Q in (19) further considers boundedness of φ ′ k : |φ ′ k (·)| ≤ Λ k , where Λ k is determined by not only set Q but also the individual subsystems F i,k , connection matrix [ω ij ], and coupling function G k , as shown in (13)- (17). Such a further consideration of φ ′ k in C Q is to measure the synchronization error caused by the heterogeneous coupling delays, see Remark 3 below.
, for all t ≥t, i ∈ N , and k ∈ K.
For later use, we denote the index set: The ith individual subsystem of (4) is defined by for Ξ ∈ R K . Such (δF i,1 , . . . , δF i,K ) depicts the difference between two adjacent subsystems.
Let (x 1 (t), . . . , x N (t)) be an arbitrary solution of system (4), with x i (t) = (x i,1 (t), . . . , x i,K (t)), i ∈ N . Set z i (t) = x i (t) − x i+1 (t), where z i (t) = (z i,1 (t), . . . , z i,K (t)), for i ∈ N − {N }. From (4), with (22), we see that (z 1 (t), . . . , z N −1 (t)) satisfies the following difference-differential system associated with system (4): where z i,k (t) = x i,k (t) − x i+1,k (t), and We shall establish the ε-approximate synchronization of system (4) by showing for some ε ≥ 0, under some conditions. The error ε off the synchronous set S naturally depends on the difference of units, represented by δF i,k , and the structure of network which is imbedded inG i,k . For the factors from the network structure, the row sums κ i of the connection matrix, defined in (7), and coupling delays τ ij , defined in (3), play important roles. Sequential contracting is an iterative argument which can be used to conclude the asymptotic behavior of z i,k (t) satisfying (23). To analyze the difference-differential equation (23), we shall recompose the terms in H i,k so that (23) can be rewritten aṡ for (i, k) ∈ A and t ≥ t 0 , where The purpose of recasting system (23) into system (26) is to exploit and utilize the dissipative structure of the coupled systems and the coupling structure contained iñ G i,k , so that a scheme can be formulated for the attraction of z i,k (t) satisfying (23) toward a neighborhood of zero, as t gets large. The dissipative property of coupled systems is certainly associated with the characters imbedded in subsystems depicted by F i,k and the coupling terms. In particular, we collect the terms involving which contributes to a contraction of z i,k (t) towards zero, if µ k > 0. In addition, e i,k contains δF i,k which depicts the difference between the ith and the (i + 1)th subsystems, σ i,k which is associated with the difference between row sums of connection matrix W , and o i,k which is associated with the inhomogeneity of coupling delays, cf. Remark 2. Moreover,w i,k collects all the remaining terms. Now, let us present how h i,k ,h i,k ,w i,k , σ i,k , and o i,k are derived, and the process to arrive (26) from (23). For the first two terms in (24), we consider the following decomposition for Ξ = (ξ 1 , . . . , ξ K ),Ξ = (ξ 1 , . . . ,ξ K ) ∈ R K , t ≥ t 0 , and (i, k) ∈ A. Such a decomposition aims to collect in f i,k certain terms involving ξ k ,ξ k , and leave the other terms in u i,k . Conditions on f i,k and u i,k will be imposed in assumption (F) below. We will illustrate such decomposition in Sections 3 and 4. ForG i,k in (25), we perform the following manipulation. We compose a matrix whose entries consist of coupling weights ω ij and row sums κ i of the connection matrix:W FromW , we further define matrixW : where T means transpose, and C is an (N − 1) × N matrix defined by ThenW is well defined and satisfies cf. [35]. In addition, it is shown in Appendix A that With (30), we shall decomposeG i,k in (25) into four parts in the following lemma.
introduced in the beginning of this section.
and κ i ,ω ij , andω ij are defined in (7), (28), and (29), respectively. Moreover, for all i ∈ N − {N } and Ξ = (ξ 1 , . . . , ξ N ) T ∈ R N . From (25) and the definition of ω ij in (28), we obtaiñ From this expression, we further decompose where the equality (33) follows from (32). We thus obtain the decomposition of G i,k in the assertion. Now, let us go back to consider H i,k defined in (24). Recall δF i,k defined in (22), f i,k and u i,k defined in (27), andg i,k , v i,k , σ i,k , and o i,k defined in Lemma 2.2. By (22), (27), and Lemma 2.2, system (23) can be rewritten as system (26): and We will see how these terms are connected to the intrinsic dynamics of the subsystems and the coupling terms in the applications in Sections 3 and 4.

Remark 2. (i)
The terms e i,k := δF i,k + σ i,k + o i,k defined in (39) collect the inhomogeneities in the coupled system (4). Indeed, from the definition in (22), δF i,k depicts the difference between two adjacent subsystems. From Lemma 2.2, σ i,k is associated with the difference between adjacent row sums of connection matrix, and o i,k collects the terms related to the heterogeneity of coupling delays. Notice that e i,k vanishes if the coupled system (4) is homogeneous, and has identical row sums and identical coupling delays in the connection (i.e., (5) and (6) hold).
(ii) Let us summarize the whole decomposition. We obtained in (24) a preliminary decomposition of H i,k into three parts: In this step, we extracted from H i,k the term δF i,k which captures the distinction between (adjacent) subsystems, as mentioned in (i).
, we further decomposed it into two parts: f i,k and u i,k , in (27). By composing matricesW andW in (28)-(29) from connection matrix W , and then applying (30), we further decomposedG i,k into four parts: In this process, the heterogeneities caused from the connection matrix and coupling delays were extracted fromG i,k into σ i,k and o i,k , respectively.
The idea of recomposing H i,k in (23) and recasting (23) into (34) originates from extending the formulations in [45,54], where identical synchronization for a coupled system comprising identical subsystems was studied. In particular, Proposition 2.1 in [54] established a decomposition of H i,k , which is in the form of the right hand side of (34) with e i,k (·) ≡ 0. By applying sequential contracting technique, this led to the convergence of z i,k (t) to zero as t → ∞, for z i,k (t) satisfying the difference-differential equations, and hence yielded the identical synchronization. The transformation (30) has been applied to diffusively coupled systems in [35]. Its application in [54] was to treat a system similar to (4), but with identical row sums for the connection matrix.
In the present investigation, taking into account distinctions between subsystems, and disparities between row sums and between coupling delays in the connection, the key idea is to develop an ingenious and effective decomposition on the associated difference-differential equations, to suitably extract and collect these inhomogeneous factors, and handle the remaining terms by the sequential contracting arguments. The extension of (30) to treat system (4) with nonidentical row sums becomes possible only after such new decomposition is derived. In the rest of this section, we shall modify and extend the sequential contracting technique to the decomposed difference-differential equations (34) to obtain the convergence of z i,k (t) to a neighborhood of zero as t → ∞, which then yields the approximate synchronization of system (4).
Next, let us introduce two assumptions. Assumption (E): For each i, j ∈ N and k ∈ K, there exists anδ k ≥ 0 such that for l ∈ K − {k}, such that the following properties hold: Functions F i,k govern the ith subsystem in (4) when there is no coupling. Hence, δ k in assumption (E) measures the difference between subsystems when they are isolated. Notably, each δF i,k defined in (22) satisfies under assumption (E). Assumption (F) imposes conditions on functions f i,k and u i,k defined in (27), where f i,k and u i,k are determined by functions F i,k . We shall illustrate how to justify these assumptions as we consider coupled FitzHugh-Nagumo equations (69) in Proposition 3.3, and coupled Lorenz equations (99) in Proposition 4.3.
Recall the definitions of σ i,k and o i,k in Lemma 2.2. Applying (43) and (44) then leads to  and (F). Assertion (iv) relies on the additional assumption (E) which measures the variation of subsystems, and on the new setting of C Q , i.e., |φ ′ k (·)| ≤ Λ k in (19). Moreover, as seen from Remark 1 and (45), the part j∈N |ω (i+1)j |)D k Λ k δ τ ofǭ i,k in assertion (iv) measures the heterogeneity caused from coupling delays between components. This explains why a new setting of C Q is needed, as compared to the one in [45,54]. (ii) Recall Remark 2 that e i,k := δF i,k + σ i,k + o i,k defined in (39) collects the inhomogeneous factors in the network system (4). As seen from (5), (6), and assumption (E),δ k , δ κ , and δ τ measure the differences between subsystems, the differences between row sums, and the variation of coupling delays, respectively. In addition, ǫ i,k defined in Proposition 2.3 decreases to zero asδ k , δ κ , and δ τ decrease to zero. This reveals that eachǭ i,k converges to zero as the coupled system (4) tends to being homogeneous, and the row sums and coupling delays in the connection tend to being identical. It will be seen in Theorem 2.6 and Remark 4 how the termsǭ i,k determine an estimate on the synchronization error. (iii) The termsβ i,k andβ i,k in Proposition 2.3 depend on the sign of κ i +ω ii . The sign of κ i is easily seen from the row sum of connection matrix, defined in (7). On the other hand, determining the sign ofω ii relies on the expression in (31) which is newly derived in Appendix A of this paper. It can be observed from (7) and (31) that, allω ii are apt to be negative if all off-diagonal entries of W are nonnegative, or the magnitudes of the negative off-diagonal entries are relatively smaller than those of the positive off-diagonal entries. Suppose that allω ii are negative. If all row sums of W are small enough, then κ i +ω ii < 0 for every i. For a nonzero connection matrix W with nonnegative off-diagonal entries and satisfying the diffusive condition (κ i = 0 for all i ∈ N ), which is a common assumption in the literature, one does have κ i +ω ii =ω ii < 0 for every i.
For later use, with the quantitiesμ i,k ,μ i,k ,β i,k ,β i,k ,μ Herein, the term η i,k is delay-dependent, since it involves the magnitude of delay τ ii .
Recall that for an arbitrary solution (x 1 (t), . . . , x N (t)) of system (4), where (23) which can be recast into system (26). In the following, let us preview the main process for showing lim sup t→∞ z i (t) ≤ ε e , for all i = 1, . . . , N −1, and some nonnegative ε e ≥ 0. This then establishes the approximate synchronization of system (4). Let us first relabel the two-dimensional indices in system (26) into one-dimensional indices, through the bijective mapping ℓ: The labeling indicates an ordering of ℓ(i, k) by the orders of i and k in succession: For later use, we define, for (i, k) ∈ A, By setting z ℓ(i,k) := z i,k = x i,k − x i+1,k , h ℓ(i,k) := h i,k ,h ℓ(i,k) :=h i,k , and w ℓ(i,k) := w i,k , we can rewrite system (26) aṡ for (i, k) ∈ A and t ≥ t 0 . For each (i, k), the equation (53) is analogous to a scalar equation: where z(t) = x(t) − y(t). In Appendix B, we introduce precisely this scalar equation (54) (relabeled as (118) therein), and summarize the convergent property of (118) in Proposition 6.1 which reveals that z(t) converges to an interval [−p,p] as t → ∞; in addition,p can be estimated. Applying such a convergent property to each z ℓ(i,k) (t), (i, k) ∈ A, we obtain the convergence of z ℓ(i,k) (t) to an interval [−p ℓ(i,k) ,p ℓ(i,k) ], with where η i,k is defined in (46). Herein, where w ℓ(i,k) is given in (53). We can derive better estimations onp ℓ(i,k) for all (i, k) ∈ A. Notice that for each (i, k), if estimates ofp ℓ(j,l) have been obtained for all (j, l) ∈Ľ i,k , we may use them to bound z ℓ(j,l) (t) which plays the role as φ j,l (θ) − φ j+1,l (θ) in (41), and hence obtain sharper estimate on |w ℓ(i,k) | max (∞), and subsequently sharper estimate onp ℓ(i,k) , by (55). We perform such a process to z 1 (t), . . . , z (N −1)K (t) successively and iteratively and obtain, for each (i, k) ∈ A, for n ≥ 1. We observe that, with p is exactly the Gauss-Seidel iteration for solving the linear system where with η i,k defined in (46),m (j,l) i,k in (47), andǭ i,k in Proposition 2.3. Notice that M is delay-dependent, since the diagonal entries η i,k involve the magnitudes of delay τ ii , as indicated in (46).
Next, for further estimations onp ℓ(i,k) , let us show how the sequential estimates {p n=0 which satisfy the iterative formula (56) is derived. This can be seen by induction in the following steps (i) and (ii).
With p is exactly the Gauss-Seidel iteration for solving the linear system (57). Recall that, for each (i, k) ∈ A, z ℓ(i,k) (t) converges to [−p ℓ(i,k) ,p ℓ(i,k) ] as t → ∞, and 0 ≤p ℓ(i,k) ≤ There is an alternative way to estimate the size of intervals [−p ℓ(i,k) ,p ℓ(i,k) ]. By employing a different sequence of upper-lower bounds in analyzing equation (54) (i.e., (118)), we can obtain a disparate estimate for the interval [−p ℓ(i,k) ,p ℓ(i,k) ] to which z ℓ(i,k) (t) converges, as stated in Proposition 6.2. Subsequently, an alternative iterative formula can be obtained, which represents the Gauss-Seidel iteration for another linear system: where e is defined in (59), η i,k ,m We introduce a delay-independent condition: By arguments similar to those for Proposition 2.4, but using Proposition 6.2 instead of Proposition 6.1, we can derive the following delay-independent result.
where DM, −LM, and −UM are the diagonal, strictly lower-triangular, and strictly upper-triangular parts ofM, respectively.
With the preparation in Lemmas 2.1, 2.2 and Propositions 2.3-2.5, the problem of approximate synchronization is converted to solving a corresponding linear algebraic system (57) or (64), which is nonhomogeneous if e = 0, and homogeneous if e = 0. If the coupled system (4) has identical subsystems, identical row sums, and identical coupling delays in the connection, expressed in (5) and (6), thenǭ i,k = 0 for all (i, k) ∈ A, and thus e = 0 and ε e =ε e = 0. In this case, the formula in (56) and the assertion in Theorem 2.6 reduce to the ones in [45], where the problem of identical synchronization was converted to solving a corresponding homogeneous linear system of algebraic equations.
This proposition follows from some energy-like estimates and iterative arguments. We sketch the main justification in Appendix C.
Notice thatq (n) are strictly decreasing with respect to n ∈ N. Thus, Q (n) with larger n provides a smaller attracting region for system (69), and thereby relaxes the condition in the synchronization criterion and yields finer estimation on the synchronization error ε e . It will be illustrated in the following Example 1 that ε e becomes smaller if considering Q * = Q (5) rather than Q * = Q (1) .
The following proposition confirms assumption (F) for system (69). Its justification is similar to Proposition 3.2 of [54] and is omitted.
Remark 5. (i) As seen from Theorem 3.4, the entries m ı in matrix M are determined by both of intrinsic dynamics of the individual FitzHugh-Nagumo neurons and the coupling terms. In particular, the terms κ i ,ω ij , andω ij are associated with the connection matrix W , as seen in (7), (28), and (29).
(ii) We can extend Remark 3.1 in [54] to make the following observations. In general, the criterion in Theorem 3.4 prefers negative (a 2 i − a i + 1)/3 + c(κ i +ω ii )L g of large magnitude and small delay τ ii , so that each −(a 2 i − a i + 1)/3 − c(κ i +ω ii )L g − τ iiβi is positive and large. In particular, if the delay in system (69) is nonzero, it requires that the coupling strength c be suitably large, not too small nor too large. When there is no delay, condition (87) reduces to the first inequality therein. If W is the zero matrix, i.e., there is no connection between neurons in system (69), then κ i = 0 andω ii = 0 for every i, and condition (87) can not hold. (iii) The first inequality in (87) implies κ i +ω ii < 0 for all i. Thus, κ i +ω ii < 0 for all i is required in the criterion of Theorem 3.4. Accordingly, as seen in Remark 3(ii), the conditions in Theorem 3.4 are more compatible with the situation that all off-diagonal entries of W are nonnegative, or the magnitudes of the negative offdiagonal entries are relatively smaller than those of the positive off-diagonal entries, and all row-sums κ i of W are small enough.
Let us demonstrate Theorem 3.4 by the following example which is a system comprising three FitzHugh-Nagumo neurons. This example aims at demonstrating the dependence of synchronization error on various heterogeneous factors of the considered system, i.e., from disparities of the subsystems in terms of parameter mismatches, of row sums of the connection matrix, and of the coupling delays. We are not aware of any such example or similar in the literature, which addresses the approximate synchronization of coupled systems with heterogeneity from connection matrix or coupling delays.
Recall that, in Proposition 3.1,q (n) , n ∈ N, are strictly decreasing with respect to n. Thus, for larger n, Q (n) provides a smaller attracting region for the dynamics of system (69), which actually yields finer estimation on the synchronization error. For example, if we take Q * = Q (5) so that assumption (D) is satisfied, then λ syn and ε e can be computed as 0393798903. The real synchronization errors can be seen in Figures 1(a)-1(d), which appear to be much smaller than the theoretical estimates. It is also interesting to see from these computations that the inhomogeneity from the row sums of connection matrix contributes the largest synchronization error. We note that these errors decrease as the magnitudes of the inhomogeneities diminish, according to our theory. Remark 6. (i) As shown in Example 1, theoretical estimation on the synchronization error depends on the estimate of the global attracting set in the phase space, and thus may be coarse. The main goal of this study is to establish approximate synchronization criterion for coupled oscillators, such as the one in Theorem 3.4. The error estimate comes with the theory to confirm that the error decreases as the inhomogeneities of the coupled system decrease. Once the criterion is met for certain parameter and delay values, the actual deviation of the asymptotical state from the synchronization set S can be seen from numerical simulations. (ii) As mentioned in the Introduction, the existing works on approximate synchronization commonly considered linear coupling without delays [18,36,59,66]. The approximate synchronization of two nonlinearly coupled systems without delay and linearly coupled systems without delay were reported in [41] and [42], respectively. In addition, it is commonly required in those works that the connection matrix satisfies diffusion condition [18,36,41,59,66], and its off-diagonal entries have the same sign [18,36,59]. Notably, the off-diagonal entries of connection matrix in our Example 1 have mixed signs and the connection matrix does not satisfy diffusion condition. To the best of our knowledge, synchronization of the system in this example can not be concluded by the methods employed in those papers, even for the case δ * κ = δ * τ = 0.
Example 2. Consider system (99) with ϑ = 1. By Proposition 4.1, (109), and Theorem 4.4, it can be computed thatR 1 =R 2 ≈ 119.4494,R 3 ≈ 157.4494, and ρ ≈ 129.4494, and functions λ ± (c) are then determined. Subsequently, we can compute to find c ϑ ≈ 63.3819. Herein, we stress that c ϑ is computed through the analytical result in Theorem 4.4, not through a numerical computation. Thus system (99) attains approximate synchronization when c > c ϑ and asymptotic synchronization, by Theorem 4.4. Figure 2 demonstrates the evolution of function Err(t) := max 1≤k≤3 {|x 1,k (t) − x 2,k (t)|}, for solutions (x i,j (t)) of (99) with ϑ = 1, and c = 64, c = 640, and c = 6400, respectively. Herein, Err(t) is referred to as the synchronization error for the corresponding solution. It appears in Figure 2 that the synchronization error decreases as the coupling strength increases, and thus strong coupling strength suppresses the synchronization error. We note that the system appears to exhibit chaotic behavior in all cases in Figure 2.
As seen from the approach in [17], an important and necessary concern for the study of asymptotic synchronization in coupled systems with parameter mismatch is finding a uniform attracting region which is independent of the coupling strength. Finding an estimate of the attracting region for coupled Lorenz equations with parameter variation is actually a challenging task, especially for those coupled through all components, such as system (99). The asymptotic synchronization of two nearly identical Lorenz equations coupled through their first components was discussed in [17], and the work concerning the existence of a uniform attracting region for Lorenz equations with parameter variation coupled through their first components can be found in [41,42]. In this section, we studied the approximate synchronization of system (99) which is coupled in three components and system (114) which is coupled in two components, and discussed the possibility for them to achieve asymptotic synchronization. We observe that a main difference between systems (99) and (114), which results in whether the asymptotic synchronization holds or not, is about the coupling configuration. The coupling is full in system (99), but there is no coupling between variables x 1,3 and x 2,3 in system (114). The consequence is that the diagonal entry in the third row of M * in (108) is β + 2c, whereas the diagonal entry in the third row ofM in (115) is β, where c is the coupling strength. It can be observed thatε ϑ,c in (116) converges to zero if this entry β is replaced by β + 2c.  5. Discussion and conclusion. We have established a framework to study approximate synchronization of coupled systems. The framework accommodates a large class of coupled systems and network systems under general coupling. The units which constitute the coupled systems can be identical or nonidentical. The nonidentical ones include possible cases of slightly different units, units with parameter mismatch, or completely disparate units. Communication delay was also taken into account. Under this setting, there naturally exists a synchronization error when the synchronization criterion is met. We have shown that the synchronization error depends on and decreases with decreasing difference between subsystems, difference between row sums of connection matrix, and difference of communication delays between components, as summarized in Remark 4. When these differences decrease to zero, approximate synchronization reduces to identical synchronization for coupled systems comprising identical subsystems, with identical row sums of connection matrix and identical communication delays. Among these inhomogeneous factors, the one from row sums of connection matrix is less obvious. In fact, identical row sums of connection matrix carries a sense of balance for total coupling weights for each unit from all units, and thus plays an important role in the coupling structure.
The present theory provides both delay-dependent and delay-independent criteria for approximate synchronization. The sufficient conditions depend on estimations of the globally attracting region for a coupled system. Finer estimates on the attracting region lead to sharper estimates on the synchronization error and weaker criterion for synchronization. In addition, in order to study asymptotic synchronization, one would expect to obtain an estimate on the attracting set, which is independent of the coupling strength, such as the coupled Lorenz systems in Example 2. However, rigorous confirmation of dissipativeness and estimation of attracting set in nonlinear systems are nontrivial mathematical tasks in general. Many of the existing results on synchronization and/or asymptotic behaviors of dynamical systems, via various mathematical treatments, also rely on such a property. As analytic estimations on globally attracting sets are difficult to obtain in general, examinations of the synchronization criteria largely rely on numerical computations. Those numerical computations are generally executed after assigning parameter values and choosing some (finitely many) initial points, and hence are not mathematically sufficient. In this paper, we were able to confirm the dissipative property, and obtained estimations of the globally attracting sets for coupled systems comprising nonidentical FitzHugh-Nagumo neurons under sigmoidal (nonlinear) coupling and coupled systems comprising parameter-mismatched Lorenz equations under linearly diffusive coupling.
The approach in this paper can be applied to establish approximate synchronization for coupled systems under quite general coupling, such as (2), and even the ones not in the form of (2) or (4). For instance, we can also establish approximate synchronization for the kinetic models on gene regulatory networks of zebrafish segmentation clocks, which are under very nonlinear coupling with transcriptional and translational delays, cf. [10,31]. The key point for such an application relies on suitable manipulation of the associated difference-differential equations. This has to be formulated according to the coupling configuration. More precisely, it hinges upon recomposing the difference-differential equations, such as (23) in this paper, into a form, such as (34), for which assertions analogous to the ones in Proposition 2.3 can be formulated. However, such formulation is quite cumbersome.
As mentioned in the Introduction, in previous studies of nonidentical synchronization, the coupled systems considered were under linear diffusive coupling, and without delays or with identical delays. The factors contributed to synchronization error were only partially uncovered. In the present investigation, via a new approach "sequential contracting", we were able to establish a framework to tackle the problem on nonidentical synchronization for coupled systems in general form. In a systematic way, we have exploited the structure behind the synchronization, including the causes of synchronization error. With this approach, some common restrictions in the literature have been lifted, as mentioned in Remark 6. Furthermore, for practical applications, all conditions in the present theory can be examined in a straightforward fashion or by numerical algebraic computation. Significantly, with this framework, the general theory for synchronization has attained a completeness and the scope on the study of synchronization is expected to be further expanded. For example, one can see how various factors, including the intrinsic terms of the individual subsystems, the coupling structure, coupling strength, and coupling delays play their roles in synchronization, such as the ones shown for coupled FitzHugh-Nagumo neurons in Theorem 3.4, Remark 5, and Example 1. Based on the present framework, we also extended our analysis to investigate asymptotic synchronization for coupled Lorenz equations and illustrated that for some coupling schemes, the synchronization error decreases as the coupling strength increases, whereas for some other coupling schemes, synchronization error can not be suppressed by increasing the coupling strength, as shown in Examples 2 and 3. This has clarified an issue which was addressed, but not settled, in the literature. From the view point of application, the factors that we found can be used to control the synchronization error in network systems. All in all, the exploration in this work has enhanced and expanded the understanding of synchronization and the dynamics of coupled systems. 6. Appendices. We outline the derivation for the entries of the transformed connection matrix in Appendix A, and summarize the properties for the scalar equation (54) in Appendix B. We estimate the globally attracting sets for the coupled FitzHugh-Nagumo system and coupled Lorenz equations comprising subsystems with mismatched parameters, in Appendices C and D, respectively.