Achieving synchronization in arrays of coupled differential systems with time-varying couplings

In this paper, we study complete synchronization of the complex dynamical networks described by linearly coupled ordinary differential equation systems (LCODEs). The coupling considered here is time-varying in both the network structure and the reaction dynamics. Inspired by our previous paper [6], the extended Hajnal diameter is introduced and used to measure the synchronization in a general differential system. Then we find that the Hajnal diameter of the linear system induced by the time-varying coupling matrix and the largest Lyapunov exponent of the synchronized system play the key roles in synchronization analysis of LCODEs with the identity inner coupling matrix. As an application, we obtain a general sufficient condition guaranteeing directed time-varying graph to reach consensus. Example with numerical simulation is provided to show the effectiveness the theoretical results.


Introduction
Complex networks have widely been used in theoretical analysis of complex systems, such as Internet, World Wide Web, communication networks, social networks. A complex dynamical network is a large set of interconnected nodes, where each node possesses a (nonlinear) dynamical system and the interaction between nodes is described as diffusion. Among them, linearly coupled ordinary differential equation systems (LCODEs) are a large class of dynamical systems with continuous time and state.
The LCODEs are usually formulated as followṡ where t ∈ R + = [0, +∞) stands for the continuous time and x i (t) ∈ R n denotes the variable state vector of the i-th node, f : R n → R n represents the node dynamic of the uncoupled system, σ ∈ R ++ = (0, +∞) denotes coupling strength, l ij ≥ 0 with i = j denotes the interaction between the two nodes and l ii = − m j =i l ij , B ∈ R n,n denotes the inner coupling matrix. The LCODEs model is widely used to describe the model in nature and engineering. For example, the authors study spike-burst neural activity and the transitions to a synchronized state using a model of linearly coupled bursting neurons in [1]; the dynamics of linearly coupled Chua circuits are studied with application to image processing, and many other cases in [2].
For decades, a large number of papers have focused on the dynamical behaviors of coupled systems [3]- [5], especially the synchronizing characteristics. The word "synchronization" comes from Greek, in this paper the concept of local complete synchronization (synchronization for simplicity) is considered (see Definition 1). For more details, we refer the readers to [6].
Synchronization of coupled systems have attracted a great deal of attention [7]- [9]. For instances, in [7], the authors considered the synchronization of a network of linearly coupled and not necessarily identical oscillators; in [8], the authors studied globally exponential synchronization for linearly coupled neural networks with time-varying delay and impulsive disturbances. Synchronization of networks with time-varying topologies was studied in [10]- [14]. For example, in [10], the authors proposed the global stability of total synchronization in networks with different topologies; in [14], the authors gave a result that the network will synchronize with the time-varying topology if the time-average is achieved sufficiently fast.
Synchronization of LCODEs has also been addressed in [15]- [17]. In [15], mathematical analysis was presented on the synchronization phenomena of LCODEs with a single coupling delay; in [16], based on geometrical analysis of the synchronization manifold, the authors proposed a novel approach to investigate the stability of the synchronization manifold of coupled oscillators; in [17], the authors proposed new conditions on synchronization of networks of linearly coupled dynamical systems with non-Lipschitz right-handsides. The great majority of research activities mention above all focused on static networks whose connectivity and coupling strengths are static. In many applications, the interaction between individuals may change dynamically. For example, communication links between agents may be unreliable due to disturbances and/or subject to communication range limitations.
In this paper, we consider synchronization of LCODEs with time-varying coupling. Similar to [15]- [17], time-varying coupling will be used to represent the interaction between individuals. In [6,13], they showed that the Lyapunov exponents of the synchronized system and the Hajanal diameter of the variational equation play key roles in the analysis of the synchronization in the discrete-time dynamical networks. In this paper, we extend these results to the continuous-time dynamical network systems. Different from [11,14], where synchronization of fast-switching systems was discussed, we focus on the framework of synchronization analysis with general temporal variation of network topologies. Additional contributions of this paper are that we explicitly show that (a) the largest projection Lyapunov exponent of a system is equal to the logarithm of the Hajanal diameter, and (b) the largest Lyapunov exponent of the transverse space is equal to the largest projection Lyapunov exponent under some proper conditions. The paper is organized as follows: in Section 2, some necessary definitions, lemmas, and hypotheses are given; in Section 3, synchronization of generalized coupled differential systems is discussed; in Section 4, criteria for the synchronization of LCODEs are obtained; in Section 5, we obtain a sufficient condition ensuring directed time-varying graph reaching consensus; in Section 6, example with numerical simulation is provided to show the effectiveness of the theoretical results; the paper is concluded in Section 7.
Notions: e n k = [0, 0, , · · · , 0, 1, 0, · · · , 0] ⊤ ∈ R n denotes the n-dimensional vector with all components zero except the k-th component 1, 1 n denotes the n-dimensional column vector with each component 1; For a set in some Euclidean space U ,Ū denotes the closure of U , U c denotes the complementary set of U , and A \ B = A ∩ B c ; For u = [u 1 , · · · , u n ] ⊤ ∈ R n , u denotes some vector norm, and for any matrix A = (a ij ) ∈ R n,m , A denotes some matrix norm induced by vector norm, for example, u 1 = n i=1 |u i | and A 1 = max j n i=1 |a ij |; For a matrix A = (a ij ) ∈ R n,m , |A| denotes a matrix with |A| = (|a ij |); For a real matrix A, A ⊤ denotes its transpose and for a complex matrix B, B * denotes its conjugate transpose; For a set in some Euclidean space W , O(W, δ) = {x : dist(x, W ) < δ}, where dist(x, W ) = inf y∈W x − y ; #J denotes the cardinality of set J; ⌊z⌋ denotes the floor function, i.e., the largest integer not more than the real number z; ⊗ denotes the Kronecker product.

