# Persistent synchronization of heterogeneous networks with time-dependent
linear diffusive coupling

Hildeberto Jardón-Kojakhmetov<EMAIL_ADDRESS>Christian
Kuehn111CK was partly supported by a Lichtenberg Professorship of the
VolkswagenStiftung. CK was also supported by a DFG Sachbeihilfe grant
444753754.222CK and IPL also acknowledge partial support of the EU within the
TiPES project funded by the European Unions Horizon 2020 research and
innovation programme under grant agreement No. 820970<EMAIL_ADDRESS>Iacopo
P. Longo333IPL was partly supported by UKRI under the grant agreement
EP/X027651/1, by the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 754462, by
MICIIN/FEDER project PID2021-125446NB-I00, by TUM International GraduateSchool
of Science and Engineering (IGSSE) and by the University of Valladolid under
project PIP-TCESC-2020.444CK and IPL also acknowledge partial support of the
EU within the TiPES project funded by the European Unions Horizon 2020
research and innovation programme under grant agreement No. 820970.
<EMAIL_ADDRESS>
###### Abstract

We study synchronization for linearly coupled temporal networks of
heterogeneous time-dependent nonlinear agents via the convergence of
attracting trajectories of each node. The results are obtained by constructing
and studying the stability of a suitable linear nonautonomous problem bounding
the evolution of the synchronization errors. Both, the case of the entire
network and only a cluster, are addressed and the persistence of the obtained
synchronization against perturbation is also discussed. Furthermore, a
sufficient condition for the existence of attracting trajectories of each node
is given. In all cases, the considered dependence on time requires only local
integrability, which is a very mild regularity assumption. Moreover, our
results mainly depend on the network structure and its properties, and achieve
synchronization up to a constant in finite time. Hence they are quite suitable
for applications. The applicability of the results is showcased via several
examples: coupled van-der-Pol/FitzHugh-Nagumo oscillators, weighted/signed
opinion dynamics, and coupled Lorenz systems.

[rug] organization=University of Groningen, Faculty of Science and Engineering
Dynamical Systems, Geometry & Mathematical Physics — Bernoulli Institute,
addressline=Nijenborgh 9, city=Groningen, postcode=9747 AG, country=The
Netherlands [tum] organization=Technical University of Munich, School of
Computation Information and Technology, Department of
Mathematics,addressline=Boltzmannstraße 3, city=Garching bei München,
postcode=85748, country=Germany [tum1] organization=Technical University of
Munich, Munich Data Science Institute,addressline=Walther-von-Dyck-Straße 10,
city=Garching bei München, postcode=85748, country=Germany

[csh] organization=Complexity Science Hub Vienna,addressline=Josefstädter
Straße 39, city=Vienna, postcode=1080, country=Austria [icl]
organization=Imperial College London, Department of Mathematics,
DynamIC,addressline=637 Huxley Building, 180 Queen’s Gate South Kensington
Campus, city=London, postcode=SW7 2AZ, country=United Kingdom

##  1 Introduction

Within the class of complex systems in nature, technology and society, an
interconnected structure is a recurrent feature [51, 64]. Dynamical systems on
networks go a long way into successfully capturing and explaining these
structures and their evolution [51]. For example, the available dynamical
theory for static networks—dynamical systems on graphs with a static
topology—allows to study problems of opinion formation [32], collective motion
[66], and epidemic dynamics [28]. In contrast to their static counterpart,
temporal networks—sometimes also referred to as time-varying networks, dynamic
networks and temporal graphs—feature a time-dependent variation of the graph
structure. It has become increasingly evident that certain types of natural
and societal interactions require a time-dependent framework [23, 24].
Examples include human proximity and communication networks [57, 63, 69],
brain networks [3], travel and transportation networks [20], distributed
computing [30], or ecological and biological networks [8] to name just a few.
However, the analysis of temporal networks is generally much more complicated
and the respective mathematical theory remains still substantially open [21,
24].

For example, the emergence of a synchronized state in a temporal network or in
some of its parts—also referred to as clusters—has been investigated
analytically and numerically mostly on a case-to-case basis [9, 27, 68];
typically under the assumption of switched dynamics—only a finite set of
coupling structures are selected over time [1, 12, 31, 33, 67]—or of fast
switching [5, 26, 50, 49, 53, 62].

More general results usually involve the study of the errors between nodes or
with respect to a synchronizing trajectory, either through a global Lyapunov
function [6, 13], a master stability function [44, 58, 45, 48, 42],
linearization along a synchronizing trajectory [40, 41], or the Hajnal
diameter [41].

In this work, we study synchronization for linearly coupled temporal networks
of $N\geq 2$ heterogeneous time-dependent agents of the form

$\dot{x}_{i}=f_{i}(t,x_{i})+\sum_{k=1}^{N}a_{ik}(t)(x_{k}-x_{i}),\quad
x_{i}=x_{i}(t)\in\mathbb{R}^{M},\,i=1,\dots,N,$ (1.1)

where $x_{i}\in\mathbb{R}^{M}$ is the state variable of node $i$ and
$f_{i}:\mathbb{R}\times\mathbb{R}^{M}\to\mathbb{R}^{M}$ its internal dynamics,
and $A:\mathbb{R}\to\mathbb{R}^{N\times N}$, $t\mapsto
A(t)=\big{(}a_{ij}(t)\big{)}$ is the generalized (or weighted) adjacency
matrix of the network. We shall assume that $A$ is a locally integrable
function allowing, for example, for discontinuous switching in the network
topology. Moreover, we emphasize that, in contrast to the majority of the
literature, we allow $a_{ij}$ to take values over the whole real line, which
implies that signed networks [59] are covered by our theory. Furthermore, the
assumptions of regularity in time required for our work are particularly mild:
the matrix function $A$ needs to be locally integrable, while the functions
$f_{i}(t,x_{i})$ are Lipschitz Carathéodory. In other words, continuity in
time is not assumed.

As for other works on synchronization, also our analysis is based on the study
of the evolution of the errors, which, however, are considered in norm
$\xi_{ij}(t)=|x_{i}(t)-x_{j}(t)|^{2}$ for $i,j=1,\dots,N$ with $i<j$. This
choice allows us to construct a suitable bounding linear system whose
stability, in terms of spectral dichotomy, is used to infer information on the
finite-time and asymptotic synchronization of the nonlinear problem. The key
advantages of our approach can be summarized as follows:

  * •

Only very mild assumptions of regularity akin to local integrability in time
are required.

  * •

The sufficient conditions for synchronization depend almost completely (and
increasingly so in case of global coupling) on the adjacency matrix
$A(t)=\big{(}a_{ij}(t))$. Moreover, such conditions are quantitative and
constructive allowing for a control-driven approach to synchronization via the
modification of a suitable part of the network.

  * •

The results are written for signed networks and for a considerably general
class of agents. On the one hand, we consider time-dependent agents, a
surprisingly rare case in the available theory for temporal networks;
synchronization results for nonautonomous agents appear for example in [11,
47] although always with static coupling. On the other hand, the nodes can
present heterogeneous dynamics and synchronization up to a constant is still
achieved, provided that the theory applies. In case the network features
identical nodes, a result of _exact synchrony_ is achieved.

  * •

The ideas used to obtain synchronyzation can be exploited to investigate its
occurrence in either the entire network or only some of its parts (clusters),
and we provide sufficient conditions for both these cases.

  * •

Being based on the roughness of the exponential dichotomy, the results of
synchrony are robust against perturbation. This fact, which is important by
itself, has also striking consequences in the case of perturbation of static
networks. For example, we show that strongly connected static networks with
positive edges always satisfy our conditions of synchronization provided that
a sufficiently large global coupling exists and therefore the achieved
synchrony is also robust.

Note also that our problem features time dependent agents. Synchronization
results for nonautonomous agents appear for example in [11, 47] although
always with static coupling. Furthermore, our main results of synchronization
are written for heterogeneous agents guaranteeing convergence up to a constant
of the individual nodes. In case the considered network features identical
nodes, a result of _exact synchrony_ is achieved.

Our work is structured as follows. In Section 2 the notation, assumptions and
some basic notions about nonautonomous linear systems are introduced. Section
3 deals with the synchronization (up to a constant) of networks of
heterogeneous non-autonomous agents which are linearly coupled via a time-
dependent network. Our approach entails the study of the synchronization
errors $\xi_{ij}(t)=|x_{i}(t)-x_{j}(t)|^{2}$ for $i,j=1,\dots,N$ with $i<j$.
Firstly we show that the maps $\xi_{ij}(t)$ for $i,j=1,\dots,N$ with $i<j$
induce a dynamical system on the positive cone of $\mathbb{R}^{N(N-1)/2}$.
Then, we introduce our fundamental assumption (H1) on the one-sided affine
upper-bound to the difference of the vector fields of any pair of uncoupled
agents in our network. The assumption (H1) is firstly introduced pointwise and
then generalized to pairs of continuous functions in Lemma 3.3, which allows
us to upper bound the evolution of any pair of trajectories of two nodes $i$
and $j$ in Lemma 3.4. Theorem 3.6 is our main result of synchronization which
uses the results above to generate a nonautonomous linear inhomogeneous system
controlling the evolution of the synchronization errors and therefore
synchronization (up to a constant) is achieved via the analysis of stability
of this problem: we check sufficient conditions for the existence of an
exponential dichotomy with projector the identity. Two additional assumptions
are made: (H2) there is a region of the phase space such that the solutions of
each node of the network with initial condition therein, remain uniformly
ultimately bounded near a globally defined bounded solution; (H3) the
heterogeneity of the nodes can be uniformly bounded on compact intervals of
time. Some corollaries address the simpler but important and possibly more
common cases that the network is subjected to a global coupling coefficient,
and that the assumption (H3) can be changed for a stronger essential
boundedness on the whole real line. Two examples complete the section to
showcase the applicability of the obtained results of synchronization: in the
first example we consider a time-dependent network of heterogeneous van der
Pol oscillators while in the second one we show how Theorem 3.6 can induce a
control-oriented approach to achieve synchronization in a ring network
featuring a contrarian node.

The idea contained in the proof of Theorem 3.6 are further explored in Section
4 to guarantee the synchronization of just a cluster (see Theorem 4.1). This
section is completed by an example of a neural network of heterogeneous
FitzHugh-Nagumo where two leading neurons induce recurrent patterns of
synchronization on their immediate but possibly shared neighbors.

In Section 5 we briefly expand on the persistence of synchronization achieved
via the previous results by using the roughness of the exponential dichotomy.
An application to the time-dependent perturbation of synchronized static
network is also presented. Incidentally, we also show that strongly connected
static networks with positive edges satisfy the conditions of synchronization
of our result provided that a sufficiently high global coupling exists. The
example of a perturbed star-network of Lorenz systems is presented at the end
of the section to showcase the appearance and persistence of synchronization.

In Section 6, we investigate sufficient conditions for the existence of local
attractors for both the uncoupled and the coupled problems. The case of
solutions that remain only uniformly ultimately bounded near a bounded
reference trajectory is also addressed. In other words, some sufficient
results guaranteeing (H2) are given.

##  2 Notation, assumptions and preliminary definitions

By the symbol $\mathbb{R}^{d}$, with $d\in\mathbb{N}$, we will denote the
$d$-dimensional Euclidean space with the usual norm $|\cdot|$. As a rule, a
singleton $\\{\xi\\}$ in $\mathbb{R}^{d}$ will be identified with the element
$\xi$ itself. For every $i=1,\dots,d$ the $i$-th component of
$\xi\in\mathbb{R}^{d}$ will be denoted by $\xi_{i}$. Moreover, the notation
$\xi\geq 0$ means that $\xi_{i}\geq 0$ for all $i=1,\dots,d$, whereas $\xi\gg
0$ means that $\xi_{i}>0$ for every $i=1,\ldots d$. The space
$\mathbb{R}^{d}_{+}$ will denote the set of points $\xi\in\mathbb{R}^{d}$ such
that $\xi\geq 0$. This notation naturally extends to vector-valued functions.
When $d=1$, we will simply write $\mathbb{R}$ instead of $\mathbb{R}^{1}$ and
thus the symbol $\mathbb{R}_{+}$ will denote the set of non-negative real
numbers. Moreover, by $B_{r}$, we denote the closed ball of $\mathbb{R}^{d}$
centered at the origin.

The symbol $\mathbb{R}^{M\times N}$, with $N,M\in\mathbb{N}$ represents the
set of matrices of dimension $M\times N$, and given $A\in\mathbb{R}^{M\times
N}$, $A^{\top}$ will denote its transpose. The space $\mathbb{R}^{M\times N}$
will be endowed with the induced norm $\|\cdot\|$ defined by
$\|A\|=\sup_{|x|=1}|Ax|$. Oftentimes, we shall write $(\mathbb{R}^{M})^{N}$ to
denote the set of $N$-tuples of vectors in $\mathbb{R}^{M}$. Although, this
set can evidently be identified with $\mathbb{R}^{M\times N}$, the notation
$(\mathbb{R}^{M})^{N}$ will be used when it is more convenient to treat each
vector in $\mathbb{R}^{M}$ separately—typically, when they correspond to the
state of a node in a network of size $N$. In this case, for every
$i=1,\dots,N$ the $i$-th vector of $x\in(\mathbb{R}^{M})^{N}$ will be denoted
by $x_{i}$—the $i$-th node of the network.

For any interval $I\subseteq\mathbb{R}$ and any $W\subset\mathbb{R}^{N}$,
$\mathcal{C}(I,W)$ will denote the space of continuous functions from $I$ to
$W$ endowed with the sup norm $\|\cdot\|_{\infty}$. Moreover, $L^{\infty}$ and
$L^{1}_{loc}$ will denote the spaces of real functions which are respectively,
essentially bounded and locally integrable on $\mathbb{R}$, i.e. belonging to
$L^{1}(I)$ for every bounded interval $I\subset\mathbb{R}$. These latter
spaces are intended as each endowed with their usual topology.

We shall work under the most general assumptions guaranteeing existence and
uniqueness of solutions for (1.1) (see for example [10]). Specifically, the
functions $t\mapsto a_{ij}(t)$ are assumed to be _locally integrable_ ,
whereas the functions $f_{i}(t,x_{i})$ are assumed to be Lipschitz
Carathéodory.

######  Definition 2.1 (Lipschitz Carathéodory functions).

A function $f\colon\mathbb{R}\times\mathbb{R}^{M}\to\mathbb{R}^{M}$ is
Lipschitz Carathéodory (in short $f\in\mathfrak{LC}$) if it satisfies
Carathéodory conditions

  1. (C1)

$f$ is Borel measurable and

  2. (C2)

for every compact set $K\subset\mathbb{R}^{M}$ there exists a real-valued
function $m^{K}\in L^{1}_{loc}$, called _$m$ -bound_ in the following, such
that for almost every $t\in\mathbb{R}$ one has $|f(t,x)|\leq m^{K}(t)$ for any
$x\in K$;

and also a Lipschitz condition

  1. (L)

for every compact set $K\subset\mathbb{R}^{M}$ there exists a real-valued
function $l^{K}\in L^{1}_{loc}$ such that for almost every $t\in\mathbb{R}$
one has $|f(t,x_{1})-f(t,x_{2})|\leq l^{K}(t)|x_{1}-x_{2}|$ for any
$x_{1},x_{2}\in K$.

Let us also recall the notion of $L^{1}_{loc}$-boundedness. A subset $S$ of
positive functions in $L^{1}_{loc}$ is bounded if for every $r>0$ the
following inequality holds

$\sup_{m\in S}\int_{-r}^{r}m(t)\,dt<\infty\,.$

In such a case we will say that $S$ is $L^{1}_{loc}$-bounded.

######  Remark 2.2.

In many cases in the applications, the functions $m^{K}$ and $l^{K}$ appearing
in (C2) and (L) are in fact constant. This is a particular case of Lipschitz
Carathéodory functions where the $m$\- and $l$-bounds are in
$L^{\infty}\subset L^{1}_{loc}$. Consequently, the theory we present still
applies. Note also that if a set of functions is uniformly bounded in
$L^{\infty}$, then it is also $L^{1}_{loc}$-bounded.

Furthermore, we will assume that the nodes are subjected to linear diffusive
coupling through a time-dependent weighted adjacency matrix
$A:\mathbb{R}\to\mathbb{R}^{N\times N}$ so that each entry $a_{ij}(\cdot)\in
L^{1}_{loc}$. More in general, we could consider the problem

$\dot{x}_{i}=f_{i}(t,x_{i})+\sum_{k=1}^{N}a_{ik}(\omega_{t})(x_{k}-x_{i}),\quad
x_{i}\in\mathbb{R}^{M},\,\omega\in\Omega,\,i=1,\dots,N,$

where $A:\Omega\to\mathbb{R}^{N\times N}$, $\omega\mapsto
A(\omega)=\big{(}a_{ij}(\omega)\big{)}$ is the generalized (or weighted)
adjacency matrix of the network with respect to a fixed $\omega$ in a complete
metric space $\Omega$. A continuous flow
$\theta:\mathbb{R}\times\Omega\to\Omega$,
$(t,\omega)\mapsto\theta(t,\omega)=\omega_{t}$ on $\Omega$ can be considered,
that, depending on the assumptions, can be completely deterministic or random:
if $\Omega=\mathbb{R}$ and $(t,s)\mapsto\theta(t,s)=t+s$, we call the network
deterministic; if $(\Omega,\mathcal{F},\mathbb{P})$ is a probability space, we
call the network random. This formalism would allow the construction of a
continuous skew-product flow [35, 34, 37] and the use of tools from
topological dynamics to investigate, for example, the propagation of
properties of synchronization. Given the complexity of the topic at hand, we
prefer to restrict ourselves to a simpler—although less powerful—formalism, to
privilege the understanding of the conditions of synchronization, the main
focus of our work. In any case, we note that, fixed a certain
$\omega\in\Omega$, a path-wise approach through the base flow $\theta$ is also
covered by our theory by simply identifying $A(\omega_{t})$ with $A(t)$ and
dropping the dependence on $\omega$.

