# Cluster Synchronization of Kuramoto Oscillators and Brain Functional
Connectivity

Rui Kato and Hideaki Ishii, Fellow, IEEE This work was supported in the part
by JSPS under Grants-in-Aid for Scientific Research (B) Grant No. JP18H01460
and for JSPS Research Fellow Grant No. JP21J21262, and by JST-Mirai program
Grant No. 18077648.R. Kato and H. Ishii are with the Department of Computer
Science, Tokyo Institute of Technology, Yokohama, 226-8502, Japan. (e-mail:
<EMAIL_ADDRESS>ishii@c.titech.ac.jp).

###### Abstract

The recent progress of functional magnetic resonance imaging techniques has
unveiled that human brains exhibit clustered correlation patterns of their
spontaneous activities. It is important to understand the mechanism of cluster
synchronization phenomena since it may reflect the underlying brain functions
and brain diseases. In this paper, we investigate cluster synchronization
conditions for networks of Kuramoto oscillators. The key analytical tool that
we use is the method of averaging, and we provide a unified framework of
stability analysis for cluster synchronization. The main results show that
cluster synchronization is achieved if (i) the inter-cluster coupling
strengths are sufficiently weak and/or (ii) the natural frequencies are
largely different among clusters. Moreover, we apply our theoretical findings
to empirical brain networks. Discussions on how to understand brain functional
connectivity and further directions to investigate neuroscientific questions
are provided.

###### Index Terms:

Brain networks, cluster synchronization, Kuramoto model, functional
connectivity, methods of averaging.

##  I Introduction

Synchronization in complex networks has attracted interests over the past
several decades (see the surveys [1, 2]). For example, synchronous generation
in power grids [3] and cooperative control of networked robots [4] are
representative applications in technological fields. On the other hand,
applications in biological fields include physiological rhythms [5] and
neuronal oscillations [6]. The phenomenon where all states are consistent
among the whole network is called _full synchronization_ or _complete
synchronization_. In contrast, it is known that human brains generate neural
activities with wide-range frequency bands, and their characteristic
statistical patterns can be observed [7]. One of the important properties
found there is that brains exhibit _cluster synchronization_ , that is, some
synchronized clusters coexist in the network.

The recent progress in brain imaging techniques allows us to noninvasively
measure the anatomical network structure in a whole brain [8, 9]. There, the
connections between selected cortical regions by nervous fibers are extracted
by, e.g., diffusion tensor imaging. On the other hand, it is known that brains
exhibit spontaneous activities even at resting state (without external
stimuli). These activities can be quantified by functional magnetic resonance
imaging (fMRI), and a characteristic correlation pattern of fMRI signals has
been found [10]. This so-called resting-state functional connectivity is
considered to serve an important role in brain functions [11]. Recently, in
the neuroscience community, functional connectivity together with anatomical
structures has been studied through graph-theoretic tools [12].

Moreover, functional connectivity has also been examined for patients subject
to brain diseases, where statistical differences between fMRI signals in
healthy brains and those in pathological brains were found. For example, the
small-world property observed in healthy brains can be disturbed by
Alzheimer’s diseases and schizophrenia [13, 14]. In epilepsy, abnormal
functional connectivities are related to the significant increase in
synchronous neural activities [15]. Hence, theoretical analysis of cluster
synchronization may be helpful for understanding the mechanism of healthy and
pathological brains [16]. On the other hand, it is increasingly important to
change the pathological state of the brain into the healthy one by external
interventions [17]. Today, electrical stimulation techniques are practical
methods to manipulate neural activities [18]. Optical stimulation methods are
also gaining attentions as a new tool to control brain states [19].
Furthermore, there is a need for developing theoretical foundations to
understand and control brain dynamics [20].

This paper considers cluster synchronization of Kuramoto oscillators with
heterogeneous natural frequencies. Here, cluster synchronization means that
the oscillators may be divided into some groups, and their phases are
eventually synchronized within the groups at possibly different frequencies.
The Kuramoto model is known as a representative mathematical model of coupled
oscillators to study collective synchronization phenomena [21]. Stability
analysis for full synchronization of Kuramoto oscillators with general network
topologies was investigated in the seminal paper [22]. After that work, many
related results have been established with several settings (see [2] and the
references therein). However, cluster synchronization of Kuramoto oscillators
has not been deeply studied except for [23, 24]. A slightly different
definition of cluster synchronization was employed in [25], where phase
cohesiveness within clusters was explored. On the other hand, some authors
considered cluster synchronization problems for general nonlinear systems. To
study cluster synchronization, graph symmetry was employed in [26, 27, 28, 29,
30, 31], and the notion known as _external equitable partitions_ was addressed
in [32, 33]. The latter works are more general than the former ones in the
sense that any symmetric network allows an external equitable partition. In
[34, 35] sufficient conditions for exponential convergence to a cluster
synchronization manifold were derived based on contraction analysis. The
effects of delays and antagonistic interactions in cluster synchronization
were explored in [36]. Moreover, control methodologies for cluster
synchronization were studied in [37, 38, 39].

Differently from the aforementioned works, our analysis is based on regarding
the intra- and inter-cluster dynamics as slow and fast subsystems,
respectively. Thus, we start with formulating the intra- and inter-cluster
dynamics in a singular perturbation form, where the smallest difference in
natural frequencies between clusters is interpreted as a perturbation
parameter. Since we are concerned with partial stability of the intra-cluster
dynamics, we then linearize the slow subsystem while the solution of the fast
subsystem is considered as a function of time. To analyze stability of the
obtained linear time-varying system, we utilize the method of averaging.
However, existing results such as one in [40] cannot be applied to our system
(see Remark 2). To overcome this difficulty, we employ nonmonotonic Lyapunov
functions to extend existing methods. We show that our framework is useful in
the study of cluster synchronization and that the time-scale property may be a
key factor in the statistical patterns of brain activities.

The contributions of this paper are summarized as follows. As a technical
contribution, we provide a stability criterion for a class of linear time-
varying systems through the method of averaging. Our result generalizes some
existing works and plays a crucial role in the subsequent analysis. The main
contribution of this paper is to provide a unified framework for stability
analysis of cluster synchronization in networks of Kuramoto oscillators.
Specifically, we show that cluster synchronization is achieved if (i) the
inter-cluster coupling strengths are sufficiently weak and/or (ii) the natural
frequencies are largely different among clusters. Compared with [24], where
these two properties are analyzed independently, we provide a tradeoff
relation between the inter-cluster connectivity and heterogeneity.
Intuitively, to achieve synchronization with clustering, large inter-cluster
heterogeneity in the frequencies is required when there are strong couplings
between clusters. Also, weak inter-cluster connectivity is required when the
oscillators have similar frequencies among clusters. The above results are
derived under the assumption that the network is partitioned into clusters
such that an external equitable partition is fulfilled. For the case without
this assumption, we explore practical stability and show that phase
cohesiveness within clusters can be made arbitrarily small by strengthening
the intra-cluster couplings.

Furthermore, we apply our theoretical findings to a brain network constructed
using empirical data. By numerical simulation, we observe that our results on
cluster synchronization are valid for the brain network. Moreover, it is shown
that cluster synchronization of Kuramoto oscillators is strongly related to
functional connectivity, which is identified by a statistical correlation
pattern of region-dependent fMRI signals in human brains. Compared with the
graph-theoretic analysis in the neuroscience community, our work is based on
dynamic properties of brains and provides a theoretical foundation to
understand functional connectivity.

The organization of this paper is described as follows. We start with
transforming the Kuramoto model into a singular perturbation form in Section
II. Based on the method of averaging, we explore a stability criterion for
linear time-varying systems in Section III, Then, we derive stability
conditions for cluster synchronization of Kuramoto oscillators in Section IV.
The case without external equitable partitions is investigated in Section V.
We consider applications to brain networks with empirical data in Section VI,
and finally, Section VII provides concluding remarks of this paper. A
preliminary version of this paper appeared as the conference paper [41]; in
the current work, we provide further extensions and the proofs of all results
along with more extensive simulations.

#### Notations

Let $\mathbb{R}$, $\mathbb{R}_{+}$, $\mathbb{Z}$, and $\mathbb{Z}_{+}$
respectively denote the set of reals, nonnegative reals, integers, and
nonnegative integers. We denote by
$\mathbb{S}^{1}\simeq\mathbb{R}/(2\pi\mathbb{Z})$ the unit circle, and by
$\mathbb{T}^{n}:=\mathbb{S}^{1}\times\cdots\times\mathbb{S}^{1}$ the
$n$-torus. The $n$-vector whose entries are all $1$ is denoted by
$\text{1}_{n}$. For a matrix $A\in\mathbb{R}^{n\times m}$,
$\operatorname{im}A$ and $\operatorname{ker}A$ stand for the image and the
kernel of $A$, respectively. The imaginary unit is denoted by $\mathrm{j}$. If
$S$ is a subspace, then $S^{\perp}$ is the orthogonal complement of $S$ with
respect to the usual inner product. A function
$\alpha\colon[0,\infty)\to[0,\infty)$ is of class $\mathcal{K}$ if it is
continuous, strictly increasing, and $\alpha(0)=0$. A function
$\beta\colon[0,\infty)\to[0,\infty)$ is of class $\mathcal{L}$ if it is
continuous, nonincreasing, and $\beta(a)\to 0$ as $a\to\infty$. We define the
distance from a point $x\in\mathbb{R}^{n}$ to a subset
$S\subset\mathbb{R}^{n}$ by $\operatorname{dist}(x,S):=\inf_{y\in S}\|x-y\|$.

##  II Problem Setting and Preliminaries

In this section we describe the formal definition of cluster synchronization
and basic assumptions. Some preliminaries that will be useful in the
subsequent sections are also provided.

###  II-A Problem Setting

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected weighted graph,
where $\mathcal{V}=\\{1,\ldots,N\\}$ is the set of nodes and
$\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ is the set of edges.
Throughout this paper, it is assumed that $\mathcal{G}$ is connected since
cluster synchronization can be trivially achieved if the network is separated.
The weighted adjacency matrix of $\mathcal{G}$ is denoted by
$A=[a_{ij}]\in\mathbb{R}^{N\times N}$, where $a_{ij}>0$ if
$(i,j)\in\mathcal{E}$ and $a_{ij}=0$ if $(i,j)\not\in\mathcal{E}$. After
selecting an enumeration and an orientation for all (undirected) edges, we
define the oriented incidence matrix of $\mathcal{G}$ by
$B=[b_{ie}]\in\mathbb{R}^{N\times|\mathcal{E}|}$, where

$\displaystyle b_{ie}:=\begin{cases}1&\text{if node $i$ is the sink node of
edge $e$},\\\ -1&\text{if node $i$ is the source node of edge $e$},\\\
0&\text{otherwise}.\end{cases}$

We also define the weight matrix
$W:=\operatorname{diag}(\\{a_{ij}\\}_{(i,j)\in\mathcal{E}})$ of $\mathcal{G}$
based on the enumeration of edges introduced above.

Throughout this paper, each node represents an oscillator and each edge
represents an interconnection of two oscillators. We consider the Kuramoto
model governed by

$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i}),\quad
i\in\\{1,\ldots,N\\},$ (1)

where $\theta_{i}\colon\mathbb{R}_{+}\to\mathbb{S}^{1}$ is the phase of the
$i$th oscillator and $\omega_{i}\in\mathbb{R}$ is its natural frequency. The
equations in (1) can be written in vector form:

$\displaystyle\dot{\theta}=\omega-BW\sin(B^{\top}\theta),$ (2)

where $\theta:=[\theta_{1}\;\cdots\;\theta_{N}]^{\top}$ and
$\omega:=[\omega_{1}\;\cdots\;\omega_{N}]^{\top}$.

If a solution $\theta(\cdot)$ of (2) satisfies
$\operatorname{dist}(\theta(t),\operatorname{span}{(\text{1}_{N})})\to 0$ as
$t\to\infty$, then the solution is said to have achieved full synchronization.
Note that this can be achieved only if
$\omega\in\operatorname{span}(\text{1}_{N})$, that is, all the natural
frequencies are identical. The subject of this paper is different from the
situations of full synchronization. We are interested in cases where the
natural frequencies are not identical and the synchronization is realized in
some groups of oscillators. Next, we introduce the notion of this phenomenon,
called cluster synchronization.

Given an integer $r\geq 2$, we denote by
$\Pi=\\{\mathcal{C}_{1},\ldots,\mathcal{C}_{r}\\}$ a _nontrivial partition_ of
$\mathcal{V}$, i.e., the following conditions hold: (i)
$\mathcal{C}_{p}\neq\emptyset$ for all $p\in\\{1,\ldots,r\\}$, (ii)
$\mathcal{C}_{p}\cap\mathcal{C}_{q}=\emptyset$ if $p\neq q$, and (iii)
$\bigcup_{p=1}^{r}\mathcal{C}_{p}=\mathcal{V}$. We call
$\mathcal{C}_{1},\ldots,\mathcal{C}_{r}$ _clusters_ of the graph
$\mathcal{G}$. It is assumed without loss of generality that every cluster
includes at least two nodes. For $p\in\\{1,\ldots,r\\}$, let
$\mathcal{G}_{p}=(\mathcal{C}_{p},\mathcal{E}_{p})$ be the subgraph of
$\mathcal{G}$ associated with cluster $\mathcal{C}_{p}$, where
$\mathcal{E}_{p}:=\\{(i,j)\in\mathcal{E}:i,j\in\mathcal{C}_{p}\\}$. In this
paper, we assume that $\mathcal{G}_{p}$ is connected for every
$p\in\\{1,\ldots,r\\}$. Also, let
$\mathcal{G}_{\mathrm{inter}}=(\mathcal{V},\mathcal{E}_{\mathrm{inter}})$ be
the subgraph of $\mathcal{G}$ containing all the inter-cluster edges, where
$\mathcal{E}_{\mathrm{inter}}:=\mathcal{E}\setminus\bigcup_{p=1}^{r}\mathcal{E}_{p}$.
The graph $\mathcal{G}$ can then be decomposed as
$\mathcal{G}=\bigcup_{p=1}^{r}\mathcal{G}_{p}\cup\mathcal{G}_{\mathrm{inter}}$.

