# Resilience of the slow component in timescale separated synchronized
oscillators

Melvyn Tyloo  Theoretical Division and Center for Nonlinear Studies (CNLS),
Los Alamos National Laboratory, Los Alamos, NM 87545, USA

###### Abstract

Physiological networks are usually made of a large number of biological
oscillators evolving on a multitude of different timescales. Phase oscillators
are particularly useful in the modelling of the synchronization dynamics of
such systems. If the coupling is strong enough compared to the heterogeneity
of the internal parameters, synchronized states might emerge where phase
oscillators start to behave coherently. Here, we focus on the case where
synchronized oscillators are divided into a fast and a slow component so that
the two subsets evolve on separated timescales. We assess the resilience of
the slow component by reducing the dynamics of the fast one and evaluating the
variance of the phase deviations when the oscillators in the two components
are subject to noise with possibly distinct correlation times. From the
general expression for the variance, we consider specific network structures
and show how the noise transmission between the fast and slow components is
affected. Interestingly, we find that oscillators that are among the most
robust when there is only a single timescale, might become the most vulnerable
when the system separates into two timescales. We also find that layered
networks seem to be insensitive to timescale separations when the noise has
homogeneous correlation time.

##  I Introduction

The synchronization dynamics of coupled phase oscillators finds numerous
applications ranging from Josephson junctions and electrical power grids to
physiological networks [1, 2, 3, 4, 5]. The collective behavior displayed by
these system is made possible by the interplay between the internal parameters
of the individual dynamical units and the interaction coupling their degrees
of freedom [6, 7, 8, 9]. Due to the nonlinear nature of the coupling together
with the complex network topology of the interaction, multiple synchronized
states might exist for the same parameters and might be visited by the system
due to perturbations or noise [10, 11, 12, 13]. Importantly, synchronization
is not always a desirable feature. For example, in electrical power grids, a
synchronous operational state ensures the good functioning and distribution of
power [14, 15, 16]. The answer is less binary in physiological systems.
Indeed, synchronization is of primal importance for some cognitive processes
in the brain ensuring an adequate level of communication between neuronal
groups [17]. Also, synchronized dynamics emerge in healthy neuronal systems
during sleep [18]. Therefore, a lack of synchronization might result in some
impairment of physiological system. However, direct connections have been
drawn between the excess of synchronization in some neuronal groups and brain
diseases [19, 20]. Therefore, the synchronization dynamics as well as its
resilience to external perturbations are topics of primal importance in order
to better understand the interplay between synchronized groups of dynamical
units.
Synchronization typically occurs in the phases or frequencies, where in the
latter case all the oscillators are in a phase-locked state with time-
independent phase differences, while in the former case all phases are the
same. Perturbations of these synchronized states can take a great variety of
forms such as external input signals injected into some internal parameters or
noisy environments [12, 21, 22, 23, 24, 25, 26, 27], interruption of the
interaction between some oscillators due to local failures [28, 29],
alteration of the dynamics of some units [30, 31]. Here, we are interested in
networks of phase oscillators in a phase-locked state where, due to damage to
a subset of oscillators or simply because of their intrinsic characteristics,
two separate timescales of the dynamics exist so that the system is divided
into a fast and a slow component. This kind of timescale separation might
occur for example in the human physiological system thanks to the wide range
of timescales reported [32]. In such a scenario, the fast oscillators adapt to
any input signal quickly compared to the ones in the slow component.
Therefore, the input signals into the oscillators belonging to the fast
component are transmitted differently to those in the slow component compared
to the case where all oscillators evolve on the same timescale. Such timescale
separation in systems of coupled phase oscillators have been used in the
modelling of power systems [33] and synchronization dynamics of Kuramoto
oscillators with attractive and repulsive couplings [34]. As a paradigmatic
model to investigate synchronization, we use Kuramoto oscillators, but the
framework presented here applies more generally to coupled dynamical systems
evolving close to a stable fixed point. We consider time-correlated noisy
inputs as in many relevant situations, dynamical systems are constantly pushed
away from their synchronized fixed point by ambient noisy conditions [35]. The
resilience of the system to such perturbations can be assessed in various
ways. One can estimate the size of the basin of attraction [36, 37], or
evaluate the amplitude of small fluctuations or the escape rate of large
fluctuations [11, 12, 26, 22, 23, 38, 39]. In this manuscript, we assess the
resilience of the slow component in the small fluctuation regime by
quantifying the phase deviations from the synchronized state. This is
important as it clarifies which features of the dynamical system make it
robust to noise when coupled oscillators evolve over multiple timescales.
Within the assumption of small fluctuation, we investigate the linear response
of the system around a stable fixed point. We first account for the timescale
separation applying Mori-Zwanzig formalism [40, 41] to the slow and fast
components. This leads to a reduced dynamics of the oscillators which is
equivalent to a Kron reduction of the Jacobian matrix [42, 43]. The latter
elucidate how the inputs in the fast component are transmitted to the slow
one. Then, solving the linear system, we calculate the variance of the phase
deviations in the slow component when time-correlated noises with distinct
typical correlation times are applied in each component. We show how the
amplitude of the excursion essentially depends on the characteristics of the
noise, as well as the system properties through the spectrum of its reduced
Jacobian. In some specific settings, we are able to further predict the
transmission of the noise from the fast to the slow component based on the
properties of the oscillators in the fast component as well as the inter-
component coupling structure. In particular, we find that some oscillators
having smaller variance when there is no timescale separation, might become
the ones with larger variance when there is a timescale separation, if there a
well connected to the fast component. Also, when the slow and fast components
are defined on a layered network, the variance is mostly insensitive to the
timescale separation when the noise correlation time is homogeneous.

In Sec. II, we introduce the model of Kuramoto oscillators with timescale
separation and apply Mori-Zwanzig formalism to obtain a reduced dynamics for
the slow component. In the same section, we then calculate the variance of the
degrees of freedom of the oscillators in the slow component subject to time-
correlated noise. In Sec. III, we numerically confirm and illustrate the
theory on various network structures. The conclusions are given in Sec. IV

##  II Timescale separation

###  II.1 Networks of phase oscillators

We are interested in the situation where, due to an external perturbation or
change in the environment, the intrinsic timescales of the individual
oscillators separate into a fast and a slow subsystem. We consider a set of
$N$ oscillators each with a compact phase degree of freedom
$\theta_{i}\in(-\pi,\pi]$ whose time-evolution is governed by the set of
coupled differential equations [7],

$\displaystyle\begin{split}d_{i}\,\dot{\theta}_{i}=\omega_{i}-\sum_{j=1}^{N}b_{ij}\,\sin(\theta_{i}-\theta_{j})+\eta_{i}\,,\end{split}$
(1)