As a rule, we shall say that a locally absolutely continuous function
$\sigma:I\subset\mathbb{R}\to\mathbb{R}^{M\times N}$,
$t\mapsto\sigma(t)=(\sigma_{1}(t),\dots,\sigma_{n}(t))$ solves (or is a
solution of) (1.1) with initial conditions at $t_{0}\in I$ given by
$\overline{x}=(\overline{x}_{1},\dots,\overline{x}_{N})\in(\mathbb{R}^{M})^{N}$,
if for every $i=1,\dots,N$, $\sigma_{i}(\cdot)$ is a solution of the integral
problem

$x_{i}(t)=\overline{x}_{i}+\int_{t_{0}}^{t}\bigg{(}f_{i}(s,x_{i}(s))+c\sum_{j=1}^{N}a_{ij}(s)(\sigma_{j}(s)-x_{i}(s))\bigg{)}\,ds,\quad
t\in I.$

We shall denote by $x(\cdot,t_{0},\overline{x}):I\to(\mathbb{R}^{M})^{N}$ such
solution.

Our objective is to identify conditions of synchronization for (1.1). Let us
clarify what we mean by synchronization in this context.

######  Definition 2.3.

We say that (1.1) synchronizes up to a constant $\mu>0$ in a synchronization
region $E\subset\mathbb{R}\times\mathbb{R}^{M}$ if all the following
properties are satisfied

  * •

there is a function $\sigma:\mathbb{R}\to(\mathbb{R}^{M})^{N}$, that is
absolutely continuous on any compact interval and maps $t\in\mathbb{R}$ to
$\sigma(t)=(\sigma_{1}(t),\dots,\sigma_{N}(t))$ such that
$(t,\sigma_{i}(t))\in E$ for all $t\in\mathbb{R}$ and all $i=1,\dots,N$, and
$\sigma(t)$ solves (1.1),

  * •

for all $(t_{0},\overline{x})\in\mathbb{R}\times(E_{t_{0}})^{N}$, where
$E_{t_{0}}=\\{x\in\mathbb{R}^{M}\mid(t_{0},x)\in E\\}$, if the absolutely
continuous function
$x(\cdot,t_{0},\overline{x}):=(x_{1}(\cdot,t_{0},\overline{x}_{1}),\dots,x_{N}(\cdot,t_{0},\overline{x}_{N}))$
solves (1.1) with $x_{i}(t_{0},t_{0},\overline{x}_{i})=\overline{x}_{i}$ for
all $i=1,\dots,N$, then $x(\cdot,t_{0},\overline{x})$ is defined for all
$t>t_{0}$ and

$\limsup_{t\to\infty}|x_{i}(t,t_{0},\overline{x}_{i})-\sigma_{i}(t)|\leq\frac{\mu}{3},\quad\text{
for all }i=1,\dots,N,$

  * •

for all $i,j=1,\dots,N$

$\limsup_{t\to\infty}|\sigma_{i}(t)-\sigma_{j}(t)|\leq\frac{\mu}{3}.$

If the previous conditions are satisfied then, the triangular inequality
guarantees that any two trajectories of two nodes $i,j$ with initial data
$\overline{x}_{i},\overline{x}_{j}\in E_{t_{0}}$ will satisfy

$\limsup_{t\to\infty}|x_{i}(t,t_{0},\overline{x}_{i})-x_{j}(t,t_{0},\overline{x}_{j})|\leq\mu.$

In particular, if $f_{i}=f$ for all $i=1,\dots,N$, then we say that (1.1)
synchronizes in a synchronization region
$E\subset\mathbb{R}\times\mathbb{R}^{M}$ if there is a solution $s(t)$ of
$\dot{x}=f(t,x)$ defined over the whole real line and with graph in $E$ such
that for all $(t_{0},\overline{x})\in\mathbb{R}\times(E_{t_{0}})^{N}$, the
solution $x(\cdot,t_{0},\overline{x})$ of (1.1), with
$x(t_{0},t_{0},\overline{x})=\overline{x}$, is defined for all $t>t_{0}$ and

$\lim_{t\to\infty}|x_{i}(t,t_{0},\overline{x}_{i})-s(t)|=0,\quad\text{ for all
}i=1,\dots,N.$

The most important results in our work require a notion of splitting of the
extended phase space of a linear time-dependent problem $\dot{y}=A(t)y$,
$y\in\mathbb{R}^{d}$, in invariant manifolds characterized by asymptotic
exponential decay either in forward or in backward time. The notions of
exponential dichotomy, dichotomy spectrum and associated splitting fulfill
this requirement. We briefly recall them here and point the interested reader
to Siegmund [60] for all the details of the locally integrable case.

######  Definition 2.4 (Exponential dichotomy and dichotomy spectrum).

Let $A:\mathbb{R}\to\mathbb{R}^{d\times d}$ be a locally integrable function
and consider the linear system

$\dot{y}=A(t)y.$ (2.1)

An invariant projector of (2.1) is a function
$P:\mathbb{R}\to\mathbb{R}^{d\times d}$ of projections $P(t)$,
$t\in\mathbb{R}$, such that

$P(t)Y(t,s)=Y(t,s)P(s),\quad\text{for all }t,s\in\mathbb{R},$

where $Y:\mathbb{R}^{2}\to\mathbb{R}^{d\times d}$ is the principal matrix
solution of (2.1) at $s\in\mathbb{R}$, i.e. $Y(\cdot,s)$ solves (2.1) with
$Y(s,s)=Id$.

The system (2.1) is said to have an _exponential dichotomy_ on $\mathbb{R}$ if
there are an invariant projector $P(\cdot)$ and constants $\gamma>0$, $K\geq
1$ such that

$\begin{split}\big{\|}Y(t,s)P(s)\big{\|}\leq Ke^{-\gamma(t-s)},\quad\text{for
all }t\geq s,&\text{ and}\\\ \big{\|}Y(t,s)\big{(}Id-P(s)\big{)}\big{\|}\leq
Ke^{\gamma(t-s)},\quad\text{for all }t\leq s,&\end{split}$

where $Id$ is the identity matrix on $\mathbb{R}^{d\times d}$. The _dichotomy
spectrum_ of (2.1) is the set

$\Sigma(A):=\\{\alpha\in\mathbb{R}\mid\dot{y}=(A(t)-\alpha\,Id)y\ \text{ has
no exponential dichotomy}\\}.$

######  Remark 2.5.

(i) Siegmund [60] showed that either $\Sigma(A)$ is empty, or it is the whole
$\mathbb{R}$, or there exists $k\in\mathbb{N}$, with $1\leq k\leq d$, such
that

$\Sigma(A)=I_{1}\cup[a_{2},b_{2}]\cup\dots\cup[a_{k-1},b_{k-1}]\cup I_{k}\,,$

where $I_{1}$ is either $[a_{1},b_{1}]$ or $(-\infty,b_{1}]$, $I_{k}$ is
either $[a_{k},b_{k}]$ or $[a_{k},\infty)$, and $a_{1}\leq b_{1}<a_{2}\leq
b_{2}<\dots\leq a_{k}\leq b_{k}$. If $A(\cdot)$ has constant entries, the
dichotomy spectrum reduces to the real parts of the eigenvalues of $A$. In
addition, a decomposition of $\mathbb{R}\times\mathbb{R}^{d}$ in spectral
manifolds holds, i.e.

$\mathbb{R}\times\mathbb{R}^{d}=\mathcal{W}_{0}\oplus\dots\oplus\mathcal{W}_{k+1}\,;$

see [60] for details. If the dichotomy spectrum of (2.1) is contained in
$(-\infty,0)$, then (2.1) admits an exponential dichotomy with invariant
projector $P(\cdot)\equiv Id$.

(ii) The notion of exponential dichotomy should be regarded as a natural
extension to the nonautonomous framework of the tools for stability analysis
of autonomous linear systems. Indeed, it is well-known that the time-dependent
eigenvalues are of little help for the investigation of stability of
nonautomous linear problems [14, 18], unless the time-dependence is periodic
(replace eigenvalues for Floquet multipliers [22]) or slow [14, 52]. On the
other hand, Lyapunov exponents successfully encapsulate the asymptotic
stability but do not necessarily provide robustness against nonlinear
perturbations, which is instead guaranteed by the roughnees of the exponential
dichotomy [14]—the Lyapunov spectrum is a subset of the dichotomy spectrum and
if the latter is a point spectrum, then it reduces to the former [56, 16]. The
downside of the notion of exponential dichotomy resides in the difficulty of
verifying it. Some explicit calculable sufficient criteria were provided by
Coppel in terms of the time-dependent eigenvalues for bounded matrix functions
[14] and by Fink in terms of row- or column-dominance of the matrix $A(t)$
[18], which we will use later in this work. Efficient numerical approximations
have also been developed using QR and singular value decomposition [15, 16,
17, 19].

##  3 Synchronization

In this section, we provide sufficient conditions for synchronization up to a
constant via comparison arguments with a suitably constructed linear system
bounding the growth of the “errors”. Given a solution
$x:\mathbb{R}\to(\mathbb{R}^{M})^{N}$ of (1.1), we shall call
_(synchronization) errors_ , the vector $\xi(t)\in\mathbb{R}^{N(N-1)/2}$
defined by

$\xi_{ij}(t)=|x_{i}(t)-x_{j}(t)|^{2},\quad i,j=1,\dots,N,\,i<j.$

Our first observation is that given a network (1.1), the synchronization
errors induce a nonautonomous dynamical system on the positive cone
$\mathbb{R}^{N(N-1)/2}_{+}$ of $\mathbb{R}^{N(N-1)/2}$.

######  Proposition 3.1.

Given a network of $N$ nodes, $N\geq 2$, in $\mathbb{R}^{M}$ as in (1.1), the
map

$x\in(\mathbb{R}^{M})^{N}\mapsto(|x_{i}-x_{j}|^{2})_{i,j=1,\dots,N,\,i<j}\in\mathbb{R}^{N(N-1)/2}$
(3.1)

induces a nonautonomous dynamical system on $\mathbb{R}^{N(N-1)/2}_{+}$.

###### Proof.

Let $(t_{0},x_{0})\in\mathbb{R}\times(\mathbb{R}^{M})^{N}$ and consider the
solution $x(\cdot,t_{0},x_{0}):\mathbb{R}\to(\mathbb{R}^{M})^{N}$ of (1.1)
with initial data $x(t_{0})=x_{0}$. Note that, by definition of solution, the
cocycle property $x(t+s,t_{0},x_{0})=x\big{(}t,s,x(s,t_{0},x_{0})\big{)}$ is
satisfied whenever all the involved terms are well-defined. Consequently, a
local cocycle is induced also on $\mathbb{R}^{N(N-1)/2}$ through the map in
(3.1). In particular, note that the obtained trajectories verify

$\begin{split}\frac{d}{dt}|x_{i}&(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0})|^{2}\\\
&=2\langle
x_{i}(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0}),f_{i}\big{(}t,x_{i}(t,t_{0},x_{0})\big{)}-f_{j}\big{(}t,x_{j}(t,t_{0},x_{0})\big{)}\rangle\\\
&\qquad+\langle
x_{i}(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0}),\sum_{k=1}^{N}a_{ik}(t)\big{(}x_{k}(t,t_{0},x_{0})-x_{i}(t,t_{0},x_{0})\big{)}\rangle\\\
&\qquad-\langle
x_{i}(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0}),\sum_{k=1}^{N}a_{jk}(t)\big{(}x_{k}(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0})\big{)}\rangle\end{split}$

for $i,j=1,\dots,N,\,i<j$. Therefore, we have that if
$|x_{i}(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0})|^{2}=0$ for some
$i,j=1,\dots,N,\,i<j$, then
$\frac{d}{dt}|x_{i}(t,t_{0},x_{0})-x_{j}(t,t_{0},x_{0})|^{2}=0$. Hence, the
positive cone $\mathbb{R}^{N(N-1)/2}_{+}$ is invariant for the induced flow. ∎

The following pair-wise property on the decoupled dynamics will be important
for our argument.

  1. (H1)

For every $i,j=1,\dots,N$ and $r>0$ there are functions
$\alpha^{r}_{ij},\beta^{r}_{ij}\in L^{1}_{loc}(\mathbb{R},\mathbb{R})$, with
$\beta^{r}_{ij}$ non-negative, such that for almost every $t\in\mathbb{R}$

$\langle
x-y,f_{i}(t,x)-f_{j}(t,y)\rangle\leq\alpha^{r}_{ij}(t)\,|x-y|^{2}+\beta^{r}_{ij}(t),\quad\text{for
all }x,y\in B_{r}.$ (3.2)

######  Remark 3.2.

Hypothesis (H1) seems somewhat technical but it is in fact not very
restrictive. In particular, note that if the network is made of identical
nodes, i.e. if $f_{i}=:f$ for all $i=1,\dots,N$, then (H1) is trivially true
with at least $\alpha^{r}_{ij}(t)=l^{r}(t)$ being the Lipschitz coefficient of
$f$ on $B_{r}$ and $\beta_{ij}^{r}(t)=0$ for all $r>0$ and $t\in\mathbb{R}$.
Indeed, using the Cauchy-Schwarz inequality and the Lipschitz continuity, we
have that

$\langle x-y,f(t,x)-f(t,y)\rangle\leq|x-y|\cdot|f(t,x)-f(t,y)|\leq
l^{r}(t)|x-y|^{2}.$

Note also that, in practical cases, one can often choose the functions
$\alpha^{r}_{ij},\beta^{r}_{ij}\in L^{1}_{loc}(\mathbb{R},\mathbb{R})$ as
constants. It goes without saying that the theory hereby developed still
applies. In fact, this simplification allows sharper results as we highlight
in some of the corollaries to our main theorems.

In any case, it is helpful to have a rough intuition of what $\alpha_{ij}$ and
$\beta_{ij}$ in Hypothesis (H1) represent. The function $\beta_{ij}^{r}(t)$
can be interpreted as a measure of how distinct the dynamics of nodes $i$ and
$j$ is; after all, $\beta_{ij}^{r}(t)\geq 0$ vanishes when $f_{i}=f_{j}$. On
the other hand $\alpha_{ij}^{r}(t)$ tells us the “tendency of the decoupled
solutions to synchronize”. In the homogeneous scalar case, if we call
$e(t)=x(t)-y(t)$, we have that hypothesis (H1) would read as
$e(t)\dot{e}(t)\leq\alpha_{ij}^{r}(t)e(t)^{2}$. If $\alpha_{ij}^{r}(t)<0$ for
all $t\in\mathbb{R}$, indeed the “decoupled error” $e(t)$ vanishes as
$t\to\infty$. This interpretation becomes also evident when we present our
main result, Theorem 3.6. If the nodes are identical, then
$\alpha_{ij}^{r}(t)=\alpha^{r}(t)<0$ for all $t\in\mathbb{R}$ implies the
dissipativity of the uncoupled problem $\dot{x}=f(t,x)$ and therefore the
existence of a forward attracting trajectory.

Next, we show that it is possible to extend (H1) to pairs of continuous
functions obtaining the same inequality almost everywhere.

######  Lemma 3.3.

Assume that (H1) holds. Then, considered two continuous functions
$\phi,\psi:I\to B_{r}$, with $I\subseteq\mathbb{R}$ and $r>0$, we have that
for almost every $t\in I$,

$\big{\langle}\phi(t)-\psi(t),f_{i}\big{(}t,\phi(t)\big{)}-f_{j}\big{(}t,\psi(t)\big{)}\big{\rangle}\leq\alpha^{r}_{ij}(t)\,|\phi(t)-\psi(t)|^{2}+\beta^{r}_{ij}(t).$
(3.3)

###### Proof.

The proof of this statement is based on the one of Lemma 3.3 in [38]. Consider
two functions $\phi,\psi:I\to B_{r}$ and let $D=\\{s_{n}\mid
n\in\mathbb{N}\\}$ be a dense subset of $I$. From (H1) we know that given
$n\in\mathbb{N}$ there is a subset $J_{n}\subset I$ of full measure such that
for every $t\in J_{n}$,

$\big{\langle}\phi(s_{n})-\psi(s_{n}),f_{i}\big{(}t,\phi(s_{n})\big{)}-f_{j}\big{(}t,\psi(s_{n})\big{)}\big{\rangle}\leq\alpha^{r}_{ij}(t)\,|\phi(s_{n})-\psi(s_{n})|^{2}+\beta^{r}_{ij}(t).$
(3.4)

Next we consider the subset $J=\bigcap_{n=1}^{\infty}J_{n}\subset I$, also of
full measure and fix $t\in J$. From the density of $D$ we find a subsequence
$(s_{n_{k}})_{k\in N}$ such that $\lim_{k\to\infty}s_{n_{k}}=t$, and from
(3.4) we deduce that for each $k\in\mathbb{N}$

$\begin{split}\big{\langle}\phi(s_{n_{k}})-\psi(s_{n_{k}}),f_{i}\big{(}t,\phi(s_{n_{k}})\big{)}&-f_{j}\big{(}t,\psi(s_{n_{k}})\big{)}\big{\rangle}\\\
&\leq\alpha^{r}_{ij}(t)\,|\phi(s_{n_{k}})-\psi(s_{n_{k}})|^{2}+\beta^{r}_{ij}(t).\end{split}$
(3.5)

Moreover, from (L) and the continuity of the functions $\phi$ and $\psi$ we
obtain

$\displaystyle\lim_{k\to\infty}f_{i}\big{(}t,\phi(s_{n_{k}})\big{)}$
$\displaystyle=f_{i}\big{(}t,\phi(t)\big{)},\quad\text{and}\quad\lim_{k\to\infty}f_{j}\big{(}t,\psi(s_{n_{k}})\big{)}=f_{j}\big{(}t,\psi(t)\big{)},$

which together with (3.5) yields (3.3) for $t\in J$, and finishes the proof. ∎

Now, we proceed to proving a differential inequality that will be the
fundamental block to subsequently construct our comparison argument.