The _cluster synchronization manifold_ associated with partition
$\Pi=\\{\mathcal{C}_{1},\ldots,\mathcal{C}_{r}\\}$ is defined by

$\displaystyle\mathcal{S}_{\Pi}:=\\{\theta\in\mathbb{T}^{N}:\theta_{i}=\theta_{j},\
i,j\in\mathcal{C}_{p},\ p\in\\{1,\ldots,r\\}\\}.$

If a solution $\theta(\cdot)$ of (2) satisfies
$\operatorname{dist}(\theta(t),\mathcal{S}_{\Pi})\to 0$ as $t\to\infty$, then
the solution is said to have achieved cluster synchronization. That is, the
phase difference of any two nodes in the same cluster converges to zero, but
that of nodes in different clusters does not necessarily converge. Stability
of the cluster synchronization manifold is formally defined below.

###### Definition 1

Given a nontrivial partition
$\Pi=\\{\mathcal{C}_{1},\ldots,\mathcal{C}_{r}\\}$ of $\mathcal{V}$, the
cluster synchronization manifold $\mathcal{S}_{\Pi}$ of (2) is said to be
_(locally) exponentially stable_ if there exist $k\geq 1$, $\lambda>0$, and
$\delta>0$ such that $\operatorname{dist}(\theta(0),\mathcal{S}_{\Pi})<\delta$
implies $\operatorname{dist}(\theta(t),\mathcal{S}_{\Pi})\leq
k\mathrm{e}^{-\lambda t}\operatorname{dist}(\theta(0),\mathcal{S}_{\Pi})$ for
all $t\in\mathbb{R}_{+}$.

Clearly, this notion requires that the cluster synchronization manifold
$\mathcal{S}_{\Pi}$ is positive invariant for (2), that is,
$\theta(0)\in\mathcal{S}_{\Pi}$ implies $\theta(t)\in\mathcal{S}_{\Pi}$ for
all $t\in\mathbb{R}_{+}$. It is known that positive invariance of
$\mathcal{S}_{\Pi}$ is ensured under the following two assumptions [32, 42]:

###### Assumption 1

For every $p\in\\{1,\ldots,r\\}$, it holds that $\omega_{i}=\omega_{j}$ for
all $i,j\in\mathcal{C}_{p}$.

###### Assumption 2

For every $p,q\in\\{1,\ldots,r\\}$ with $p\neq q$, it holds that
$\sum_{k\in\mathcal{C}_{q}}a_{ik}=\sum_{k\in\mathcal{C}_{q}}a_{jk}$ for all
$i,j\in\mathcal{C}_{p}$.

Assumption 1 restricts the natural frequencies of oscillators in each cluster
to be identical. This is a natural requirement since cluster synchronization
implies that the oscillators rotate coherently within clusters. Assumption 2
states that any two nodes in the same cluster have equal weight sums with
respect to other clusters. Notice that this assumption depends only on inter-
cluster edge weights. It can thus be considered as a generalization of the
trivial case where there is no edge between clusters; in such a case, it is
clear that under Assumption 1, cluster synchronization can be realized.
However, the condition seems somewhat restrictive in that the equality in
Assumption 2 must hold strictly. We will relax this aspect in Section V. The
following proposition states that these assumptions are necessary for
analyzing stability of the cluster synchronization manifold [42].

###### Proposition 1

The cluster synchronization manifold is positive invariant for (2) if
Assumptions 1 and 2 hold. These are also necessary when, for every
$p,q\in\\{1,\ldots,r\\}$ with $p\neq q$,
$\max_{i\in\mathcal{C}_{p}}\dot{\theta}_{i}(t)\neq\max_{j\in\mathcal{C}_{q}}\dot{\theta}_{j}(t)$
holds for all $t\in\mathbb{R}_{+}$.

###### Remark 1

A graph partition satisfying the condition in Assumption 2 is called an
_external equitable partition_ , which was introduced in [43]. Such a graph
partition is known as a sufficient condition for the positive invariance of
cluster synchronization manifolds for many classes of systems [32]. In [42],
the authors showed that an external equitable partition is necessary and
sufficient for the invariance in the Kuramoto model when the oscillators’
frequencies are largely different between clusters as stated in Proposition 1.
We also remark that if the network size is large, there are a huge number of
partitions that belong to the class of external equitable partitions. The
trivial case is where every single node is grouped into a single cluster. For
algorithms to find the external equitable partitions with the minimal
cardinality, we refer to [44, 45].

###  II-B Preliminaries

To analyze synchronization, it is useful to consider relative phases instead
of phases themselves. Thus, we transform the Kuramoto model (2) into the
dynamics of phase differences. We then show that those dynamics can be written
by differential equations in a singular perturbation form.

Recall that we have decomposed the graph $\mathcal{G}$ as
$\mathcal{G}=\bigcup_{p=1}^{r}\mathcal{G}_{p}\cup\mathcal{G}_{\mathrm{inter}}$.
With appropriate ordering of the edge indices, we can write the incidence
matrix $B$ and the weight matrix $W$ in the following forms:

$\displaystyle
B=\begin{bmatrix}B_{\mathrm{intra}}&B_{\mathrm{inter}}\end{bmatrix},\quad
W=\begin{bmatrix}W_{\mathrm{intra}}&0\\\ 0&W_{\mathrm{inter}}\end{bmatrix},$

where $B_{\mathrm{intra}}$ and $B_{\mathrm{inter}}$ are the incidence matrices
of $\bigcup_{p=1}^{r}\mathcal{G}_{p}$ and $\mathcal{G}_{\mathrm{inter}}$,
respectively, and $W_{\mathrm{intra}}$ and $W_{\mathrm{inter}}$ are the weight
matrices of $\bigcup_{p=1}^{r}\mathcal{G}_{p}$ and
$\mathcal{G}_{\mathrm{inter}}$, respectively.

As with [24], we consider a spanning tree
$\mathcal{T}=(\mathcal{V},\tilde{\mathcal{E}})$ of $\mathcal{G}$, where
$\tilde{\mathcal{E}}\subset\mathcal{E}$ is the set of selected edges. Let
$\tilde{B}\in\mathbb{R}^{N\times(N-1)}$ be the incidence matrix of
$\mathcal{T}$. Without loss of generality, we can write
$\tilde{B}=[\tilde{B}_{\mathrm{intra}}\;\tilde{B}_{\mathrm{inter}}]$, where
$\tilde{B}_{\mathrm{intra}}$ and $\tilde{B}_{\mathrm{inter}}$ are associated
with the edges in $\bigcup_{p=1}^{r}\mathcal{E}_{p}$ and
$\mathcal{E}_{\mathrm{inter}}$, respectively. We denote by $n$ and $m$ the
dimensions of the range spaces of $\tilde{B}_{\mathrm{intra}}^{\top}$ and
$\tilde{B}_{\mathrm{inter}}^{\top}$, respectively. Note that $n+m=N-1$.

Now, we would like to calculate a matrix
$R\in\mathbb{R}^{|\mathcal{E}|\times(N-1)}$ such that
$B^{\top}=R\tilde{B}^{\top}$. Under an appropriate partition of $R$, we have
the relation

$\displaystyle\begin{bmatrix}B_{\mathrm{intra}}^{\top}\\\
B_{\mathrm{inter}}^{\top}\end{bmatrix}=R\begin{bmatrix}\tilde{B}_{\mathrm{intra}}^{\top}\\\
\tilde{B}_{\mathrm{inter}}^{\top}\end{bmatrix}=\begin{bmatrix}R_{1}&R_{2}\\\
R_{3}&R_{4}\end{bmatrix}\begin{bmatrix}\tilde{B}_{\mathrm{intra}}^{\top}\\\
\tilde{B}_{\mathrm{inter}}^{\top}\end{bmatrix}.$ (3)

Since $\tilde{B}$ has full column rank, we have
$R=B^{\top}(\tilde{B}^{\top})^{\dagger}=B^{\top}(\tilde{B}^{\dagger})^{\top}$.
The problem is thus to calculate the Moore–Penrose inverse
$\tilde{B}^{\dagger}$ of the column-wise partitioned matrix $\tilde{B}$. To
this end, we borrow a useful lemma from [46].

###### Lemma 1

Let $M_{1}\in\mathbb{R}^{m\times n_{1}}$ and $M_{2}\in\mathbb{R}^{m\times
n_{2}}$. Define the orthogonal projectors $P_{1},P_{2}\in\mathbb{R}^{m\times
m}$ by $P_{i}:=I_{m}-M_{i}M_{i}^{\dagger}$ for $i\in\\{1,2\\}$. Then, we have

$\displaystyle\begin{bmatrix}M_{1}&M_{2}\end{bmatrix}^{\dagger}=\begin{bmatrix}(P_{2}M_{1})^{\dagger}\\\
(P_{1}M_{2})^{\dagger}\end{bmatrix}$

if and only if $\operatorname{im}M_{1}\cap\operatorname{im}M_{2}=\\{0\\}$.

We define the two matrices
$P_{\mathrm{intra}}:=I_{N}-\tilde{B}_{\mathrm{intra}}\tilde{B}_{\mathrm{intra}}^{\dagger}$
and
$P_{\mathrm{inter}}:=I_{N}-\tilde{B}_{\mathrm{inter}}\tilde{B}_{\mathrm{inter}}^{\dagger}$.
Note that these are orthogonal projectors onto
$\operatorname{ker}\tilde{B}_{\mathrm{intra}}^{\top}$ and
$\operatorname{ker}\tilde{B}_{\mathrm{inter}}^{\top}$, respectively. Since
$\operatorname{im}\tilde{B}_{\mathrm{intra}}\cap\operatorname{im}\tilde{B}_{\mathrm{inter}}=\\{0\\}$,
Lemma 1 provides

$\displaystyle\tilde{B}^{\dagger}=\begin{bmatrix}\tilde{B}_{\mathrm{intra}}&\tilde{B}_{\mathrm{inter}}\end{bmatrix}^{\dagger}=\begin{bmatrix}(P_{\mathrm{inter}}\tilde{B}_{\mathrm{intra}})^{\dagger}\\\
(P_{\mathrm{intra}}\tilde{B}_{\mathrm{inter}})^{\dagger}\end{bmatrix}.$

It thus follows that

$\displaystyle R$ $\displaystyle=\begin{bmatrix}B_{\mathrm{intra}}^{\top}\\\
B_{\mathrm{inter}}^{\top}\end{bmatrix}\begin{bmatrix}(\tilde{B}_{\mathrm{intra}}^{\top}P_{\mathrm{inter}})^{\dagger}&(\tilde{B}_{\mathrm{inter}}^{\top}P_{\mathrm{intra}})^{\dagger}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}B_{\mathrm{intra}}^{\top}(\tilde{B}_{\mathrm{intra}}^{\top}P_{\mathrm{inter}})^{\dagger}&B_{\mathrm{intra}}^{\top}(\tilde{B}_{\mathrm{inter}}^{\top}P_{\mathrm{intra}})^{\dagger}\\\
B_{\mathrm{inter}}^{\top}(\tilde{B}_{\mathrm{intra}}^{\top}P_{\mathrm{inter}})^{\dagger}&B_{\mathrm{inter}}^{\top}(\tilde{B}_{\mathrm{inter}}^{\top}P_{\mathrm{intra}})^{\dagger}\end{bmatrix}.$

We notice that
$\operatorname{ker}B_{\mathrm{intra}}^{\top}=\operatorname{ker}\tilde{B}_{\mathrm{intra}}^{\top}$,
which follows from the fact that $B_{\mathrm{intra}}^{\top}\text{1}_{n}=0$ and
$\tilde{B}_{\mathrm{intra}}^{\top}\text{1}_{n}=0$. Hence, we have
$B_{\mathrm{intra}}^{\top}(\tilde{B}_{\mathrm{inter}}^{\top}P_{\mathrm{intra}})^{\dagger}=0$
since
$\operatorname{im}P_{\mathrm{intra}}=\operatorname{ker}B_{\mathrm{intra}}^{\top}$.
The four matrices in (3) are finally given by

$\displaystyle R_{1}$
$\displaystyle=B_{\mathrm{intra}}^{\top}(\tilde{B}_{\mathrm{intra}}^{\top}P_{\mathrm{inter}})^{\dagger},$
$\displaystyle\quad R_{2}$ $\displaystyle=0,$ $\displaystyle R_{3}$
$\displaystyle=B_{\mathrm{inter}}^{\top}(\tilde{B}_{\mathrm{intra}}^{\top}P_{\mathrm{inter}})^{\dagger},$
$\displaystyle\quad R_{4}$
$\displaystyle=B_{\mathrm{inter}}^{\top}(\tilde{B}_{\mathrm{inter}}^{\top}P_{\mathrm{intra}})^{\dagger}.$

Define the vectors of intra- and inter-cluster phase differences as
$x:=\tilde{B}_{\mathrm{intra}}^{\top}\theta$ and
$z:=\tilde{B}_{\mathrm{inter}}^{\top}\theta$, respectively. For the Kuramoto
model (2), we choose a perturbation parameter $\varepsilon$ as the reciprocal
of the smallest natural frequency difference between clusters. Without loss of
generality, we assume $\varepsilon$ to be positive, namely,

$\displaystyle\varepsilon:=\frac{1}{\min_{(i,j)\in\mathcal{E}_{\mathrm{inter}}}|\omega_{i}-\omega_{j}|}.$
(4)

Then, under Assumption 1, we obtain the state equations in the following
singular perturbation form:

$\displaystyle\dot{x}$ $\displaystyle=f(x,z),$ (5)
$\displaystyle\varepsilon\dot{z}$ $\displaystyle=\eta+\varepsilon g(x,z),$ (6)

where $\eta:=\varepsilon\tilde{B}_{\mathrm{inter}}^{\top}\omega$ is a constant
vector representing the ratio of natural frequency differences with respect to
$\varepsilon$, and