for $i=1,...\,N$ . The natural frequency of the $i$th oscillator is denoted
$\omega_{i}$ , the structure of the coupling network is given by elements
$b_{ij}$ of the adjacency matrix [44]. Ambient noise is modelled at the $i$th
oscillator by $\eta_{i}$ and is taken as a time-correlated noise, uncorrelated
in space, i.e.
$\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\eta_{0,i}^{2}\delta_{ij}\exp[-|t-t^{\prime}|/\tau_{i}]$
, where $\tau_{i}$ is the typical correlation time of the noise at the $i$th
oscillator. The non-negative parameters $d_{i}$’s define the individual
timescale of each oscillator. Removing the noise term, Eq. (1) may have
multiple stable fixed points of the dynamics which essentially depend on the
coupling topology and strength as well as the distribution of natural
frequencies. Below, we assume that such a stable fixed point
$\\{\theta_{i}^{*}\\}$ exists and that the noise term is small enough such
that the dynamics remains inside the initial basin of attraction. In the
present scenario, we assume that we have two sets of oscillators that we
denote $\mathcal{F}$ and $\mathcal{S}$ , respectively with $N_{\mathcal{F}}$
and $N_{\mathcal{S}}$ oscillators, such that

$\displaystyle d_{i}=\begin{cases}\overline{d}&i\in\mathcal{F}\\\
\underline{d}&i\in\mathcal{S}\end{cases}$ (2)

with $\overline{d}\ll\underline{d}$ . The latter means that oscillators
belonging to $\mathcal{S}$ have a much slower intrinsic timescale than those
belonging to $\mathcal{F}$ . In the following, we focus on the dynamics of the
oscillators in the slow component. Within the assumption of timescale
separation, one can rewrite Eq. (1) as,

$\displaystyle\begin{split}\,\dot{\theta}_{i}&=\omega_{i}-\sum_{j=1}^{N}b_{ij}\,\sin(\theta_{i}-\theta_{j})+\eta_{i}\,,\,i\in\mathcal{S}\\\
\epsilon\,\dot{\theta}_{i}&=\omega_{i}-\sum_{j=1}^{N}b_{ij}\,\sin(\theta_{i}-\theta_{j})+\eta_{i}\,,\,i\in\mathcal{F}\,,\end{split}$
(3)

where we defined $\underline{d}/\overline{d}=\epsilon^{-1}$ and, without loss
of generality, set $\underline{d}$ to unity. In the limit $\epsilon\rightarrow
0$ , the oscillators within $\mathcal{F}$ instantaneously adapt their phases.
In the next section, we consider the dynamical system Eqs. (3) in the vicinity
of a stable fixed point and perform a singular perturbation analysis using
Mori-Zwanzig formalism.

###  II.2 Near-equilibrium and reduced dynamics

Even though we consider Kuramoto oscillators, the following approach applies
in general to coupled dynamical systems that have a stable fixed point around
which they evolve and where linearization is valid. To analyze the resilience
of the slow component, we consider the dynamics of the system close to a fixed
point $\\{\theta_{i}^{*}\\}$ . In particular, we are interested in the time-
evolution of the phase deviations $x_{i}(t)=\theta_{i}(t)-\theta_{i}^{*}$ for
$i\in\mathcal{S}$ and $y_{i}(t)=\theta_{i}(t)-\theta_{i}^{*}$ for
$i\in\mathcal{F}$ whose dynamics at the first order reads,

$\displaystyle\begin{bmatrix}\dot{\bf x}\\\ \epsilon\,\dot{\bf
y}\end{bmatrix}=\begin{bmatrix}{\bf J}_{\mathcal{S}\mathcal{S}}&{\bf
J}_{\mathcal{S}\mathcal{F}}\\\ {\bf J}_{\mathcal{F}\mathcal{S}}&{\bf
J}_{\mathcal{F}\mathcal{F}}\end{bmatrix}\begin{bmatrix}{\bf x}\\\ {\bf
y}\end{bmatrix}+\begin{bmatrix}{\bm{\eta}}_{\mathcal{S}}\\\
{\bm{\eta}}_{\mathcal{F}}\end{bmatrix}$ (4)

where we defined the matrix

$\displaystyle
J_{ij}=\begin{cases}b_{ij}\cos(\theta_{i}^{*}-\theta_{j}^{*})&i\neq j\\\
-\sum_{k=1}^{N}b_{ik}\cos(\theta_{i}^{*}-\theta_{k}^{*})&i=j\,,\end{cases}$
(5)

which is the Jacobian of the system and is a Laplacian matrix when phase
differences are between $-\frac{\pi}{2}$ and $\frac{\pi}{2}$ . Using Mori-
Zwanzig formalism [40, 41] with $\bf x$ and $\bf y$ being respectively the
resolved and unresolved variables (see [45] for an introduction), one can
express the first row of Eq. (4) as,

$\displaystyle\dot{x}_{i}$
$\displaystyle=\sum_{j=1}^{N_{\mathcal{S}}}{J_{\mathcal{S}\mathcal{S}}}_{ij}\,x_{j}+{\eta_{\mathcal{S}}}_{i}+\sum_{\alpha=1}^{N_{\mathcal{F}}}\int_{0}^{t}\epsilon^{-1}e^{\nu_{\alpha}\epsilon^{-1}(t-t^{\prime})}$
(6)
$\displaystyle\times\sum_{k,m=1}^{N_{\mathcal{F}}}\sum_{l=1}^{N_{\mathcal{S}}}{J_{\mathcal{FS}}}_{kl}x_{l}(t^{\prime})w_{\alpha,k}{J_{\mathcal{SF}}}_{im}w_{\alpha,m}\,{\rm
d}t^{\prime}$ $\displaystyle+$
$\displaystyle\sum_{\alpha=1}^{N_{\mathcal{F}}}\int_{0}^{t}\epsilon^{-1}e^{\nu_{\alpha}\epsilon^{-1}(t-t^{\prime})}\sum_{k,m=1}^{N_{\mathcal{F}}}{\eta_{\mathcal{F}}}_{k}(t^{\prime})w_{\alpha,k}{J_{\mathcal{SF}}}_{im}w_{\alpha,m}\,{\rm
d}t^{\prime}\,,$

where we denoted ${\bf w}_{\alpha}$ the eigenvectors of ${\bf
J}_{\mathcal{F}\mathcal{F}}$ with corresponding eigenvalues $\nu_{\alpha}<0$ .
In Eq. (6) , the first two terms on the right-hand side are Markovian while
the other ones have memory. We are interested in the time-evolution of the
slow components ${\bf x}$ when there is a timescale separation, i.e.
$\epsilon\rightarrow 0$ . Taking the latter limit in Eq. (6) yields in a
matrix form,