Preliminaries
In this section we will give some necessary definitions, lemmas, and hypotheses. Consider the following general coupled differential system: with initial state where t 0 ∈ R + denotes the initial time, t ∈ R + denotes the continuous time and x i (t) = [x i 1 (t), · · · , x i n (t)] ∈ R n denotes the variable state of the i-th node, i = 1, 2, · · · , m.
For the functions f i : R nm × R + → R n , i = 1, 2, · · · , m, we make following assumption: Assumption A1 (a) There exists a function f : R n → R n such that f i (s, s, · · · , s, t) = f (s) for all i = 1, 2, · · · , m, s ∈ R n , and t ≥ 0 is uniformly locally Lipschitz continuous: there exists a locally bounded function K(x, y) such that We say a function g(y) : R q → R p is locally bounded if for any compact set K ⊂ R q , there exists M > 0 such that g(y) ≤ M holds for all y ∈ K.
The first item of assumption A1 ensures that the diagonal synchronization manifold is an invariant manifold for (2). If x 1 (t) = x 2 (t) = · · · = x m (t) = s(t) ∈ R n is the synchronized state, then the synchronized state s(t) satisfiesṡ (t) = f (s(t)). ( Since f (·) is C 1 -smooth, then s(t) can be denoted by the corresponding continuous semi-flow s(t) = ϑ (t) s 0 of the intrinsic system (3). For ϑ (t) , we make following assumption: Assumption A2 The system (3) has an asymptotically stable attractor: there exists a compact set A ⊂ R n such that (a) A is invariant through the system (3), i.e., ϑ (t) A ⊂ A for all t ≥ 0; (b) There exists an open bounded neighborhood U of A such that t≥0 ϑ (t)Ū = A; (c) A is topologically transitive, i.e., there exists s 0 ∈ A such that ω(s 0 ), the ω limit set of the trajectory ϑ (t) s 0 , is equal to A [3].
Definition 1 Let A m denote the Cartesian product A × · · · × A (m times). Local complete synchronization (synchronization for simplicity) is defined in the sense that the set is an asymptotically stable attractor in R nm . That is, for the coupled dynamical system (2), differences between components converge to zero if the initial states are picked sufficiently near S A m , i.e., if the components are all close to the attractor A and if their differences are sufficiently small.
Next we give some lemmas which will be used later, and the proofs can be seen in Appendices.
Lemma 2 Under the assumptions A1 and A2, there exists a compact neighborhood W of A such that We have the following variational equation near the synchronized state s(t).
Or in matrix form: where DF t (s(t)) denotes the Jacobin matrix DF t (ŝ(t)) for simplicity. From [18], we can give the results on the existence, uniqueness, and continuous dependence of the equations (2) and (4).
Lemma 3 Under the assumption A2, each of the differential equations (2) and (4) has a unique solution which is continuously dependent on the initial condition.
Thus, the solution of the linear system (4) can be written in matrix form.
where u k (t, t 0 , s 0 ) denotes the k-th column and is the solution of the following Cauchy problem: Immediately, according to Lemma 3, we can conclude that the solution of the following Cauchy problem can be written as δx(t) = U (t, t 0 , s 0 )δx 0 . We define the time varying Jacobin matrix DF t by the following way: with s(t 0 ) = s 0 , where 2 R nm,nm is the collection of all the subsets of R nm,nm .
Definition 3 For a time varying system denoted by DF, we can define its Hajnal diameter of the variational system (4) as follows: where for a R nm,nm matrix in block matrix form: U = (U ij ) m i,j=1 with U ij ∈ R n,n , its Hajnal diameter is defined as follows: Lemma 4 (Grounwell-Beesack's inequality) If function v(t) satisfies the following condition: where b(t) ≥ 0 and a(t) are some measurable functions, then we have Based on the assumption A1, for the solution matrix U , we have the following lemma: Lemma 5 Under the assumption A1, we have 1. m j=1 U ij (t, t 0 , s 0 ) =Ȗ (t, t 0 , s 0 ), whereȖ (t, t 0 , s 0 ) denotes the solution matrix of the following Cauchy problem: 2. For any given t ≥ 0 and the compact set W given in lemma 2, U (t + t 0 , t 0 , s 0 ) is bounded for all t 0 ≥ 0 and s 0 ∈ W and equi-continuous with respect to s 0 ∈ W .
Let P = (P ij ) m i,j=1 be a R nm,nm matrix with P ij ∈ R n,n satisfying: (a) P i1 = 1 √ m P 0 for some orthogonal matrix P 0 ∈ R n,n and all i = 1, 2, · · · , m; (b) P is also an orthogonal matrix in R nm,nm . We also write P and its inverse P −1 = P ⊤ in the form where P 1 = 1 √ m 1 m ⊗ P 0 and P 2 ∈ R nm,n(m−1) . According to Lemma 5, we have Since P ⊤ 2 P 1 = 0 which implies that each row of P ⊤ 2 is located in the subspace orthogonal to the subspace {1 m ⊗ ξ, ξ ∈ R n }, we can conclude that P ⊤ 2 U (t, t 0 , s 0 )P 1 = 0. Then, we have , α(t, t 0 , s 0 ) ∈ R n,n(m−1) denotes a appropriate matrix we omit its accurate expression. One can see thatŨ (t, t 0 , s 0 ) is the solution matrix of the following linear differential system.

Definition 4
We define the following linear differential system by the projection variational system of (4) along the directions P 2 : where D P F t (s(t)) = P ⊤ 2 DF t (s(t))P 2 .
Definition 5 For any time varying variational system DF : R + × R n → 2 R nm,nm , we define the Lyapunov exponent of the variational system (4) as follows: where u ∈ R nm and s(t 0 ) = s 0 .
Similarly, we can define the projection Lyapunov exponents by the following projection timevarying variation: Then, we have the following lemma: Remark 1 From Lemma 6, we can see that the largest projection Lyapunov exponent is independent of the choice of matrix P .
Consider the time-varying driven by some metric dynamical system (MDS) (Ω, B, P, ̺ (t) ), where Ω is the compact state space, B is the σ-algebra, P is the probability measure, and ̺ (t) is a continuous semi-flow. Then, the variational equation (4) is independent of the initial time t 0 and can be rewritten as follows: In this case, we denote the solution matrix, the projection solution matrix, and the solution matrix on the synchronization space by U (t, s 0 , ω 0 ),Ũ (t, s 0 , ω 0 ), andȖ (t, s 0 , ω 0 ), respectively. For simplicity, we write them as U (t),Ũ (t), andȖ (t), respectively. Also, we write the Lyapunov exponents and the projection Lyapunov exponent as follows: We add the following assumption: The following are involving linea differential systems. For more details, we refer the readers to [19]. For a continuous scalar function u(t), we denote its Lyapunov exponent by The following properties will be used later: , where c k , k = 1, 2, · · · , n, are constants; For the following linear differential system: where x(t) ∈ R n , a transformation x(t) = L(t)y(t) is said to be a Lyapunov transformation if L(t) satisfies : It can be seen that the class of Lyapunov transformations form a group and the linear system for y(t) should beẏ where Then, we say system (10) is a reducible system of system (9). We define the adjoint system of (9) bẏ If letting V (t) be the fundamental matrix of (9), then [V −1 (t)] * is the fundamental matrix of (11). Thus, we say the system (9) is a regular system if the adjoint systems (9) and (11) have convergent Lyapunov exponent series: {α 1 , · · · , α n } and {β 1 , · · · , β n }, respectively, which satisfy α i + β i = 0 for i = 1, 2, · · · , n; or its reducible system (10) is also regular.

General synchronization analysis
In this section we provide a methodology based on the previous theoretical analysis to judge whether a general differential system can be synchronized or not.
Theorem 1 Suppose that W ∈ R n is the compact subset given in Lemma 2, and assumptions A1 and A2 are satisfied. If then the coupled system (2) is synchronized.
Proof: The main techniques of the proof come from [3,6] with some modifications. Let ϑ (t) be the semi-flow of the uncoupled system (3). By the condition (12), there exist d sat- According to the equi-continuity of U (t 0 + t(s 0 ), t 0 , s 0 ), there exists δ > 0 such that for any s According to the compactness of W , there exists a finite positive number set T = {t 1 , t 2 , · · · , t v } with t j ≥ T 1 for all j = 1, 2, · · · , v such that for any s 0 ∈ W , there exists t j ∈ T such that diam(U (t 0 + t j , t 0 , s 0 )) < 1 3 for all t 0 ≥ 0. Let x(t) be the collective states {x 1 (t), · · · , x m (t)} which is the solution of the coupled system (2) with initial condition x i (t 0 ) = x i 0 , i = 1, 2, · · · , m. And, let s(t) be the solution of the synchronization state equation where ξ ij kl (t) ∈ R mn , i, j = 1, 2, · · · , m, k, l = 1, 2, · · · , n, are obtained by the mean value principle of the differential functions. Letting DF t (ξ(t)) = ( , we can write the equations above in matrix form: and denote its solution matrix byÛ (t + t 0 , t 0 , Then, for any t > 0 there exists K 2 > 0 such that DF t+t 0 (ξ(t + t 0 )) ≤ K 2 for all t ∈ T and t 0 ≥ 0 according to the 3-th item of the assumption A1. Then, we have Picking α sufficiently small such that for each x 0 ∈ W α , there exists t ∈ T such that m j=1 ∆x j (t+ t 0 ) < r 0 2 and diam(Û (t + t 0 , t 0 , x 0 )) < 1 2 for all t 0 ≥ 0. Thus, we are to prove synchronization step by step.
For any x 0 ∈ W α , there exists t ′ = t(x 0 ) ∈ T such that Then, re-initiated with time t ′ + t 0 and condition x(t ′ + t 0 ), continuing with the phase above, we can obtain that lim t→∞ max i,j x i (t) − x j (t) = 0. Namely, the coupled system (2) is synchronized. Furthermore, from the proof, we can conclude that the convergence is exponential with rate O(δ t ) where δ = sup s 0 ∈W diam(DF t , s 0 ), and uniform with respect to t 0 ≥ 0 and x 0 ∈ W α . This completes the proof.
Remark 2 According to the assumption A2 that attractor A is asymptotically stable and the properties of the compact neighbor W given in Lemma 2, we can conclude that the quantity is independent on the choice of W .
If the time-variation is driven by some MDS {Ω, B, PP, ̺ (t) } and there exists a metric dynamical system {W × Ω, F, P, π (t) }, where F is the product σ− algebra on W × Ω, P is the probability measure, and π (t) (s 0 , ω) = (θ (t) s 0 , ̺ (t) ω). From Theorem 1, we have Corollary 1 Suppose that the conditions in Lemma7 are satisfied, W × Ω is compact in the topology defined in this MDS, the semi-flow π (t) is continuous, and on W × Ω the Jacobian be the Lyapunov exponents of this MDS with multiplicity and {σ i } n i=1 correspond the synchronization space. If where Erg π (W × Ω) denotes the ergodic probability measure set supported in the MDS {W × Ω, F, P, π (t) }, then the coupled system (2) is synchronized.

Synchronization of LCODEs with identity inner coupling matrix and time-varying couplings
In this section we study synchronization in linearly coupled ordinary differential equation systems (LCODEs) with time-varying couplings. Considering the following LCODEs with identity inner coupling matrix: where x i (t) ∈ R n denotes the state variable of the i-th node, f (·) : R n → R n is a differential map, σ ∈ R + denotes coupling strength, and l ij (t) denotes the coupling coefficient from node j to i at time t, for all i = j, are supposed to satisfy the following assumption. Here, we highlight that the inner coupling matrix is the identity matrix. Assumption A4 (a) l ij (t) ≥ 0, i = j, are measurable and l ii (t) = − m j=1,j =i l ij (t); (b) There exists M 1 > 0 such that |l ij (t)| ≤ M 1 for all i, j = 1, 2, · · · , m.
In the following we omit coupling strength σ since it is just a positive scalar. Similarly, we can define the Hajnal diameter of the following linear system: Let V (t) = (v ij (t)) m i,j=1 be the fundamental solution matrix of the system (14). Then, its solution matrix can be written as V (t, t 0 ) = V (t)V (t 0 ) −1 . Thus, the Hajnal diameter of the system (14) can be defined as follows: If log(diam(L)) + µ < 0, then the LCODEs (13) is synchronized.

By Theorem 1, we have
Proof: Considering the variational equation of (13): LetȖ (t, t 0 , s 0 ) be the solution matrix of the synchronized state system (6) and V (t, t 0 ) = (v ij (t, t 0 )) m i,j=1 be the solution matrix of the linear system (14). We can see that V (t, t 0 ) ⊗ U (t, t 0 , s 0 ) is the solution matrix of the variational system (15). Then, This implies that the Hajnal diameter of the variational system (15) is less than e µ diam(L). This completes the proof according to Theorem 1. For the linear system (14), we firstly have From Lemmas 6 and 7, we have Let I be the set consisting of all compact time intervals in [0, +∞) and G be the the set consisting of all graph with vertex set N = {1, 2, · · · , m}. Define where G(I, δ) = {N, E} is a graph with vertex set N and its edge set E is defined as follows: there exists an edge from vertex j to vertex i if and only if t 2 t 1 l ij (τ )dτ > δ. Namely, we say that there is a δ-edge from vertex j to i across I = [t 1 , t 2 ].

Definition 6
We say that the LCODEs (13) has a δ-spanning tree across the time interval I if the corresponding graph G(I, δ) has a spanning tree.

For a stochastic matrix
Then, we can also define that V is δ-scrambling if η(V ) > δ.
Theorem 3 Suppose the assumption A4 is satisfied. diam(L) < 1 if and only if there exist δ > 0 and T > 0 such that the LCODEs (13) has a δ-spanning tree across any T -length time interval.
Remark 3 Different from [14], we don't need to assume that L(t) has zero column sums and the time-average is achieved sufficiently fast.
Before proving this theorem, we need the following lemma: Lemma 9 If the LCODEs (13) has a δ-spanning tree across any T -length time interval, then there exist δ 1 > 0 and T 1 > 0 such that V (t, t 0 ) is δ 1 -scrambling for any T 1 -length time interval.

Consensus analysis of multi-agent system with directed timevarying graphs
If let n = 1, f ≡ 0 and σ = 1 in system (13), then we havė In this case, if the assumption A4 is satisfied, then the synchronization analysis of system (16) become another important research field named consensus problems.
Definition 7 We say the differential system (16) reaches consensus if for any In graph view, the coefficients matrix of (16) L(t) = (l ij (t)) ∈ R m,m is equal to the negative graph Laplacian associated with the digraph G(t) at time t, where G(t) = (V, E(t), A(t)) is a weighted digraph (or directed graph) with m vertices, the set of nodes V = {v 1 , · · · , v m }, set of edges E(t) ⊆ V × V, and the weighted adjacency matrix A(t) = (a ij (t)) with nonnegative adjacency elements a ij (t). An edge of G(t) is denoted by e ij (t) = (v i , v j ) ∈ E(t) if there is an directed edge from vertex i to vertex j at time t. The adjacency elements associated with the edges of the graph are positive, i.e., e ij (t) ∈ E(t) ⇐⇒ a ij (t) > 0, for all i, j ∈ {1, 2, · · · , m}. It is assumed that a ii (t) = 0 for all i ∈ {1, 2, · · · , m}. The in-degree and out-degree of node v i at time t are, respectively, defined as follows: The degree matrix of digraph G(t) at time t is defined as D(t) = diag(deg out (v 1 (t)), · · · , deg out (v m (t))). The graph Laplacian associated with the digraph G(t) at time t is defined as Let G(I, δ) defined as before. We say that the digraph G(t) has a δ-spanning tree across the time interval I if G(I, δ) has a spanning.
Theorem 4 Suppose the assumption A4 is satisfied. The system (16) reaches consensus if and only if there exist δ > 0 and T > 0 such that the corresponding digraph G(t) has a δ-spanning tree across any T -length time interval.
Proof: Since f ≡ 0, we have µ = 0 in Theorem 2. This completes the proof according to Theorem 3 and Theorem 2.  In this section, a numerical example is given to demonstrate the effectiveness of the presented results on synchronization of LCODEs with time-varying couplings. The Lyapunov exponents are be computed numerically. By this way, we can verify the the synchronization criterion and analyze synchronization numerically. We use the Rössler system [14,24] as the node dynamics where a = 0.165, b = 0.2, and c = 10. Figure 1 shows the dynamical behaviors of the Rössler system (18) with random initial value in [0, 1] that includes a chaotic attractor [14,24]. The network with time-varying topology we used here is NW small-world network with a time-varying coupling which was introduced as the blinking model in [11,25]. The time-varying network model algorithm is presented as follows. We divide the time axis into intervals of length τ , in each interval: (a) Begin with a nearest neighbor coupled network consisting of m nodes arranged in a ring, where each node i is adjacent to its 2k-nearest neighbor nodes; (b) Add a connection between each pair of nodes with probability p, which usually is a random number between [0, 0.1]; For more details, we refer the readers to [11]. Figure 2 shows the time-varying structure of shortcut connections in the blinking model with m = 50 and k = 3.
In this example, the parameters are set as m = 50, k = 3, τ = 1, and p = 0.04. Then small-world network can be generated with the coupling graph Laplacian L(G(t)) = −L(t). The dynamical network system can be describe as following: Let e(t) = max 1≤i<j≤50 x i (t) − x j (t) denotes the maximum distance between nodes at time t. Let E = T +R T e(t)dt, for some sufficiently large T > 0 and R > 0. Let H = µ + ς defined in Corollary 4. Figure 3 shows convergence of the maximum distance between nodes during the topology evolution with different coupling strength σ. It can be seen from Figure 3 that the dynamical network system (19) can be synchronized with σ = 0.4 and σ = 0.5.
We pick the evolution time length to be 200. Let T = 190, R = 10. And choose initial state randomly from the interval [0, 1]. Figure 4 shows the variation of E and H with respect to the coupling strength σ. It can be seen that the parameter (coupling strength σ) region where H is negative coincides with that of synchronization , i.e., where E is near zero. This verified the theoretical result (Corollary 4). In addition, we find that σ ≈ 0.38 is the threshold for synchronizing the coupled systems in this case.

Conclusions
In this paper, we present a theoretical framework for synchronization analysis of general coupled differential dynamical systems. The extended Hajnal diameter is introduced to measure the synchronization. The coupling between nodes is time-varying in both the network structure and the reaction dynamics. Inspired by the approaches in [6,13], we show that the Hajnal diameter of the linear system induced by the time-varying coupling matrix and the largest Lyapunov exponent of the synchronized system play the key roles in synchronization analysis of LCODEs. These results extend synchronization analysis of discrete-time network in [6] to continuous-time case. As an application, we obtain a sufficient condition ensuring directed time-varying graph reaching consensus which is very general, and the way we get this result is different from [23]. An example of numerical simulation is provided to show the effectiveness the theoretical results. Additional contributions of this paper are that we explicitly show that the largest projection Lyapunov exponent, the Hajanal diameter and the largest Lyapunov exponent of the transverse space equal in coupled differential systems (see lemma 6 and lemma 7), which was proved in [6] for couple discrete-time systems. according to the 3-th item of the assumption A1. Write the solution of (5) δx(t) = U (t, t 0 , s 0 )δx 0 as Then, According to Lemma 4, we have δx(t + t 0 ) ≤ δx 0 + K t 0 δx 0 e (t−τ )K dτ = e Kt δx 0 . This implies that U (t + t 0 , t 0 , s 0 ) ≤ e Kt for all s 0 ∈ W and t 0 ≥ 0.
By induction, we can conclude that χ[v jk (t)] ≤ σ k for all j > k. For j < k, χ[v jk (t)] = −∞ due to the lower-triangularity of the matrixV (t).
Proof of Lemma 8: Since L(t) satisfies assumption A 4 , if the initial condition is u(t 0 ) = 1 m , then the solution must be u(t) = 1 m , which implies that each row sum of V (t, t 0 ) is one. Then, we will prove all elements in V (t, t 0 ) are nonnegative. Consider the ith column of V (t, t 0 ) denoted by V i (t, t 0 ) which can be regarded as the solution of the following equation σl i 0 j (u j (t) − u i 0 (t)) ≥ 0. This implies that min i=1,2,··· ,m u i (t) is always nondecreasing for all t ≥ t 0 . Therefore, u i (t) ≥ 0 holds for all i = 1, 2, · · · , m and t ≥ t 0 . We can conclude that V (t, t 0 ) is a stochastic matrix. The proof is completed. Proof of Lemma 9: Consider the following Cauchy problem: σl ij (t)u j (t) u i (t 0 ) = 1 i = k 0 otherwise , i = 1, 2, · · · , m.
Noting thatu k (t) ≥ σl kk u k , we have u k (t) ≥ e −M 1 (t−t 0 ) . For each i = k, since u i (t) ≥ 0 for all i = 1, 2, · · · , m and t ≥ t 0 , we have So, if there exists a δ-edge from vertex j to i across [t 0 , t 0 + T ], then we have v ij (t 0 + T, t 0 ) ≥ e −M 1 T δ. Let δ 2 = min{e −M 1 T , e −M 1 T δ}. We can see that V (t, t 0 ) has a δ 2 spanning tree across any T -length time interval. Therefore, according to [26,27], there exist δ 1 > 0 and T 1 = (m−1)T such that V (t, t 0 ) is δ 1 scrambling across any T 1 -length time interval. The Lemma is proved.