$\displaystyle f(x,z)$
$\displaystyle:=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{intra}}W_{\mathrm{intra}}\sin(R_{1}x)$
$\displaystyle\quad{}\quad{}-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\sin(R_{3}x+R_{4}z),$
$\displaystyle g(x,z)$
$\displaystyle:=-\tilde{B}_{\mathrm{inter}}^{\top}B_{\mathrm{intra}}W_{\mathrm{intra}}\sin(R_{1}x)$
$\displaystyle\quad{}\quad{}-\tilde{B}_{\mathrm{inter}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\sin(R_{3}x+R_{4}z).$

The equations (5) and (6), respectively, describe the intra- and inter-cluster
dynamics.

We now introduce the concept of partial stability defined in Chap. 4 of [47].

###### Definition 2

Consider the dynamical system described by

$\displaystyle\dot{x}_{1}=f_{1}(x_{1},x_{2}),\quad\dot{x}_{2}=f_{2}(x_{1},x_{2}).$

Suppose that $f_{1}(0,x_{2})=0$ for all $x_{2}$. Then, the zero solution
$x_{1}(t)\equiv 0$ is said to be _partially exponentially stable uniformly in
$x_{2}$_ if there exist $k\geq 1$, $\lambda>0$, and $\delta>0$ such that
$\|x_{1}(0)\|<\delta$ implies $\|x_{1}(t)\|\leq k\mathrm{e}^{-\lambda
t}\|x_{1}(0)\|$ for all $t\in\mathbb{R}_{+}$, and this statement holds
uniformly with respect to the initial state $x_{2}(0)$.

Recall that under Assumptions 1 and 2, the cluster synchronization manifold
$\mathcal{S}_{\Pi}$ is positive invariant for (2). This implies that
$f(0,z)=0$ for all $z\in\mathbb{R}^{m}$. For the Kuramoto model (2) and its
associated state equations (5) and (6), exponential stability of
$\mathcal{S}_{\Pi}$ in the sense of Definition 1 is equivalent to partial
exponential stability of $x(t)\equiv 0$ in the sense of Definition 2.
Therefore, in what follows, we will explore partial stability of the dynamical
system described by (5) and (6) to derive conditions for cluster
synchronization.

##  III Stability Theory via Averaging

In this section we present a sufficient condition for stability of a class of
linear time-varying systems. We use an averaging principle to show that there
is a nonmonotonic Lyapunov function for the given linear time-varying system
characterized by a small perturbation parameter. The result in this section
will play a crucial role for deriving stability conditions for cluster
synchronization in Section IV.

First, we give a general result involving nonmonotonic Lyapunov functions (see
[48] for details). Briefly, different from the requirement of $V$ to be a
Lyapunov function, $V(x(t))$ must decrease only at certain sampled instants
along solutions $x(t)$ and time intervals where it increases may be present.

###### Proposition 2

Consider the nonlinear system described by

$\displaystyle\dot{x}=f(t,x),$ (7)

where $f\colon\mathbb{R}_{+}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a
continuous function. Suppose that $f(t,0)=0$ for all $t\in\mathbb{R}_{+}$ and
that the solutions of (7) for $t\in\mathbb{R}_{+}$ exist and are unique in a
neighborhood of $x=0$. Let $V\colon\mathbb{R}^{n}\to\mathbb{R}$ be a
continuous function such that the following properties hold:

  1. 1.

There exist $c_{1},c_{2}>0$ such that

$\displaystyle c_{1}\|x\|^{2}\leq V(x)\leq c_{2}\|x\|^{2}$

for all $x\in\mathbb{R}^{n}$.

  2. 2.

For some time sequence $\\{t_{k}\\}_{k\in\mathbb{Z}_{+}}$ with the properties
that $t_{0}=0$, $t_{k}<t_{k+1}$ for all $k\in\mathbb{Z}_{+}$, and
$t_{k}\to\infty$ as $k\to\infty$, there exist $c_{3}>0$ and $\delta>0$ such
that for any solution $x\colon\mathbb{R}_{+}\to\mathbb{R}^{n}$ of (7)
satisfying $\|x(0)\|<\delta$, it holds

$\displaystyle\frac{V(x(t_{k+1}))-V(x(t_{k}))}{t_{k+1}-t_{k}}\leq-
c_{3}\|x(t_{k})\|^{2}$

for all $k\in\mathbb{Z}_{+}$.

Then, the zero solution $x(t)\equiv 0$ is exponentially stable.

Now, consider the system of linear differential equations

$\displaystyle\dot{x}=\varepsilon A(t,\varepsilon)x,$ (8)

where $A\colon\mathbb{R}_{+}\times\mathbb{R}_{+}\to\mathbb{R}^{n\times n}$ is
continuous and bounded in its arguments and $\varepsilon>0$ is a small
parameter. Assume that the time average of $A(t,\varepsilon)$ exists for
$\varepsilon=0$. That is, the matrix

$\displaystyle
A_{\mathrm{av}}:=\lim_{T\to\infty}\frac{1}{T}\int_{t}^{t+T}A(s,0)\,\mathrm{d}s$
(9)

can be defined uniformly in $t$. Also, it is assumed that there exist a class
$\mathcal{L}$ function $\beta\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ and a
class $\mathcal{K}$ function $\alpha\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$
such that

$\displaystyle\left\|\frac{1}{T}\int_{t}^{t+T}[A(s,\varepsilon)-A_{\mathrm{av}}]\,\mathrm{d}s\right\|\leq\beta(T)+\alpha(\varepsilon)$
(10)

for all $t\in\mathbb{R}_{+}$.

The following theorem can be established using Proposition 2.

###### Theorem 1

Suppose that the matrix $A_{\mathrm{av}}$ in (9) is a Hurwitz matrix and that
the inequality (10) holds. Then, there exists $\varepsilon^{*}>0$ such that if
$\varepsilon<\varepsilon^{*}$, then the zero solution $x(t)\equiv 0$ of (8) is
exponentially stable.

###### Proof:

See Appendix -A. ∎

###### Remark 2

Here, we give some technical comments on Theorem 1 in comparison with existing
results. By the method of averaging, a Lyapunov function for the average
system can be employed to prove stability of the original system (see Chap. 10
in [40]). To apply the method in [40], it is required that $A(t,\varepsilon)$
in (8) is Lipschitz with respect to $\varepsilon$ uniformly in $t$, i.e.,
there exists $L>0$ such that for every $\varepsilon_{1},\varepsilon_{2}>0$
sufficiently small, we have $\|A(t,\varepsilon_{1})-A(t,\varepsilon_{2})\|\leq
L|\varepsilon_{1}-\varepsilon_{2}|$ uniformly in $t$. In the present case,
however, the matrix $A(t,\varepsilon)$ in (8) does not satisfy the Lipschitz
condition. Thus, we cannot use the conventional method in [40] to prove
stability of our system. To overcome this difficulty, we employed nonmonotonic
Lyapunov functions and applied Proposition 2 to prove stability. The use of
nonmonotonic Lyapunov functions was also considered for linear time-varying
systems in [49, 50]. In those papers, the authors addressed the case where
$A(t,\varepsilon)$ is independent of $\varepsilon$. Thus, Theorem 1 can be
viewed as an extension of the results in [49, 50]. On the other hand, the
authors of [51] considered the case where the uniform average in (9) does not
exist, and the time-dependent average matrix was used in stability analysis.
In contrast, we consider the case where the uniform average in (9) can be
defined for $\varepsilon=0$, but it may not necessarily be defined for
$\varepsilon>0$.

###### Remark 3

We note that in Proposition 2, the positive-definite function $V$ need not be
differentiable. That is, one has the advantage that a nonsmooth Lyapunov
function candidate can be used in stability analysis. This property may be
useful to study synchronization problems. For related results using nonsmooth
Lyapunov functions, we refer to, e.g., [52].

##  IV Cluster Synchronization Conditions for Kuramoto Oscillators

In this section we provide our main results. Through the analytical method
developed in the previous section, we derive sufficient conditions for
exponential stability of a cluster synchronization manifold associated with
the Kuramoto model. This section addresses two problems. First, we consider
cluster synchronization with two or more clusters. Then, we restrict our
attention to the case of two clusters and carry out further analysis.
Numerical examples are also presented to demonstrate our results.

###  IV-A Multi-Cluster Case

Recall that the problem under consideration is to analyze partial stability of
the nonlinear system described by (5) and (6). First, we linearize the intra-
cluster dynamics (5) at $x=0$, which corresponds to the cluster
synchronization manifold, while $z$ is interpreted as a function of time. To
do so, let us introduce a new time variable $\tau:=t/\varepsilon$. In this
time scale, the equations (5) and (6) are written as

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}\tau}$ $\displaystyle=\varepsilon
f(x,z),$ (11) $\displaystyle\frac{\mathrm{d}z}{\mathrm{d}\tau}$
$\displaystyle=\eta+\varepsilon g(x,z).$ (12)

For a fixed $\varepsilon$, let $\zeta(\cdot,\varepsilon)$ denote the solution
of the initial value problem for (12) under $x=0$:

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}\tau}=\eta+\varepsilon g(0,z),\quad
z(0)=z_{0},$ (13)

where we have omitted the dependence on $z_{0}$ because it does not affect the
results. Note that when $\varepsilon=0$, the solution to (13) is given by
$\zeta(\tau,0)=z_{0}+\eta\tau$. Linearization of $f(x,z)$ around $x=0$ and
$z=\zeta(\tau,\varepsilon)$ leads us to the Jacobian

$\displaystyle J(\tau,\varepsilon)$ $\displaystyle:=\left.\frac{\partial
f(x,z)}{\partial x}\right|_{(x,z)=(0,\zeta(\tau,\varepsilon))}$
$\displaystyle=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{intra}}W_{\mathrm{intra}}R_{1}$
$\displaystyle\quad{}-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\operatorname{diag}[R_{4}\cos(\zeta(\tau,\varepsilon))]R_{3}.$
(14)

The linearized system is thus described by

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}\tau}=\varepsilon
J(\tau,\varepsilon)x.$ (15)

The following proposition motivates us to analyze stability based on the
linearization principle.

###### Proposition 3

Consider the nonlinear system (5) and (6). The zero solution $x(t)\equiv 0$ is
partially exponentially stable uniformly in $z$ if and only if the linear
time-varying system (15) is (uniformly) exponentially stable.

###### Proof:

Interpreting $z$ as a function of time, we can use the standard stability
analysis techniques (e.g., Chap. 4 of [40]) to conclude the proof.
Specifically, the sufficiency part follows from Lyapunov’s indirect method and
the necessity part follows from the converse Lyapunov theorem. ∎

We are ready to derive cluster synchronization conditions. Particularly, we
analyze stability of (15) based on Theorem 1. From (4), we note that the
smaller $\varepsilon$ is, the larger the inter-cluster natural frequency
differences are. The following result provides sufficient conditions for
cluster synchronization. Its proof is given in the next subsection.

###### Theorem 2

Let $\mathcal{S}_{\Pi}$ be the cluster synchronization manifold associated
with the Kuramoto model (2), where
$\Pi=\\{\mathcal{C}_{1},\ldots,\mathcal{C}_{r}\\}$ is a nontrivial partition
of $\mathcal{V}$ such that the subgraphs
$\mathcal{G}_{1},\ldots,\mathcal{G}_{r}$ are all connected. Suppose that
Assumptions 1 and 2 hold. Then, the following statements hold:

  1. 1.

There exists $a^{*}>0$ such that if

$\displaystyle a_{ij}<a^{*},\quad(i,j)\in\mathcal{E}_{\mathrm{inter}},$

then $\mathcal{S}_{\Pi}$ is exponentially stable.

  2. 2.

There exists $\omega^{*}>0$ such that if

$\displaystyle|\omega_{i}-\omega_{j}|>\omega^{*},\quad(i,j)\in\mathcal{E}_{\mathrm{inter}},$

then $\mathcal{S}_{\Pi}$ is exponentially stable.

Theorem 2 indicates that cluster synchronization can be achieved if (i) the
inter-cluster coupling strengths are sufficiently weak and/or (ii) the natural
frequencies are sufficiently different between clusters. Our approach provides
a unified framework for stability analysis in that the two conditions 1) and
2) in the above theorem are simultaneously derived. We must highlight that
this is achieved based on the use of a single nonmonotonic Lyapunov function.
This allows us to characterize a certain tradeoff between the two conditions.
We must note that these two stability conditions are similar to those derived
in [24] (Theorems 3.2 and 3.3). There, different methods are applied to obtain
the two conditions separately regarding the inter-cluster edge weights and the
natural frequency differences.

Furthermore, our result shows that there is a critical value of the natural
frequency differences between clusters, above which the cluster
synchronization manifold is exponentially stable. This has an advantage since
Theorem 3.3 in [24] requires that $|\omega_{i}-\omega_{j}|\to\infty$ for all
$(i,j)\in\mathcal{E}_{\mathrm{inter}}$. The main difference in the analysis
lies in the fact that our result does not depend on a particular Lyapunov
function. Instead, we found a nonmonotonic Lyapunov function for the cluster
synchronization manifold to prove the exponential stability via Proposition 2.
We will see in Section IV-B that a specific value of the natural frequency
difference for which cluster synchronization is realized can be found.

###### Remark 4

The analysis in this section is based on the so-called two-timescale
averaging, that is, averaging is applied to a singularly perturbed system
[53]. However, our goal was to show partial stability of the slow dynamics and
stability of the fast dynamics is not involved. A similar setting was
considered in [54], where the fast subsystem is restricted to a scalar system
and the fast state is interpreted as a time variable. In contrast, the fast
subsystem given in (6) has vector states, and hence, the result in [54] is not
applicable.

###### Remark 5