$\displaystyle\dot{\bf x}$ $\displaystyle={\bf
J}_{\mathcal{S}\mathcal{S}}\,{\bf x}-{\bf J}_{\mathcal{S}\mathcal{F}}\,{\bf
J}^{-1}_{\mathcal{F}\mathcal{F}}\,{\bf J}_{\mathcal{F}\mathcal{S}}\,{\bf
x}+{\bm{\eta}}_{\mathcal{S}}-{\bf J}_{\mathcal{S}\mathcal{F}}\,{\bf
J}^{-1}_{\mathcal{F}\mathcal{F}}\,{\bm{\eta}}_{\mathcal{F}}$ (7)
$\displaystyle={\bf J}_{\rm red}\,{\bf x}+{\bm{\xi}}\,,$

where in the second line we defined the reduced Jacobian ${\bf J}_{\rm
red}={\bf J}_{\mathcal{S}\mathcal{S}}-{\bf J}_{\mathcal{S}\mathcal{F}}\,{\bf
J}^{-1}_{\mathcal{F}\mathcal{F}}\,{\bf J}_{\mathcal{F}\mathcal{S}}$ , and
denoted the noise term as ${\bm{\xi}}={\bm{\eta}}_{\mathcal{S}}-{\bf
J}_{\mathcal{S}\mathcal{F}}\,{\bf
J}^{-1}_{\mathcal{F}\mathcal{F}}\,{\bm{\eta}}_{\mathcal{F}}$ . It is
interesting to note that the reduced dynamics given by Eq. (7) can be obtained
by a Kron reduction [42, 43] of the fast component of the system (see [46] for
an example). The dynamics of the slow component is then governed by Eq. (7)
where the effective noise at the $i$th oscillator is a combination of the
noise at the $i$th oscillator combined with a superposition of the noise
inputs at oscillators belonging to the fast component. For undirected coupling
as we consider in the following, one has that ${\bf
J}_{\mathcal{F}\mathcal{S}}={\bf J}_{\mathcal{S}\mathcal{F}}^{\top}$ . The
linear system Eq. (7) can be solved by expanding over the eigenmodes of ${\bf
J}_{\rm red}$ denoted ${\bf u}_{\alpha}$ , with corresponding eigenvalues
$\lambda_{\alpha}$ , $\alpha=1,...N_{\mathcal{S}}=|\mathcal{S}|$ . As ${\bf
J}_{\rm red}$ is also the negative of a Laplacian matrix, one has that
$0>\lambda_{1}\geq...\geq\lambda_{N_{\mathcal{S}}}$ with
$u_{1,i}=1/\sqrt{N_{\mathcal{S}}}$ . The general solution to Eq. (7) is given
by,

$\displaystyle
x_{i}(t)=\sum_{\alpha=2}^{N_{\mathcal{S}}}\int_{0}^{t}\exp[\lambda_{\alpha}(t-t^{\prime})]{\bf
u}_{\alpha}\cdot{\bm{\xi}}(t^{\prime})\,{\rm d}t^{\prime}u_{\alpha,i}\,,$ (8)

for $i=2,...N_{\mathcal{S}}$ and where we omitted the first mode in the sum as
any perturbation along it does not modify the system . This expression can be
used directly to calculate the moments of the phase deviations.

Figure 1: Connectivity where oscillators in the fast component are
disconnected. The fast oscillators are shown in red while the slow oscillators
are shown in blue.

###  II.3 Fluctuations from the synchronized state

Various characteristics of the response can be used to determine the
resilience of the coupled oscillators. When subject to stochastic inputs, a
natural choice is to evaluate the magnitude of the deviations from the
synchronized fixed point by calculating the variance of the phase deviations.
The variance of the phase deviations in the slow component is calculated from
Eq. (8) and reads in the long time limit,

$\displaystyle\langle x_{i}^{2}\rangle$
$\displaystyle=\eta_{0,\mathcal{S}}^{2}\sum_{\alpha=2}^{N_{\mathcal{S}}}\frac{u_{\alpha,i}^{2}}{\lambda_{\alpha}(\lambda_{\alpha}-\tau_{\mathcal{S}}^{-1})}$
(9) $\displaystyle+\eta_{0,\mathcal{F}}^{2}$
$\displaystyle\sum_{\alpha,\beta=2}^{N_{\mathcal{S}}}\frac{(\lambda_{\alpha}+\lambda_{\beta}-2\tau_{\mathcal{F}}^{-1})\Gamma_{\alpha\beta}\,u_{\alpha,i}u_{\beta,i}}{(\lambda_{\alpha}+\lambda_{\beta})(\tau_{\mathcal{F}}^{-1}-\lambda_{\alpha})(\tau_{\mathcal{F}}^{-1}-\lambda_{\beta})}\,,$

with the scalar $\Gamma_{\alpha\beta}={\bf u}^{\top}_{\alpha}{\bf
J}_{\mathcal{S}\mathcal{F}}{\bf J}_{\mathcal{F}\mathcal{F}}^{-2}{\bf
J}_{\mathcal{F}\mathcal{S}}{\bf u}_{\beta}$ . In Eq. (9), we set the standard
deviation of the ambient noise in the slow and fast components respectively,
to $=\eta_{0,\mathcal{S}}$ and $\eta_{0,\mathcal{F}}$ . We also set distinct
homogeneous correlation times for the noise in each component as
$\tau_{i}=\tau_{\mathcal{S}}$ for $i\in\mathcal{S}$ and
$\tau_{i}=\tau_{\mathcal{F}}$ for $i\in\mathcal{F}$ . While the contribution
to the variance from the additive noise in the slow component is essentially
given by the position of the oscillators on the slowest eigenmodes, the effect
of the noise coming from the fast component involves combinations of
eigenmodes. The precise combination depends on the effective reduced dynamics
through $\Gamma_{\alpha\beta}$ . The shortest timescale in the system is set
by the oscillators belonging to $\mathcal{F}$ . However, by tuning the
correlation time of the noise $\tau=\tau_{\mathcal{F}}=\tau_{\mathcal{S}}$ ,
one can investigate the regimes where the $\lambda_{N_{\mathcal{S}}}\tau\ll 1$
and $\lambda_{2}\tau\gg 1$ . Indeed, in the limit where the noise correlation
time is shorter than the timescales of the slow component and the same in both
components, the variance becomes,

$\displaystyle\langle x_{i}^{2}\rangle$
$\displaystyle=\eta_{0,\mathcal{S}}^{2}\,\tau\sum_{\alpha=2}^{N_{\mathcal{S}}}\frac{u_{\alpha,i}^{2}}{(-\lambda_{\alpha})}+\eta_{0,\mathcal{F}}^{2}\,\tau\sum_{\alpha,\beta=2}^{N_{\mathcal{S}}}\frac{2\,\Gamma_{\alpha\beta}\,u_{\alpha,i}u_{\beta,i}}{(-\lambda_{\alpha}-\lambda_{\beta})}\,.$
(10)