######  Lemma 3.4.

Assume that (H1) holds and that for every $i=1,\dots,N$ there is a bounded
absolutely continuous function $\sigma_{i}:\mathbb{R}\to\mathbb{R}^{M}$, such
that $\sigma(t)=(\sigma_{1}(t),\dots,\sigma_{N}(t))$ solves (1.1). Then, for
every $i,j=1,\dots,N$, the following inequality holds,

$\begin{split}\frac{1}{2}\frac{d}{dt}&\,|\sigma_{i}(t)-\sigma_{j}(t)|^{2}\leq\delta_{ij}(t)\,|\sigma_{i}(t)-\sigma_{j}(t)|^{2}+\beta^{\rho}_{ij}(t)+\\\\[-2.0pt]
&\qquad\quad+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}(t)-a_{ik}(t)\big{)}\,\big{(}|\sigma_{i}(t)-\sigma_{k}(t)|^{2}-|\sigma_{j}(t)-\sigma_{k}(t)|^{2}\big{)}\end{split}$

where

$\delta_{ij}(t)=\alpha^{\rho}_{ij}(t)-\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq i,j\end{subarray}}^{N}\big{(}a_{ik}(t)+a_{jk}(t)\big{)}\Big{)}.$

###### Proof.

By definition, we have that

$\begin{split}\frac{1}{2}\frac{d}{dt}|\sigma_{i}(t)-\sigma_{j}(t)|^{2}&=\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),f_{i}\big{(}t,\sigma_{i}(t)\big{)}-f_{j}\big{(}t,\sigma_{j}(t)\big{)}\big{\rangle}+\\\
&\qquad\qquad+\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),\sum_{k=1}^{N}a_{ik}(t)\big{(}\sigma_{k}(t)-\sigma_{i}(t)\big{)}\big{\rangle}+\\\
&\qquad\qquad-\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),\sum_{k=1}^{N}a_{jk}(t)\big{(}\sigma_{k}(t)-\sigma_{j}(t)\big{)}\big{\rangle}.\end{split}$

From Lemma 3.3 and the boundedness of $\sigma(t)$, we immediately have that
for almost every $t\in\mathbb{R}$,

$\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),f_{i}\big{(}t,\sigma_{i}(t)\big{)}-f_{j}\big{(}t,\sigma_{j}(t)\big{)}\big{\rangle}\leq\alpha^{\rho}_{ij}(t)\,|\sigma_{i}(t)-\sigma_{j}(t)|^{2}+\beta^{\rho}_{ij}(t),$

for some $\rho>0$. On the other hand, recalling that $2\langle
x,y\rangle=|x|^{2}+|y|^{2}-|y-x|^{2}$, the following chain of inequalities
holds true

$\begin{split}&\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),\,\sum_{k=1}^{N}\big{[}a_{ik}(t)\big{(}\sigma_{k}(t)-\sigma_{i}(t)\big{)}-a_{jk}(t)\big{(}\sigma_{k}(t)-\sigma_{j}(t)\big{)}\big{]}\big{\rangle}\\\
&\qquad=-\big{(}a_{ij}(t)+a_{ji}(t)\big{)}\,|\sigma_{i}(t)-\sigma_{j}(t)|^{2}-\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}a_{ik}(t)\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),\,\sigma_{i}(t)-\sigma_{k}(t)\big{\rangle}\\\
&\qquad\qquad\quad-\sum_{\begin{subarray}{c}k=1\\\ k\neq
i\end{subarray}}^{N}a_{jk}(t)\big{\langle}\sigma_{i}(t)-\sigma_{j}(t),\,\sigma_{k}(t)-\sigma_{j}(t)\big{\rangle}\\\
&\qquad=-\big{(}a_{ij}(t)+a_{ji}(t)\big{)}\,|\sigma_{i}(t)-\sigma_{j}(t)|^{2}+\\\
&\qquad\qquad\quad-\,\sum_{\begin{subarray}{c}k=1\\\ k\neq
j\end{subarray}}^{N}\frac{a_{ik}(t)}{2}\big{(}|\sigma_{i}(t)-\sigma_{j}(t)|^{2}+|\sigma_{i}(t)-\sigma_{k}(t)|^{2}-|\sigma_{k}(t)-\sigma_{j}(t)|^{2}\big{)}+\\\
&\qquad\qquad\quad-\sum_{\begin{subarray}{c}k=1\\\ k\neq
i\end{subarray}}^{N}\frac{a_{jk}(t)}{2}\,\big{(}|\sigma_{i}(t)-\sigma_{j}(t)|^{2}+|\sigma_{k}(t)-\sigma_{j}(t)|^{2}-|\sigma_{i}(t)-\sigma_{k}(t)|^{2}\big{)}\\\
&\qquad=-\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{N}\big{(}a_{ik}(t)+a_{jk}(t)\big{)}\Big{)}\,|\sigma_{i}(t)-\sigma_{j}(t)|^{2}\\\
&\qquad\qquad\quad+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}(t)-a_{ik}(t)\big{)}\,\big{(}|\sigma_{i}(t)-\sigma_{k}(t)|^{2}-|\sigma_{j}(t)-\sigma_{k}(t)|^{2}\big{)}\end{split}$

Gathering all the previous formulas we obtain the sought-for inequality. ∎

In order to discuss the question of synchronization for systems like (1.1), we
still need two more properties besides (H1). First of all, we shall assume
that for every $i=1,\dots,N$, there is a reference trajectory
$\sigma_{i}:\mathbb{R}\to\mathbb{R}^{M}$, and a tubular neighborhood around
it, towards which the dynamics of the node $i$ converges as time increases. If
the trajectories $\sigma_{i}(\cdot)$ are in fact (locally) attractive, then we
shall talk of synchronization of attracting trajectories. In order to
rigorously present our assumption, let us recall the notion of nonautonomous
set.

######  Definition 3.5.

A nonautonomous set is a subset of the extended phase space
$\mathcal{U}\subset\mathbb{R}\times\mathbb{R}^{M}$. A $t$-fiber of
$\mathcal{U}$ is defined as
$\mathcal{U}_{t}=\\{x\in\mathbb{R}^{M}\mid(t,x)\in\mathcal{U}\\}$. A
nonautonomous set is called forward invariant if
$x(t,t_{0},\mathcal{U}_{t_{0}})\subseteq\mathcal{U}_{t}$ for all $t>t_{0}$. In
general, $\mathcal{U}$ is said to have a topological property (such as
compactness or closedness) if each fiber of $\mathcal{U}$ has this property.
We shall say that this property is uniform if it is uniformly true for every
fiber.

Our assumption on local uniform ultimate boundedness of solutions reads as,

  1. (H2)

consider (1.1) and assume that there are
$\mathcal{U}\in\mathbb{R}\times\mathbb{R}^{M}$, $\mu\geq 0$, and for every
$i=1,\dots,N$ there is a bounded absolutely continuous function
$\sigma_{i}:\mathbb{R}\to\mathbb{R}^{M}$ such that
$\sigma(t)=(\sigma_{1}(t),\dots,\sigma_{N}(t))^{\top}\in\mathcal{U}_{t}$
solves (1.1), $\sup_{i=1,\dots,N}\|\sigma_{i}(\cdot)\|_{\infty}<\rho$ for some
$\rho>0$, and if additionally
$x(\cdot,t_{0},\overline{x})=(x_{1}(\cdot,t_{0},\overline{x}_{1}),\dots,x_{N}(\cdot,t_{0},\overline{x}_{N}))^{\top}$
solves (1.1) with initial conditions
$x_{i}(t_{0})=\overline{x}_{i}\in\mathcal{U}_{t_{0}}$, then
$x(\cdot,t_{0},\overline{x})$ is defined for all $t\geq t_{0}$ and

$\limsup_{t\to\infty}|x_{i}(t,t_{0},\overline{x}_{i})-\sigma_{i}(t)|\leq\frac{\mu}{3},\quad\text{for
all }i=1,\dots,N.$ (3.6)

In Section 6, we provide sufficient conditions for (H2) to hold true. The
assumption (H2) can be alternatively read as the existence of an inflowing
invariant open ball of diameter $\mu>0$ for the networked system. Indeed, if
this is true a globally defined bounded solution of (1.1) can be obtained as
pullback limit of any initial condition on the boundary of such ball.

Furthermore, we shall assume that heterogeneity of the nodes can be uniformly
bounded on compact intervals of time, i.e.

  1. (H3)

the set $\\{\beta_{ij}^{\rho}(\cdot+t)\in L^{1}_{loc}\mid
t\in\mathbb{R},\,i,j=1,\dots,N\\}$ is $L^{1}_{loc}$-bounded, and call

$0\leq\mu_{1}:=\sup_{\tau\in\mathbb{R}}\int_{\tau}^{\tau+1}|\beta(s)|\,ds<\infty,$

where $\beta(t)=(2\beta^{\rho}_{ij}(t))_{i,j=1,\dots,N,\,i<j}^{\top}$.

###  3.1 Synchronization up to a constant of the entire network

The following theorem is our main result of synchronization up to a constant
for heterogeneous time-dependent linearly coupled networks. Note that, due to
assumption (H2), any function $x(t)=(x_{1}(t),\dots,x_{N}(t))$ which solves
(1.1) with initial data in $\mathcal{U}$ also satisfies (3.6). Hence, in order
to investigate the synchronization of the system (1.1) for initial conditions
of the nodes in $\mathcal{U}$, it is sufficient to study the asymptotic
behavior of $|\sigma_{i}(t)-\sigma_{j}(t)|$ for all $i,j=1,\dots,N$.

######  Theorem 3.6.

Assume that (H1), (H2) and (H3) hold and fix any $M>\mu_{1}$. If there is
$t_{0}\in\mathbb{R}$ such that for all $i,j=1,\dots,N$, with $i<j$, and for
almost every $t>t_{0}$,

$\begin{split}\delta_{ij}(t):=\alpha^{\rho}_{ij}(t)-\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq i,j\end{subarray}}^{N}\big{(}a_{jk}(t)+a_{ik}(t)\big{)}\Big{)}<0\\\
\end{split}$ (3.7)

and additionally

$\overline{\gamma}:=\inf_{t\in\mathbb{R}}\min_{\begin{subarray}{c}i,j=1,\dots,N,\\\
i\neq
j\end{subarray}}\Big{\\{}\underbrace{2|\delta_{ij}(t)|-\sum_{\begin{subarray}{c}k=1\\\
k\neq
i,j\end{subarray}}^{N}\big{|}a_{jk}(t)-a_{ik}(t)\big{|}}_{=:\gamma_{ij}(t)}\Big{\\}}>-\log\left(1-\frac{\mu_{1}}{M}\right),$
(3.8)

then, for every $\varepsilon>0$ there is
$T(\varepsilon)=\frac{1}{\gamma}\ln\left(4\rho^{2}/\varepsilon\right)>0$ such
that for all $i,j=1,\dots,N$

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\varepsilon+M,\qquad\text{for
}t-t_{0}>T(\varepsilon).$

In other words, (1.1) synchronizes up to a constant in finite time.

###### Proof.

Firstly, let us introduce the continuous function
$\eta:\mathbb{R}\to\mathbb{R}_{+}$ defined by

$\eta(a)=\begin{cases}0&\text{if }a\leq 0,\\\ a&\text{if }a>0,\end{cases}$

which will be used multiple times within this proof. Consider the
synchronization errors given by the vector of $N(N-1)/2$ components

$\xi(t)=\big{(}\xi_{ij}(t)\big{)}_{i,j=1,\dots,N,\,i<j}^{\top},\qquad\text{where}\qquad\xi_{ij}(t)=|\sigma_{i}(t)-\sigma_{j}(t)|^{2}.$

Thanks to Lemma 3.4, we have that $\dot{\xi}(t)$ is an _under-function_ with
respect to the initial value problem $u^{\prime}=E(t)u+\beta(t)$,
$u(t_{0})=\xi(t_{0})$, i.e.,

$\dot{\xi}(t)\leq E(t)\xi(t)+\beta(t),\quad\text{for all
}t>t_{0},\,t_{0}\in\mathbb{R},$ (3.9)

where $\beta(t)=(2\beta^{\rho}_{ij}(t))_{i,j=1,\dots,N,\,i<j}^{\top}$ (see
(H1)) and $E(t)$ is the time-dependent square matrix defined row-wise as
follows: let us label each of the $N(N-1)/2$ rows by
$e^{ij}\in\mathbb{R}^{N(N-1)/2}$, where $i,j=1,\dots,N,\,i<j$, i.e.

$e^{ij}(t)=\big{(}e^{ij}_{lk}(t)\big{)}_{l,k=1,\dots,N,\,l<k};$

then, fixed $i,j=1,\dots,N,\,i<j$, we have that for $l,k=1,\dots,N,\,l<k$,

$e^{ij}_{lk}(t)=\begin{cases}2\delta_{ij}(t)&\text{if }l=i\text{ and
}k=j,\\\\[4.30554pt] \eta\big{(}a_{jk}(t)-a_{ik}(t)\big{)}&\text{if }l=i\text{
and }k\neq j,\\\\[4.30554pt] \eta\big{(}a_{ik}(t)-a_{jk}(t)\big{)}&\text{if
}l=j,\\\\[4.30554pt] 0&\text{otherwise}.\end{cases}$ (3.10)

Note that the entries of $E$ outside the diagonal are either positive or null.
In other words, the coefficients obtained through the inequality in Lemma 3.4,
have been further bounded from above in the sense that only the terms with
positive sign are left. Moreover, due the inequality in (3.7), the entries on
the diagonal, that is $e^{ij}_{ij}(t)=2\delta_{ij}(t)$, are strictly negative
and the matrix is row-dominant [18, Definition 7.10] for $t\geq t_{0}$ .
Consequently, the linear homogeneous problem $\dot{u}=E(t)u$ has dichotomy
spectrum contained in $(-\infty,0)$ for $t\geq t_{0}$ (and thus it admits an
exponential dichotomy with projector the identity on $[t_{0},\infty)$) [18,
Theorem 7.16)]. In fact a stronger property of exponential decay holds, i.e.,
denoted by $U(t,s)$ the principal matrix solution of $\dot{u}=E(t)u$ at
$s\in\mathbb{R}$, it holds that

$\|U(t,s)\|\leq e^{-\overline{\gamma}(t-s)},\quad\text{for all }t\geq s\geq
t_{0}.$ (3.11)

Incidentally, this means that the solution of the non-homeogeneous linear
problem $u^{\prime}=E(t)u+\beta(t)$, $u(t_{0})=\xi(t_{0})$ which is given by
the variation of constants formula

$u\big{(}t,t_{0},\xi(t_{0})\big{)}=U(t,t_{0})\xi(t_{0})+\int_{t_{0}}^{t}U(t,s)\beta(s)\,ds,$
(3.12)

where the integral is understood component-wise, is defined for all $t\geq
t_{0}$ thanks to (3.11).

We shall prove that (3.9) implies $\xi(t)\leq
u\big{(}t,t_{0},\xi(t_{0})\big{)}$ for all $t\geq t_{0}$. Considering
$\varepsilon>0$ and
$\overline{\varepsilon}=\varepsilon(1,\dots,1)\in\mathbb{R}^{N(N-1)/2}$, note
that from the continuity with respect to initial data, it is enough to prove
that $\xi(t)\ll
u\big{(}t,t_{0},\xi(t_{0})+\overline{\varepsilon}\big{)}=:u(t,\varepsilon)$
for all $t\geq t_{0}$. Assume, on the contrary, that there is a first time
$t_{1}>t_{0}$ for which the equality holds for some component. By simplicity
of notation, let the first be such a component. Then,

$\xi_{ij}(t)<u_{ij}(t,\varepsilon)\quad\text{for
}t\in[t_{0},t_{1})\qquad\text{and}\qquad\xi_{12}(t_{1})=u_{12}(t_{1},\varepsilon).$
(3.13)

Denote by $g_{\xi},g_{u}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$, the
Carathéodory functions defined by

$\begin{split}g_{\xi}(t,v)=2\delta_{12}(t)v&+\sum_{\begin{subarray}{c}k=1\\\
k\neq
1,2\end{subarray}}^{N}\Big{[}\eta\big{(}a_{2k}(t)-a_{1k}(t)\big{)}\xi_{1k}(t)+\eta\big{(}a_{1k}(t)-a_{2k}(t)\big{)}\xi_{2k}(t)\Big{]}\end{split}$
$\begin{split}g_{u}(t,v)=2\delta_{12}(t)v&+\\!\sum_{\begin{subarray}{c}k=1\\\
k\neq
1,2\end{subarray}}^{N}\\!\Big{[}\eta\big{(}a_{2k}(t)-a_{1k}(t)\big{)}u_{1k}(t,\varepsilon)+\eta\big{(}a_{1k}(t)-a_{2k}(t)\big{)}u_{2k}(t,\varepsilon)\Big{]}\end{split}$

and consider the scalar Carathéodory differential problems
$\dot{v}=g_{\xi}(t,v)$, $\dot{v}=g_{u}(t,v)$. Due to Lemma 3.4 and the
assumptions in (3.13), the following scalar differential inequalities holds
true,

$\dot{\xi}_{12}(t)\leq g_{\xi}\big{(}t,\xi_{12}(t)\big{)}\leq
g_{u}\big{(}t,\xi_{12}(t)\big{)},\quad\text{for all }t\in[t_{0},t_{1}],$

and the comparison theorem for Carathéodory scalar differential equations (see
Olech and Opial [43]) yields

$\xi_{12}(t)\leq
v\big{(}t,t_{0},\xi_{12}(t_{0})\big{)}<v\big{(}t,t_{0},\xi_{12}(t_{0})+\varepsilon\big{)}=u_{12}(t,\varepsilon)$