Synchronization in networks with periodically switching topologies was
considered in [55], and the nonperiodic case was explored in [56]. In those
papers, it is allowed that the graph is disconnected at certain time instants,
but it must be connected in some sense. By contrast, our setting is more
general since the network dynamics becomes unstable in some time intervals,
that is, the pointwise Jacobian in (15) may have unstable eigenvalues. We have
utilized the inequality (10) to show that fast oscillations generated from the
inter-cluster dynamics support cluster synchronization. Moreover, we note that
our results can also address chimera states considered in [57, 58], where
synchronous and asynchronous groups coexist in the network. This is because
each oscillator in the asynchronous group can be regarded as an isolated
cluster.

###  IV-B Proof of Theorem 2

Consider the linearized system (15). We define the uniform average of
$J(\tau,\varepsilon)$ by

$\displaystyle J_{\mathrm{av}}$
$\displaystyle:=\lim_{T\to\infty}\frac{1}{T}\int_{\tau}^{\tau+T}J(s,0)\,\mathrm{d}s$
$\displaystyle=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{intra}}W_{\mathrm{intra}}R_{1}.$
(16)

The above matrix is well defined because each element of $\cos(\zeta(\tau,0))$
in (14) is periodic in $\tau$. If these elements have a common period, then
the vector-valued function $\cos(\zeta(\tau,0))$ is periodic in $\tau$, and
otherwise it is almost periodic in $\tau$. Clearly, $J_{\mathrm{av}}$ is
determined by only the intra-cluster edge weights, and thus, it can be written
as $J_{\mathrm{av}}=\operatorname{blk-diag}(J_{1},\ldots,J_{r})$. Note that
$J_{i}$ is a Hurwitz matrix if and only if $\mathcal{G}_{i}$ is connected (see
Lemma 3.1 in [24]). Therefore, the average Jacobian $J_{\mathrm{av}}$ is a
Hurwitz matrix.

The following lemma characterizes the relation between the original and the
average Jacobians.

###### Lemma 2

There exists $\gamma>0$ such that

  1. (i)

$\gamma=O(\|W_{\mathrm{inter}}\|_{\infty})$ and

  2. (ii)

for all $\varepsilon>0$ and all $T>0$, the following inequality holds
uniformly in $\tau\in\mathbb{R}_{+}$:

$\displaystyle\left\|\frac{1}{T}\int_{\tau}^{\tau+T}[J(s,\varepsilon)-J_{\mathrm{av}}]\,\mathrm{d}s\right\|\leq\gamma\left(\frac{1}{T}+\varepsilon\right).$

###### Proof:

See Appendix -B. ∎

It can be observed from Lemma 2 that the time average of the difference
between $J(\tau,\varepsilon)$ and $J_{\mathrm{av}}$ is small enough if one of
the following holds:

  1. 1.

$\gamma$ is small;

  2. 2.

$T$ is large and $\varepsilon$ is small.

To satisfy the condition in (10), we can choose a class $\mathcal{L}$ function
$\beta(T)=\gamma/T$ and a class $\mathcal{K}$ function
$\alpha(\varepsilon)=\gamma\varepsilon$.

Let $\Phi(\tau+T,\tau)$ be the state-transition matrix of (15), and let
$\Phi_{\mathrm{av}}(\tau+T,\tau):=\mathrm{e}^{\varepsilon TJ_{\mathrm{av}}}$.
Then, we consider
$H(\tau+T,\tau):=\Phi(\tau+T,\tau)-\Phi_{\mathrm{av}}(\tau+T,\tau)$. By the
Peano–Baker series of the transition matrices, we have

$\displaystyle H(\tau+T,\tau)$
$\displaystyle=\varepsilon\int_{\tau}^{\tau+T}[J(s,\varepsilon)-J_{\mathrm{av}}]\,\mathrm{d}s$
$\displaystyle\quad{}+\sum_{i=2}^{\infty}\varepsilon^{i}\int_{\tau}^{\tau+T}J(s_{1},\varepsilon)\cdots\int_{t}^{s_{i-1}}J(s_{i},\varepsilon)\,\mathrm{d}s_{i}\cdots\mathrm{d}s_{1}$
$\displaystyle\quad{}-\sum_{i=2}^{\infty}\frac{(\varepsilon
T)^{i}J_{\mathrm{av}}^{i}}{i!}.$

Thus, it can be obtained that

$\displaystyle\|H(\tau+T,\tau)\|$
$\displaystyle\leq\varepsilon\left\|\int_{\tau}^{\tau+T}[J(s,\varepsilon)-J_{\mathrm{av}}]\,\mathrm{d}s\right\|$
$\displaystyle\quad{}+2\sum_{i=2}^{\infty}\frac{(\rho\varepsilon T)^{i}}{i!},$

where
$\rho:=\sup_{\tau\in\mathbb{R}_{+},\varepsilon>0}\|J(\tau,\varepsilon)\|$.
Notice from (14) that $\rho$ is finite. Then, Lemma 2 implies that

$\displaystyle\|H(\tau+T,\tau)\|\leq\gamma\varepsilon
T\left(\frac{1}{T}+\varepsilon\right)+2(\mathrm{e}^{\rho\varepsilon
T}-1-\rho\varepsilon T).$

We note that $J_{\mathrm{av}}$ has only the negative of the nonzero
eigenvalues of the graph Laplacians associated with the subgraphs
$\mathcal{G}_{1},\ldots,\mathcal{G}_{r}$. This implies that
$\|\Phi_{\mathrm{av}}(\tau+T,\tau)\|\leq\mathrm{e}^{-\lambda_{2}\varepsilon
T}$, where $\lambda_{2}$ is the smallest algebraic connectivity of
$\mathcal{G}_{1},\ldots,\mathcal{G}_{r}$. We thus obtain

$\displaystyle\|\Phi(\tau+T,\tau)\|$
$\displaystyle\leq\|\Phi_{\mathrm{av}}(\tau+T,\tau)\|+\|H(\tau+T,\tau)\|$
$\displaystyle\leq\mathrm{e}^{-\lambda_{2}\varepsilon T}+\gamma\varepsilon
T\left(\frac{1}{T}+\varepsilon\right)$
$\displaystyle\quad{}+2(\mathrm{e}^{\rho\varepsilon T}-1-\rho\varepsilon T)$
$\displaystyle=:\kappa(\gamma,\varepsilon,T).$ (17)

The problem is to show $\kappa(\gamma,\varepsilon,T)<1$ for some $T>0$ if
$\gamma$ or $\varepsilon$ is sufficiently small. To find such a constant $T$,
we first let $\varepsilon$ be fixed. At $\gamma=0$, we have
$\kappa(0,\varepsilon,T)=\mathrm{e}^{-\lambda_{2}\varepsilon
T}+2(\mathrm{e}^{\rho\varepsilon T}-1-\rho\varepsilon T)$. Because
$\rho<\infty$, there exists $T>0$ such that $\kappa(0,\varepsilon,T)<1$. This
means by continuity that $\kappa(\gamma,\varepsilon,T)<1$ is satisfied for
sufficiently small $\gamma>0$, that is, small
$\|W_{\mathrm{inter}}\|_{\infty}$ by Lemma 2 (i); by the definition of
$W_{\mathrm{inter}}$, this corresponds to taking the edge weights of
$\mathcal{G}_{\mathrm{inter}}$ to be small. Similarly, if we let $\gamma$ be
fixed, the same conclusion can be obtained for sufficiently small
$\varepsilon>0$, whose procedure has been shown in the proof of Theorem 1.
Notice that a specific value of $\varepsilon$ such that
$\kappa(\gamma,\varepsilon,T)<1$ for a suitably chosen $T$ can be found
numerically. By definition of $\varepsilon$ in (4), we conclude that the
natural frequencies are required to be sufficiently different among clusters.
The proof is now complete.

###  IV-C Two-Cluster Case

In the previous subsections we have not restricted the number of clusters.
However, the obtained stability conditions are somewhat conservative. We here
consider cluster synchronization with two clusters ($r=2$) and derive a more
specific stability condition, which may be less conservative.

For $i\in\mathcal{C}_{1}$ and $j\in\mathcal{C}_{2}$, define

$\displaystyle\bar{\omega}:=|\omega_{i}-\omega_{j}|,\quad\bar{a}:=\sum_{k\in\mathcal{C}_{2}}a_{ik}+\sum_{k\in\mathcal{C}_{1}}a_{jk}.$

We restrict our attention to the case where $\bar{\omega}>\bar{a}$. In that
case, the solutions to (12) converge to a periodic orbit with the period

$\displaystyle
T_{2}=\frac{2\pi\bar{\omega}}{\sqrt{\bar{\omega}^{2}-\bar{a}^{2}}}$ (18)

(see Lemma 3.4 in [24]). If $\bar{\omega}<\bar{a}$, the solutions to (12)
converge to constant values. This corresponds to frequency synchronization as
considered in [22].

Given a nontrivial partition $\Pi=\\{\mathcal{C}_{1},\mathcal{C}_{2}\\}$, one
can assume without loss of generality that
$R_{4}=\text{1}_{|\mathcal{E}_{\mathrm{inter}}|}$. We note that
$\varepsilon=1/\bar{\omega}$. Similarly to the previous subsection, we obtain
the linearized system

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}\tau}=\varepsilon
J_{2}(\tau,\varepsilon)x,$ (19)

where
$J_{2}(\tau,\varepsilon):=J_{\mathrm{intra}}+J_{\mathrm{inter}}\cos(\zeta(\tau,\varepsilon))$.
Here, we have defined
$J_{\mathrm{intra}}:=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{intra}}W_{\mathrm{intra}}R_{1}$
and
$J_{\mathrm{inter}}:=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}R_{3}$.
We remark that the Jacobian cannot be transformed into the above form in the
multi-cluster case.

The following result provides an explicit stability condition. In particular,
this condition includes that in Theorem 3.6 in [24], where
$J_{\mathrm{intra}}$ is assumed to be of the form $J_{\mathrm{intra}}=\alpha
I$ for some $\alpha<0$.

###### Theorem 3

Let $\mathcal{S}_{\Pi}$ be the cluster synchronization manifold associated
with the Kuramoto model (2), where $\Pi=\\{\mathcal{C}_{1},\mathcal{C}_{2}\\}$
is a nontrivial partition of $\mathcal{V}$ such that the subgraphs
$\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are connected. Suppose that
Assumptions 1 and 2 hold and that $\bar{\omega}>\bar{a}$. If the matrices
$J_{\mathrm{intra}}$ and $J_{\mathrm{inter}}$ commute, i.e., if
$J_{\mathrm{intra}}J_{\mathrm{inter}}=J_{\mathrm{inter}}J_{\mathrm{intra}}$,
then $\mathcal{S}_{\Pi}$ is exponentially stable.

###### Proof:

The inter-cluster dynamics under $x=0$ is governed by

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}\tau}=1-\varepsilon\bar{a}\sin(z).$

Since $\varepsilon\bar{a}<1$, the solution is periodic with the period $T_{2}$
in (18) and has zero average. This implies that the uniform average is given
as

$\displaystyle
J_{2,\mathrm{av}}:=\frac{1}{T_{2}}\int_{\tau}^{\tau+T_{2}}J_{2}(s,0)\,\mathrm{d}s=J_{\mathrm{intra}}.$

Notice that for the two cluster case, the inequality in Lemma 2 is translated
into the following equality:

$\displaystyle\frac{1}{T_{2}}\int_{\tau}^{\tau+T_{2}}[J_{2}(s,\varepsilon)-J_{2,\mathrm{av}}]\,\mathrm{d}s=0.$
(20)

Since $J_{\mathrm{intra}}$ and $J_{\mathrm{inter}}$ commute, we have

$\displaystyle
J_{2}(\tau,\varepsilon)\int_{\sigma}^{\tau}J_{2}(s,\varepsilon)\,\mathrm{d}s=\int_{\sigma}^{\tau}J_{2}(s,\varepsilon)\,\mathrm{d}sJ_{2}(\tau,\varepsilon)$

for all $\sigma,\tau\in\mathbb{R}_{+}$. In the proof of Theorem 2, we employed
the Peano–Baker series of state transition matrices. This takes a complex form
in general, leading us to a conservative stability condition. Here, by the
above property, the Peano–Baker series becomes simpler in its form (see Chap.
4 of [59]). Let $\Phi_{2}(\tau+T_{2},\tau)$ denote the state transition matrix
of (19), and let
$\Phi_{2,\mathrm{av}}(\tau+T_{2},\tau):=\mathrm{e}^{\varepsilon
T_{2}J_{2,\mathrm{av}}}$. Then, we have

$\displaystyle\Phi_{2}(\tau+T_{2},\tau)-\Phi_{2,\mathrm{av}}(\tau+T_{2},\tau)$
$\displaystyle=\sum_{i=2}^{\infty}\frac{\varepsilon^{i}}{i!}\left(\int_{\tau}^{\tau+T_{2}}[J_{2}(s,\varepsilon)-J_{2,\mathrm{av}}]\,\mathrm{d}s\right)^{i}.$

Due to (20), it follows that
$\Phi_{2}(\tau+T_{2},\tau)-\Phi_{2,\mathrm{av}}(\tau+T_{2},\tau)=0$, which
implies

$\displaystyle\Phi_{2}(\tau+T_{2},\tau)=\Phi_{2,\mathrm{av}}(\tau+T_{2},\tau)=\mathrm{e}^{\varepsilon
T_{2}J_{2,\mathrm{av}}}.$

Since $J_{2,\mathrm{av}}$ is a Hurwitz matrix,
$\|\Phi_{2}(\tau+T_{2},\tau)\|<1$ can be established. Therefore, we conclude
the proof. ∎

###### Remark 6

Two-cluster synchronization can be considered as a nonlinear version of
_bipartite consensus_ , which is studied in the presence of antagonistic
interactions [60, 61]. The Kuramoto model with antagonistic interactions is
considered in [62], where the underlying network is represented by a signed
graph. In this paper we do not consider signed graphs explicitly, but our
results can be extended to that setting by utilizing a signed external
equitable partition condition [32] in place of Assumption 2.

###  IV-D Numerical Examples

We present two numerical examples. The first one is related to Theorem 2 and
the second one demonstrates Theorem 3.

Figure 1: Network partitioned into three clusters Figure 2: Time responses of
phases of nodes in the three clusters (Example 1)