In the other limit where the noise correlation time is the longest timescale,
one has,

$\displaystyle\langle x_{i}^{2}\rangle$
$\displaystyle=\eta_{0,\mathcal{S}}^{2}\,\sum_{\alpha=2}^{N_{\mathcal{S}}}\frac{u_{\alpha,i}^{2}}{\lambda_{\alpha}^{2}}+\eta_{0,\mathcal{F}}^{2}\,\sum_{\alpha,\beta=2}^{N_{\mathcal{S}}}\frac{\Gamma_{\alpha\beta}\,u_{\alpha,i}u_{\beta,i}}{\lambda_{\alpha}\lambda_{\beta}}\,.$
(11)

Comparing the two limiting cases Eqs. (10) and (11) , one remarks that in both
variances, a significant contribution might come from the slowest eigenmodes.
Note also that Eq. (9) is more generally valid in the case where
$\tau_{\mathcal{F}}$ and $\tau_{\mathcal{S}}$ are different.

In the Appendix A , we give the variance of the phases when there is no
timescale separation.

To obtain more insights into the contribution from the fast component, let us
consider specific situations in the strong coupling limit.

###  II.4 Strong coupling limit

In the strong coupling limit, one has $|\theta_{i}^{*}-\theta_{j}^{*}|\ll 1$
$\forall i,j$ , so that one can approximate the Jacobian Eq. (5) as,

$\displaystyle J_{ij}=\begin{cases}b_{ij}&i\neq j\\\
-\sum_{k=1}^{N}b_{ik}&i=j\,.\end{cases}$ (12)

Within this coupling limit and some other assumptions that are specified
below, one can further consider network structures that give more insights
about Eq. (9) .

Figure 2: Connectivity where oscillators in the fast component with inter-
component connections are connected to all the oscillators in the slow
component. The fast oscillators are shown in red while the slow oscillators
are shown in blue. Figure 3: Connectivity where a single oscillator in the
slow component is connected to all the oscillators in the fast one. The fast
oscillators are shown in red while the slow oscillators are shown in blue.

####  II.4.1 Disconnected oscillators in the fast component

In the simple scenario where only a single oscillator $l$ belongs to the fast
component while all the others are in the slow one, ${\bf
J}_{\mathcal{F}\mathcal{F}}$ is a scalar such that
${J}_{\mathcal{F}\mathcal{F}}^{-2}=k_{l}^{-2}$ is the inverse of the squared
weighted degree of the fast oscillator indexed by $l$ . Therefore, one has
$\Gamma_{\alpha\beta}=\left(\sum_{j\in\mathcal{N}(l)}u_{\alpha,j}b_{lj}\right)\left(\sum_{k\in\mathcal{N}(l)}u_{\beta,k}b_{lk}\right)/k_{l}^{2}$
, where $\mathcal{N}(l)$ is the set of oscillators in the slow component
connected to the fast oscillator $l$ . The contribution from the second term
in Eq. (9) therefore crucially depends on the location of the oscillators on
the slowest eigenmodes of the reduced Jacobian. This situation easily
generalizes to the case of multiple oscillators in the fast component that are
not connected as shown in Fig. 1 . One then has
$\Gamma_{\alpha\beta}=\sum_{l\in\mathcal{F}}\left(\sum_{j\in\mathcal{N}(l)}u_{\alpha,j}b_{lj}\right)\left(\sum_{k\in\mathcal{N}(l)}u_{\beta,k}b_{lk}\right)/k_{l}^{2}$
, where we took the sum over all the oscillators in the fast component.

####  II.4.2 All-to-all coupling from fast to slow component

When the oscillators in the fast components that have inter-component
connections are connected to all the oscillators in the slow component with
homogeneous coupling, i.e. $b_{lj}=b_{l}>0$ $\forall
j\in\mathcal{N}(l)=\mathcal{S}$ , with $l\in\mathcal{F}$ (see Fig. 2), the
second term in Eq. (9) vanishes. Indeed, in such a situation, the columns of
the matrix ${\bf J}_{\mathcal{S}\mathcal{F}}={\bf
J}_{\mathcal{F}\mathcal{S}}^{\top}$ corresponding to the inter-component
coupling are full of ones and as ${\bf u}_{1}\cdot{\bf
u}_{\alpha}=\delta_{1\alpha}$ by orthogonality of the eigenmodes, one has that
$\Gamma_{\alpha\beta}=0$ for all $\alpha,\beta=1,...,N_{\mathcal{S}}$ .
Intuitively, if the signal from one oscillator in the fast component is
transmitted with the same strength to all the oscillators in the slow one,
then it will result in an overall phase shift without affecting the fixed
point. However, if the coupling is not homogeneous between the slow and fast
components, the variance will be different from zero.

####  II.4.3 All-to-all from slow to fast component

In the opposite case where only a single oscillator in the slow component is
homogeneously connected to all the oscillators in the fast one (see Fig. 3),
the intra-component structure of the coupling within the fast component does
not influence the propagation of the noise. Indeed, if $b_{lj}=b_{l}>0$
$\forall j\in\mathcal{M}(l)=\mathcal{F}$ where here $\mathcal{M}(l)$ is the
set of oscillators in the fast component connected to the $l$th oscillator in
the slow component, one has
$\Gamma_{\alpha\beta}=b_{l}^{2}u_{\alpha,l}u_{\beta,l}N_{\mathcal{F}}$ .
Therefore, only the size of the fast component and not the intra-component
network structure of the oscillators influence the variance in the slow one.
This can be generalized to the case where multiple oscillators in slow
component are connected to all the ones in the fast component. One then has,
$\Gamma_{\alpha\beta}=\sum_{l\in\mathcal{S}}b_{l}^{2}u_{\alpha,l}u_{\beta,l}m^{-2}$
with $m$ the number of oscillators in the slow component that have all-to-all
connections to the fast component. Here again, the structure of the coupling
within the fast component does not influence the variance.

####  II.4.4 Layered networks

An interesting case arises when oscillators are connected on layered networks
so that the fast component is on one layer, the slow one on another layer, and
the two layers are connected together. In the specific scenario where each
fast oscillator is connected to a single distinct oscillator in the slow
component and the number of units in the layers is the same, one has that
${\bf J}_{\mathcal{S}\mathcal{F}}$ is a diagonal matrix (up to a permutation
of the oscillators indices). This is illustrated in Fig. 4 . If one further
assumes that the inter-layer coupling is homogeneous, i.e. ${\bf
J}_{\mathcal{S}\mathcal{F}}=\tilde{b}\,\mathbb{I}$ is a multiple of the
identity matrix, one has that the eigenbases of ${\bf J}_{\rm red}$ and ${\bf
J}_{\mathcal{S}\mathcal{S}}={\bf J}_{\mathcal{F}\mathcal{F}}$ satisfy the
following relations,