for every $t\in[t_{0},t_{1}]$ and in particular for $t=t_{1}$. However, this
contradicts (3.13). Hence, it must be that such $t_{1}\in\mathbb{R}$ does not
exist and $\xi(t)\ll u\big{(}t,t_{0},\xi(t_{0})+\overline{\varepsilon}\big{)}$
for all $t\geq t_{0}$. In turn, the arbitrariness on $\varepsilon>0$ gives us
the sought for ordering of vector solutions $\xi(t)\leq
u\big{(}t,t_{0},\xi(t_{0})\big{)}$ for all $t\geq t_{0}$. Then, from (3.12) we
immediately obtain that

$|\xi(t)|\leq\|U(t,t_{0})\|\,|\xi(t_{0})|+\int_{t_{0}}^{t}\|U(t,s)\|\,|\beta(s)|\,ds.$
(3.14)

Concerning the first term, we have that $\|U(t,t_{0})\|\,|\xi(t_{0})|\leq
4\rho^{2}e^{-\overline{\gamma}(t-t_{0})}$, thanks to (3.11) and (H2).
Therefore, for any $\varepsilon>0$, it holds

$\|U(t,t_{0})\|\,|\xi(t_{0})|<\varepsilon,\qquad\text{whenever
}t-t_{0}>\frac{1}{\overline{\gamma}}\ln\left(\frac{4\rho^{2}}{\varepsilon}\right).$
(3.15)

We shall, thus, analyze the second term, the integral
$\int_{t_{0}}^{t}\|U(t,s)\|\,|\beta(s)|\,ds$. Note that, thanks to (H3),
$\mu_{1}:=\sup_{\tau\in\mathbb{R}}\int_{\tau}^{\tau+1}|\beta(s)|\,ds<\infty$.
Then, we have that

$\begin{split}\int_{t_{0}}^{t}\|U(t,s)\|\,|\beta(s)|\,ds&\leq\int_{t_{0}}^{t}|\beta(s)|e^{-\overline{\gamma}(t-s)}\,ds=\int_{0}^{t-t_{0}}\\!\\!\\!|\beta(t-u)|e^{-\overline{\gamma}u}\,du\\\
&\leq\int_{0}^{\infty}\\!\\!|\beta(t-u)|e^{-\overline{\gamma}u}\,du\leq\sum_{n=0}^{\infty}\int_{n}^{n+1}\\!\\!\\!\\!|\beta(t-u)|e^{-\overline{\gamma}u}\,du\\\
&\leq\sum_{n=0}^{\infty}e^{-\overline{\gamma}n}\int_{n}^{n+1}\\!\\!\\!|\beta(t-u)|\,du=\frac{\mu_{1}}{1-e^{-\overline{\gamma}}}.\end{split}$
(3.16)

Furthermore, since by assumption $\overline{\gamma}>-\log(1-\mu_{1}/M)$, then

$\frac{\mu_{1}}{1-e^{-\overline{\gamma}}}<M.$

This inequality, together with (3.14) and (3.15) gives the sought-for result.
∎

######  Remark 3.7.

A closer look at the proof of Theorem 3.6 shows that the fundamental step was
to show that, under the given assumptions, the homogeneous linear system
$\dot{u}=E(t)u$ has dichotomy spectrum contained in $(-\infty,0)$. Yet, the
property of row-dominance is only a sufficient condition for the existence of
an exponential dichotomy with projector the identity. Therefore, it is worth
noting that, although we privileged row-dominance in order to give a set of
easily computable inequalities, other sufficient conditions may be just as
effective in guaranteeing synchronization. This also includes the weaker
requirement that the Lyapunov spectrum is contained in $(-\infty,0)$.
Numerical methods to approximate the Lyapunov and the dichotomy spectrum under
the assumption of integral separation can be found in [15, 16, 17, 19]. See
also Remark 2.5 for further details.

We notice that in our main Theorem 3.6, the synchronization error might be
largely influenced by the heterogeneity of the nodes. In the next corollary we
address such an issue provided that the network has a global coupling
strength.

######  Corollary 3.8.

Suppose that Theorem 3.6 holds for (1.1). Then, the same Theorem holds for the
system with global coupling strength

$\dot{x}_{i}=f_{i}(t,x_{i})+c\sum_{k=1}^{N}a_{ik}(t)(x_{k}-x_{i}),\qquad c\geq
1$ (3.17)

but with a synchronization error

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\varepsilon+\frac{1}{c}M.$ (3.18)

###### Proof.

The proof is a slight adaptation of the proof of Theorem 3.6, therefore, we
refer to the notation used there. First, let $\alpha_{ij}\leq 0$ in (H1). In
such a case, it suffices to consider $\alpha_{ij}=0$. It follows as in the
proof of Theorem 3.6 that the error dynamics $\dot{\xi}(t)$ is an under
function with respect to the initial value problem
$u^{\prime}=cE(t)u+\beta(t)$. In turn, this initial value problem is smoothly
equivalent, via time rescaling, to $u^{\prime}=E(t)u+\frac{1}{c}\beta(t)$. The
rest of the proof follows by the same arguments as in the proof of Theorem 3.6
but with $\beta(t)$ replaced by $\frac{1}{c}\beta(t)$. Next, we consider the
case when $\alpha_{ij}>0$. Since we assume that Theorem 3.6 holds for $c=1$,
it necessarily follows that
$\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}(t)+a_{ik}(t)\big{)}\Big{)}>0$. Let

$\begin{split}\bar{\delta}_{ij}&\coloneqq\alpha_{ij}-c\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq i,j\end{subarray}}^{N}\big{(}a_{jk}(t)+a_{ik}(t)\big{)}\Big{)}\\\
&=c\bigg{(}\frac{\alpha_{ij}}{c}-\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}(t)+a_{ik}(t)\big{)}\Big{)}\bigg{)}\leq
c\delta_{ij}.\end{split}$ (3.19)

By the previous inequality, we now have the same arguments as above for the
error dynamics $\dot{\xi}_{ij}(t)$. ∎

Analogous results to Corollary 3.8 are known, for example, for the static case
with $a_{ij}\geq 0$ [46] provided that the internal dynamics are dissipative.
Nevertheless, our result also covers the time-varying case and allows for the
weights to be negative.

Let us next present a corollary that might be of additional use in practical
cases. Assume that the function $\beta$ in (H3) is in fact essentially
bounded. Then a stronger result is available that highlights the role of the
coupling strength in the network. Particularly, a result of sharp
synchronization appears when $\|\beta(\cdot)\|_{L^{\infty}}=0$.

######  Corollary 3.9.

Under the assumptions of Theorem 3.6, if additionally $\beta(\cdot)\in
L^{\infty}$, then, for every $\varepsilon>0$ there is
$T(\varepsilon)=\frac{1}{\gamma}\ln\left(4\rho^{2}/\varepsilon\right)>0$ such
that for all $i,j=1,\dots,N$

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\varepsilon+\frac{\|\beta(\cdot)\|_{L^{\infty}}}{\overline{\gamma}},\qquad\text{for
}t-t_{0}>T(\varepsilon).$

Moreover, if the considered system has a global coupling coefficient $c>0$ as
in (3.17), then the previous inequality reads as

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\varepsilon+\frac{\|\beta(\cdot)\|_{L^{\infty}}}{c\,\overline{\gamma}},\qquad\text{for
}t-t_{0}>T(\varepsilon).$

If $f_{i}=f$ for all $i=1,\dots,N$, and (H2) holds with $\mu=0$, then (1.1)
synchronizes.

###### Proof.

The result is an easy consequence of Theorem 3.6 once (3.16) is changed for

$\begin{split}\int_{t_{0}}^{t}\|U(t,s)\|\,|\beta(s)|\,ds&\leq\|\beta(\cdot)\|_{L^{\infty}}\int_{0}^{t-t_{0}}e^{-\overline{\gamma}u}\,du=\frac{\|\beta(\cdot)\|_{L^{\infty}}}{\overline{\gamma}}\big{(}1-e^{-\overline{\gamma}(t-t_{0})}\big{)}.\end{split}$

Then, the case of global coupling is inherited from the proof of Corollary
3.8. Finally, the result of sharp synchronization for a network of identical
agents follows immediately by taking $\beta_{ij}(t)=0$ for all
$t\in\mathbb{R}$ and all $i,j=1,\dots,N$ (see also Remark 3.2). ∎

Our last corollary is concerned with the strengthening of (H1) from a local to
a global property, i.e. when the functions $\alpha^{r}_{ij},\beta^{r}_{ij}\in
L^{1}_{loc}(\mathbb{R},\mathbb{R})$ do not depend on $r>0$. Note that this is
for example the case when identical nodes with a global Lipschitz constant are
chosen.

  1. (H1enumi)

Assume that there are functions $\alpha_{ij},\beta_{ij}\in
L^{1}_{loc}(\mathbb{R},\mathbb{R})$, with $\beta_{ij}$ non-negative, such that
(H1) holds with $\alpha^{r}_{ij}=\alpha_{ij}$ and $\beta^{r}_{ij}=\beta_{ij}$
for every $r>0$ and every $i,j=1,\dots,N$.

It is immediateto check that an analogous version of Lemma 3.3 holds true also
for (H1enumi).

If (H1enumi) is in force, then the assumption (H2) can be weakened in the
sense that boundedness for the solution $\sigma$ in (H2) becomes dispensable.

  1. (H2enumi)

Assume that (H2) holds, except the solution
$\sigma(t)=(\sigma_{1}(t),\dots,\sigma_{N}(t))^{\top}\in\mathcal{U}_{t}$
(1.1), is not necessarily bounded.

######  Corollary 3.10.

Assume that (H1enumi), (H2enumi) and (H3) hold. Then with the respective
additional assumptions, the results of Theorem 3.6, and Corollaries 3.8 and
3.9 still hold true.

###### Proof.

The results are immediate once one notices that the boundedness of $\sigma(t)$
in (H2) was only needed to fix $\rho>0$ so that (H1) could be used in Lemma
3.4 for the suitable $\alpha^{\rho}_{ij},\beta^{\rho}_{ij}\in
L^{1}_{loc}(\mathbb{R},\mathbb{R})$. The same argument of Lemma 3.4 can now be
repeated to construct a differential inequality holding for any (possibly not
bounded) solution of (1.1). ∎

To highlight some the results of this section, let us present a couple of
examples.

######  Example 3.11 (Synchronization of heterogeneous van der Pol oscillators
with time-dependent perturbation parameter).

Let us consider a network of $N\geq 2$ heterogeneous van der Pol oscillators.
Each oscillator has internal dynamics given by

$\begin{split}u_{i}^{\prime}&=v_{i}+b_{i}u_{i}-\displaystyle\frac{u_{i}^{3}}{3}\\\
v_{i}^{\prime}&=-\varepsilon_{i}(t)u_{i}.\end{split}$ (3.20)

Notice that since the individual oscillators are heterogeneous, both in
amplitude and phase, we expect accordingly that $\alpha_{ij}>0$ and
$\beta_{ij}>0$ in (H1). Nevertheless, in this example, we shall show that a
strong enough coupling allows us to synchronize the oscillators. This
observation goes in hand with (H2) since for a weak enough coupling we indeed
expect the existence of a large enough inflowing invariant ball for the
coupled dynamics. Finally, since the individual nodes are all oscillators,
they differ on bounded quantities, hence (H3) holds.

For this example, we randomly choose the parameters $b_{i}\in(\frac{1}{2},1)$
from a uniform distribution, which influence the amplitude of each oscillator.
Furthermore, we let
$\varepsilon_{i}(t)=\varepsilon_{0}\left(1+\frac{1}{2}\sin(\omega_{i}t)\right)$,
with $\omega_{i}>0$ also picked uniformly at random in the interval
$\omega_{i}\in(1,2)$, and $0<\varepsilon_{0}\ll 1$. Naturally, we can choose
any other time-varying behavior of $\varepsilon(t)$ as long as it is positive
but sufficiently small, i.e., $0<\varepsilon(t)\ll 1$ for almost all $t\geq
t_{0}$. Regarding the network, we choose a piece-wise constant adjacency
matrix $A(t)=[a_{ij}(t)]$ where $a_{ij}\in\left\\{0,1\right\\}$ updates
randomly every $\Delta t$ units of time (we do ensure that the underlying
graph is always connected). Finally, we also assume that the interconnection
occurs on the fast timescales, meaning that the model we consider reads as

$\begin{split}u_{i}^{\prime}&=v_{i}+b_{i}u_{i}-\displaystyle\frac{u_{i}^{3}}{3}+\frac{c}{N}\sum_{j=1}^{N}a_{ij}(t)(u_{j}-u_{i})\\\
v_{i}^{\prime}&=-\varepsilon_{i}(t)u_{i}+\frac{c}{N}\sum_{j=1}^{N}a_{ij}(t)(v_{j}-v_{i}).\end{split}$
(3.21)

We notice that since the weights $a_{ij}$ are nonnegative and the dynamics of
the decoupled nodes are ultimately bounded, one can guarantee that Theorem
3.6, and especially Corollary 3.8, hold for $c$ large enough. In Figure 1 we
show a corresponding simulation for $N=5$ and $\Delta t=50$, while in Figure 2
we show the results for a similar setup, but with $N=100$ and $\Delta t=5$
(for practicality we only show the error dynamics for the latter one). In both
cases we compare the effect of the global coupling $c$, as described in
Corollary 3.8 and verify that the synchronization error decreases as $c$
increases. The shown pictures are representative among ten different
simulations, all of which show a similar qualitative behavior.

Figure 1: A temporal network of $N=5$ heterogeneous van der Pol oscillators,
and update time $\Delta t=50$, as described above. Each column corresponds to
a particular value of $c$ in (3.21). In the first row, we plot the projection
to the $(u,v)$-plane of the solutions. Since (3.21) is nonautonomous, the
apparent intersections are just due to the projection. The second and third
rows show the maximum of pairwise errors for each component. According to
Corollary 3.8, there is a large enough $c$ leading to synchronization. This is
verified as one goes from left to right in the plots. We notice that the
attractor, for example in the right most column, does not appear regular due
to the time-dependent random switching of the network topology. Figure 2:
Error plots analogous to those in Figure 1, corresponding to a simulation with
$N=100$ and $\Delta t=5$. Observe that, again, as $c$ increases the
synchronization error decreases, as guaranteed by Corollary 3.8.

In the next example, we exploit Theorem 3.6 to achieve synchronization in the
case when some weights in the network can be negative. We also briefly argue
how our synchronization results on a static network may persist under small
enough time-varying perturbations of the weights. The latter case is further
detailed in section 5.1.

######  Example 3.12 (Compensating contrarians for consensus on a ring
network).

Consensus dynamics is very important in several fields of science [4, 61, 65].
Let us consider a ring network with $N$ nodes as sketched in Figure 3.

$x_{1}$$x_{2}$$x_{3}$$x_{i-2}$$x_{i-1}$$x_{i}$$x_{i+1}$$x_{i+2}$$x_{N-1}$$x_{N}$
Figure 3: A ring network where, _except for the contrarian node $x_{1}$_, all
nodes interact in a bidirectional way with its nearest $2$-neighbors. The
contrarian node $x_{1}$ has a negative influence on its $2$-neighbors, while
the neighbors do not influence $x_{1}$. In this example we assume that all the
conformist weights (black arrows) are $1$, while the contrarian weights (red
arrows) are $-a$, $a>0$. In this setup, _consensus in not achieved_. We use
the results of Theorem 3.6 to determine what influence $x_{2}$ must have on
$x_{1}$ to reach consensus.

In this example, each node is scalar $x_{i}\in\mathbb{R}$ and interacts with
its nearest $2$-neighbors. We assume that the network dynamics are governed by
the well-known consensus protocol $\dot{x}=-Lx$, where $L$ denotes the
(signed) Laplacian. Component-wise, the dynamics of each node is determined by
$\dot{x}_{i}=\sum_{j=i-2}^{i+2}a_{ij}(x_{j}-x_{i})$ with $j\in\mathbb{Z}\text{
mod }N$, $a_{ii}=0$, and where $a_{ij}\neq 0$ denotes a connection from node
$j$ towards node $i$. From now on, we stick to the previously described
$2$-nearest neighbor topology.

Motivated by [25], we identify $a_{ij}>0$ with a “conformist” influence, while
$a_{ij}<0$ is referred to as a “contrarian” influence. Let us consider that
there is one contrarian node. Without loss of generality, let node $x_{1}$ be
the contrarian represented as $a_{i1}<0$ for $i=1,2,N-1,N$. Moreover, we
assume for now that neighboring nodes to $x_{1}$ do not influence it, that is
$a_{1j}=0$ for $j=1,2,N-1,N$, and that the rest of the nodes are conformists,
i.e., the remaining nonzero weights $a_{ij}$, are positive. We recall that in
the case where all the $a_{ij}$’s are positive, the dynamics of the consensus
protocol lead to convergence of the solutions to some finite value. In the
case of a signed Laplacian, where negative weights are allowed, determining
the stability of the protocol is considerably harder. In particular, in the
above setting consensus does not hold. In this example we are going to use
Theorem 3.6 to determine how strongly would node $x_{2}$ need to influence
node $x_{1}$ (that is choose $a_{12}$) to overcome $x_{1}$’s negative
influence so that the protocol can reach consensus.

For simplicity, let $a_{i1}=-a$, $a>0$, and $a_{ij}=1$ for $i,j\geq 2$. With
the notation of Theorem 3.6, we have that $\alpha_{ij}=0$ and $\mu_{1}=0$.
Moreover, notice that for all $i,j\in\left\\{4,\ldots,N-2\right\\}$ one has
$\delta_{ij}<0$. Accounting for the symmetry $\delta_{ij}=\delta_{ji}$ we have

$\begin{split}\delta_{12}&=-\left(a_{12}+\frac{3}{2}-a\right)\\\
\delta_{13}=\delta_{1N}=\delta_{1(N-1)}&=-\left(\frac{a_{12}}{2}+\frac{3}{2}-a\right)\\\
\delta_{23}=\delta_{2N}&=-\left(\frac{5}{2}-a\right).\end{split}$ (3.22)

Since $a_{12}$ has no influence on $\delta_{23}$, we impose the further
constraint that $a<\frac{5}{2}$. On the other hand, regarding $\gamma$ let