###### Example 1

Consider the network shown in Fig. 1, which is partitioned into three clusters
as indicated. The adjacency matrix of this network is given by

$\displaystyle
A=\left[\begin{array}[]{cc|ccc|ccc}0&a_{1}&b_{1}&b_{1}/2&0&0&0&0\\\
a_{1}&0&0&b_{1}/2&b_{1}&0&0&0\\\ \hline\cr b_{1}&0&0&a_{2}&0&b_{2}&0&0\\\
b_{1}/2&b_{1}/2&a_{2}&0&a_{2}&0&b_{2}&0\\\ 0&b_{1}&0&a_{2}&0&0&0&b_{2}\\\
\hline\cr 0&0&b_{2}&0&0&0&a_{3}&a_{3}\\\ 0&0&0&b_{2}&0&a_{3}&0&a_{3}\\\
0&0&0&0&b_{2}&a_{3}&a_{3}&0\end{array}\right].$

Notice that in this case, Assumption 2 is satisfied for any parameter value.
We set all the edge weights to be $1$, that is,
$a_{1}=a_{2}=a_{3}=b_{1}=b_{2}=1$. The natural frequencies are set as
$\omega_{i}=0$ for $i\in\mathcal{C}_{1}$, $\omega_{i}=5$ for
$i\in\mathcal{C}_{2}$, and $\omega_{i}=10$ for $i\in\mathcal{C}_{3}$. We
present in Fig. 2 the phase evolutions of the oscillators, where the
oscillators are synchronized within clusters and the inter-cluster
oscillations can be observed. Theoretically, we obtain the tradeoff curve for
which $\kappa(\gamma,\varepsilon,T)$ in the proof of Theorem 2 is equal to one
as shown in Fig. 3. The cluster synchronization manifold is guaranteed to be
exponentially stable for any pair $(\gamma,\varepsilon)$ in the left-bottom
area.

Figure 3: Tradeoff curve for connectivity and heterogeneity between clusters

###### Example 2

Consider the network shown in Fig. 4, where the network is partitioned into
two clusters. The adjacency matrix is given by

$\displaystyle A=\left[\begin{array}[]{ccc|ccc}0&a_{1}&0&0&0&b\\\
a_{1}&0&a_{1}&0&b&0\\\ 0&a_{1}&0&b&0&0\\\ \hline\cr 0&0&b&0&a_{2}&0\\\
0&b&0&a_{2}&0&a_{2}\\\ b&0&0&0&a_{2}&0\end{array}\right].$

In this case we have the linearized system (19) with

$\displaystyle J_{\mathrm{intra}}=\begin{bmatrix}-2a_{1}&a_{1}&0&0\\\
a_{1}&-2a_{1}&0&0\\\ 0&0&-2a_{2}&a_{2}\\\ 0&0&a_{2}&-2a_{2}\end{bmatrix},$

which is a Hurwitz matrix, and

$\displaystyle J_{\mathrm{inter}}=\begin{bmatrix}-b&0&0&-b\\\ 0&-b&-b&0\\\
0&-b&-b&0\\\ -b&0&0&-b\end{bmatrix}.$

One can observe that $J_{\mathrm{intra}}$ and $J_{\mathrm{inter}}$ commute if
$a_{1}=a_{2}$ or $b=0$. Thus, we set $a_{1}=a_{2}=1$ and $b=1$. The simulation
result with $\omega_{i}=5$ for $i\in\mathcal{C}_{1}$ and $\omega_{i}=1$ for
$i\in\mathcal{C}_{2}$ is shown in Fig. 5 (top). In this setting, Theorem 3 can
provide a less conservative stability condition compared with Theorem 2, which
requires that $\bar{\omega}>125$ for achieving cluster synchronization. Next,
we would like to present an example where cluster synchronization does not
occur. We change the edge weights within cluster $\mathcal{C}_{1}$ to a much
smaller value $a_{1}=0.01$. Then, $J_{\mathrm{intra}}$ and
$J_{\mathrm{inter}}$ do not commute, and the simulation result with the same
natural frequencies is presented in Fig. 5 (bottom). This indicates that the
cluster synchronization manifold is not necessarily stable.

Figure 4: Network with two clusters Figure 5: Time responses of phases of
nodes in the two clusters (Example 2, stable case (top) and unstable case
(bottom))

##  V Practical Stability without External Equitable Partitions

In Section IV we have studied cluster synchronization under the assumption of
external equitable partitions. This hypothesis may be restrictive to apply to
real-world networks. Here, we remove Assumption 2 and explore the question of
whether it is possible to achieve arbitrarily small partial phase cohesiveness
within clusters by tuning network parameters. To do so, we consider practical
stability under perturbation from an external equitable partition. In
particular, we are concerned with the stability concept defined below.

###### Definition 3

The cluster synchronization manifold $\mathcal{S}_{\Pi}$ of (2) is
_$(\delta,\rho)$ -practically stable_ if for given $\delta>0$ and $\rho>0$,
there exists $T>0$ such that
$\operatorname{dist}(\theta(0),\mathcal{S}_{\Pi})<\delta$ implies
$\operatorname{dist}(\theta(t),\mathcal{S}_{\Pi})\leq\rho$ for all $t\geq T$.

Similar definitions of practical stability under perturbations can be found in
[63, 64]. To proceed, we impose the following assumption instead of Assumption
2.

###### Assumption 3

There exists $K\geq 0$ such that for every $p,q\in\\{1,\ldots,r\\}$ with
$p\neq q$, it holds that $|\sum_{k\in\mathcal{C}_{q}}(a_{ik}-a_{jk})|\leq K$
for all $i,j\in\mathcal{C}_{p}$.

We remark that as long as the edge weights are finite, there always exists $K$
such that the condition in Assumption 3 is satisfied. Also, we note that
Assumption 2 corresponds to the case where $K=0$, and hence, $K$ represents
how large the variation from an external equitable partition is.

The following result states that the distance between the solutions of (2) and
the cluster synchronization manifold can be made arbitrarily close by
strengthening the connectivity within clusters.

###### Theorem 4

Let $\mathcal{S}_{\Pi}$ be the cluster synchronization manifold associated
with the Kuramoto model (2), where $\Pi$ is a nontrivial partition of
$\mathcal{V}$ such that the subgraphs $\mathcal{G}_{1},\ldots,\mathcal{G}_{r}$
are all connected. Suppose that Assumptions 1 and 3 hold. Then, for any
$\rho>0$, there exist $\delta>0$ and $a^{\prime}>0$ depending on $K$ such that
if

$\displaystyle a_{ij}>a^{\prime}\quad(i,j)\in\mathcal{E}_{p},\quad
p\in\\{1,\ldots,r\\},$

then $\mathcal{S}_{\Pi}$ is $(\delta,\rho)$-practically stable.

###### Proof:

Instead of (5), we consider the nonlinear system $\dot{x}=h(t,x)$, where $h$
is a function such that $h(t,x)=f(x,z(t))$ with $z(t)$ satisfying (6). This
can be considered as the perturbed system

$\displaystyle\dot{x}=h_{0}(t,x)+p(t,x),$ (21)

where $h_{0}(t,x)$ is a vector field associated with the external equitable
partition condition and $p(t,x)$ is a bounded perturbation term. In
particular, these are given as

$\displaystyle h_{0}(t,x)$
$\displaystyle:=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{intra}}W_{\mathrm{intra}}\sin(R_{1}x)$
$\displaystyle\quad{}-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W^{\prime}_{\mathrm{inter}}\sin(R_{3}x+R_{4}z(t)),$
$\displaystyle p(t,x)$
$\displaystyle:=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}(W_{\mathrm{inter}}-W^{\prime}_{\mathrm{inter}})$
$\displaystyle\quad{}\times\sin(R_{3}x+R_{4}z(t)),$

where $W^{\prime}_{\mathrm{inter}}$ is a certain weight matrix with respect to
inter-cluster edges for which the condition in Assumption 3 holds with $K=0$.
Also, for all $t\in\mathbb{R}_{+}$ and all $x\in\mathbb{R}^{n}$, we have

$\displaystyle\|p(t,x)\|\leq\|\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}\|\|W_{\mathrm{inter}}-W^{\prime}_{\mathrm{inter}}\|_{\infty}\leq
cK,$

where $c:=\|\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}\|$.

Assume that the zero solution $x^{\prime}(t)\equiv 0$ to the unperturbed
system $\dot{x}^{\prime}=h_{0}(t,x^{\prime})$ is exponentially stable, i.e.,
there exists $k>0$, $\lambda>0$, and $\delta^{\prime}>0$ such that
$\|x(t_{0})\|<\delta^{\prime}$ implies $\|x(t)\|\leq
k\mathrm{e}^{-\lambda(t-t_{0})}\|x(t_{0})\|$ for all $t\geq t_{0}$. We note
that $\lambda$ can be made arbitrarily large by strengthening the algebraic
connectivities of the subgraphs $\mathcal{G}_{1},\ldots,\mathcal{G}_{r}$. This
can be seen from the result in Section IV-A. According to the converse
Lyapunov theorem, the unperturbed system allows a Lyapunov function
$V\colon\mathbb{R}_{+}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$. Let
$c_{1},c_{2},c_{3},c_{4}>0$ be constants such that for all
$t\in\mathbb{R}_{+}$ and all $x\in\mathbb{R}^{n}$,

$\displaystyle c_{1}\|x\|^{2}\leq V(t,x)\leq c_{2}\|x\|^{2},$
$\displaystyle\frac{\partial V(t,x)}{\partial t}+\frac{\partial
V(t,x)}{\partial x}h_{0}(t,x)\leq-c_{3}\|x\|^{2},$
$\displaystyle\left\|\frac{\partial V(t,x)}{\partial x}\right\|\leq
c_{4}\|x\|$

whenever $\|x\|<d$ for some $d>0$. In what follows, we show that the solutions
to the perturbed system (21) converge to a certain small residual set.

Assume that $cK<(c_{3}/c_{4})\sqrt{c_{1}/c_{2}}d$. From Lemma 9.2 in [40], any
solution to (21) such that $\|x(0)\|<\sqrt{c_{1}/c_{2}}d$ satisfies the
following inequality with some $T>0$:

$\displaystyle\|x(t)\|<\frac{c_{4}}{c_{3}}\sqrt{\frac{c_{2}}{c_{1}}}cK,\quad
t\geq T.$ (22)

By choosing $c_{1},c_{2},c_{3},c_{4}$ as in Theorem 4.14 in [40], we can
observe that the larger $\lambda$ is, the smaller $c_{4}/c_{3}$ is, while
$c_{2}/c_{1}$ remains the same. Thus, by strengthening the algebraic
connectivities of $\mathcal{G}_{1},\ldots,\mathcal{G}_{r}$, we can make the
right-hand side in (22) small arbitrarily. For the same reason, the hypothesis
that $cK<(c_{3}/c_{4})\sqrt{c_{1}/c_{2}}d$ is allowed. Finally, we verify that
for any $\rho>0$, $(\delta,\rho)$-practical stability is realizable for some
$\delta>0$. Given $\rho>0$, we make $\lambda$ large such that
$(c_{4}/c_{3})\sqrt{c_{2}/c_{1}}cK=\rho$. Then, by defining
$\delta:=\sqrt{c_{1}/c_{2}}d$, we conclude that if $\|x(0)\|<\delta$, then
$x(t)$ eventually converges to the set
$\\{x\in\mathbb{R}^{n}:\|x\|\leq\rho\\}$. This immediately results in the
$(\delta,\rho)$-practical stability of the cluster synchronization manifold
$\mathcal{S}_{\Pi}$. Therefore, the proof is now complete. ∎

Theorem 4 indicates that even if Assumption 2 is not satisfied, clustered
patterns of synchronization can arise. Furthermore, we can design the network
parameters $a_{ij}$ to realize a desired partial phase cohesiveness within
clusters, which is characterized by $\rho$ in the above theorem. We note that
this theorem is helpful to proceed from analysis to control of synchronization
patterns.

###### Remark 7

Here, we discuss differences between our result and those in [25], where
partial phase cohesiveness within clusters is considered with a slightly
different setting. In this paper we consider the perturbation of the edge
weights from an external equitable partition. This can be seen as a
generalization of the case of completely separated clusters, i.e., the case
where $a_{ij}=0$ for all $(i,j)\in\mathcal{E}_{\mathrm{inter}}$. Therefore,
combining the result in this section with those in Section IV, we can conclude
that if the natural frequencies are largely different among clusters, the
inter-cluster couplings do not necessarily have to be small for realizing
partial synchronization patterns. This fact cannot be obtained from the
results in [25].

##  VI Applications to Brain Networks

This section focuses on applications of the cluster synchronization conditions
derived so far to brain networks. We first explain the construction of human
brain networks, followed by the exploration of synchronization patterns in
neural activities. However, neural activities of whole brains cannot be
measured directly, and fMRI measures blood flow-dependent signals. Thus, we
then consider hemodynamic responses [65] and the relation between
synchronization patterns of neural activities and functional connectivity
obtained from fMRI signals.

###  VI-A Construction of Brain Networks from Empirical Data

Mapping of human brain networks with noninvasive methods has been an active
research area in the neuroscience community. Here, we consider the brain
networks identified in [9], whose data including their weighted adjacency
matrices can be found at the USC Multimodal Connectivity Database [66]. In
[9], the authors used diffusion spectrum imaging to extract the network
structure with 998 cortical subregions of equal sizes for 5 human subjects.
Those subregions are classified into anatomically determined 66 cortical
regions.

Following [67, 68], we process the empirical data to construct a
representative network for computer simulation. We first binarize all the 5
weighted adjacency matrices and then average them into the 66 cortical
regions, that is, sum up all the weights in each region and divide it by the
number of the target regions. After that, we average the 66-region networks of
5 human subjects and normalize it so that the weights belong to the range
$[0,10]$. The resulting representative network is shown in Fig. 6, where the
image was generated by BrainNet Viewer [69] and only 30% of edges are
visualized. We also present in Fig. 7 the obtained adjacency matrix. Note that
this network is nonsymmetric but is almost symmetric.