$\displaystyle{\bf J}_{\rm red}{\bf u}_{\alpha}=\lambda_{\alpha}{\bf
u}_{\alpha}$ $\displaystyle=({\bf
J}_{\mathcal{S}\mathcal{S}}-\tilde{b}^{2}{{\bf
J}_{\mathcal{S}\mathcal{S}}}^{-1}){\bf u}_{\alpha}\,,$ (13)
$\displaystyle({\bf J}_{\mathcal{S}\mathcal{S}}-\tilde{b}^{2}{{\bf
J}_{\mathcal{S}\mathcal{S}}}^{-1}){\bf v}_{\beta}$
$\displaystyle=(\mu_{\beta}-\tilde{b}^{2}\mu_{\beta}^{-1}){\bf v}_{\beta}={\bf
J}_{\rm red}{\bf v}_{\beta}\,,$ (14)

with $\alpha,\beta=1,...,N_{\mathcal{S}}$ and where the eigenmodes of ${\bf
J}_{\mathcal{S}\mathcal{S}}$ are denoted ${\bf v}_{\beta}$ with corresponding
eigenvalue $\mu_{\beta}$ . Noticing that
$(\mu_{\beta}-\tilde{b}^{2}\,\mu_{\beta}^{-1})$ is a monotonically increasing
function of $\mu_{\beta}$ , one has a one-to-one correspondence between the
eigenmodes of ${\bf J}_{\rm red}$ and ${\bf J}_{\mathcal{S}\mathcal{S}}$ such
that ${\bf u}_{\alpha}={\bf v}_{\alpha}$ and
$\lambda_{\alpha}=(\mu_{\alpha}-\tilde{b}^{2}\,\mu_{\alpha}^{-1})$ for
$\alpha=1,...,N_{\mathcal{S}}$ . Given the specific structure of the coupling,
one further has that $\mu_{\beta}=\gamma_{\beta}-\tilde{b}$ where
$\gamma_{\beta}$ are the eigenvalues of the Jacobian in each of the layer when
removing the inter-layer coupling. Assuming that the noise amplitudes as well
as the correlation times are the same in both components, one can rewrite Eq.
(9) as,

$\displaystyle\langle x_{i}^{2}\rangle$
$\displaystyle=\eta_{0}^{2}\sum_{\alpha=2}^{N_{\mathcal{S}}}\frac{v_{\alpha,i}^{2}(\mu_{\alpha}^{2}+\tilde{b}^{2})}{(\mu_{\alpha}^{2}-\tilde{b}^{2})(\mu_{\alpha}^{2}-\tilde{b}^{2}-\mu_{\alpha}\tau^{-1})}\,.$
(15)

One can also calculate the variance when there is no timescale separation by
remarking that the eigenmodes of the full Jacobian Eq. (5) are given by $[{\bf
v}_{\alpha},\pm{\bf v}_{\alpha}]^{\top}/\sqrt{2}$ with corresponding
eigenvalues $(\mu_{\alpha}\pm\tilde{b})$ for
$\alpha=1,...,(N_{\mathcal{S}}+N_{\mathcal{F}})$ . Using Eq. (17) , the
variance then reads,

$\displaystyle\langle x_{i}^{2}\rangle$
$\displaystyle=\eta_{0}^{2}\sum_{\alpha=2}^{N_{\mathcal{S}}}\frac{v_{\alpha,i}^{2}(-\mu_{\alpha}\tau^{-1}+\mu_{\alpha}^{2}+\tilde{b}^{2})}{(\mu_{\alpha}^{2}-\tilde{b}^{2})(\tau^{-2}-2\mu_{\alpha}\tau^{-1}-\tilde{b}^{2}+\mu_{\alpha}^{2})}$
(16)
$\displaystyle+\frac{\eta_{0}^{2}(N_{\mathcal{S}}+N_{\mathcal{F}})^{-1}}{2\tilde{b}(\tau^{-1}+2\tilde{b})}\,.$

While Eqs. (15), (16) are different in general, when the correlation time of
the noise becomes short, the two variances only differ by a constant given by
the second term in the right-hand side of Eq. (16) .

Figure 4: Connectivity in layers where the fast oscillators are shown in red
in the top layer while the slow oscillators are shown in blue in the bottom
layer. Each oscillator in one of the layer can be connected only to a single
oscillator in the other layer.

In the following section, we illustrate and confirm numerically the results
discussed so far.

Figure 5: Sum of the variance of the phase deviations at every oscillator in
the slow component as a function of the coupling strength $K$ . The variance
is obtained by time-evolving Eq. (1) for a Watts-Strogatz network of 50
oscillators with heterogeneous natural frequencies (${\rm var}[\omega]=0.01$)
where the ratio of the damping parameters in the slow and fast component is
$\overline{d}/\underline{d}=0.1$ . The correlation times of the noise are
fixed to
$\tau_{\mathcal{S}}/\underline{d}=\tau_{\mathcal{F}}/\underline{d}=0.05$ .
Each dot is obtained by time-averaging the variance over a single simulation
of the dynamics while the solid curve gives the theory that is the sum of the
variances in Eq. (9) . The noise amplitudes are
$\eta_{0,\mathcal{S}}/\underline{d}=\eta_{0,\mathcal{F}}/\underline{d}=0.05$ .
The black dashed and dotted-dashed lines give respectively, $\sim K^{-1}$ and
$\sim K^{-2}$ .

##  III Numerical Results

We first confirm numerically Eq. (9) in Fig. 5 where we consider a Watts-
Strogatz network [44] with $m=4$ initial nearest neighbors and a rewiring
probability $p_{\rm rewiring}=0.1$ , of size
$N_{\mathcal{S}}+N_{\mathcal{F}}=50$ , with respectively $N_{\mathcal{S}}=35$
and $N_{\mathcal{F}}=15$ oscillators in each component. The damping parameters
are taken homogeneous in each component and satisfy
$\overline{d}/\underline{d}=0.1$ in order to effectively simulate a timescale
separated system . We fix the correlation times of the noise
$\tau_{\mathcal{S}}/\underline{d}=\tau_{\mathcal{F}}/\underline{d}=0.05$ and
tune the coupling weights of the network as $b_{ij}\rightarrow K\,b_{ij}$
$\forall i,j=1,...,50$ in order to visit the different regimes. As predicted
by the theory Eq. (9) and its two limiting cases Eqs. (10), (11) , when the
noise correlation time is shorter than the system’s timescales, one has that
$\sum_{i\in\mathcal{S}}\langle x_{i}^{2}\rangle\sim 1/K$ while in the case
where the noise is much slower than the system’s timescales one has
$\sum_{i\in\mathcal{S}}\langle x_{i}^{2}\rangle\sim 1/K^{2}$ . Moreover, the
theory given by the solid line matches the numerical results given by the
dots.