$\gamma_{ij}=2|\delta_{ij}|-\sum_{k\neq i,j}|a_{jk}-a_{ik}|.$ (3.23)

So, accounting for the symmetry mentioned above, we look to satisfy the
inequalities

$\begin{split}\gamma_{12}&=2\left|a_{12}+\frac{3}{2}-a\right|-3>0\\\
\gamma_{13}&=2\left|\frac{a_{12}}{2}+\frac{3}{2}-a\right|-|1-a_{12}|-2>0\\\
\gamma_{23}&=2\left|\frac{5}{2}-a\right|-2>0,\end{split}$ (3.24)

where, as for the $\delta_{ij}$’s above, the rest of the $\gamma_{ij}$’s are
all nonnegative. We notice that $\gamma_{12}>0$ can always be satisfied by an
appropriate choice of $a_{12}$ and that since $a_{12}$ has no influence on
$\gamma_{23}$, we have the further restriction $a<\frac{3}{2}$. However
$\gamma_{13}>0$ further imposes that $a<1$. With this restriction,
$\delta_{ij}<0$ for any $a_{12}>0$. So, from Theorem 3.6, we can conclude that
for $0<a<1$ the given consensus protocol achieves consensus provided that the
contrarian influence is compensated by $a_{12}>a$. This last inequality is
obtained from the requirement $\gamma_{12}>0$.

It is evident from the analysis performed in this example, that the weights
$a_{ij}$ can be time-varying. For more general details see section 5.1.
Moreover, the contrarian influence does not have to be homogeneous, the result
above will hold as long as all of them have modulus less than $1$. In Figure 4
we show a couple of simulations that verify our arguments with the conformist
weights fixed to $1$, the contrarian weights given by
$a_{i1}=a_{i1}(t)=-\frac{1}{2}+\frac{1}{2}\sin(\omega_{i}t)$ with uniformly
distributed random frequencies, and the conformist compensation $a_{12}=0$, on
the left, and $a_{12}=1$ on the right.

Figure 4: Simulations of the consensus protocol $\dot{x}=-Lx$, for a ring
network as in Figure 3. For all simulations, the initial conditions are
randomly chosen, the conformist weights are set to $1$, and the contrarian
weights are $a_{i1}=-\frac{1}{2}+\frac{1}{2}\sin(\omega_{i}t)$, $i=2,3,N-1,N$,
for some randomly chosen frequency $\omega_{i}>0$. In the left we show the
dynamics where no compensation is implemented, i.e. $a_{12}=0$. One can
clearly notice that not only the overall dynamics is unstable (due to the
negative weights), but there is no consensus at all. On the right we verify
that from Theorem 3.6 we have deduced that, for example, a weight $a_{12}=1$
leads to consensus.

##  4 On the synchronization of a cluster

In practical cases, the synchronization of the whole network can be sometimes
excessive and one is rather interested in the behaviour of a limited portion
of the network. Hereby, we show how our arguments can be used to analyze also
such case.

######  Theorem 4.1.

Let $2<n\leq N$ integer and consider an ordered set $J$ of $n$ nodes from
(1.1) identified by their indices, i.e.
$J=\\{j_{1},\dots,j_{n}\\}\subset\\{1,\dots,N\\}$, with $j_{m}<j_{m+1}$ for
all $m=1,\dots,n-1$. Assume that (H2) holds for (1.1), and that (H1) and (H3)
are true for all $i,j\in J$.

Furthermore, assume that $\\{a_{jk}(\cdot+t)-a_{ik}(\cdot+t)\in
L^{1}_{loc}\mid t\in\mathbb{R},\,i,j\in J,k\notin J\\}$ is
$L^{1}_{loc}$-bounded, and call

$0\leq\mu_{2}:=2\rho^{2}\max_{\begin{subarray}{c}i,j\in J,\\\ k\notin
J\end{subarray}}\sup_{\tau\in\mathbb{R}}\int_{\tau}^{\tau+1}|a_{jk}(s)-a_{ik}(s)|\,ds<\infty.$

Fixed any $M>\mu_{1}+\mu_{2}(N-n)\sqrt{2n(n-1)}$, if there is
$t_{0}\in\mathbb{R}$ such that for all $i,j\in J$ with $i<j$ and almost every
$t>t_{0}$,

$\begin{split}\delta_{ij}(t):=\alpha^{\rho}_{ij}(t)-\Big{(}a_{ij}(t)+a_{ji}(t)+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq i,j\end{subarray}}^{N}\big{(}a_{jk}(t)+a_{ik}(t)\big{)}\Big{)}<0\\\
\end{split}$

and

$\begin{split}\overline{\gamma}_{J}:=\inf_{t\in\mathbb{R}}\min_{\begin{subarray}{c}i,j=1,\dots,N,\\\
i\neq j\end{subarray}}\Big{\\{}2|\delta_{ij}(t)|&-\sum_{k\in
J\setminus\\{i,j\\}}\big{|}a_{jk}(t)-a_{ik}(t)\big{|}\Big{\\}}\\\
&>-\log\left(1-\frac{\mu_{1}+\mu_{2}(N-n)\sqrt{2n(n-1)}}{M}\right),\end{split}$

then for every $\varepsilon>0$ there is
$T(J,\varepsilon)=\frac{1}{\overline{\gamma}_{J}}\ln\left(4\rho^{2}/\varepsilon\right)>0$
such that

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\varepsilon+M,\qquad\text{for
}t-t_{0}>T(J,\varepsilon).$

In other words, $J$ synchronizes up to a constant in finite time.

###### Proof.

For the sake of simplicity, we shall assume that $J=\\{1,\dots,n\\}$. Then, we
can proceed akin to the proof of Theorem 3.6 but now considering the vector
$\zeta(t):=\big{(}|\sigma_{i}(t)-\sigma_{j}(t)|^{2}\big{)}_{i,j\in
J,\,i<j}^{\top}$. Thanks to Lemma 3.4, we have that $\dot{\zeta}(t)$ is an
_under-function_ with respect to the initial value problem
$u^{\prime}=C(t)u+\nu_{J}(t)$, $u(t_{0})=\zeta(t_{0})$, i.e.,

$\dot{\zeta}(t)\leq C(t)\zeta(t)+\nu_{J}(t),\quad\text{for all
}t>t_{0},\,t_{0}\in\mathbb{R},$

where $\nu_{J}(t)=(\nu_{ij}(t))_{i,j\in J,\,i<j}^{\top}$ is defined by

$\nu_{ij}(t)=\eta\left(2\beta^{\rho}_{ij}(t)+\sum_{k\notin
J}\big{(}a_{jk}(t)-a_{ik}(t)\big{)}\,\big{(}\xi_{ik}(t)-\xi_{jk}(t)\big{)}\right),$
(4.1)

and $C(t)$ is the time-dependent matrix constructed exactly as the matrix
$E(t)$ in the proof of Theorem 3.6 where now $n$ appears in place of $N$
everywhere. Now, note that, under the given assumptions, $\dot{u}=C(t)u$
admits an exponential dichotomy on $[t_{0},\infty)$ with projector the
identity and exponential rate of convergence $\gamma_{J}>0$; as for Theorem
3.6, the linear homogeneous system $\dot{u}=C(t)u$ is row-dominant for $t\geq
t_{0}$ [18, Theorem 7.16)]. Then, with analogous reasoning, we obtain that for
all $t\geq t_{0}$,

$\begin{split}\big{(}|\sigma_{i}(t)-&\sigma_{j}(t)|^{2}\big{)}_{i,j=1,\dots,n,\,i<j}^{\top}\\\
\leq&U(t,t_{0})\big{(}|\sigma_{i}(t_{0})-\sigma_{j}(t_{0})|^{2}\big{)}_{i,j=1,\dots,n,\,i<j}^{\top}+\int_{t_{0}}^{t}U(t,s)\nu(s)\,ds,\end{split}$

where $U(t,t_{0})$ is the principal matrix solution at $t_{0}$ of
$\dot{u}=C(t)u$. Then, fixed $\varepsilon>0$ and reasoning as in the proof of
Theorem 3.6, but recalling that $\nu$ is defined by (4.1), we come to the
analogous chain of inequalities as (3.16), that is,

$\begin{split}\int_{t_{0}}^{t}&\|U(t,s)\|\,\Big{(}|\beta(s)|+\sqrt{\frac{n(n-1)}{2}}\max_{i,j\in
J}\sum_{k\notin
J}|a_{jk}(s)-a_{ik}(s)|\,|\xi_{ik}(t)-\xi_{jk}(t)|\Big{)}\,ds\\\
&\leq\int_{0}^{t-t_{0}}\\!\\!\Big{(}|\beta(t-u)|+4\rho^{2}\sqrt{\frac{n(n-1)}{2}}\max_{i,j\in
J}\sum_{k\notin
J}|a_{jk}(t-u)-a_{ik}(t-u)|\Big{)}e^{-\overline{\gamma}_{J}u}\,du\\\
&\leq\sum_{n=0}^{\infty}e^{-\overline{\gamma}_{J}n}\big{(}\mu_{1}+\mu_{2}(N-n)\sqrt{2n(n-1)}\big{)}=\frac{\mu_{1}+\mu_{2}(N-n)\sqrt{2n(n-1)}}{1-e^{-\overline{\gamma}_{J}}}.\end{split}$

Therefore, one has for every $i,j=1,\dots,n$,

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}\leq\varepsilon+\frac{\mu_{1}+\mu_{2}(N-n)\sqrt{2n(n-1)}}{1-e^{-\overline{\gamma}_{J}}}<\varepsilon+M,$

whenever
$t-t_{0}>T(J,\varepsilon)=\frac{1}{\overline{\gamma}_{J}}\ln\left(4\rho^{2}/\varepsilon\right)>0$,
which concludes the proof. ∎

Also in the case of cluster synchronization, a sharper result can be obtained
by considering a stronger assumption than the uniform
$L^{1}_{loc}$-boundedness in the statement of Theorem 4.1, i.e. boundedness in
$L^{\infty}$.

######  Corollary 4.2.

Under the assumptions of Theorem 4.1, if additionally $|\beta(\cdot)|\in
L^{\infty}$, and also $a_{jk}(\cdot),a_{ik}(\cdot)\in L^{\infty}$ for all
$i,j\in J$ and $k\notin J$, then, for every $\varepsilon>0$ there is
$T(\varepsilon)=\frac{1}{\gamma_{J}}\ln\left(4\rho^{2}/\varepsilon\right)>0$
such that for all $i,j=1,\dots,N$

$|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\varepsilon+\frac{1}{\overline{\gamma}_{J}}\Big{(}\|\beta(\cdot)\|_{L^{\infty}}+2\rho^{2}(N-n)\sqrt{2n(n-1)}\max_{\begin{subarray}{c}i,j\in
J,\\\ k\notin
J\end{subarray}}\|a_{jk}(\cdot)-a_{ik}(\cdot)\|_{L^{\infty}}\Big{)},$

for $t-t_{0}>T(\varepsilon)$.

###### Proof.

The result is an easy consequence of Theorem 4.1 reasoning as for Corollary
3.9. ∎

The next example highlights how the up-to-a-constant synchronization achieved
in finite time through Theorem 4.1, can be used to produce recurrent patterns
of cluster synchronization alternated by intervals with no synchrony.

######  Example 4.3 (Clustering in a FitzHugh-Nagumo network).

To showcase Theorem 4.1, let us consider a network of heterogeneous FitzHugh-
Nagumo neurons. The $i$-th neuron has dynamics given by

$\begin{split}\dot{x}_{i}&=c_{i}x_{i}-x_{i}^{3}-y_{i}+I_{i}\\\
\dot{y}_{i}&=\varepsilon(x_{i}+a_{i}-b_{i}y_{i}),\end{split}\qquad\qquad
i=1,\ldots,N.$ (4.2)

In this model $x_{i}$ represents the $i$-th membrane’s voltage and $y_{i}$ the
$i$-th recovery variable [55, 54]. Regarding the parameters, we have that
$c_{i}>0$ modulates the amplitude of oscillations, $I_{i}$ accounts for the
stimulus current, and $\varepsilon$ stands for the difference in timescales
between the voltage and the recovery variables. The parameters $a_{i}>0$ and
$b_{i}>0$, together with $I_{i}$, determine whether the neuron is in
excitatory or in refractory mode. For this example we let the neurons be
slightly heterogoeneous, and unless otherwise stated, we randomly assign
parameter values according to the following: $c_{i}\in[0.75,1]$,
$a_{i}\in[-0.3,\,0.3]$, $b\in[0.1,2]$, $I_{i}\in[0,.01]$ and
$\varepsilon=0.05$. We notice that with these parameters, isolated neurons may
or may not oscillate. We moreover have similar arguments as in example 3.11
regarding (H1), (H2), and (H3).

We setup the network as follows: an underlying connected, unweighted, and
directed graph of $N$ nodes is randomly generated by selecting an adjacency
matrix $A=[a_{ij}]_{i,j=1,\ldots,N}$ with $a_{ij}\in\left\\{0,1\right\\}$,
from a (discrete) uniform distribution. Below we shall specify a time-varying
change in the weights of this underlying network, but we emphasize that no new
edges are created. Then, two nodes are chosen at random. Let $(x_{l},y_{l})$
and $(x_{k},y_{k})$ be such neurons; we set their parameters to
$(c_{l},I_{l},a_{l},b_{l})=(0.5,0.1,0.3,1.4)$ and
$(c_{k},I_{k},a_{k},b_{k})=(0.75,0.15,0.3,1.4)$, which ensures that, in
isolation, the $l$ and $k$ neurons are oscillating. Furthermore, we call _the
neighbors of neuron $l$ (resp. $k$)_ the nodes $i$ for which there is a
directed edge $a_{il}$ from $l$ to $i$ (resp. $a_{ik}$ from $k$ to $i$). The
set of neighbors of neuron $l$ (resp. $k$) is denoted by $\mathcal{N}_{l}$
(resp. $\mathcal{N}_{k}$).

Figure 5: A representative simulation for $N=15$ nodes. As described in the
main text, within the black/red shaded time intervals, neighboring neurons
synchronize with the dashed black/red neuron. This is because during such
time-frames, the conditions of Theorem 4.1 hold. Indeed one can particularly
observe that, since $k=12$ is a neighbor of $l=9$, the $k$-th neuron (dashed
red) has larger amplitude according to its own parameters during the red time-
frames, but synchronizes with the smaller amplitude oscillator $l=9$ (dashed
black) during the black time-frames. We further notice that during the overlap
of the time-frames, some trajectories also seem to synchronize. This, however,
is not characterized in the example, and may very well depend on further
connectivity properties. Nevertheless, notice that all neighbors of $k=12$,
except for $1$ and $15$, are also neighbors of $l$. During the overlap of the
time-frames we hence see a common cluster that does not include the
aforementioned neighbors.

Next, the (nonzero) weights of the network are set to
$a_{ij}=\frac{1}{100}$555where $A=[a_{ij}]$ is the adjacency matrix of the
previously randomly generated graph, and we identify a weighted edge from $j$
to $i$ with $a_{ij}$. If there is no edge from node $j$ to node $i$, then
$a_{ij}=0\,\forall t$., and the weights of the outgoing edges of the $l$ and
$k$ neurons are updated periodically as follows:

$a_{il}(t)=\begin{cases}\bar{a},&\sin(\omega_{l}t)\geq 0\\\
\frac{1}{100},&\sin(\omega_{l}t)<0\end{cases},\qquad
a_{ik}(t)=\begin{cases}\bar{a},&\sin(\omega_{k}t)\geq 0\\\
\frac{1}{100},&\sin(\omega_{k}t)<0,\end{cases}$ (4.3)

for some positive frequencies $\omega_{l}$ and $\omega_{k}$. This example has
been setup so that whenever $a_{il}=\bar{a}$ is sufficiently large and
$a_{ik}=\frac{1}{100}$ (or viceversa), Theorem 4.1 holds for $J$ being the
neuron $l$ together with its neighbors $\mathcal{N}_{l}$ (or viceversa).

Figure 6: Neurons that are not neighbors of the $l$ nor the $k$ neurons.

In figure 5 we present a representative simulation for $N=15$ and $\bar{a}=3$.
The frequencies $\omega_{l}$ and $\omega_{k}$ have been chosen so that during
the time interval with a black (resp. red) background, the cluster is formed
by the $l$ (resp. $k$) neuron, shown as the black (resp. red) dashed curve,
and its neighbors. So, indeed notice that along the “black interval” the $l$
neuron and its neighbors form a cluster while outside such intervals no
synchronization seems to ensue (the same for the red intervals and the $k$
neuron). Outside the black and red intervals, the network is weakly connected
with $a_{ij}=\frac{1}{100}$, except for the first $50$ time-units where, for
comparison, the network is disconnected. The black and red intervals overlap
differently because for this simulation we chose $\omega_{l}$ and $\omega_{k}$
incommensurate with each other. In figure 6 we show the nodes that are not
neighbors of either the $l$ or the $k$ neurons.

##  5 Persistence of synchronization

One of the interesting features of Theorems 3.6 and 4.1 is the property of
inherent robustness of the achieved synchronization against perturbation of
both, the dynamics of the individual nodes, and of the adjacency matrix. The
underlying reason is the roughness of the exponential dichotomy guaranteeing
synchronization in Theorems 3.6 and 4.1. In this sense, the following result
reminds in spirit the one in [47], although only static networks are therein
considered. In this brief subsection, we aim to make these relations more
explicit. For simplicity of notation, we shall treat the case of the entire
network, although the same ideas can be applied also to the synchronization of
clusters as we briefly highlight in Example 5.4.

######  Theorem 5.1.

Let $f:\mathbb{R}\times\mathbb{R}^{M}\to\mathbb{R}^{M}$ be a $\mathfrak{LC}$
function and $A:\mathbb{R}\to\mathbb{R}^{N\times N}$ be locally integrable,
and consider the network