Figure 6: Image of the structural brain network Figure 7: Adjacency matrix of
the structural brain network

###  VI-B Cluster Synchronization of Neural Activities

In general, functional connectivity is similar to the underlying structural
connectivity. However, they are not consistent and there are some regions
between which functional connectivity is strong but structural connectivity is
weak [67]. This fact can be explained by cluster synchronization phenomena as
explored in Section IV.

We briefly introduce a mathematical model for neural activities in brain
networks. Because of the complexity of human brains, the model is required to
represent dynamical behaviors at an appropriate spatial scale. Here, we are
interested in macroscopic behaviors of the populations of neurons. The
macroscopic dynamics of brains were investigated in [70, 71], and the so-
called Wilson–Cowan model was proposed. This model consists of the populations
of neurons with excitatory and inhibitory types, and can represent
oscillations in the phase space of the mean numbers of activated excitatory
and inhibitory neurons. Such a model was used to analyze functional
connectivity in [72, 67].

Here, we adopt the simpler Kuramoto model, which was employed in [68, 73]. The
derivation of a multiplexed Kuramoto model from the Wilson–Cowan model was
investigated in [74]. In functional connectivity, it is known that delays in
the communication among oscillators also play a certain role [72]. We,
however, ignore the time delays and focus on the effect of heterogeneity of
the natural frequencies. Thus, the model considered here is the one in (2).
Using the adjacency matrix obtained above, we simulate the phase evolution of
the Kuramoto oscillators with various choices of parameters.

We consider three clusters in the brain network, each consisting of 22 nodes
as represented in Figs. 6 and 7. To measure the synchronization level in each
cluster, we introduce the order parameters defined by
$r_{p}(t)=|(1/|\mathcal{C}_{p}|)\sum_{i\in\mathcal{C}_{p}}\mathrm{e}^{\,\mathrm{j}\theta_{i}(t)}|$
for $p\in\\{1,2,3\\}$. In geometric interpretations, the order parameters
represent the distance from the origin to the mean position of the oscillators
on the unit circle. That is, $r_{p}(t)=1$ when all the oscillators in the
cluster $\mathcal{C}_{p}$ are located at the same point, and $r_{p}(t)=0$ when
they are well spread with equal distances. We also define the global order
parameter to quantify the synchronization level among the whole network as
$r(t)=|(1/N)\sum_{i=1}^{N}\mathrm{e}^{\,\mathrm{j}\theta_{i}(t)}|$.
Furthermore, we calculate the Pearson correlation of the neural activities
among cortical regions, where the neural signal in the $i$th region is written
as $z_{i}(t)=\sin(\theta_{i}(t))$.

We simulate the neural activities with the Kuramoto model for $100$ seconds.
Then, after removing the data of the first $40$ seconds, we compute the
Pearson correlation coefficients among all pairs of regions. In the
simulations, we produce the oscillators’ frequencies $\omega_{i}/2\pi$ by
Gaussian distributions with mean $\mu$ and standard deviation $\sigma$. For
the first simulation, we set $(\mu,\sigma)=(60,0.5)$, and the order parameters
as well as the correlation matrix are presented in Fig. 8. One can observe in
the left figures that the synchronization levels of the individual clusters
(in the top plot) are high but the global synchronization level (in the bottom
plot) is also high. This indicates that specific correlation patterns of
neural activities do not appear as we can also confirm from the right figure.

Next, we demonstrate that the heterogeneity in the natural frequencies plays a
key role for the appearance of synchronization patterns. We produce the
natural frequencies in clusters $\mathcal{C}_{1}$, $\mathcal{C}_{2}$, and
$\mathcal{C}_{3}$ by Gaussian distributions with $(\mu,\sigma)=(50,0.5)$,
$(\mu,\sigma)=(60,0.5)$, and $(\mu,\sigma)=(70,0.5)$, respectively. The
simulation result is presented in Fig. 9, and it can be seen that the order
parameters within clusters are moderately high while the global order
parameter is averagely low. Hence, a correlation pattern can be observed. This
result exhibits that the heterogeneity in the natural frequencies is important
for the functional patterns observed in human brains and is also consistent
with the conclusions discussed in [75].

Finally, we aim to realize stronger synchronization levels within clusters.
From Theorem 4, this can be realized by strengthening the intra-cluster
coupling strengths. Thus, we double the intra-cluster coupling strengths while
the natural frequencies are kept the same as in the second simulation. The
obtained simulation result is shown in Fig. 10. In this case, the order
parameters within clusters are very high while the global order parameter is
averagely low.

Figure 8: Synchronization levels for natural frequencies under Gaussian
distributions with $(\mu,\sigma)=(60,0.5)$. Left: Order parameters for
individual clusters (top) and the whole system (bottom). Right: Correlation
matrix.

Figure 9: Synchronization levels for natural frequencies under Gaussian
distributions with $(\mu,\sigma)=(50,0.5)$ for cluster $\mathcal{C}_{1}$,
$(\mu,\sigma)=(60,0.5)$ for cluster $\mathcal{C}_{2}$, and
$(\mu,\sigma)=(70,0.5)$ for cluster $\mathcal{C}_{3}$. Left: Order parameters
for individual clusters (top) and the whole system (bottom). Right:
Correlation matrix.

Figure 10: Synchronization levels with doubled intra-cluster coupling
strengths. Left: Order parameters for individual clusters (top) and the whole
system (bottom). Right: Correlation matrix.

###  VI-C Hemodynamic Responses and Functional Connectivity

As already stated in the introduction, the fMRI techniques allow us to measure
the underlying brain activities. These are detected as slow fluctuations
(about $0.1$ Hz) of the blood oxygen level-dependent (BOLD) signals, which are
related to the blood volume and the deoxyhemoglobin content of the cortical
regions. This is based on the fact that the brain activity at a region
consumes oxygen in the blood flowing around that region. Functional
connectivity is then defined as the correlation level of the BOLD signals
among different cortical regions. In the last part of this section, we discuss
the relation of cluster synchronization with functional connectivity.

Following [65], we describe a mathematical model of hemodynamic responses. The
model consists of two parts. The first part represents a linear relation
between the neural activity and the regional cerebral blood flow (rCBF). In
particular, the neural activity $z_{i}$ causes the increase of the
vasodilatory signal $s_{i}$, which is integrated to obtain the blood inflow
$f_{i}$. The signal $s_{i}$ is regulated by the self feedback and by the blood
flow feedback. The second part links rCBF and fMRI BOLD signals with nonlinear
relation and is called the Balloon-Windkessel model. In this model, the
difference between the inflow $f_{i}$ and the outflow $f^{\mathrm{out}}_{i}$
results in the increase of the regional blood volume $v_{i}$. The
deoxyhemoglobin content $q_{i}$ changes depending on the blood flow delivery
and the concentration itself. The BOLD signal $y_{i}$ is then obtained,
depending nonlinearly on $v_{i}$ and $q_{i}$. In particular, the hemodynamic
model is described by

$\displaystyle\dot{s}_{i}$ $\displaystyle=z_{i}-\kappa s_{i}-\gamma(f_{i}-1),$
(23) $\displaystyle\dot{f}_{i}$ $\displaystyle=s_{i},$ (24)
$\displaystyle\tau\dot{v}_{i}$
$\displaystyle=f_{i}-f^{\mathrm{out}}_{i}(v_{i}),$ (25)
$\displaystyle\tau\dot{q}_{i}$
$\displaystyle=f_{i}\frac{E(f_{i})}{E_{0}}-f^{\mathrm{out}}_{i}\frac{q_{i}}{v_{i}},$
(26)

where $f^{\mathrm{out}}_{i}(v_{i})=v_{i}^{1/\alpha}$ and
$E(f_{i})=1-(1-E_{0})^{1/f_{i}}$. The BOLD signal is calculated as

$\displaystyle
y_{i}=V_{0}(k_{1}(1-q_{i})+k_{2}(1-q_{i}/v_{i})+k_{3}(1-v_{i})),$ (27)

where $V_{0}=0.02$, $k_{1}=7E_{0}$, $k_{2}=2$, and $k_{3}=2E_{0}-0.2$. We show
in Table I the meanings of the parameters appearing in (23)–(26) and their
values estimated in [65].

TABLE I: Parameters in the hemodynamic model Symbol | Meaning | Value
---|---|---
$\kappa$ | Decay rate of signals | 0.65
$\gamma$ | Elimination rate of signals | 0.41
$\tau$ | Time constant for blood flow | 0.98
$\alpha$ | Stiffness exponent | 0.33
$E_{0}$ | Resting oxygen extraction fraction | 0.34

It is known that a characteristic correlation pattern of the BOLD signals
appears in the low frequency range less than $0.1$ Hz [76]. This low-frequency
synchronization can be explained as follows. The state equations of the
hemodynamic model in (23)–(26) have an equilibrium point at
$[s_{i}\;f_{i}\;v_{i}\;q_{i}]^{\top}=[0\;1\;1\;1]^{\top}$. Since each signal
fluctuates around the equilibrium with small magnitudes, we linearize the
hemodynamic model. In what follows, we use $s_{i}$, $f_{i}$, $v_{i}$, and
$q_{i}$ to represent the deviation from this point with abuse of notation. The
linearized model is described by the series connection of the following two
systems:

$\displaystyle\begin{bmatrix}\dot{s}_{i}\\\ \dot{f}_{i}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}-\kappa&-\gamma\\\
1&0\end{bmatrix}\begin{bmatrix}s_{i}\\\ f_{i}\end{bmatrix}+\begin{bmatrix}1\\\
0\end{bmatrix}z_{i},$ $\displaystyle\eta_{i}$
$\displaystyle=\begin{bmatrix}0&1\end{bmatrix}\begin{bmatrix}s_{i}\\\
f_{i}\end{bmatrix}$

and

$\displaystyle\begin{bmatrix}\dot{v}_{i}\\\ \dot{q}_{i}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}-1/(\tau\alpha)&0\\\
-(1-\alpha)/(\tau\alpha)&-1/\tau\end{bmatrix}\begin{bmatrix}v_{i}\\\
q_{i}\end{bmatrix}+\begin{bmatrix}1/\tau\\\
\beta(E_{0})\end{bmatrix}\eta_{i},$ $\displaystyle y_{i}$
$\displaystyle=V_{0}\begin{bmatrix}k_{2}-k_{3}&-k_{1}-k_{2}\end{bmatrix}\begin{bmatrix}v_{i}\\\
q_{i}\end{bmatrix},$

where $\beta(E_{0}):=[E_{0}+(1-E_{0})\ln(1-E_{0})]/(\tau E_{0})$. Both are
stable systems, and the hemodynamic model can be interpreted as a low-pass
filter [68]. This implies that synchronization at relatively high frequency
can be detected as the correlation of the BOLD signals at low frequency. In
particular, by a frequency analysis, it can be observed that a neural signal
$z_{i}(t)$ with frequency around $60$ Hz is identified as a BOLD signal
$y_{i}(t)$ with peak frequency around $0.1$ Hz.

We compute the BOLD signals with the hemodynamic model for $100$ seconds. The
natural frequencies of the Kuramoto oscillators in each cluster are
respectively generated by Gaussian distributions with means $50$, $60$, and
$70$, and standard deviation $0.5$, which are the same setting as the second
simulation in the previous subsection. We then compute the Pearson correlation
coefficients after removing the data of the first $40$ seconds. We present in
Fig. 11 the obtained functional connectivity map. From this figure, we see
that a pattern of functional connectivity can be observed similarly to the
correlation pattern of the neural activities examined in the previous
subsection.

Figure 11: Correlation matrix for fMRI BOLD signals

We here give some discussions. This paper has been concerned with
synchronization patterns in oscillator networks. In particular, weak inter-
cluster connections and large heterogeneity in the natural frequencies among
clusters are key factors of cluster synchronization. These factors also play
certain roles in human brain networks as demonstrated above. Thus, it would be
of interest to analyze functional brain networks from these viewpoints and
moreover see if healthy and pathological states of brains may be characterized
based on this criterion. Compared with the graph-theoretic analysis [12] and
numerical studies [68] in the neuroscience community, this paper has provided
a theoretical analysis based on dynamical models of oscillators to explain the
robust patterns in functional connectivity of brains. There are some possible
future directions of this work. One is the brain intervention for cluster
synchronization as considered in [16]. Furthermore, the analysis of remote
synchronization is important since it is known that some brain regions located
remotely can be synchronized [77].

##  VII Concluding Remarks

In this paper we have considered cluster synchronization of Kuramoto
oscillators. Through a stability criterion based on averaging, we have shown
that cluster synchronization is achieved if (i) the inter-cluster couplings
are sufficiently weak and/or (ii) the natural frequencies are sufficiently
different among clusters. Practical stability has been investigated to remove
the hypothesis on network partitions. We have also examined applications to
empirical brain networks and have discussed the relation between the cluster
synchronization problem for the Kuramoto model and the problem of functional
connectivity in brain networks.

###  -A Proof of Theorem 1

The average system associated with (8) is given by

$\displaystyle\dot{\xi}=\varepsilon A_{\mathrm{av}}\xi,$ (28)

where $A_{\mathrm{av}}$ is defined in (9). Since $A_{\mathrm{av}}$ is a
Hurwitz matrix, there exists a positive-definite matrix
$P\in\mathbb{R}^{n\times n}$ such that
$A_{\mathrm{av}}^{\top}P+PA_{\mathrm{av}}\prec 0$. We define the positive-
definite function $V(x):=x^{\top}Px$. This is clearly a Lyapunov function for
(28), and thus, there exist $c_{1},c_{2},c_{3}>0$ such that for all
$x\in\mathbb{R}^{n}$,

$\displaystyle c_{1}\|x\|^{2}\leq V(x)\leq c_{2}\|x\|^{2},$ (29)
$\displaystyle\frac{\partial V(x)}{\partial x}A_{\mathrm{av}}x\leq-
c_{3}\|x\|^{2}.$ (30)