Figure 6: Variance of the phase deviations at every oscillator in the slow
component when noise is only applied on the fast component. The variance is
obtained by time-evolving Eq. (1) for a modified Erdős-Rényi network of 13
oscillators where 3 oscillators in the fast component are connected to all the
oscillators in the slow component. The ratio of the damping parameters in the
slow and fast component is $\overline{d}/\underline{d}=0.1$ . The correlation
time of the noise in the fast component is fixed to
$\tau_{\mathcal{F}}/\underline{d}=0.05$ for the orange and green data points
and $\tau_{\mathcal{F}}/\underline{d}=10$ for the red and blue ones. The
natural frequencies are such that ${\rm var}[\omega]=0.01$ for the green and
blue points and ${\rm var}[\omega]=0.05$ for the red and orange ones. Each dot
is obtained by time-averaging the variance over a single simulation of the
dynamics while the solid curve gives the theory Eq. (9) . The noise amplitude
is $\eta_{0,\mathcal{F}}/\underline{d}=0.1$ . Figure 7: Variance of the phase
deviations at every oscillator when there no timescale separation [panels (a),
(c) and orange and green data points in panel (e)] and when there is a
timescale separation [panels (b), (d) and red and blue data points in panel
(e)], for a modified Watts-Strogatz network of
$N_{\mathcal{S}}+N_{\mathcal{F}}=30$ where the yellow node (with largest
variance) in (b) is connected to all the oscillators in the fast component.
The oscillators depicted by grey squares are in the fast component while all
the others are in the slow component. The variance is numerically obtained by
time-averaging over a single realization of the dynamics Eq. (1) . (e)
Comparison between the theory Eq. (9) given by the solid lines and the
numerical simulations given by the dots when there is no timescale separation
(in orange and green) and when there is a timescale separation (in red and
blue). The ratio of the damping parameters in the slow and fast component is
$\overline{d}/\underline{d}=0.02$ . The correlation times of the noise are
fixed to
$\tau_{\mathcal{F}}/\underline{d}=\tau_{\mathcal{F}}/\underline{d}=0.05$
(orange and red) and
$\tau_{\mathcal{F}}/\underline{d}=\tau_{\mathcal{F}}/\underline{d}=50$ (green
and blue). The natural frequencies are such that ${\rm var}[\omega]=0.03$ .
The noise amplitude is $\eta_{0,\mathcal{F}}/\underline{d}=0.2$ . Figure 8:
Variance of the phase deviations at every oscillator when there no timescale
separation [panels (a), (c) and orange and green data points in panel (e)] and
when there is a timescale separation [panels (b), (d) and red and blue data
points in panel (e)], for a modified Watts-Strogatz network of
$N_{\mathcal{S}}+N_{\mathcal{F}}=30$ where the yellow node (with largest
variance) in (b) is connected to all the oscillators in the fast component.
The oscillators depicted by grey squares are in the fast component while all
the others are in the slow component. The variance is numerically obtained by
time-averaging over a single realization of the dynamics Eq. (1) . (e)
Comparison between the theory Eq. (9) given by the solid lines and the
numerical simulations given by the dots when there is no timescale separation
(in orange and green) and when there is a timescale separation (in red and
blue). The ratio of the damping parameters in the slow and fast component is
$\overline{d}/\underline{d}=0.02$ . The correlation times of the noise are
fixed to $\tau_{\mathcal{F}}/\underline{d}=0.05$
$\tau_{\mathcal{F}}/\underline{d}=50$ (orange and red) and
$\tau_{\mathcal{F}}/\underline{d}=50$ ,
$\tau_{\mathcal{F}}/\underline{d}=0.05$ (green and blue). The natural
frequencies are such that ${\rm var}[\omega]=0.03$ . The noise amplitude is
$\eta_{0,\mathcal{F}}/\underline{d}=0.2$ . Figure 9: Variance of the phase
deviations at every oscillator in the slow component when noise is only
applied on both the slow and the fast components. The variance is obtained by
time-evolving Eq. (1) for a layered network made of two copies of an Erdős-
Rényi network of 15 oscillators. One layer is the fast component while the
other is the slow one. Pairs of corresponding oscillators in the two layers
have a single inter-component connection. The correlation time of the noise in
both components is fixed to $\tau/\underline{d}=0.05$ for the red and orange
data points, $\tau/\underline{d}=50$ for the green and blue ones, and
$\tau/\underline{d}=2$ for the purple and pink ones. The natural frequencies
are identical and vanishing for all oscillators. Each dot is obtained by time-
averaging the variance over a single simulation of the dynamics while the
solid curve gives the theory Eq. (9) when there is no timescale separation
($\overline{d}/\underline{d}=1$) for the orange, green and purple data points,
and when there is a timescale separation ($\overline{d}/\underline{d}=0.05$)
for the red, blue and pink ones. The noise amplitude is
$\eta_{0,\mathcal{F}}/\underline{d}=0.1$ .

Next, in Sec. II.4 , we discuss particular examples in the strong coupling
limit. We numerically investigate these examples also going beyond the strong
coupling limit.

###  III.1 All-to-all coupling from fast to slow component

We consider the situation described in Sec. II.4.2 where oscillators in the
fast component that have some connection to the slow one are connected to all
of them. The numerical results are shown on Fig. 6 where this particular
setting has been simulated for a modified Erdős-Rényi network of
$N_{\mathcal{S}}+N_{\mathcal{F}}=13$ oscillators with various correlation time
of the noise and different heterogeneity levels in the natural frequencies.
Here, only the fast component is subjected to noise, while the slow component
is noiseless. One observes that increasing the heterogeneity in the natural
frequencies of the oscillators induces larger variances for the phase
deviations. Not shown on Fig. 6 is the homogeneous case of oscillators with
identical natural frequencies, for which the variance vanishes as predicted in
Sec. II.4.2 . While the heterogeneity increases the noise transmission from
the fast to the slow component, one observes that the amplitude of the
deviations is still rather small in Fig. 5 . Besides ${\bf
u}^{\top}_{\alpha}{\bf J}_{\mathcal{S}\mathcal{F}}$ being small, this is
because the oscillators in the fast component with inter-component connections
have relatively large degrees which directly reduce the noise transmission in
$\Gamma_{\alpha\beta}={\bf u}^{\top}_{\alpha}{\bf
J}_{\mathcal{S}\mathcal{F}}{\bf J}_{\mathcal{F}\mathcal{F}}^{-2}{\bf
J}_{\mathcal{F}\mathcal{S}}{\bf u}_{\beta}$ .