$\dot{x}_{i}=f(t,x_{i})+\sum_{k=1}^{N}a_{ik}(t)(x_{k}-x_{i}),\quad
x_{i}\in\mathbb{R}^{M},\,i=1,\dots,N,$

Moreover, assume that the assumptions of Theorem 3.6 are satisfied with
$\mu=0$ in (H2), and thus sharp synchronization of the entire network is
achieved. The following statements are true.

  * •

For every $\delta>0$, such that if
$f_{i}:\mathbb{R}\times\mathbb{R}^{M}\to\mathbb{R}^{M}$, for $i=1,\dots,N$,
are $\mathfrak{LC}$ functions and

$\sup_{\begin{subarray}{c}x\in\mathbb{R}^{M},\,t\in\mathbb{R}\end{subarray}}|f(t,x)-f_{i}(t,x)|<\delta,$

the perturbed network (1.1) synchronizes up to a constant
$4\rho\delta\sqrt{N(N-1)/2}/\overline{\gamma}$, provided that condition (H2)
is still satisfied with $\mu=0$.

  * •

If $B:\mathbb{R}\to\mathbb{R}^{N\times N}$, defined by
$B(t)=(b_{ij}(t))_{i,j=1,\dots,N}$ is locally integrable and
$\sup_{t\in\mathbb{R}_{+}}|B(t)|<\overline{\gamma}/4$, then the perturbed
network

$\dot{x}_{i}=f(t,x_{i})+\sum_{k=1}^{N}\big{(}a_{ik}(t)+b_{ik}(t)\big{)}(x_{k}-x_{i}),\quad
x_{i}\in\mathbb{R}^{M},\,i=1,\dots,N,$

achieves sharp synchronization.

###### Proof.

In order to prove the first statement, note that for all $t\in\mathbb{R}$,

$\begin{split}\langle x-y,&f_{i}(t,x)-f_{j}(t,y)\rangle=\langle
x-y,f_{i}(t,x)-f(t,x)\rangle+\\\ &\qquad\qquad\qquad+\langle
x-y,f(t,x)-f(t,y)\rangle+\langle x-y,f(t,y)-f_{j}(t,y)\rangle\\\ &\leq
l^{r}(t)|x-y|^{2}+2r\delta,\qquad\text{for all }x,y\in B_{r}.\end{split}$

For the previous chain of inequality we have used Remark 3.2, Cauchy-Schwarz
inequality and the assumption of $f_{i},f_{j}\in\mathfrak{LC}$. Therefore (H1)
and (H3) hold true with $\beta_{ij}(t)=2r\delta$ for all $t\in\mathbb{R}$ and
all $i,j=1,\dots,N$, while (H2) is satisfied by assumption. Hence, Corollary
3.9 applies. Noting that
$\|\beta(\cdot)\|_{L^{\infty}}=4\rho\delta\sqrt{N(N-1)/2}$, one has that for
all $i,j=1,\dots,N$,

$\lim_{t\to\infty}|\sigma_{i}(t)-\sigma_{j}(t)|^{2}<\frac{4\rho\delta\sqrt{N(N-1)/2}}{\overline{\gamma}}.$

The second statement is a direct consequence of the roughness of the
exponential dichotomy and it is obtained applying [14, Proposition 4.1]. ∎

###  5.1 An application to time-dependent perturbations of static networks

It is well-known that strongly diffusely coupled, static networks of identical
nodes locally synchronize provided that the global coupling overcomes a
certain threshold, see e.g., [46, 47]. The case of nonidentical nodes, under
further constraints, has also been considered, see for example [70]. In this
section we discuss how our theory relates to this classic result both in terms
of synchronization and time-dependent perturbation. Furthermore, we showcase
such relations by means of an example on a star network at the end of the
section.

Firstly, we show consistency of the two theories: the inequalities of Theorem
3.6 are always verified by static strongly connected networks of identical
nodes (and as a matter of fact even more general static networks), provided
that the global coupling overcomes a certain threshold.

######  Corollary 5.2.

Consider a static connected network of identical nodes satisfying (H2) and
(H3),

$\dot{x}_{i}=f(x_{i})+c\sum_{k=1}^{N}a_{ik}(x_{k}-x_{i}),x_{i}\in\mathbb{R}^{M},\,i=1,\dots,N,$

with $a_{ij}\in\mathbb{R}$ for all $i,j=1,\dots,N$ and global coupling
strength $c>0$. Assume, furthermore, that
$2(a_{ij}+a_{ji})+\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}\Big{(}a_{jk}+a_{ik}-\big{|}a_{jk}-a_{ik}\big{|}\Big{)}>0$
for every $i,j=1,\dots,N$—e.g. this is true if the graph is strongly connected
and $a_{ij}\geq 0$ for all $i,j=1,\dots,N$. Then, there is $\overline{c}>0$
such that for $c>\overline{c}$ the assumptions of Theorem 3.6 are satisfied
and the network achieves synchronization.

###### Proof.

The assumptions on the edges of the network imply that
$a_{ij}+a_{ji}+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}+a_{ik}\big{)}>0$ for all $i,j=1,\dots,N$.
Hence, there is $c_{1}>0$ such that for any pair $i,j=1,\dots,N$ with $i<j$,

$\delta_{ij}=l^{\rho}-c\Big{(}a_{ij}+a_{ji}+\frac{1}{2}\sum_{\begin{subarray}{c}k=1\\\
k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}+a_{ik}\big{)}\Big{)}<0,\qquad\text{for
}c>c_{1},$

where $l^{\rho}$ is the Lipschitz coefficient for $f$ on $B_{\rho}$ and the
same notation of Theorem 3.6 has been used. Since $\beta_{ij}$ can be taken
equal to zero for identical nodes (see Remark 3.2), we also have that for
$c>c_{1}$,

$2|\delta_{ij}|-c\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\big{|}a_{jk}-a_{ik}\big{|}=2c(a_{ij}+a_{ji})-2l^{\rho}+c\sum_{\begin{subarray}{c}k=1\\\
k\neq
i,j\end{subarray}}^{N}\big{(}a_{jk}+a_{ik}-\big{|}a_{jk}-a_{ik}\big{|}\big{)}.$

Again, the assumptions on the network topology imply that a $\overline{c}\geq
c_{1}$ exists such that both inequalities of Theorem 3.6 are satisfied for
$c>\overline{c}$ and in such a case the system synchronizes. ∎

Of further interest, is the fact that the synchronization achieved beyond the
coupling threshold presented in the previous corollary is robust against small
perturbations of the type discussed in Theorem 5.1.

######  Corollary 5.3.

Consider the static network of identical nodes,

$\dot{x}_{i}=f(x_{i})+c\sum_{k=1}^{N}a_{ik}(x_{k}-x_{i}),x_{i}\in\mathbb{R}^{M},\,i=1,\dots,N,$

with coupling coefficients $a_{ik}\geq 0$ for all $i,k=1,\dots,N$ and global
coupling strength $c$ greater than the coupling threshold $\overline{c}$ in
Corollary 5.2. Then, any sufficiently small (possibly time-dependent)
perturbation in the sense of Theorem 5.1 will also produce synchronization.

###### Proof.

The result is a direct consequence of Theorem 5.1 and Corollary 5.2. ∎

######  Example 5.4 (Synchronization of chaotic oscillators on star networks).

Although dynamics on undirected star networks can be studied, for example, by
spectral methods at least in the static case, we can also use Theorem 3.6 and
the Corollaries 5.2 and 5.3 to find conditions that lead to synchronization on
time-varying directed ones.

We consider the following setup: the central node (the hub), we call it
$x_{1}$, has directed outgoing edges with weights $a_{i1}=a\in\mathbb{R}$,
$i=2,\ldots,N$, with the rest of the network, while the remainder of the nodes
(the leaves, or satellites) have corresponding weight $a_{1j}=b\in\mathbb{R}$,
$j=2,\ldots,N$; every other weight $a_{ij}$ is zero due to the considered
graph structure—a schematic representation is shown in Figure 7. We shall
hereby assume that all the nodes have identical dynamics (hence
$\beta^{r}_{ij}=0$ for all $r>0$ and all $i,j=1,\dots,N$), and that a global
attractor exists and $\alpha^{r}_{ij}(t)<\alpha\in\mathbb{R}$ for all
$i,j=1,\dots,N$ in a sufficiently big ball of radius $r>0$ containing the
global attractor. Notice, however, that if $\alpha_{ij}(t)\leq 0$ for all
$i,j=1,\ldots,N$ and all $t\in\mathbb{R}$, then it suffices to let $\alpha=0$.

$x_{1}$$x_{N}$$x_{2}$$x_{3}$$x_{4}$$x_{5}$$x_{6}$$x_{7}$$a$$b$ Figure 7: A
representation of the star network setup used in Example 5.4. Node $x_{1}$ is
set as the hub and it has directed outgoing edges with positive weight
$a_{i1}=a>0$, $i=2,\ldots,N$. Every other node $x_{j}$, for $j=2,\ldots,N$ has
only one directed outgoing edge towards $x_{1}$, which has weight
$a_{1j}=b<0$.

Under the described setting, one finds that, for any pair of satellites, i.e.
$1<i<j$,

$0<2(a_{ij}+a_{ji})+\sum_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\Big{(}a_{jk}+a_{ik}-\big{|}a_{jk}-a_{ik}\big{|}\Big{)}=2a\quad\Leftrightarrow\quad
a>0,$ (5.1)

whereas, for the hub ($i=1$) and any other satellite ($j>1$),

$\begin{split}0<2(a_{1j}+a_{j1})&+\sum_{\begin{subarray}{c}k=1\\\ k\neq
1,j\end{subarray}}^{N}\Big{(}a_{jk}+a_{1k}-\big{|}a_{jk}-a_{1k}\big{|}\Big{)}=2(b+a)+(N-2)(b-|b|),\end{split}$
(5.2)

and the previous inequality is satisfied if and only if