For a fixed $t$, we consider the time interval $[t,t+T]$, where $T>0$ is
determined later. Let $\Phi(t+T,t)$ denote the state transition matrix of (8),
i.e., $x(t+T)=\Phi(t+T,t)x(t)$. The Peano–Baker series of $\Phi(t+T,t)$ is
given by

$\displaystyle\Phi(t+T,t)$
$\displaystyle=I+\varepsilon\int_{t}^{t+T}A(s_{1},\varepsilon)\,\mathrm{d}s_{1}$
$\displaystyle\quad{}+\varepsilon^{2}\int_{t}^{t+T}A(s_{1},\varepsilon)\int_{t}^{s_{1}}A(s_{2},\varepsilon)\,\mathrm{d}s_{2}\mathrm{d}s_{1}+\cdots$
$\displaystyle=I+\varepsilon\int_{t}^{t+T}A(s_{1},\varepsilon)\,\mathrm{d}s_{1}$
$\displaystyle\quad{}+\sum_{i=2}^{\infty}\varepsilon^{i}\int_{t}^{t+T}A(s_{1},\varepsilon)\cdots\int_{t}^{s_{i-1}}A(s_{i},\varepsilon)\,\mathrm{d}s_{i}\cdots\mathrm{d}s_{1}.$
(31)

Let $\Phi_{\mathrm{av}}(t+T,t)$ denote the state transition matrix of (28),
i.e., $\xi(t+T)=\Phi_{\mathrm{av}}(t+T,t)\xi(t)$. This can be calculated as

$\displaystyle\Phi_{\mathrm{av}}(t+T,t)$
$\displaystyle=\mathrm{e}^{\varepsilon TA_{\mathrm{av}}}$
$\displaystyle=I+\varepsilon
TA_{\mathrm{av}}+\sum_{i=2}^{\infty}\frac{(\varepsilon
T)^{i}A_{\mathrm{av}}^{i}}{i!}.$ (32)

We then consider the difference
$H(t+T,t):=\Phi(t+T,t)-\Phi_{\mathrm{av}}(t+T,t)$ between these two transition
matrices. From (31) and (32), we have

$\displaystyle H(t+T,t)$
$\displaystyle=\varepsilon\int_{t}^{t+T}[A(s,\varepsilon)-A_{\mathrm{av}}]\,\mathrm{d}s$
$\displaystyle\quad{}+\sum_{i=2}^{\infty}\varepsilon^{i}\int_{t}^{t+T}A(s_{1},\varepsilon)\cdots\int_{t}^{s_{i-1}}A(s_{i},\varepsilon)\,\mathrm{d}s_{i}\cdots\mathrm{d}s_{1}$
$\displaystyle\quad{}-\sum_{i=2}^{\infty}\frac{(\varepsilon
T)^{i}A_{\mathrm{av}}^{i}}{i!}.$

It follows that

$\displaystyle\|H(t+T,t)\|$
$\displaystyle=\varepsilon\left\|\int_{t}^{t+T}[A(s,\varepsilon)-A_{\mathrm{av}}]\,\mathrm{d}s\right\|$
$\displaystyle\quad{}+\sum_{i=2}^{\infty}\varepsilon^{i}\int_{t}^{t+T}\|A(s_{1},\varepsilon)\|\cdots\int_{t}^{s_{i-1}}\|A(s_{i},\varepsilon)\|\,\mathrm{d}s_{i}\cdots\mathrm{d}s_{1}$
$\displaystyle\quad{}-\sum_{i=2}^{\infty}\frac{(\varepsilon
T)^{i}\|A_{\mathrm{av}}\|^{i}}{i!}$ $\displaystyle\leq\varepsilon
T[\beta(T)+\alpha(\varepsilon)]+2(\mathrm{e}^{\rho\varepsilon
T}-1-\rho\varepsilon T),$

where $\rho:=\sup_{t\in\mathbb{R}_{+},\varepsilon>0}\|A(t,\varepsilon)\|$ and
we have used the inequality (10) in the last line. Let us define

$\displaystyle h(\varepsilon,T):=\varepsilon
T[\beta(T)+\alpha(\varepsilon)]+2(\mathrm{e}^{\rho\varepsilon
T}-1-\rho\varepsilon T).$ (33)

Again, we consider the Lyapunov function candidate $V(x)$. Along the solutions
to (8), we have

$\displaystyle V(x(t+T))-V(x(t))$
$\displaystyle=x^{\top}(t)[\Phi^{\top}(t+T,t)P\Phi(t+T,t)-P]x(t)$
$\displaystyle=x^{\top}(t)[\Phi_{\mathrm{av}}^{\top}(t+T,t)P\Phi_{\mathrm{av}}(t+T,t)-P$
$\displaystyle\quad{}+\Phi_{\mathrm{av}}^{\top}(t+T,t)PH(t+T,t)$
$\displaystyle\quad{}+H^{\top}(t+T,t)P\Phi_{\mathrm{av}}(t+T,t)$
$\displaystyle\quad{}+H^{\top}(t+T,t)PH(t+T,t)]x(t).$

For the first two terms in the far right-hand side, it can be guaranteed from
(29) and (30) that

$\displaystyle
x^{\top}(t)[\Phi_{\mathrm{av}}^{\top}(t+T,t)P\Phi_{\mathrm{av}}(t+T,t)-P]x(t)$
$\displaystyle\leq c_{2}(\mathrm{e}^{-\frac{c_{3}}{c_{2}}\varepsilon
T}-1)\|x(t)\|^{2}.$

Note that $\|P\|\leq c_{2}$ and
$\|\Phi_{\mathrm{av}}(t+T,t)\|\leq\sqrt{c_{2}/c_{1}}\mathrm{e}^{-\frac{c_{3}}{2c_{2}}\varepsilon
T}$. Hence, by using $h(\varepsilon,T)$ defined in (33), we obtain

$\displaystyle V(x(t+T))-V(x(t))$ $\displaystyle\leq
c_{2}\left[(\mathrm{e}^{-\frac{c_{3}}{c_{2}}\varepsilon
T}-1)+2\sqrt{\frac{c_{2}}{c_{1}}}\mathrm{e}^{-\frac{c_{3}}{2c_{2}}\varepsilon
T}h(\varepsilon,T)+h^{2}(\varepsilon,T)\right]$
$\displaystyle\quad{}\times\|x(t)\|^{2}$
$\displaystyle=:\tilde{h}(\varepsilon,T)\|x(t)\|^{2}.$ (34)

It can be verified that $\tilde{h}(0,T)=0$. Moreover, we have

$\displaystyle\limsup_{\varepsilon\to
0^{+}}\frac{\tilde{h}(\varepsilon,T)-\tilde{h}(0,T)}{\varepsilon}=-T\left(c_{3}-2c_{2}\sqrt{\frac{c_{2}}{c_{1}}}\beta(T)\right).$

Because $\beta$ is a strictly decreasing function that converges to zero, the
positive constant $T$ can be chosen such that
$c_{3}-2c_{2}\sqrt{c_{2}/c_{1}}\beta(T)>0$. This implies that for sufficiently
small $\varepsilon>0$, we have $\tilde{h}(\varepsilon,T)<0$. It follows from
(34) that $V(x(t+T))-V(x(t))<0$ whenever $x(t)\neq 0$.

Notice that the above procedure can be performed for every
$t\in\mathbb{R}_{+}$. Hence, there always exists a time sequence
$\\{t_{k}\\}_{k\in\mathbb{Z}_{+}}$ satisfying the hypotheses in Proposition 2.
Therefore, we conclude that there exists $\varepsilon^{*}>0$ such that if
$\varepsilon<\varepsilon^{*}$, the zero solution $x(t)\equiv 0$ of (8) is
exponentially stable. ∎

###  -B Proof of Lemma 2

From (14) and (16) we have

$\displaystyle\int_{\tau}^{\tau+T}[J(s,\varepsilon)-J_{\mathrm{av}}]\,\mathrm{d}s$
$\displaystyle=-\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\int_{\tau}^{\tau+T}\operatorname{diag}[R_{4}\cos(\zeta(s,\varepsilon))]\,\mathrm{d}sR_{3}.$

Because each row of $R_{4}$ contains single $1$ or $-1$ and other elements are
zeros, it holds that $\|R_{4}x\|_{\infty}=\|x\|_{\infty}$ for all
$x\in\mathbb{R}^{m}$. It can be verified that

$\displaystyle\left\|\int_{\tau}^{\tau+T}[J(s,\varepsilon)-J_{\mathrm{av}}]\,\mathrm{d}s\right\|$
$\displaystyle\leq\|\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\|\|R_{3}\|\left\|\int_{\tau}^{\tau+T}\cos(\zeta(s,\varepsilon))\,\mathrm{d}s\right\|_{\infty},$

where we have used the fact that the spectrum norm of the diagonal matrix in
the integral part is equal to its induced $\infty$-norm. Note that the
solution $\zeta$ to (13) can be written as

$\displaystyle\zeta(\tau,\varepsilon)$ $\displaystyle=z_{0}+\eta\tau$
$\displaystyle\quad{}-\varepsilon\tilde{B}_{\mathrm{inter}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}R_{4}\int_{0}^{\tau}\sin(\zeta(s,\varepsilon))\,\mathrm{d}s,$

where we note that $R_{4}$ can be moved outside of the argument of
$\sin(\cdot)$ by the same reason as mentioned above. For notational
simplicity, let
$\Psi:=\tilde{B}_{\mathrm{inter}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}R_{4}$.
Then, we consider

$\displaystyle\mathrm{e}^{\,\mathrm{j}\zeta(s,\varepsilon)}=\mathrm{e}^{\,\mathrm{j}(z_{0}+\eta
s)}\odot\mathrm{e}^{-\mathrm{j}\varepsilon\Psi\int_{0}^{s}\sin(\zeta(r,\varepsilon))\,\mathrm{d}r},$

where the exponential function of a vector is to be taken elementwise such as
$\mathrm{e}^{x}=[\mathrm{e}^{x_{1}}\,\cdots\,\mathrm{e}^{x_{n}}]^{\top}$ and
$\odot$ denotes the Hadamard (elementwise) product. By defining
$\psi(s):=[\operatorname{diag}(\eta)]^{-1}\mathrm{e}^{\,\mathrm{j}(z_{0}+\eta
s)}$, we have

$\displaystyle\int_{\tau}^{\tau+T}\mathrm{e}^{\,\mathrm{j}(z_{0}+\eta
s)}\odot\mathrm{e}^{-\mathrm{j}\varepsilon\Psi\int_{0}^{s}\sin(\zeta(r,\varepsilon))\,\mathrm{d}r}\,\mathrm{d}s$
$\displaystyle=\left[-\mathrm{j}\psi(s)\odot\mathrm{e}^{-\mathrm{j}\varepsilon\Psi\int_{0}^{s}\sin(\zeta(r,\varepsilon))\,\mathrm{d}r}\right]_{\tau}^{\tau+T}$
$\displaystyle\quad{}+\varepsilon\int_{\tau}^{\tau+T}\psi(s)\odot\Psi\sin(\zeta(s,\varepsilon))\odot\mathrm{e}^{-\mathrm{j}\varepsilon\Psi\int_{0}^{s}\sin(\zeta(r,\varepsilon))\,\mathrm{d}r}\,\mathrm{d}s.$

Noting that $\|\psi(s)\|_{\infty}\leq 1$, we obtain

$\displaystyle\left\|\int_{\tau}^{\tau+T}\mathrm{e}^{\,\mathrm{j}\zeta(s,\varepsilon)}\,\mathrm{d}s\right\|_{\infty}\leq
2+\varepsilon T\|\Psi\|_{\infty}.$

Similarly, we have

$\displaystyle\left\|\int_{\tau}^{\tau+T}\mathrm{e}^{-\mathrm{j}\zeta(s,\varepsilon)}\,\mathrm{d}s\right\|_{\infty}\leq
2+\varepsilon T\|\Psi\|_{\infty}.$

Combining the above two inequalities with the equality
$\cos(\zeta(s,\varepsilon))=[\mathrm{e}^{\,\mathrm{j}\zeta(s,\varepsilon)}+\mathrm{e}^{-\mathrm{j}\zeta(s,\varepsilon)}]/2$
yields

$\displaystyle\left\|\frac{1}{T}\int_{\tau}^{\tau+T}\cos(\zeta(s,\varepsilon))\,\mathrm{d}s\right\|_{\infty}\leq\frac{2}{T}+\varepsilon\|\Psi\|_{\infty}.$

Therefore, we obtain

$\displaystyle\left\|\frac{1}{T}\int_{\tau}^{\tau+T}[J(s,\varepsilon)-J_{\mathrm{av}}]\,\mathrm{d}s\right\|$
$\displaystyle\leq\|\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\|\|R_{3}\|\left(\frac{2}{T}+\varepsilon\|\Psi\|_{\infty}\right)$
$\displaystyle\leq\gamma\left(\frac{1}{T}+\varepsilon\right),$

where
$\gamma:=\|\tilde{B}_{\mathrm{intra}}^{\top}B_{\mathrm{inter}}W_{\mathrm{inter}}\|\|R_{3}\|\max\\{2,\|\Psi\|_{\infty}\\}$.
It is clear that $\gamma=O(\|W_{\mathrm{inter}}\|_{\infty})$. The proof is now
complete. ∎