###  III.2 All-to-all from slow to fast component

In the other situation where some oscillators in the slow components are
connected to a large fraction of the oscillators in the fast component, we
showed in Sec. II.4.3 that their variance is more important than oscillators
with fewer or no connection to the fast component. This result is particularly
interesting and intriguing, as in the regular situation where there is no
timescale separation, oscillators that have a larger number of connections to
the other elements typically have a smaller variance [47]. Indeed, this is
first illustrated in Fig. 7 where the variance of each oscillator is given by
the color map when there is no timescale separation in panels (a) and (c), and
when there is a timescale separation in panels (b) and (d) within the
synchronized oscillators. In both correlation time limits (with
$\tau_{\mathcal{S}}=\tau_{\mathcal{F}}$), one observes that oscillators
belonging to $\mathcal{S}$ having larger variances in Fig. 7(b) and (d) are
among the ones with smaller variances in panels (a) and (c). Comparing Fig.
7(c) and (d), one also observes that the timescale separation modifies the
variance of the oscillators far from the fast component. The theory and the
numerical simulations for the two systems with and without timescale
separation are confirmed in Fig. 7(e) . We then move to the case where the
noise correlation times are different in each component. In Fig. 8 , we show
both limits $\lambda_{N_{\mathcal{S}}}\tau_{\mathcal{F}}\ll 1$ and
$\lambda_{2}\tau_{\mathcal{S}}\gg 1$ [panels (a), (b) and orange and red data
points in panel (e)], and $\lambda_{N_{\mathcal{S}}}\tau_{\mathcal{F}}\gg 1$
and $\lambda_{2}\tau_{\mathcal{S}}\ll 1$ [panels (c), (d) and green and blue
data points in panel (e)]. Similar amplification of the fluctuations as in the
previous case are observed. However, comparing Fig. 8(c) and (d), one remarks
that some oscillators well connected to the fast component keep rather small
variances while others having fewer connections become more vulnerable.

###  III.3 Layered networks

Here, we check the theory when the system is defined on a layered network,
i.e. the slow and fast components each corresponds to one layer. We consider
the specific setting where both layers have the same network connectivity and
the inter-layer coupling is made through single connections between
corresponding oscillators in each component (such structure are sometimes
called multiplex [44]). In Fig. 9 , the numerical simulations for the variance
(dots) match the theory Eqs. (9) and (17) (solid lines) for various
correlation times of the noise [homogeneous in panel (a), i.e.
$\tau_{\mathcal{F}}=\tau_{\mathcal{S}}$ and inhomogeneous in panel (b)]. The
red, blue and pink data points correspond to the situation where there is a
timescale separation between the two components, while for the orange, green
and purple ones, there is no timescale separation. Interestingly, one observes
that the two different situations produce similar variances for the phases. As
predicted in Sec. II.4.4 , in the limit where $\tau$ is the longest timescale
in the system, the variances in the two situations only differ by a constant,
which is small when comparing the blue and green data points in Fig. 9(a) .
However, in Fig. 9(b) where the noise correlation times are distinct in each
component, one sees that the variances with and without timescale separation
differ.

##  IV Conclusion

Physiological systems are composed of a multitude of synchronized dynamical
units evolving on various timescales. It is therefore relevant to investigate
how these different timescales impact the synchronization dynamics of
networked phase oscillators. Here, we considered networks of synchronized
phase oscillators where a timescale separation divides the units into a slow
and a fast component. Using Mori-Zwanzig formalism, we derived a reduced
dynamical system describing the time-evolution of the slow component which we
used to assess the resilience of the slow component by calculating the
variance of the phase deviations. We obtained a closed-form expression for the
variance of each oscillator as a function of the eigenmodes of the reduced
Jacobian. Interestingly, noise propagation from the fast to the slow component
essentially depends on the mixing of the different eigenmodes. The precise
mixing is given by the inter- and intra-component coupling structures. In
particular, we showed that oscillators that have a small variance when there
is no timescale separation, might have a strongly amplified variance when
there is a timescale separation and they have numerous connections to the fast
component. Also, for homogeneous long correlation time of the noise, we found
that when the fast and slow component are connected over a layered structure,
the variance of the oscillators is mostly insensitive to a timescale
separation.
The theory presented here highlights the importance of timescales to assess
the resilience of coupled phase oscillators. Some oscillators that might be
the most robust within one ratio of the timescales, might become the most
fragile ones for another ratio (see Figs. 7 and 8).
While the results of this manuscript were obtained for Kuramoto oscillators,
they apply more generally to coupled dynamical system evolving close to a
stable fixed so that the linear approximation is valid.
Future research should consider more than one timescale separation, evaluate
the resilience of the fast component, and go beyond the small fluctuation
framework presented here.

## Acknowledgments

We thank Francesco Caravelli for useful discussions. This work has been
supported by the Laboratory Directed Research and Development program of Los
Alamos National Laboratory under project numbers 20220797PRD2 and 20220774ER
and by U.S. DOE/OE as part of the DOE Advanced Sensor and Data Analytics
Program.

##  Appendix A No timescale separation

One can also calculate the variance of the oscillators belonging to
$\mathcal{S}$ (and also $\mathcal{F}$) when there is no timescale separation.
Assuming as previously that the oscillators in $\mathcal{S}$ and $\mathcal{F}$
are subject to noises with correlation times $\tau_{\mathcal{S}}$ and
$\tau_{\mathcal{F}}$ respectively, and that the standard deviations are
homogeneous, one has,

$\displaystyle\langle x_{i}^{2}\rangle$
$\displaystyle=\eta_{0}^{2}\sum_{\alpha,\beta=2}^{N_{\mathcal{S}}+N_{\mathcal{F}}}\sum_{j\in\mathcal{S}}\frac{(\gamma_{\alpha}+\gamma_{\beta}+-2\tau_{\mathcal{S}}^{-1})q_{\alpha,j}q_{\beta,j}q_{\alpha,i}q_{\beta,i}}{(\gamma_{\alpha}+\gamma_{\beta})(\tau_{\mathcal{S}}^{-1}-\gamma_{\alpha})(\tau_{\mathcal{S}}^{-1}-\gamma_{\beta})}$
(17)
$\displaystyle+\eta_{0}^{2}\sum_{\alpha,\beta=2}^{N_{\mathcal{S}}+N_{\mathcal{F}}}\sum_{j\in\mathcal{F}}\frac{(\gamma_{\alpha}+\gamma_{\beta}-2\tau_{\mathcal{F}}^{-1})q_{\alpha,j}q_{\beta,j}q_{\alpha,i}q_{\beta,i}}{(\gamma_{\alpha}+\gamma_{\beta})(\tau_{\mathcal{F}}^{-1}-\gamma_{\alpha})(\tau_{\mathcal{F}}^{-1}-\gamma_{\beta})}\,,$