$\text{either}\quad\leavevmode\hbox to16.14pt{\vbox
to16.14pt{\pgfpicture\makeatletter\hbox{\hskip 8.07158pt\lower-8.07158pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{7.87158pt}{0.0pt}\pgfsys@curveto{7.87158pt}{4.3474pt}{4.3474pt}{7.87158pt}{0.0pt}{7.87158pt}\pgfsys@curveto{-4.3474pt}{7.87158pt}{-7.87158pt}{4.3474pt}{-7.87158pt}{0.0pt}\pgfsys@curveto{-7.87158pt}{-4.3474pt}{-4.3474pt}{-7.87158pt}{0.0pt}{-7.87158pt}\pgfsys@curveto{4.3474pt}{-7.87158pt}{7.87158pt}{-4.3474pt}{7.87158pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{-3.41666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{A}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\begin{cases}b<0,\\\
a>-b(N-1),\end{cases}\qquad\text{or}\qquad\leavevmode\hbox to15.9pt{\vbox
to15.9pt{\pgfpicture\makeatletter\hbox{\hskip 7.9522pt\lower-7.9522pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{7.7522pt}{0.0pt}\pgfsys@curveto{7.7522pt}{4.28146pt}{4.28146pt}{7.7522pt}{0.0pt}{7.7522pt}\pgfsys@curveto{-4.28146pt}{7.7522pt}{-7.7522pt}{4.28146pt}{-7.7522pt}{0.0pt}\pgfsys@curveto{-7.7522pt}{-4.28146pt}{-4.28146pt}{-7.7522pt}{0.0pt}{-7.7522pt}\pgfsys@curveto{4.28146pt}{-7.7522pt}{7.7522pt}{-4.28146pt}{7.7522pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.54167pt}{-3.41666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{B}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\begin{cases}b\geq
0,\\\ a\geq-b.\end{cases}$

Since we assume $a>0$ (due to (5.1)), then B holds always true, while in case
A ($b<0$), some care needs to be taken in choosing $a$ depending on the size
of the network and the modulus of $b$. It is also straightforward to extend
the previous arguments to the time-varying case $a_{ij}=a_{ij}(t)$, provided
that the $a_{ij}$’s are bounded since it would suffice that for almost all
$t\geq t_{0}$,

$\begin{split}a\leq\min_{i=2,\ldots,N}\left\\{a_{i1}(t)\right\\}\,\qquad\text{and}\qquad
b\leq\min_{j=2,\ldots,N}\left\\{a_{1j}(t)\right\\}.\end{split}$ (5.3)

In conclusion, on a star network and provided that (H1)-(H3) hold, no matter
the (time-varying, possibly negative) influence of the satellites, the hub can
always induce synchronization if $a>0$ is sufficiently big. Moreover the
synchronization error can be made as small as desired if a global coupling
$c>0$ intervenes—see Corollaries 3.8 and 5.2.

To verify our arguments, we consider a directed star network where the nodes
are Lorenz systems, that is, the internal dynamics of each node is given by

$\begin{split}\dot{x}_{i}&=\sigma(y_{i}-x_{i})\\\ \dot{y}_{i}&=x_{i}(\rho-
z_{i})-y_{i}\\\ \dot{z}_{i}&=x_{i}y_{i}-\beta z_{i},\end{split}$ (5.4)

where, unless otherwise is stated, the parameters are $\sigma=10$, $\rho=28$,
$\beta=\frac{8}{3}$. We focus on the less immediate case of negative edges and
set the influence of the leaves to the hub to $a_{1j}=b=-1$, $j=2,\ldots,N$.
Whenever a simulation depends on a random choice of parameters, we always show
a representative among dozens of simulations. All simulations have been
performed in Matlab using ODE45 with initial conditions randomly chosen near
the origin. For the simulations we shall use an error

$\hat{e}=\sqrt{e_{x}^{2}+e_{y}^{2}+e_{z}^{2}}$ (5.5)

defined by

$e_{\zeta}=e_{\zeta}(t)=\max_{ij}|\zeta_{i}(t)-\zeta_{j}(t)|,\qquad\zeta=x,y,z,\;i=1,\ldots,N,$
(5.6)

corresponding to the maximum of pairwise errors at each time $t$.

First, in Figure 8 we show simulations for $N=5$ nodes, for three different
scenarios, see the description in the caption. The common feature is that
following the analysis presented above, synchronization of the network can be
achieved.

Figure 8: Simulations of three distinct numerical experiments for a star
network of $N=5$ Lorenz oscillators. For convenience we show only the time
series of the $z$-coordinate on the left and the corresponding error $\hat{e}$
(5.5) on the right. In the first plot, we let the network be disconnected for
the first 10 time units. At $t=10$ we connect the network with $a_{i1}=a=5$
and $a_{1j}=b=-1$ (see A above), and global coupling $c=2$. With these
parameters, Corollary 5.2 guarantees synchronization, as verified in the plot
for $t\geq 10$. To showcase the persistence of synchronization, the second
plot shows a similar setting as the first, but with all weights perturbed as
$a_{ij}\mapsto a_{ij}(1+\frac{1}{10}\sin(\omega_{ij}t))$ for some randomly
chosen frequencies $\omega_{ij}\in(\pi,2\pi)$. Note that for the considered
static example $1/10<\min_{i,j=1,\dots,N}\overline{\gamma}_{ij}/4$ for all
$t>10$. Therefore, Corollary 5.3 applies. Finally, the third plot shows the
case where $a=a(t)=4+3\tanh\left(\frac{1}{5}(t-10)\right)$. Hence this plot
shows the transition from (5.2) not holding $(a<4)$ to where it does. Notice
that since (5.1) holds, we observe that first the satellites tend to form a
cluster, to later transition to full synchronization.

Finally, we consider a large network of heterogeneous Lorenz systems. For this
we let $N=200$, and all the leaves have randomly picked parameters
$\sigma_{i}\in(\sigma-1,\sigma+1)$, $\rho_{i}\in(\rho-1,\rho+1)$,
$\beta_{i}\in(\beta-1,\beta+1)$. Similar to the second experiment in Figure 8,
we let $a_{1j}=b\left(1+\frac{1}{10}\sin(\omega_{1j}t)\right)$, $b=-1$, and
$a_{j1}=-2(N-1)b\left(1+\frac{1}{10}\sin(\omega_{j1}t)\right)$; for some
randomly set frequencies $\omega_{ij}\in(\pi,2\pi)$, and $c=1$. For
presentation purposes we prefer to show in Figure 9 only the corresponding
error (5.5).

Figure 9:  Synchronization error (5.5) for a temporal star network with
$N=200$ heterogeneous Lorenz oscillators. In this figure, for the first $10$
time units, the network is disconnected. Afterwards, for $t\geq 10$ the
network is connected as described above so that Theorem 3.6 holds, leading to
synchronization.

##  6 Existence of local attractors

In this final section, we provide a sufficient condition for the existence of
an attracting trajectory for each equation of (1.1) so that (H2) is satisfied.
The treatment largely owes to the theory developed by Caraballo et al. [11]. A
substantial generalization intervenes in the sense that now all the considered
inequalities involve locally integrable functions in place of constants and
more than two coupled systems are considered. This is in line with the rest of
the treatment in our work. Our fundamental inspiration for this type of
generalization is the work by Longo et al. [37].

###  6.1 The uncoupled problem

Given $f\in\mathfrak{LC}$, we shall consider the following assumptions,

  1. (SL1)

(_local one-sided Lipschitz continuity_) there exists a nonempty forward
invariant uniformly bounded nonautonomous set
$\mathcal{U}\subset\mathbb{R}\times\mathbb{R}^{M}$ and a function $l\in
L^{1}_{loc}$ such that for almost every $t\in\mathbb{R}$,
$B_{r}\subset\mathcal{U}_{t}$ and,

$\qquad 2\langle x_{1}-x_{2},f(t,x_{1})-f(t,x_{2})\rangle\leq
l(t)|x_{1}-x_{2}|^{2},\ \ \text{for all }x_{1},x_{2}\in\mathcal{U}_{t}.$ (6.1)

  1. (SL2)

(_dissipativity_) given (SL1), there are constants $K\geq 1$ and $\gamma>0$
such that

$\exp\left(\int_{t_{0}}^{t}l(s)\,ds\right)\leq
Ke^{-\gamma(t-t_{0})},\quad\text{for any }t_{0}\leq t.$

For practical reasons, we shall assume that the origin, denoted by the vector
$\mathbf{0}\in\mathbb{R}^{M}$, belongs to $\mathcal{U}_{t}$ for almost every
$t\in\mathbb{R}$. This hypothesis is not restrictive: let us assume that there
is a set $S\subset\mathbb{R}$ with positive measure for which
$\mathbf{0}\notin\mathcal{U}_{t}$ for all $t\in S$. Since $\mathcal{U}$ is
nonempty, uniformly bounded and forward invariant there is at least one entire
solution $\zeta(t)$ whose graph is contained in $\mathcal{U}$. Then, the time
dependent change of variables $(t,y)=(t,x-\zeta(t))$ returns a new vector
field $\widetilde{f}$ and a new forward invariant uniformly bounded
nonautonomous set $\widetilde{\mathcal{U}}$ for which
$\mathbf{0}\in\widetilde{\mathcal{U}}_{t}$ for almost all $t\in\mathbb{R}$.
This fact will be important in the proof of some of the following results.

Next, we show that condition (SL1) can be extended to pairs of continuous
functions and the inequality (6.1) keeps holding almost everywhere.

######  Proposition 6.1.

Let $f\in\mathfrak{LC}$ and assume that (SL1) holds. If $I\subset\mathbb{R}$
is an interval and $\phi,\psi\in C(I,\mathbb{R}^{M})$ with
$\phi(t),\psi(t)\in\mathcal{U}_{t}$ for almost every $t\in I$, then

$2\langle\phi(t)-\psi(t),f\big{(}t,\phi(t)\big{)}-f\big{(}t,\psi(t)\big{)}\rangle\leq
l(t)|\phi(t)-\psi(t)|^{2}$

for almost every $t\in I$.

###### Proof.

A proof of this statement can be obtained reasoning as for Lemma 3.3. ∎

Note that (SL1) guarantees that each pair of solutions $x(t),y(t)$ of the
uncoupled problem $\dot{x}=f(t,x)$ with initial conditions in $\mathcal{U}$
will converge in forward time. Indeed, one immediately has that

$\frac{d}{dt}|x(t)-y(t)|^{2}\leq l(t)|x(t)-y(t)|^{2}.$

Hence, using (SL2)we have that for any $t,t_{0}\in\mathbb{R}$, with $t\geq
t_{0}$, and $x_{0},y_{0}\in\mathcal{U}_{t_{0}}$,

$|x(t)-y(t)|^{2}\leq Ke^{-\gamma(t-t_{0})}|x_{0}-y_{0}|^{2}.$

In order to understand towards what these solutions are converging, we firstly
have to show that (SL1) implies a more standard dissipative condition.

######  Proposition 6.2.

Let $f\in\mathfrak{LC}$ and assume that (SL1) holds. Then, $f$ is
asymptotically dissipative, that is, for almost every $t\in\mathbb{R}$, it
satisfies

$2\langle x,f(t,x)\rangle\leq\alpha(t)|x|^{2}+\beta(t),\quad\text{for all
}x\in\mathcal{U}_{t},$ (6.2)

where $\alpha(\cdot),\beta(\cdot)\in L^{1}_{loc}$ and there is
$0<\overline{\gamma}<\gamma$ such that

$\exp\left(\int_{t_{0}}^{t}\alpha(s)\,ds\right)\leq
Ke^{-\overline{\gamma}(t-t_{0})},\quad\text{for any }t_{0}\leq t.$ (6.3)

###### Proof.

Firstly notice that

$\begin{split}\langle
x,y\rangle&=\frac{1}{2}(|x+y|^{2}-|x|^{2}-|y|^{2})\leq\frac{1}{2}(|x|^{2}+|y|^{2}+2|x||y|-|x|^{2}-|y|^{2})\\\
&=\frac{1}{2}\big{(}|x|^{2}+|y|^{2}-(|x|-|y|)^{2}\big{)}\leq|x|^{2}+|y|^{2}.\end{split}$

Fix $0<\varepsilon<\gamma/2$ and $x\in\mathcal{U}_{t}$. Since (SL1) holds, we
have that

$\begin{split}2\langle x,f(t,x)\rangle&\leq 2\langle
x,f(t,\mathbf{0})\rangle+l(t)|x|^{2}\leq
2\varepsilon|x|^{2}+\frac{2}{\varepsilon}|f(t,\mathbf{0})|^{2}+l(t)|x|^{2}\\\
&=\big{(}2\varepsilon+l(t)\big{)}|x|^{2}+\frac{2}{\varepsilon}|f(t,\mathbf{0})|^{2}.\end{split}$

Denoted $\alpha(t):=2\varepsilon+l(t)$ and
$\beta(t):=2/\varepsilon|f(t,\mathbf{0})|^{2}$, note that we have

$\exp\left(\int_{t_{0}}^{t}\alpha(s)\,ds\right)=e^{2\varepsilon(t-t_{0})}\exp\left(\int_{t_{0}}^{t}l(s)\,ds\right)\leq
Ke^{(2\varepsilon-\gamma)(t-t_{0})},\quad\text{for any }t_{0}\leq t,$

which concludes the proof with $\overline{\gamma}=\gamma-2\varepsilon$. ∎

It is now possible to prove that any uncoupled system $\dot{x}=f(t,x)$
satisfying (SL1) admits a bounded local pullback and forward attracting
trajectory, provided that the set,

$\\{\beta_{t}(\cdot)\\}_{t\in\mathbb{R}}=\\{|f(t+\cdot,\mathbf{0})|^{2}\mid
t\in\mathbb{R}\\},$

is $L^{1}_{loc}$-bounded.

######  Proposition 6.3.

Let $f\in\mathfrak{LC}$, satisfying (SL1) and (SL2) Moreover, let
$\alpha,\beta\in L^{1}_{loc}$ be the functions provided by Proposition 6.2. If
$\\{\beta_{t}(\cdot)\\}_{t\in\mathbb{R}}$ is $L^{1}_{loc}$-bounded, then
$\dot{x}=f(t,x)$ has a bounded local pullback attractor made of a single
globally defined trajectory which is also forward attracting.

###### Proof.

Thanks to Proposition 6.2 we have that for any $t\in\mathbb{R}$, $s>0$ and
$x_{0}\in\mathcal{U}_{t-s}$,

$|x(t,t-s,x_{0})|^{2}\leq
Ke^{-\overline{\gamma}s}|x_{0}|^{2}+K\int_{t-s}^{t}\beta(u)e^{-\overline{\gamma}(t-u)}\,du.$

Reasoning as for (3.16), and taking $x_{0}\in\mathcal{U}_{t-s}$, one obtains
that

$|x(t,t-s,x_{0})|^{2}\leq
Ke^{-\overline{\gamma}s}r^{2}+\frac{K\mu}{1-e^{-\overline{\gamma}}},$

where
$\mu:=\sup_{\tau\in\mathbb{R}}\int_{\tau}^{\tau+1}|\beta(s)|\,ds<\infty$, and
$\mathcal{U}_{s}\subset B_{r}$ for some $r>0$ by assumption. Therefore, the
ball of radius $\varepsilon+K\mu/(1-e^{-\overline{\gamma}})$ for
$\varepsilon>0$ chosen as small as desired, pullback absorbs in finite time
the fiber $\mathcal{U}_{t-s}$. This fact implies that there is a universe of
attraction contained in $B_{r}$ that is pullback absorbed into the ball of
radius $\varepsilon+K\mu/(1-e^{-\overline{\gamma}})$. Therefore, a unique
bounded pullback attractor exists [29, Theorem 3.27]. On the other hand,
Proposition 6.1 and the remark thereafter imply forward convergence of all the
trajectories of $\dot{x}=f(t,x)$. Therefore, one immediately has that all the
sections of the pullback attractor are singleton sets. This fact concludes the
proof. ∎

If instead of (SL1) and (SL2) the weaker conditions (6.2) and (6.3) are
assumed, a weaker result of uniform ultimate boundedness of solution is
obtained.

######  Corollary 6.4.

Let $f\in\mathfrak{LC}$, satisfying (6.2) and (6.3). If
$\\{\beta_{t}(\cdot)\\}_{t\in\mathbb{R}}$ is $L^{1}_{loc}$-bounded, then
$\dot{x}=f(t,x)$ has a bounded local pullback attractor and the rest of
trajectores starting in $\mathcal{U}$ are uniformly ultimately bounded.

###### Proof.

The existence of a unique bounded local pullback attractor is achieved
following the same initial steps in the proof of Proposition 6.3. Now consider
$t>t_{0}$. Thanks to Proposition 6.2 and taking $x_{0}\in\mathcal{U}_{t_{0}}$,
one obtains that

$|x(t,t_{0},x_{0})|^{2}\leq
Ke^{-\overline{\gamma}(t-t_{0})}|x_{0}|^{2}+K\int_{t_{0}}^{t}\beta(u)e^{-\overline{\gamma}(t-u)}\,du.$

Hence, reasoning as for (3.16), and taking $x_{0}\in\mathcal{U}_{t_{0}}$, one
obtains that

$|x(t,t_{0},x_{0})|^{2}\leq
Ke^{-\overline{\gamma}(t-t_{0})}r^{2}+\frac{K\mu}{1-e^{-\overline{\gamma}}},$

and for any $\varepsilon>0$ there is a time $T(\varepsilon)>0$ such that if
$t-t_{0}>T(\varepsilon)$ then

$|x(t,t_{0},x_{0})|^{2}\leq\varepsilon+\frac{K\mu}{1-e^{-\overline{\gamma}}},$

which completes the proof. ∎

###  6.2 The coupled problem

Now we turn our attention to the coupled problem (1.1). The natural question
is if, despite the coupling, each node still has a pullback and forward
attracting trajectory. An important role in our argument shall be played by
the $L^{1}_{loc}$ linear operator defined for any $t\in\mathbb{R}$ by

$\mathbb{R}^{M}\ni u\mapsto(2A(t)-Id_{N\times N}L(t))u$

where $L(t)=(l^{\rho}_{i}(t))^{\top}_{i=1,\dots,N}$, and for each
$i=1,\dots,N$, $l_{i}(\cdot)$ is the $L^{1}_{loc}$ function provided by (SL1).
As for the uncoupled case, we firstly deal with the forward attractivity.

######  Lemma 6.5.

Consider $f_{i}\in\mathfrak{LC}$ for $i=1,\dots,N$ satisfying (SL1), and the
networked system (1.1) where, for almost every $t\in\mathbb{R}$,
$a_{ij}(t)\geq 0$ for all $i,j=1,\dots,N$, with $i\neq j$ and $a_{ii}(t)=0$
for all $i=1,\dots,N$. Then, considered any pair of absolutely continuous
functions $y(t)$ and $z(t)$ respectively solving (1.1) with initial data
$\overline{y},\overline{z}\in\mathcal{U}^{i}_{t_{0}}$ for some
$t_{0}\in\mathbb{R}$, it holds that

$\big{(}|y_{i}(t)-z_{i}(t)|^{2}\big{)}^{\top}_{i=1,\dots,N}\leq
U(t,t_{0})(|\overline{y}_{i}-\overline{z}_{i}|^{2})^{\top}_{i=1,\dots,N}$

where $U(t,t_{0})$ is the principal matrix solution of
$\dot{u}=(2A(t)-Id_{N\times N}L(t))u$ at $t_{0}\in\mathbb{R}$.

###### Proof.

Note that for every $i=1,\dots,N$ and using the Cauchy-Schwarz inequality,

$\begin{split}\frac{d}{dt}|y_{i}(t)&-z_{i}(t)|^{2}=2\big{\langle}y_{i}(t)-z_{i}(t),f_{i}\big{(}t,y_{i}(t)\big{)}-f_{i}\big{(}t,z_{i}(t)\big{)}\big{\rangle}+\\\
&\qquad\qquad\,+2\big{\langle}y_{i}(t)-z_{i}(t),\sum_{k=1}^{N}a_{ik}(t)\big{(}y_{k}(t)-y_{i}(t)-z_{k}(t)+z_{i}(t)\big{)}\big{\rangle}\\\
&\leq
l_{i}(t)|y_{i}(t)-z_{i}(t)|^{2}-2\sum_{k=1}^{N}a_{ik}(t)\big{\langle}y_{i}(t)-z_{i}(t),y_{i}(t)-z_{i}(t)\big{\rangle}+\\\
&\qquad\qquad\,+2\sum_{k=1}^{N}a_{ik}(t)\big{\langle}y_{i}(t)-z_{i}(t),y_{k}(t)-z_{k}(t)\big{\rangle}\\\
&\leq
l_{i}(t)|y_{i}(t)-z_{i}(t)|^{2}-2\sum_{k=1}^{N}a_{ik}(t)|y_{i}(t)-z_{i}(t)|^{2}+\\\
&\qquad\qquad\,+2\sum_{k=1}^{N}a_{ik}(t)\big{[}|y_{i}(t)-z_{i}(t)|^{2}+|y_{k}(t)-z_{k}(t)|^{2}\big{]}\\\
&=l_{i}(t)|y_{i}(t)-z_{i}(t)|^{2}+2\sum_{k=1}^{N}a_{ik}(t)|y_{k}(t)-z_{k}(t)|^{2}.\end{split}$

Therefore, we have that for $t_{0}\leq t$ where it is well-defined, the vector
$(|y_{i}(t)-z_{i}(t)|^{2})^{\top}_{i=1,\dots,N}$ is an _under-function_ with
respect to the initial value problem $u^{\prime}=(2A(t)-Id_{N\times N}L(t))u$,
$u(s)=|y(s)-z(s)|^{2}$, where $L(t)=(l^{\rho}_{i}(t))^{\top}_{i=1,\dots,N}$.
Therefore, we can reason as for the proof of Theorem 3.6 to obtain that for
all $t>t_{0}$,

$\big{(}|y_{i}(t)-z_{i}(t)|^{2}\big{)}^{\top}_{i=1,\dots,N}\leq
U(t,t_{0})(|y_{i}(t_{0})-z_{i}(t_{0})|^{2})^{\top}_{i=1,\dots,N}$

where $U(t,s)$ is the principal matrix solution of $\dot{u}=(2A(t)-Id_{N\times
N}L(t))u$ at $s\in\mathbb{R}$. Incidentally, this shows that $y(t)$ and $z(t)$
are defined for all $t\geq t_{0}$. ∎

Lemma 6.5 allows to immediately obtain a sufficient condition for forward
convergence of the trajectories of each node. If the linear differential
problem $\dot{u}=2(A(t)-Id_{N\times N}L(t))u$ has dichotomy spectrum contained
in $(-\infty,0)$, then each node has a bounded forward attracting trajectory.

######  Theorem 6.6.

Consider $f_{i}\in\mathfrak{LC}$ for $i=1,\dots,N$ and assume that they all
satisfy (SL1), each within a forward invariant uniformly bounded nonautonomous
set $\mathcal{U}^{i}\subset\mathbb{R}\times\mathbb{R}^{M}$. Moreover, assume
that $\dot{u}=2(A(t)-Id_{N\times N}L(t))u$ has dichotomy spectrum contained in
$(-\infty,0)$. Then, there is an absolutely continuous function
$\sigma:\mathbb{R}\to\mathbb{R}^{M\times N}$,
$t\mapsto\sigma(t)=\big{(}\sigma_{i}(t)\big{)}_{i=1,\dots,N}$ that solves
(1.1) and such that, for every $i=1,\dots,N$, $\sigma_{i}(\cdot)$ is pullback
attracting for

$\dot{x}_{i}=f_{i}(t,x_{i})+\sum_{j=1}^{N}a_{ij}(t)(\sigma_{j}(t)-x_{i}),$

and if $y(t)=\big{(}y_{i}(t)\big{)}_{i=1,\dots,N}$ solves (1.1) with
$y_{i}(t_{0})\in\mathcal{U}^{i}_{t_{0}}$ for some $t_{0}\in\mathbb{R}$, then
for every $i=1,\dots,N$,

$\lim_{t\to\infty}|y_{i}(t)-\sigma_{i}(t)|=0.$

###### Proof.

Since $\dot{u}=2(A(t)-Id_{N\times N}L(t))u$ has dichotomy spectrum contained
in $(-\infty,0)$, we have that there are $K\geq 1$ and $\gamma>0$, such that

$|U(t,t_{0})|\leq Ke^{-\gamma(t-t_{0})}\quad\text{for }t_{0}\leq t,$

where $U(t,t_{0})$ is the principal matrix solution of
$\dot{u}=(2A(t)-Id_{N\times N}L(t))u$ at $t_{0}\in\mathbb{R}$. Therefore,
thanks to Lemma 6.5, we have that for any $t_{0}\in\mathbb{R}$ and
$\overline{y}=(\overline{y}_{i})_{i=1,\dots,N},\overline{z}=(\overline{z}_{i})_{i=1,\dots,N}\in\mathbb{R}^{M}$
with $\overline{y}_{i},\overline{z}_{i}\in\mathcal{U}^{i}_{t_{0}}$ for all
$i=1,\dots,N$, the solutions
$y(t,t_{0},\overline{y})=\big{(}y_{i}(t,t_{0},\overline{y}_{i})\big{)}_{i=1,\dots,N}$
and
$z(t,t_{0},\overline{z})=\big{(}z_{i}(t,t_{0},\overline{z}_{i})\big{)}_{i=1,\dots,N}$
of (1.1) with initial data $y(t_{0})=\overline{y}$ and
$z(t_{0})=\overline{z}$, respectively, are defined for all $t\geq t_{0}$ and
in particular,

$\big{(}|y_{i}(t,t_{0},\overline{y}_{i})-z_{i}(t,t_{0},\overline{z}_{i})|^{2}\big{)}^{\top}_{i=1,\dots,N}\leq
Ke^{-\gamma(t-t_{0})}(|\overline{y}_{i}-\overline{z}_{i}|^{2})^{\top}_{i=1,\dots,N},$

for any $t>t_{0}$. Particularly, since each nonautonomous set
$\mathcal{U}^{i}$, $i=1,\dots,N$, is uniformly bounded, then there is $r>0$
such that $\mathcal{U}^{i}_{t}\subset B_{\sqrt{r}/2}$ for all $t\in\mathbb{R}$
and $i=1,\dots,N$. Then, fixed $\varepsilon>0$, and considered
$s\geq\log(rK/\varepsilon)/\gamma=:T(r,\varepsilon)$, one has that for any
$t\in\mathbb{R}$,

$|y_{i}(t,t-s,\overline{y}_{i})-z_{i}(t,t-s,\overline{z}_{i})|^{2}\leq\varepsilon\quad\text{for
$s\geq T(r,\varepsilon)$}.$

Consequently, the function $\zeta(t)=(0,\dots,0)^{\top}\in\mathbb{R}^{M}$
pullback and forward attracts all the trajectories
$\big{(}|y_{i}(t,t-s,\overline{y}_{i})-z_{i}(t,t-s,\overline{z}_{i})|^{2}\big{)}^{\top}_{i=1,\dots,N}$
of the dynamical system induced by Proposition 3.1. In particular, this
implies that there is a globally defined and absolutely continuous function
$\sigma:\mathbb{R}\to\mathbb{R}^{M\times N}$ that solves (1.1) and satisfies
the thesis. ∎

As we have recalled in Remark 2.5, row-dominance is a sufficient condition for
the existence of an exponential dichotomy. Therefore, the stronger
dissipativity condition than (SL1) guarantees the persistence of an attractor
also for the coupled network as we show in the next corollary. Particularly,
we shall henceforth assume that

  1. (SL1enumi)

for every $i=1,\dots,N$ there exists a function $l_{i}\in L^{1}_{loc}$ such
that, (6.1) is verified and

$\sup_{\begin{subarray}{c}t\in\mathbb{R},\\\
i=1,\dots,N\end{subarray}}l_{i}(t)<0\quad\text{and}\quad\gamma=\sup_{\begin{subarray}{c}t\in\mathbb{R},\\\
i=1,\dots,N\end{subarray}}\bigg{\\{}|l_{i}(t)|-2\sum_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}^{N}a_{ik}(t)\bigg{\\}}>0.$

######  Corollary 6.7.

Consider the networked system (1.1) where, for almost every $t\in\mathbb{R}$,
$a_{ij}(t)\geq 0$ for all $i,j=1,\dots,N$, with $i\neq j$ and $a_{ii}(t)=0$
for all $i=1,\dots,N$, and assume further that $f_{i}\in\mathfrak{LC}$
satisfies (SL1enumi) for all $i=1,\dots,N$. Then, the thesis of Theorem 6.6
holds true.

###### Proof.

The assumption (SL1enumi) implies that the linear problem
$\dot{u}=(2A(t)-Id_{N\times N}L(t))u$ is row dominant and that its dichotomy
spectrum is contained in $(-\infty,0)$. Therefore the thesis is a direct
consequence of Theorem 6.6. ∎

##  7 Conclusion and Discussion

In this paper we have considered general linearly and diffusely coupled
networks. First, the nodes themselves, with internal dynamics $f_{i}(t,x_{i})$
are of the Lipschitz Carathéodoy class, meaning that, in particular,
continuity in time is not even required. Second, the interaction between nodes
is also quite general. In essence we consider temporal interconnections,
described by $a_{ij}(t)\in\mathbb{R}$, that are locally integrable. This
includes (but it is not restricted to), for example, piece-wise continuous
changes in the topology of the network. Our main goal has been to provide
quantitative conditions under which synchronization of the nodes is achieved.
One of the main highlights of the results we provide is that they mainly
depend on the network structure, i.e., on the $a_{ij}$’s. This offers
important advantages compared to some other criteria that can, for example,
depend on spectral properties, and at the same time allows for a control
approach to synchronization. Among the results we have presented we emphasize
the synchronization (up to a constant) and the synchronization of clusters, of
temporal networks. We have presented some further results, for example for
networks with global couplings, where the synchronization error can be made
arbitrarily small. Another striking feature of the developed sufficient
results of synchronization is their robustness against perturbation due to the
roughness of the exponential dichotomy on which they are based.

A limitation of the presented theory is that we require that all the
components of (higher dimensional) nodes are connected in a uniform way—no
inner-coupling matrix is considered; although it might be possible to also
include this case by using weighted inner products for our proofs in analogy
with the idea of the Mahalanobis distance, which is often employed in
statistics in such an inner coupling between variables. In summary, as we
exemplified through the paper, several quite general and practically relevant
situations for synchronization can be covered by our theory. As future work,
one could attempt to improve the presented results by allowing that only some
of the components of the nodes are to be connected. This is reminiscent of
under-actuated control laws. Another interesting extension would include
nonlinear interaction functions, as well as fully adaptive networks where the
network topology depends on the node dynamics and vice versa node dynamics
depends on the topology [7].

Finally, we notice that since a topological theory for the construction of
continuous flows for Carathéodory differential equations exists for ordinary
[2, 34, 35, 37], delay [36, 38] and parabolic [39] differential problems, not
only the ideas here exposed could be extended to these contexts but also
additional results of propagation of synchronization in the hull could be
explored.

## References

  * [1] R. Amritkar and C.-K. Hu, Synchronized state of coupled dynamics on time-varying networks, Chaos: An Interdisciplinary Journal of Nonlinear Science, 16 (2006), p. 015117.
  * [2] Z. Artstein, Topological dynamics of an ordinary differential equation, Journal of Differential Equations, 23 (1977), pp. 216–223.
  * [3] D. S. Bassett, N. F. Wymbs, M. P. Rombach, M. A. Porter, P. J. Mucha, and S. T. Grafton, Task-based core-periphery organization of human brain dynamics, PLoS computational biology, 9 (2013), p. e1003171.
  * [4] L. Becchetti, A. Clementi, and E. Natale, Consensus dynamics: An overview, ACM SIGACT News, 51 (2020), pp. 58–104.
  * [5] I. V. Belykh, V. N. Belykh, and M. Hasler, Blinking model and synchronization in small-world networks with a time-varying coupling, Physica D: Nonlinear Phenomena, 195 (2004), pp. 188–206.
  * [6] V. N. Belykh, I. V. Belykh, and M. Hasler, Connection graph stability method for synchronized coupled chaotic systems, Physica D: nonlinear phenomena, 195 (2004), pp. 159–187.
  * [7] R. Berner, T. Gross, C. Kuehn, J. Kurths, and S. Yanchuk, Adaptive dynamical networks, arXiv:2304.05652, (2023), pp. 1–.
  * [8] B. Blonder, T. W. Wey, A. Dornhaus, R. James, and A. Sih, Temporal dynamics and network analysis, Methods in Ecology and Evolution, 3 (2012), pp. 958–972.
  * [9] S. Boccaletti, D.-U. Hwang, M. Chavez, A. Amann, J. Kurths, and L. M. Pecora, Synchronization in dynamical networks: Evolution along commutative graphs, Physical Review E, 74 (2006), p. 016102.
  * [10] A. Bressan and B. Piccoli, Introduction to the mathematical theory of control, vol. 1, American institute of mathematical sciences Springfield, 2007\.
  * [11] T. Caraballo, P. E. Kloeden, and A. Neuenkirch, Synchronization of systems with multiplicative noise, Stochastics and Dynamics, 8 (2008), pp. 139–154.
  * [12] M. Cenk Eser, E. S. Medeiros, M. Riza, and A. Zakharova, Edges of inter-layer synchronization in multilayer networks with time-switching links, Chaos: An Interdisciplinary Journal of Nonlinear Science, 31 (2021), p. 103119.
  * [13] M. Chen, Synchronization in time-varying networks: a matrix measure approach, Physical Review E, 76 (2007), p. 016104.
  * [14] W. A. Coppel, Dichotomies in stability theory, vol. 629, Springer, 2006\.
  * [15] L. Dieci and E. S. Van Vleck, Lyapunov and other spectra: a survey, Collected lectures on the preservation of stability under discretization (Fort Collins, CO, 2001), 197 (2002), p. 218.
  * [16] , Lyapunov spectral intervals: theory and computation, SIAM Journal on Numerical Analysis, 40 (2002), pp. 516–542.
  * [17] , Lyapunov and Sacker–Sell spectral intervals, Journal of dynamics and differential equations, 19 (2007), pp. 265–293.
  * [18] A. M. Fink, Almost periodic differential equations, vol. 377, 1974.
  * [19] G. Froyland, T. Hüls, G. P. Morriss, and T. M. Watson, Computing covariant lyapunov vectors, Oseledets vectors, and dichotomy projectors: A comparative numerical study, Physica D: Nonlinear Phenomena, 247 (2013), pp. 18–39.
  * [20] R. Gallotti and M. Barthelemy, The multilayer temporal network of public transport in Great Britain, Scientific data, 2 (2015), pp. 1–8.
  * [21] D. Ghosh, M. Frasca, A. Rizzo, S. Majhi, S. Rakshit, K. Alfaro-Bittner, and S. Boccaletti, The synchronized dynamics of time-varying networks, Physics Reports, 949 (2022), pp. 1–63.
  * [22] J. K. Hale, Ordinary Differential Equations, Wiley-Interscience, New York, 1969.
  * [23] P. Holme, Modern temporal network theory: a colloquium, The European Physical Journal B, 88 (2015), pp. 1–30.
  * [24] P. Holme and J. Saramäki, Temporal networks, Physics reports, 519 (2012), pp. 97–125.
  * [25] H. Hong and S. H. Strogatz, Conformists and contrarians in a Kuramoto model with identical natural frequencies, Phys. Rev. E, 84 (2011), p. 046202.
  * [26] R. Jeter and I. Belykh, Synchronization in on-off stochastic networks: Windows of opportunity, IEEE Transactions on Circuits and Systems I: Regular Papers, 62 (2015), pp. 1260–1269.
  * [27] R. Jeter, M. Porfiri, and I. Belykh, Dynamics and control of stochastically switching networks: beyond fast switching, Temporal Network Theory, (2019), pp. 269–304.
  * [28] I. Kiss, J. Miller, and P. Simon, Mathematics of Epidemics on Networks: From Exact to Approximate Models, Springer, 2017.
  * [29] P. E. Kloeden and M. Rasmussen, Nonautonomous dynamical systems, no. 176, American Mathematical Soc., 2011.
  * [30] F. Kuhn, N. Lynch, and R. Oshman, Distributed computation in dynamic networks, in Proceedings of the forty-second ACM symposium on Theory of computing, 2010, pp. 513–522.
  * [31] D. Liberzon and A. S. Morse, Basic problems in stability and design of switched systems, IEEE control systems magazine, 19 (1999), pp. 59–70.
  * [32] T. Liggett, Stochastic Interacting Systems: Contact, Voter and Exclusion Processes, Springer, 2013.
  * [33] B. Liu, W. Lu, and T. Chen, Synchronization in complex networks with stochastically switching coupling structures, IEEE Transactions on Automatic Control, 57 (2011), pp. 754–760.
  * [34] I. P. Longo, Topologies of continuity for Carathéodory differential equations with applications in non-autonomous dynamics, Phd Thesis, Universidad de Valladolid, 2018.
  * [35] I. P. Longo, S. Novo, and R. Obaya, Topologies of ${L}^{p}_{loc}$ type for Carathéodory functions with applications in non-autonomous differential equations, Journal of Differential Equations, 263 (2017), pp. 7187–7220.
  * [36] , Topologies of continuity for Carathéodory delay differential equations with applications in non-autonomous dynamics, Discrete and Continuous Dynamical Systems, 39 (2019), pp. 5491–5520.
  * [37] , Weak topologies for Carathéodory differential equations: continuous dependence, exponential dichotomy and attractors, Journal of Dynamics and Differential Equations, 31 (2019), pp. 1617–1651.
  * [38] , Monotone skew-product semiflows for Carathéodory differential equations and applications, Journal of Dynamics and Differential Equations, (2021), pp. 1–37.
  * [39] I. P. Longo, R. Obaya, and A. M. Sanz, Topologies of continuity for Carathéodory parabolic PDEs from a dynamical perspective, Discrete and Continuous Dynamical Systems, 43 (2023), pp. 2213–2240.
  * [40] J. Lu and G. Chen, A time-varying complex dynamical network model and its controlled synchronization criteria, IEEE Transactions on Automatic Control, 50 (2005), pp. 841–846.
  * [41] W. Lu, F. M. Atay, and J. Jost, Synchronization of discrete-time dynamical networks with time-varying couplings, SIAM Journal on Mathematical Analysis, 39 (2008), pp. 1231–1259.
  * [42] R. Mulas, C. Kuehn, and J. Jost, Coupled dynamics on hypergraphs: master stability of steady states and synchronization, Phys. Rev. E, 101 (2020), p. 062313.
  * [43] C. Olech and Z. Opial, Sur une inégalité différentielle, in Annales Polonici Mathematici, vol. 3, 1960, pp. 247–254.
  * [44] H. Othmer and L. Scriven, Non-linear aspects of dynamic pattern in cellular networks, J. Theor. Biol., 43 (1974), pp. 83–112.
  * [45] L. Pecora and T. Carroll, Master stability functions for synchronized coupled systems, Phys. Rev. Lett., 80 (1998), pp. 2109–2112.
  * [46] T. Pereira, Stability of synchronized motion in complex networks, arXiv preprint arXiv:1112.2297, (2011).
  * [47] T. Pereira, J. Eldering, M. Rasmussen, and A. Veneziani, Towards a general theory for coupling functions allowing persistent synchronisation, arXiv preprint arXiv:1304.7679, (2013).
  * [48] M. Porfiri, A master stability function for stochastically coupled chaotic maps, Europhysics Letters, 96 (2011), p. 40014.
  * [49] M. Porfiri, D. J. Stilwell, and E. M. Bollt, Synchronization in random weighted directed networks, IEEE Transactions on Circuits and Systems I: Regular Papers, 55 (2008), pp. 3170–3177.
  * [50] M. Porfiri, D. J. Stilwell, E. M. Bollt, and J. D. Skufca, Random talk: Random walk and synchronizability in a moving neighborhood network, Physica D: Nonlinear Phenomena, 224 (2006), pp. 102–113.
  * [51] M. A. Porter and J. P. Gleeson, Dynamical systems on networks, Frontiers in Applied Dynamical Systems: Reviews and Tutorials, 4 (2016).
  * [52] C. Pötzsche, Exponential dichotomies of linear dynamic equations on measure chains under slowly varying coefficients, Journal of Mathematical Analysis and Applications, 289 (2004), pp. 317–335.
  * [53] S. Rakshit, B. K. Bera, E. M. Bollt, and D. Ghosh, Intralayer synchronization in evolving multiplex hypernetworks: Analytical approach, SIAM Journal on Applied Dynamical Systems, 19 (2020), pp. 918–963.
  * [54] J. Rauch and J. A. Smoller, Qualitative theory of the FitzHugh-Nagumo equations, (1978).
  * [55] C. Rocsoreanu, A. Georgescu, and N. Giurgiteanu, The FitzHugh-Nagumo model: bifurcation and dynamics, vol. 10, Springer Science & Business Media, 2012.
  * [56] R. J. Sacker and G. R. Sell, A spectral theory for linear differential systems, Journal of Differential Equations, 27 (1978), pp. 320–358.
  * [57] J. Saramäki and E. Moro, From seconds to months: an overview of multi-scale dynamics of mobile telephone calls, The European Physical Journal B, 88 (2015), pp. 1–10.
  * [58] L. Segel and S. Levin, Application of nonlinear stability theory to the study of the effects of diffusion on predator-prey interactions, AIP Conf. Proceed., 27 (1976), p. 123.
  * [59] G. Shi, C. Altafini, and J. Baras, Dynamics over signed networks, SIAM Rev., 61 (2019), pp. 229–257.
  * [60] S. Siegmund, Dichotomy spectrum for nonautonomous differential equations, Journal of Dynamics and Differential Equations, 14 (2002), pp. 243–258.
  * [61] D. P. Spanos, R. Olfati-Saber, and R. M. Murray, Dynamic consensus on mobile networks, in IFAC world congress, 2005, pp. 1–6.
  * [62] D. J. Stilwell, E. M. Bollt, and D. G. Roberson, Sufficient conditions for fast switching synchronization in time-varying network topologies, SIAM Journal on Applied Dynamical Systems, 5 (2006), pp. 140–156.
  * [63] A. Stopczynski, V. Sekara, P. Sapiezynski, A. Cuttone, M. M. Madsen, J. E. Larsen, and S. Lehmann, Measuring large-scale social networks with high resolution, PloS one, 9 (2014), p. e95978.
  * [64] S. H. Strogatz, Exploring complex networks, nature, 410 (2001), pp. 268–276.
  * [65] A. Tahbaz-Salehi and A. Jadbabaie, On consensus over random networks, in 44th Annual Allerton Conference, Citeseer, 2006.
  * [66] T. Vicsek and A. Zafiris, Collective motion, Phys. Rep., (2012), pp. 71–140.
  * [67] G. Wen, W. Yu, G. Hu, J. Cao, and X. Yu, Pinning synchronization of directed networks with switching topologies: A multiple lyapunov functions approach, IEEE Transactions on Neural Networks and Learning Systems, 26 (2015), pp. 3239–3250.
  * [68] Y. Zhang and S. H. Strogatz, Designing temporal networks that synchronize under resource constraints, Nature communications, 12 (2021), p. 3273.
  * [69] Y.-Q. Zhang, X. Li, J. Xu, and A. V. Vasilakos, Human interactive patterns in temporal networks, IEEE Transactions on Systems, Man, and Cybernetics: Systems, 45 (2014), pp. 214–222.
  * [70] J. Zhao, D. J. Hill, and T. Liu, Synchronization of dynamical networks with nonidentical nodes: Criteria and control, IEEE Transactions on Circuits and Systems I: Regular Papers, 58 (2010), pp. 584–594.