## References

  * [1] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou.  Synchronization in complex networks.  Phys. Rep., 469(3):93–153, 2008.
  * [2] F. Dörfler and F. Bullo.  Synchronization in complex networks of phase oscillators: A survey.  Automatica, 50(6):1539–1564, 2014.
  * [3] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa.  Spontaneous synchrony in power-grid networks.  Nat. Phys., 9(3):191–197, 2013.
  * [4] B. A. Francis and M. Maggiore.  Flocking and Rendezvous in Distributed Robotics.  Springer, 2016.
  * [5] L. Glass.  Synchronization and rhythmic processes in physiology.  Nature, 410:277–284, 2001.
  * [6] G. Buzsáki and A. Draguhn.  Neuronal oscillations in cortical networks.  Science, 304(5679):1926–1929, 2004.
  * [7] C. J. Honey, R. Kötter, M. Breakspear, and O. Sporns.  Network structure of cerebral cortex shapes functional connectivity on multiple time scales.  Proc. Natl. Acad. Sci. USA, 104(24):10240–10245, 2007.
  * [8] Y. Iturria-Medina, R. C. Sotero, E. J. Canales-Rodríguez, Y. Alemán-Gómez, and L. Melie-Garcıa.  Studying the human brain anatomical network via diffusion-weighted MRI and graph theory.  NeuroImage, 40(3):1064–1076, 2008.
  * [9] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey, V. J. Wedeen, and O. Sporns.  Mapping the structural core of human cerebral cortex.  PLoS Biol., 6(7):e159, 2008.
  * [10] M. H. Lee, C. D. Smyser, and J. S. Shimony.  Resting-state fMRI: A review of methods and clinical applications.  Am. J. Neuroradiol., 34(10):1866–1872, 2013.
  * [11] M. Fox and M. Raichle.  Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging.  Nat. Rev. Neurosci., 8:700–711, 2007.
  * [12] E. Bullmore and O Sporns.  Complex brain networks: Graph theoretical analysis of structural and functional systems.  Nat. Rev. Neurosci., 10(3):186–198, 2009.
  * [13] C. J. Stam, B. F. Jones, G. Nolte, M. Breakspear, and P. Scheltens.  Small-world networks and functional connectivity in Alzheimer’s disease.  Cereb. Cortex, 17(1):92–99, 2007.
  * [14] M.-E. Lynall, D. S. Bassett, R. Kerwin, P. J. McKenna, M. Kitzbichler, U. Muller, and E. Bullmore.  Functional connectivity and brain networks in schizophrenia.  J. Neurosci., 30(28):9477–9487, 2010.
  * [15] S.C. Ponten, F. Bartolomei, and C. J. Stam.  Small-world networks and epilepsy: Graph theoretical analysis of intracerebrally recorded mesial temporal lobe seizures.  Clin. Neurophysiol., 118(4):918–927, 2007.
  * [16] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti.  A framework to control functional connectivity in the human brain.  In Proc. IEEE Conf. on Decision and Control, pages 4697–4704, 2019\.
  * [17] S. F. Muldoon, F. Pasqualetti, S. Gu, M. Cieslak, S. T. Grafton, J. M. Vettel, and D. S. Bassett.  Stimulation-based control of dynamic brain networks.  PLoS Comput. Biol., 12(9):1–23, 2016.
  * [18] S. F. Cogan.  Neural stimulation and recording electrodes.  Annu. Rev. Biomed. Eng., 10(1):275–309, 2008.
  * [19] O. Yizhar, L. E. Fenno, T. J. Davidson, M. Mogri, and K. Deisseroth.  Optogenetics in neural systems.  Neuron, 71(1):9–34, 2011.
  * [20] E. Tang and D. S. Bassett.  Colloquium: Control of dynamics in brain networks.  Rev. Mod. Phys., 90(3):031003, 2018.
  * [21] Y. Kuramoto.  Chemical Oscillations, Waves, and Turbulence.  Springer, 1984.
  * [22] A. Jadbabaie, N. Motee, and M. Barahona.  On the stability of the Kuramoto model of coupled nonlinear oscillators.  In Proc. American Control Conf., pages 4296–4301, 2004.
  * [23] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti.  Exact and approximate stability conditions for cluster synchronization of Kuramoto oscillators.  In Proc. American Control Conf., pages 205–210, 2019.
  * [24] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti.  Stability conditions for cluster synchronization in networks of heterogeneous Kuramoto oscillators.  IEEE Trans. Control Netw. Syst., 7(1):302–314, 2020.
  * [25] Y. Qin, Y. Kawano, O. Portoles, and M. Cao.  Partial phase cohesiveness in networks of networks of Kuramoto oscillators.  IEEE Trans. Autom. Control, 2021, to appear.
  * [26] A. Pogromsky, G. Santoboni, and H. Nijmeijer.  Partial synchronization: From symmetry towards stability.  Physica D, 172(1–4):65–87, 2002.
  * [27] I. Stewart, M. Golubitsky, and M. Pivato.  Symmetry groupoids and patterns of synchrony in coupled cell networks.  SIAM J. Appl. Dyn. Syst., 2(4):609–646, 2003.
  * [28] V. N. Belykh, G. V. Osipov, V. S. Petrov, J. A. K. Suykens, and J. Vandewalle.  Cluster synchronization in oscillatory networks.  Chaos, 18(3):037106, 2008.
  * [29] W. Lu, B. Liu, and T. Chen.  Cluster synchronization in networks of coupled nonidentical dynamical systems.  Chaos, 20(1):013120, 2010.
  * [30] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy.  Cluster synchronization and isolated desynchronization in complex networks with symmetries.  Nat. Commun., 5:4079, 2014.
  * [31] F. Sorrentino, L. M. Pecora, A. M. Hagerstrom, T. E. Murphy, and R. Roy.  Complete characterization of the stability of cluster synchronization in complex dynamical networks.  Sci. Adv., 2(4):e1501737, 2016.
  * [32] M. T. Schaub, N. O’Clery, Y. N. Billeh, J.-C. Delvenne, R. Lambiotte, and M. Barahona.  Graph partitions and cluster synchronization in networks of oscillators.  Chaos, 26(9):094821, 2016.
  * [33] L. V. Gambuzza and M. Frasca.  A criterion for stability of cluster synchronization in networks with external equitable partitions.  Automatica, 100:212–218, 2019.
  * [34] Q.-C. Pham and J.-J. E. Slotine.  Stable concurrent synchronization in dynamic system networks.  Neural Netw., 20(1):62–77, 2007.
  * [35] Z. Aminzare, B. Dey, E. N. Davison, and N. E. Leonard.  Cluster synchronization of diffusively coupled nonlinear systems: A contraction-based approach.  J. Nonlinear Sci., 2018.
  * [36] W. Xia and M. Cao.  Clustering in diffusively coupled networks.  Automatica, 47(11):2395–2405, 2011.
  * [37] W. Wu, W. Zhou, and T. Chen.  Cluster synchronization of linearly coupled complex networks under pinning control.  IEEE Trans. Circuits Syst. I, Reg. Papers, 56(4):829–839, 2009\.
  * [38] X. B. Lu and B. Z. Qin.  Adaptive cluster synchronization in complex dynamical networks.  Phys. Lett. A, 373(40):3650–3658, 2009.
  * [39] X. Liu and T. Chen.  Cluster synchronization in directed networks via intermittent pinning control.  IEEE Trans. Neural Netw., 22(7):1009–1020, 2011.
  * [40] H. K. Khalil.  Nonlinear Systems.  Prentice Hall, 3rd edition, 2002.
  * [41] R. Kato and H. Ishii.  Averaging and cluster synchronization of Kuramoto oscillators.  In Proc. European Control Conf., 2021.
  * [42] L. Tiberi, C. Favaretto, M. Innocenti, D. S. Bassett, and F. Pasqualetti.  Synchronization patterns in networks of Kuramoto oscillators: A geometric approach for analysis and control.  In Proc. IEEE Conf. on Decision and Control, pages 481–486, 2017\.
  * [43] D. M. Cardoso, C. Delorme, and P. Rama.  Laplacian eigenvectors and eigenvalues and almost equitable partitions.  Eur. J. Comb., 28(3):665–673, 2007.
  * [44] N. O’Clery, Y. Yuan, G.-B. Stan, and M. Barahona.  Observability and coarse graining of consensus dynamics through the external equitable partition.  Phys. Rev. E, 88(4):042805, 2013.
  * [45] S. Zhang, M. Cao, and M. K. Camlibel.  Upper and lower bounds for controllable subspaces of networks of diffusively coupled agents.  IEEE Trans. Autom. Control, 59(3):745–750, 2014.
  * [46] J. K. Baksalary and O. M. Baksalary.  Particular formulae for the Moore–Penrose inverse of a columnwise partitioned matrix.  Linear Algebra Its Appl., 421(1):16–23, 2007.
  * [47] W. M. Haddad and V. S. Chellaboina.  Nonlinear Dynamical Systems and Control: A Lyapunov-Based Approach.  Princeton University Press, 2008.
  * [48] A. N. Michel, L. Hou, and D. Liu.  Stability of Dynamical Systems.  Birkhäuser, 2015.
  * [49] D. Aeyels and J. Peuteman.  A new asymptotic stability criterion for nonlinear time-variant differential equations.  IEEE Trans. Autom. Control, 43(7):968–971, 1998.
  * [50] D. Aeyels and J. Peuteman.  On exponential stability of nonlinear time-varying differential equations.  Automatica, 35(6):1091–1100, 1999.
  * [51] R. L. Kosut, B. D. O. Anderson, and I. M. Y. Mareels.  Stability theory for adaptive systems: Method of averaging and persistency of excitation.  IEEE Trans. Autom. Control, 32(1):26–34, 1987.
  * [52] Z. Lin, B. Francis, and M. Maggiore.  State agreement for continuous-time coupled nonlinear systems.  SIAM J. Control Optim., 46(1):288–307, 2007.
  * [53] J. K. Hale.  Ordinary Differential Equations.  Robert E. Krieger, 2nd edition, 1980.
  * [54] Y. Qin, Y. Kawano, B. D. O. Anderson, and M. Cao.  Partial exponential stability analysis of slow-fast systems via periodic averaging.  IEEE Trans. Autom. Control, 2021, to appear.
  * [55] D. J. Stilwell, E. M. Bollt, and D. G. Roberson.  Sufficient conditions for fast switching synchronization in time-varying network topologies.  SIAM J. Appl. Dyn. Syst., 5(1):140–156, 2006.
  * [56] M. Arcak.  Passivity as a design tool for group coordination.  IEEE Trans. Autom. Control, 52(8):1380–1390, 2007.
  * [57] D. M. Abrams and S. H. Strogatz.  Chimera states for coupled oscillators.  Phys. Rev. Lett., 93(17):174102, 2004.
  * [58] Y. S. Cho, T. Nishikawa, and A. E. Motter.  Stable chimeras and independently synchronizable clusters.  Phys. Rev. Lett., 119(8):084101, 2017.
  * [59] W. J. Rugh.  Linear System Theory.  Prentice Hall, 2nd edition, 1996.
  * [60] C. Altafini.  Consensus problems on networks with antagonistic interactions.  IEEE Trans. Autom. Control, 58(4):935–946, 2013.
  * [61] A. V. Proskurnikov, A. S. Matveev, and M. Cao.  Opinion dynamics in social networks with hostile camps: Consensus vs. polarization.  IEEE Trans. Autom. Control, 61(6):1524–1536, 2016.
  * [62] R. Delabays, P. Jacquod, and F. Dörfler.  The Kuramoto model on oriented and signed graphs.  SIAM J. Appl. Dyn. Syst., 18(1):458–480, 2019.
  * [63] J. M. Montenbruck, M. Bürger, and F. Allgöwer.  Practical synchronization with diffusive couplings.  Automatica, 53:235–243, 2015.
  * [64] C. Lageman and Z. Sun.  Consensus on spheres: Convergence analysis and perturbation theory.  In Proc. IEEE Conf. on Decision and Control, pages 19–24, 2016\.
  * [65] K. J. Friston, A. Mechelli, R. Turner, and C. J. Price.  Nonlinear responses in fMRI: The Balloon model, Volterra kernels, and other hemodynamics.  NeuroImage, 12(4):466–477, 2000.
  * [66] J. Brown, J. Rudie, A. Bandrowski, J. VanHorn, and S. Bookheimer.  The UCLA multimodal connectivity database: A web-based platform for brain connectivity matrix sharing and analysis.  Front. Neuroinform., 6:28, 2012.
  * [67] C. J. Honey, O. Sporns, L. Cammoun, X. Gigandet, J. P. Thiran, R. Meuli, and P. Hagmann.  Predicting human resting-state functional connectivity from structural connectivity.  Proc. Natl. Acad. Sci. USA, 106(6):2035–2040, 2009.
  * [68] J. Cabral, E. Hugues, O. Sporns, and G. Deco.  Role of local network oscillations in resting-state functional connectivity.  NeuroImage, 57(1):130–139, 2011.
  * [69] M. Xia, J. Wang, and Y. He.  BrainNet Viewer: A network visualization tool for human brain connectomics.  PLoS One, 8(7):e68910, 2013.
  * [70] S.-I. Amari.  Characteristics of randomly connected threshold-element networks and network systems.  Proc. IEEE, 59(1):35–47, 1971.
  * [71] H. R. Wilson and J. D. Cowan.  Excitatory and inhibitory interactions in localized populations of model neurons.  Biophys. J., 12(1):1–24, 1972.
  * [72] G. Deco, V. Jirsa, A. R. McIntosh, O. Sporns, and R. Kötter.  Key role of coupling, delay, and noise in resting brain fluctuations.  Proc. Natl. Acad. Sci. USA, 106(25):10302–10307, 2009.
  * [73] W. H. Lee and S. Frangou.  Linking functional connectivity and dynamic properties of resting-state networks.  Sci. Rep., 7:16610, 2017.
  * [74] M. Sadilek and S. Thurner.  Physiologically motivated multiplex kuramoto model describes phase diagram of cortical activity.  Sci. Rep., 5:10015, 2015.
  * [75] A. Ponce-Alvarez, G. Deco, P. Hagmann, G. L. Romani, D. Mantini, and M. Corbetta.  Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity.  PLoS Comput. Biol., 11(2):1–23, 2015.
  * [76] B. Biswal, F. Z. Yetkin, V. M. Haughton, and J. S. Hyde.  Functional connectivity in the motor cortex of resting human brain using echo-planar MRI.  Magn. Reson. Med., 34(4):537–541, 1995.
  * [77] V. Nicosia, M. Valencia, M. Chavez, A. Díaz-Guilera, and V. Latora.  Remote synchronization reveals network symmetries and functional modules.  Phys. Rev. Lett., 110(17):174102, 2013.