where ${\bf q}_{\alpha}$ are the eigenvectors of the Jacobian Eq. (5) , with
corresponding eigenvalues $\gamma_{\alpha}$ .

## References

  * Wiesenfeld  _et al._ [1998] K. Wiesenfeld, P. Colet, and S. H. Strogatz, Phys. Rev. E 57, 1563 (1998).
  * Pikovsky  _et al._ [2003] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_ (Cambridge university press, 2003).
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
  * Strogatz [2014] S. H. Strogatz, Westview Press, Boulder (2014).
  * Stiefel and Ermentrout [2016] K. M. Stiefel and G. B. Ermentrout, Journal of neurophysiology 116, 2950 (2016).
  * Winfree [1967] A. T. Winfree, Journal of theoretical biology 16, 15 (1967).
  * Kuramoto [1975] Y. Kuramoto, in  _Lecture Notes in Physics 39, International Symposium on Mathematical Problems in Theoretical Physics_, edited by H. Araki (Springer, Berlin, 1975).
  * Kuramoto [1984] Y. Kuramoto,  _Chemical Oscillations, Waves and Turbulence_ , Springer Series in Synergetics, Vol. 19 (Springer Berlin Heidelberg, 1984).
  * Strogatz [2000] S. H. Strogatz, Physica D 143, 1 (2000).
  * Kramers [1940] H. Kramers, Physica 7, 284 (1940).
  * Dykman [1990] M. I. Dykman, Phys. Rev. A 42, 2020 (1990).
  * DeVille [2012] L. DeVille, Nonlinearity 25, 1473 (2012).
  * Rodrigues  _et al._ [2016] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, Physics Reports 610, 1 (2016).
  * Machowski  _et al._ [2008] J. Machowski, J. W. Bialek, and J. R. Bumby,  _Power System Dynamics_ , 2nd ed. (Wiley, Chichester, U.K, 2008).
  * Blaabjerg  _et al._ [2006] F. Blaabjerg, R. Teodorescu, M. Liserre, and A. V. Timbus, IEEE Transactions on industrial electronics 53, 1398 (2006).
  * Dörfler  _et al._ [2013] F. Dörfler, M. Chertkov, and F. Bullo, Proc. Natl. Acad. Sci. 110, 2005 (2013).
  * Fries [2005] P. Fries, Trends in cognitive sciences 9, 474 (2005).
  * González  _et al._ [2023] J. González, M. Cavelli, A. B. Tort, P. Torterolo, and N. Rubido, Plos one 18, e0290146 (2023).
  * Uhlhaas and Singer [2006] P. J. Uhlhaas and W. Singer, neuron 52, 155 (2006).
  * Popovych and Tass [2014] O. V. Popovych and P. A. Tass, Frontiers in neurology 5, 268 (2014).
  * Hindes and Myers [2015] J. Hindes and C. R. Myers, Chaos 25 (2015).
  * Schäfer  _et al._ [2017] B. Schäfer, M. Matthiae, X. Zhang, M. Rohden, M. Timme, and D. Witthaut, Phys. Rev. E 95, 060203(R) (2017).
  * Hindes and Schwartz [2018] J. Hindes and I. B. Schwartz, Chaos 28, 071106 (2018).
  * Ronellenfitsch  _et al._ [2018] H. Ronellenfitsch, J. Dunkel, and M. Wilczek, Physical Review Letters 121, 208301 (2018).
  * Tyloo  _et al._ [2018] M. Tyloo, T. Coletta, and P. Jacquod, Phys. Rev. Lett. 120, 084101 (2018).
  * Halekotte and Feudel [2020] L. Halekotte and U. Feudel, Scientific reports 10, 11783 (2020).
  * Tyloo [2022a] M. Tyloo, Journal of Physics: Complexity 3, 03LT01 (2022a).
  * Soltan  _et al._ [2017] S. Soltan, D. Mazauric, and G. Zussman, IEEE Transactions on Control of Network Systems 4, 288 (2017).
  * Delabays  _et al._ [2022] R. Delabays, L. Pagnier, and M. Tyloo, IFAC-PapersOnLine 55, 270 (2022), 9th IFAC Conference on Networked Systems NECSYS 2022.
  * Wang and Wang [2019] Z. Wang and Y. Wang, IEEE Transactions on Control of Network Systems 6, 338 (2019).
  * Tyloo [2023] M. Tyloo, Journal of Physics: Complexity 4, 045005 (2023).
  * Gao  _et al._ [2020] R. Gao, R. L. van den Brink, T. Pfeffer, and B. Voytek, Elife 9, e61277 (2020).
  * Kokotovic  _et al._ [1980] P. V. Kokotovic, J. J. Allemong, J. R. Winkelman, and J. H. Chow, Automatica 16, 23 (1980).
  * Kirillov  _et al._ [2020] S. Y. Kirillov, V. Klinshov, and V. I. Nekorkin, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2020).
  * van Kampen [1976] N. G. van Kampen, Phys. Rep. 24, 171 (1976).
  * Wiley  _et al._ [2006] D. A. Wiley, S. H. Strogatz, and M. Girvan, Chaos 16, 015103 (2006).
  * Menck  _et al._ [2013] P. J. Menck, J. Heitzig, N. Marwan, and J. Kurths, Nat. Phys. 9, 89 (2013).
  * Tyloo [2022b] M. Tyloo, Chaos 32 (2022b).
  * Hindes  _et al._ [2023] J. Hindes, I. B. Schwartz, and M. Tyloo, Chaos 33 (2023).
  * Mori [1965] H. Mori, Progress of theoretical physics 33, 423 (1965).
  * Zwanzig [1973] R. Zwanzig, Journal of Statistical Physics 9, 215 (1973).
  * Kron [1939] G. Kron,  _Tensor Analysis of Networks_ (Wiley, Cambridge, England, 1939).
  * Dörfler and Bullo [2012] F. Dörfler and F. Bullo, IEEE Transactions on Circuits and Systems I: Regular Papers 60, 150 (2012).
  * Newman [2018] M. Newman,  _Networks_ (Oxford University Press, 2018).
  * Caravelli and Lin [2023] F. Caravelli and Y. T. Lin, arXiv preprint arXiv:2308.13653 (2023).
  * Tyloo  _et al._ [2023] M. Tyloo, M. Vuffray, and A. Y. Lokhov, arXiv preprint arXiv:2310.00458 (2023).
  * Tyloo  _et al._ [2019] M. Tyloo, L. Pagnier, and P. Jacquod, Sci. Adv. 5, eaaw8359 (2019).