A Framework to Control Functional Connectivity in the Human Brain

TL;DR

This paper introduces a control framework for brain functional connectivity using interconnected oscillators, enabling the regulation of neural synchronization patterns to restore healthy brain dynamics.

Contribution

It presents a novel, minimally invasive method to manipulate neural oscillator parameters for desired connectivity patterns, validated through realistic brain models.

Findings

Framework effectively controls synchronization patterns

Robust to parameter mismatches and inaccuracies

Validated with realistic neurovascular simulations

Abstract

In this paper, we propose a framework to control brain-wide functional connectivity by selectively acting on the brain's structure and parameters. Functional connectivity, which measures the degree of correlation between neural activities in different brain regions, can be used to distinguish between healthy and certain diseased brain dynamics and, possibly, as a control parameter to restore healthy functions. In this work, we use a collection of interconnected Kuramoto oscillators to model oscillatory neural activity, and show that functional connectivity is essentially regulated by the degree of synchronization between different clusters of oscillators. Then, we propose a minimally invasive method to correct the oscillators' interconnections and frequencies to enforce arbitrary and stable synchronization patterns among the oscillators and, consequently, a desired pattern of functional connectivity. Additionally, we show that our synchronization-based framework is robust to parameter mismatches and numerical inaccuracies, and validate it using a realistic neurovascular model to simulate neural activity and functional connectivity in the human brain.

Equations

$${\dot{\theta }}_i={\omega }_i+\sum _{j\ne i}{a}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$${\mathcal{S}}_{\mathcal{P}}=\{\theta \in {\mathbb{T}}^n:{\theta }_i={\theta }_j\text{ for all }i,j\in {\mathcal{P}}_k,k=1,\dots ,m\}.$$

$$v_i^{\mathsf{T}}=[\underset{{\sum }_{j=1}^{i-1}|{\mathcal{P}}_j|}{\underbrace{0,\dots ,0}},\underset{|{\mathcal{P}}_i|}{\underbrace{1,\dots ,1}},\underset{{\sum }_{j=1+1}^n|{\mathcal{P}}_j|}{\underbrace{0,\dots ,0}}],$$

$${\xi }_{k\ell }=\{\begin{array}{ll}{\nu }_{k\ell }\overline{\sigma }(G_k(i{\omega }^{(k\ell )})),& \text{if }\overline{\lambda }(J_k)\le \overline{\lambda }(J_{\ell }),\\ {\nu }_{k\ell }\frac{\overline{\sigma }(G_k(0))}{\overline{\sigma }(G_{\ell }(0))}\overline{\sigma }(G_{\ell }(i{\omega }^{(k\ell )})),& \text{if }\overline{\lambda }(J_{\ell })<\overline{\lambda }(J_k),\end{array}$$

$$\displaystyle\min\limits_{\Delta,\mu}\$$

$${\bar{V}}_{\mathcal{P}}^{\mathsf{T}}(\bar{A}+\Delta )V_{\mathcal{P}}=0,$$

$$B_{\text{span}}^{\mathsf{T}}(\omega +\mu )=0,$$

$$H^{\text{c}}\circ \Delta =0$$

$$\bar{A}+\Delta \ge 0,$$

$$\omega +\mu \ge 0,$$

$${\lambda }_{\text{max}}(\Xi (A+\Delta ,\omega +\mu ))<1,$$

$$\displaystyle\min\limits_{\Delta}\$$

$${\bar{V}}_{\mathcal{P}}^{\mathsf{T}}(\bar{A}+\Delta )V_{\mathcal{P}}=0,$$

$$H^{\text{c}}\circ \Delta =0,$$

$$\bar{A}+\Delta \ge 0.$$

$$Y_k$$

$$T_{k+1}$$

$$Z_{k+1}$$

$$Q_{k+1}$$

$${\Pi }_{{\mathcal{Z}}_2}(W)$$

$$=W-{\bar{V}}_{\mathcal{P}}{\bar{V}}_{\mathcal{P}}^{\mathsf{T}}W{V}_{\mathcal{P}}{V}_{\mathcal{P}}^{\mathsf{T}},$$

$$Z^{\ast }$$

$$=\bar{A}+\arg \underset{\begin{array}{c}\Delta \text{ s.t. }\text{eqref}eq:controlinvariancea,\text{eqref}eq:controlinvarianceb,\text{eqref}eq:controlinvariancec\end{array}}{\min }\Vert \Delta {\Vert }_{\text{F}}^2,$$

$$\sum _{k\in {\mathcal{P}}_{\ell }}{a}_{ik}-a_{jk}=0$$

$$\displaystyle\min\limits_{\mu}\$$

$$B_{\text{span}}^{\mathsf{T}}(\omega +\mu )=0,$$

$$\omega +\mu \ge 0,$$

$${\lambda }_{\text{max}}(\Xi (A,\omega +\mu ))<1.$$

$${\omega }_{\text{av}}$$

$$=[\underset{|{\mathcal{P}}_1|}{\underbrace{{\omega }_{\text{av}}^{(1)},\dots ,{\omega }_{\text{av}}^{(1)}}},\dots ,\underset{|{\mathcal{P}}_m|}{\underbrace{{\omega }_{\text{av}}^{(m)},\dots ,{\omega }_{\text{av}}^{(m)}}}{]}^{\mathsf{T}}.$$

$${\dot{\theta }}_i={\tilde{\omega }}_i+\sum _{j\ne i}{\tilde{a}}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$${\dot{\theta }}_j-{\dot{\theta }}_i={\omega }_j+\delta {\omega }_i-{\omega }_i-\delta {\omega }_j$$

$$+\sum _{z=1}^n\left[(a_{jz}+\delta a_{jz})\sin ({\theta }_z-{\theta }_j)-(a_{iz}+\delta a_{iz})\sin ({\theta }_z-{\theta }_i)\right]$$

$$={\omega }_j-{\omega }_i+\sum _{z=1}^n\left[a_{jz}\sin ({\theta }_z-{\theta }_j)-a_{iz}\sin ({\theta }_z-{\theta }_i)\right]+{\delta }_{ij},$$

$$|{\theta }_j(t)-{\theta }_i(t)|\le c\gamma \forall t\ge T,$$

$$V(x_{\text{intra}})=\sum _{k=1}^m{d}_k{x}_{\text{intra}}^{(k)\mathsf{T}}{P}_k{x}_{\text{intra}}^{(k)},$$

$$\dot{V}$$

$$\le -\frac{1}{2}{\lambda }_{\text{min}}(DS+S^{\mathsf{T}}D){\Vert x_{\text{intra}}\Vert }^2,$$

$${\dot{s}}_i=z_i-{\kappa }_i{s}_i-{\gamma }_i(f_i-1),{\dot{f}}_i=s_i,$$

$$\tau {\dot{\mu }}_i=f_i-{\mu }_i^{1/\alpha },{\tau }_i{\dot{q}}_i=f_{i}E(f_i,{\rho }_i)/{\rho }_i-{\mu }_i^{1/\alpha }{q}_i/{\mu }_i.$$

$$y_i=V_0(k_1(1-q_i)+k_2(1-q_i/v_i)+k_3(1-{\mu }_i)),$$

Full text

I Introduction

The structural (i.e., matrix of anatomical connections between brain regions) and functional (i.e., matrix of correlation coefficients between the activity of brain regions) connectivity of the brain vary across healthy individuals and those affected by neurological or psychiatric disorders, and can be used as biomarkers to detect or predict pathological conditions. While structural connectivity changes rather slowly over time and can be measured accurately via diffusion imaging techniques [1], functional connectivity depends on the instantaneous neural activity and is affected, for instance, by the tasks being performed and external stimuli [2]. Today, common measures of functional connectivity rely on resting-state functional magnetic resonance imaging (rs-fMRI) timeseries to quantify the level of correlated activity between brain regions. The relationships between structural and functional connectivity have recently received considerable attention [3, 4], and the tantalizing idea of controlling functional states by leveraging or modifying brain structure has given birth to a new, thrilling, field of research [5, 6, 7].

In this paper, we leverage the connection between structural and functional connectivity, and propose a framework to control functional connectivity by selectively modifying structural connectivity and the regions’ intrinsic frequencies (see Fig. 1). In particular, building on prior studies [8, 9], we model the brain’s neural activity as the phases of a collection of interconnected Kuramoto oscillators, and postulate that the level of functional connection between two regions is proportional to the level of synchronization between the phases of the oscillators associated with the two regions. Then, we derive conditions and methods to tune the oscillators’ interconnection weights and natural frequencies so as to enforce arbitrary synchronization patterns and, consequently, brain-wide functional connectivity. We remark that the control mechanisms used in our framework are biologically plausible. For instance, changes in the spontaneous neural activity (i.e., oscillators’ frequencies) are typical of the brain, involve natural modifications in regional metabolism of the neurons, and can alternatively be induced by a number of non-invasive stimulation techniques [10]. Changes to the structural interconnections (i.e., oscillators’ interconnections), instead, can arise from different chemical or electrical mechanisms including, at the microscale, Hebbian plasticity [11] and short-term synaptic facilitation [12].

Related work. The discovery of oscillatory or rhythmic brain activity dates back almost a century. Yet, control-theoretic studies that exhaust the oscillatory nature of brain states have been sparse and of relatively recent date. Some authors focus on localized desynchronization of neural activity [13, 14, 15], which is desirable in individuals affected by epilepsy or Parkinson’s disease, and others use synchronization phenomena to describe cognitive and functional brain states [16, 17, 18]. To the best of our knowledge, a framework to control the pattern of brain-wide functional connectivity is still missing, and is proposed for the first time in this paper.

At the core of our framework to model and control functional connectivity is the concept of cluster synchronization in a network of oscillators, where groups of oscillators behave cohesively but independently from other clusters. For the case of oscillators with Kuramoto dynamics as used in this work, [19, 20] explore approximate notions of cluster synchronization in simplified configurations, while [21] provides exact invariance conditions for arbitrary cluster synchronization manifolds. Our recent work introduces rigorous [22] and approximate [23] stability conditions for cluster synchronization, which are also used here. Compared to the above references, this paper focuses on the control of cluster synchronization, rather than on its enabling conditions.

Paper contribution. The contributions of this paper are twofold. On the technical side, we formulate and solve a network optimization problem to enforce stable cluster synchronization among interconnected Kuramoto oscillators (Section III). We provide a two-step procedure to compute the smallest (as measured by the Frobenius norm) perturbation of the network weights and the oscillators’ natural frequencies so as to achieve a desired and arbitrary synchronization pattern. Notably, the proposed algorithm allows for the modification of only a selected subset of the network parameters, as typically constrained in applications. We also prove that cluster synchronization is robust to parameter mismatches and numerical inaccuracies, which complements the theoretical derivations in [22, 23], and strengthen the applicability of our control methods to work in practice.

On the application side, this work contains the first mathematically rigorous and neurologically plausible framework to control functional connectivity in the brain, and takes a significant step to fill the gap between empirical studies on oscillatory neural activity [8, 24, 9] and the recent technical body of work inspired by neural synchronization [19, 21, 22, 23]. In Section V, we apply our control technique to an empirically-reconstructed structural brain network, and validate our results by computing the correlation of resting-state fMRI signals obtained through a realistic hemodynamic model. As a minor contribution, our work extends [9] by allowing heterogeneous Kuramoto dynamics.

Mathematical notation. The sets $\mathbb{R}_{>0}$, $\mathbb{S}^{1}$ and $\mathbb{T}^{n}$ denote the positive real numbers, the unit circle, and the $n$-dimensional torus, respectively. We represent the vector of all ones with $\mathds{1}$. The Frobenius and $\ell_{2}$ norms are denoted as $\|\cdot\|_{\text{F}}$ and $\|\cdot\|$, respectively, and $A\circ B$ is the Hadamard product between matrices $A$ and $B$. A (block-)diagonal matrix is denoted by $\mathrm{(blk)diag}(\cdot)$. We let $\operatorname{i}=\sqrt{-1}$. Let $A\geq 0$ represent an element-wise inequality on the entries of $A$, $A^{+}$ the element-wise nonnegative part of $A$, and $A\succ 0$ a positive definite matrix $A$. We let $\lambda_{i}(A)$ and $\sigma_{i}(A)$ denote the $i$-th eigenvalue and the $i$-th singular value of $A\in\mathbb{R}^{n\times n}$, respectively, and $\lambda_{\text{max}}(A)=\max_{i}|\lambda_{i}(A)|$ and $\lambda_{\text{min}}(A)=\min_{i}|\lambda_{i}(A)|$. Finally, we let $\overline{\lambda}(A)=\frac{1}{n}\sum_{i}\lambda_{i}(A)$ and $\overline{\sigma}(A)=\frac{1}{n}\sum_{i}\sigma_{i}(A)$.

II Problem setup and preliminary notions

The aim of this work is to control network parameters so that groups of brain regions exhibit a high degree of functional connectivity. In this context, functional interactions are defined as the pairwise correlation between hemodynamic signals recorded in two brain regions. One model used to simulate such hemodynamic signals is described by a set of nonlinear differential equations [25] that can be approximated in the frequency domain as a linear low-pass filter [9]. Because the only input to such hemodynamic model is the oscillatory neural activity, the formation of strongly (functionally) connected brain regions can be promoted by controlling the synchronization level of their neural dynamics. We follow [9] to model such neural dynamics with a sparse network of heterogeneous Kuramoto oscillators that are connected to each other according to the anatomical architecture of the human brain, more specifically known as white matter tracts.111We assume that at each node of a structural brain network there exists a community of excitatory and inhibitory neurons whose dynamical state is in a regime of self-sustained oscillation. In other words, the neurons’ firing rates delineate a limit cycle, and their dynamics can be approximated by a single variable, which is the angle (or phase) on this cycle. Ultimately, the problem of generating desired patterns of functional connectivity reduces to the one of controlling cluster synchronization in a network of heterogeneous Kuramoto oscillators.

To be precise, let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a weighted digraph, where $\mathcal{V}=\{1,\dots,n\}$ and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ represent the oscillators, or nodes, and their interconnection edges, respectively. The $i$-th oscillator’s dynamics reads as:

[TABLE]

where $\omega_{i}\in\mathbb{R}_{>0}$ denotes the natural frequency of the $i$-th oscillator, $\theta_{i}\in\mathbb{S}^{1}$ is its phase, $a_{ij}\in\mathbb{R}_{>0}$ is the weight of the edge $(j,i)\in\mathcal{E}$, with $a_{ij}=0$ when $(j,i)\not\in\mathcal{E}$, and $A=[a_{ij}]$ is the weighted adjacency matrix of $\mathcal{G}$.

To characterize synchronized trajectories among subsets of oscillators, let $\mathcal{P}=\{\mathcal{P}_{1},\dots,\mathcal{P}_{m}\}$ be a nontrivial partition of $\mathcal{V}$, where each cluster contains at least two oscillators and its graph is strongly connected.222As the brain is densely connected [26], this assumption is not restrictive. We say that a network exhibits cluster synchronization when the oscillators can be partitioned so that the phases of the oscillators in each cluster evolve identically. Formally, we define the cluster synchronization manifold associated with the partition $\mathcal{P}$ as

[TABLE]

Then, the network is cluster-synchronized with partition $\mathcal{P}$ when the phases of the oscillators belong to $\mathcal{S}_{\mathcal{P}}$ at all times. Without loss of generality, the oscillators are labeled so that $\mathcal{P}_{k}=\{\sum_{\ell=1}^{k-1}|\mathcal{P}_{\ell}|+1,\dots,\sum_{\ell=1}^{k}|\mathcal{P}_{\ell}|\}$, where $|\mathcal{P}_{\ell}|$ denotes the cardinality of the set $\mathcal{P}_{\ell}$.

Because our control framework leverages conditions for the invariance and stability of the cluster synchronization manifold to modify the network weights and oscillators’ natural frequencies, we briefly recall useful preliminary results that have recently been established in [21, 22, 23]. Specifically, given a desired network partition $\mathcal{P}=\{\mathcal{P}_{1},\dots,\mathcal{P}_{m}\}$, invariance of $\mathcal{S}_{\mathcal{P}}$ is guaranteed by the following conditions:

• (C1)

The natural frequencies satisfy $\omega_{i}=\omega_{j}$ for every $i,j\in\mathcal{P}_{k}$ and $k~{}\in~{}\{1,\dots,m\}$. Equivalently, $B_{\text{span}}^{\mathsf{T}}\omega=0$,

where $B_{\text{span}}\in\mathbb{R}^{|\mathcal{V}|\times|\bigcup_{k}\mathcal{E}_{\text{span},k}|}$ is the incidence matrix of $\bigcup_{k=1}^{m}\mathcal{T}_{k}$, with $\mathcal{T}_{k}=(\mathcal{P}_{k},\mathcal{E}_{\text{span},k})$ being a spanning tree of the digraph $\mathcal{G}_{k}$ of the isolated cluster $\mathcal{P}_{k}$;

• (C2)

The network weights satisfy $\bar{V}_{\mathcal{P}}^{\mathsf{T}}\bar{A}V_{\mathcal{P}}=0$,

where $V_{\mathcal{P}}\in\mathbb{R}^{n\times m}$ is the characteristic matrix of the network defined as $V_{\mathcal{P}}=\left[{v_{1}}/{\|v_{1}\|},\dots,{v_{m}}/{\|v_{m}\|}\right]$, with

[TABLE]

$\bar{V}_{\mathcal{P}}\in\mathbb{R}^{n\times(n-m)}$ is an orthonormal basis of the orthogonal subspace to the image of $V_{\mathcal{P}}$, and $\bar{A}=A-A\circ V_{\mathcal{P}}V_{\mathcal{P}}^{\mathsf{T}}$ is the matrix of inter-cluster connections only (see also [21]).

We assume that the isolated clusters are locally stable:

• (A$1$)

The dynamics (1), with $a_{ij}=0$ when $i,j$ belong to different clusters, converges exponentially fast to $\mathcal{S}_{\mathcal{P}}$.

Notice that Assumption (A$1$) is satisfied when $\mathcal{G}_{k}$ has symmetric weights and condition (C$1$) holds [22, Lemma 3.1][27, Theorem 5.1]. In our case, (A$1$) is not restrictive because structural brain networks are typically symmetric [5].333In a general case, one can ensure that Assumption (A$1$) is satisfied simply by pairing the control mechanism developed in the next session with an independent one that makes intra-cluster connections symmetric.

Let $\omega^{(k\ell)}$ denote the natural frequency difference between any two nodes in disjoint clusters $\mathcal{P}_{k}$ and $\mathcal{P}_{\ell}$. If (C$1$) and (C$2$) hold, then a tight approximate condition for $\mathcal{S}_{\mathcal{P}}$ to be locally exponentially stable is [23]:

• (C3)

The natural frequencies and the network weights are such that $\lambda_{\text{max}}(\Xi(A,\omega))<1$, with $\Xi=[\xi_{k\ell}]$ and

[TABLE]

where $G_{k}(s)=(sI-J_{k})^{-1}$, $J_{k}$ is the Hurwitz stable Jacobian matrix of the intra-cluster phase difference dynamics and $\nu_{k\ell}$ is a function of the inter-cluster weights. Due to space constraints, we refer the interested reader to [21, 23] for a detailed discussion on conditions (C1), (C2), and (C3).

III Control of cluster synchronization

In this section, we propose a control mechanism to obtain a prescribed and robust configuration of synchronized oscillatory patterns. Towards this aim, we consider a network $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and an arbitrary partition $\mathcal{P}=\{\mathcal{P}_{1},\dots,\mathcal{P}_{m}\}$ of $\mathcal{V}$. The proposed control technique is minimally invasive in the sense that it looks for the smallest correction (in the Frobenius norm sense) of inter-cluster network weights and oscillators’ natural frequencies that renders the cluster synchronization manifold $\mathcal{S}_{\mathcal{P}}$ invariant and locally stable. In practice, a modification of the network parameters will require either the exploitation of neural plasticity or localized surgical intervention for the modification of the network weights and structure, and pharmacological or electromagnetic influence for the refinement of the brain regions’ natural frequencies. In mathematical terms, the approach is encoded into solving the following constrained minimization problem:

[TABLE]

where $\Delta$ is the correction of the network matrix, $\mu$ is the correction of the natural frequencies vector, and the $(i,j)$-th entry of $\Delta$ is zero if $i$, $j$ belong to the same partition $\mathcal{P}_{k}$, $k\in\{1,\dots,m\}$. Further, $H$ is the [math]-$1$ adjacency matrix of $\mathcal{H}=(\mathcal{V},\mathcal{E}_{\mathcal{H}})$, which is the graph encoding the set of edges $\mathcal{E}_{\mathcal{H}}\subseteq\mathcal{E}$ that is allowed to be modified, and $H^{\text{c}}=\mathds{1}\mathds{1}^{\mathsf{T}}-H$. That is, the $(i,j)$-th entry of a solution $\Delta^{*}$ to problem (3) is zero when the corresponding $(i,j)$-th entry of $H$ is zero. The optimization problem (3) is illustrated in Fig. 2.

Constraints (3a) and (3b) are equivalent to conditions (C$2$) and (C$1$), respectively, for the invariance of $\mathcal{S}_{\mathcal{P}}$. Constraint (3c) restricts the corrective action to a subset of all the possible interconnections, in affinity with the practical limitations of localized interventions. Constraints (3d) and (3e) are due to biological compatibility and require the inter-cluster weights of the perturbed network and oscillators’ natural frequencies to be nonnegative. Finally, Constraint (3f) corresponds to (C3) and guarantees the (local) stability of $\mathcal{S}_{\mathcal{P}}$. In particular, the latter constraint makes the above problem non-convex and, therefore, potentially intractable from a numerical viewpoint. To overcome this issue, we next propose a suboptimal, yet numerically more tractable, control strategy. Specifically, we decouple (3) into two simpler subproblems. The first one solves for the smallest correction of inter-cluster weights satisfying (3a), (3c), and (3d), whereas the second one solves for the smallest correction of the oscillators’ natural frequencies satisfying (3b), (3e) and (3f).

III-A Inter-cluster structural control for invariance of $\mathcal{S}_{\mathcal{P}}$

We first address the problem of computing the smallest correction of inter-cluster weights such that constraints (3a), (3c), and (3d) are satisfied. Specifically, we focus on the following minimization problem:

[TABLE]

The optimization problem (4) is convex and, when feasible, it can be efficiently solved by means of standard optimization techniques. Feasibility of (4) depends on the constraint graph $\mathcal{H}$ (see Remark 1). In what follows, we present a simple and efficient projection-based algorithm to solve this problem.

Theorem III.1

(Smallest sparse inter-cluster correction)* Assume that the problem (4) is feasible, and consider the matrix sequence $\{Z_{k}\}_{k\geq 0}$ generated via the following iterative procedure:*

[TABLE]

where $Z_{0}=\bar{A}$, and $T_{0}=Q_{0}=0$. Then, the sequence $\{Z_{k}\}_{k\geq 0}$ converges to a matrix $Z^{*}$, and a minimizer of (4) subject to (4a), (4b), and (4c), has the form $\Delta^{*}=Z^{*}-\bar{A}$.

Proof:

Let $\Pi_{\mathcal{Z}}(W)=\arg\min_{Z\in\mathcal{Z}}\|Z-W\|_{\text{F}}^{2}$ denote the projection (in the Frobenius norm sense) of $W$ onto a convex set $\mathcal{Z}$, and define the closed convex sets $\mathcal{Z}_{1}=\{Z\in\mathbb{R}^{n\times n}\;:\;Z\geq 0\text{ and }H^{\text{c}}\circ Z=\bar{A}\}$ and $\mathcal{Z}_{2}=\{Z\in\mathbb{R}^{n\times n}\;:\;\bar{V}_{\mathcal{P}}^{\mathsf{T}}ZV_{\mathcal{P}}=0\}$. Note that $\Pi_{\mathcal{Z}_{1}}(W)=H\circ W^{+}+H^{\text{c}}\circ\bar{A}$ and, by Lemma .1 in the Appendix,

[TABLE]

for any $W$. Hence, the sequence $\{Z_{k}\}_{k\geq 0}$ generated by (5) coincides with the sequence generated by Dykstra’s projection algorithm [29] applied to the projections onto $\mathcal{Z}_{1}$ and $\mathcal{Z}_{2}$. Since the problem (4) is feasible, $\mathcal{Z}_{1}\cap\mathcal{Z}_{2}\neq\emptyset$, and the latter sequence converges to a matrix $Z^{*}=\Pi_{\mathcal{Z}_{1}\cap\mathcal{Z}_{2}}(\bar{A})=\arg\min_{Z\in\mathcal{Z}_{1}\cap\mathcal{Z}_{2}}\|Z-\bar{A}\|_{\text{F}}^{2}$ [29]. Finally,

[TABLE]

and the statement follows. ∎

Remark 1

(Sufficient condition for the feasibility of (4))* Recall from [21] that condition (C$2$) is equivalent to:*

[TABLE]

for every $i,j\in\mathcal{P}_{z}$ and for all $z,\ell\in\{1,\dots,m\}$, with $z\neq\ell$. Notice that, if for every $i\in\mathcal{P}_{z}$ there exists at least one $(k,i)\in\mathcal{E}_{\mathcal{H}}$, $k\in\mathcal{P}_{\ell}$, a solution to (6) can always be found and problem (4) is feasible. $\square$

III-B Frequency tuning for local stability of $\mathcal{S}_{\mathcal{P}}$

We now turn to the problem of computing the smallest correction of natural frequencies such that constraints (3b), (3e), and (3f) are satisfied. That is,

[TABLE]

Theorem III.2

(Feasibility of problem (7))* There always exists a correction $\mu$ satisfying (7a), (7b), and (7c).*

Proof:

Consider the vector $\mu=[\mu_{1},\dots,\mu_{n}]^{\mathsf{T}}$. Note that we can find some $\mu_{i}$ such that (i) $\omega_{i}+\mu_{i}=\omega_{j}+\mu_{j}>0$ for all $i,j\in\mathcal{P}_{k}$, $k\in\{1,\dots,m\}$, and (ii) $|\omega_{i}+\mu_{i}-(\omega_{j}+\mu_{j})|>\eta$, for all $i\in\mathcal{P}_{k}$, $j\in\mathcal{P}_{\ell}$, $k,\ell\in\{1,\dots,m\}$, $k\neq\ell$, and $\eta>0$ arbitrarily large. From (i), $\mu$ satisfies (7a) and (7b). Further, since each nonzero entry of $\Xi(A,\omega+\mu)$ in (2) behaves as a low-pass filter, by fact (ii) $\lambda_{\text{max}}(\Xi(A,\omega+\mu))$ can be made arbitrarily small. This implies that there always exists a vector $\mu$ satisfying (7c) and concludes the proof. ∎

An optimal solution to (7) is typically difficult to compute, because of the eigenvalue constraint (7c). However, several heuristics can be used to compute a suboptimal correction in (7). For instance, we next outline an effective procedure to find a suboptimal solution to (7). Let $\omega_{\text{av}}^{(k)}=\frac{1}{|\mathcal{P}_{k}|}\sum_{i\in\mathcal{P}_{k}}\omega_{i}$ denote the average frequency within each cluster, and let

[TABLE]

Further, define the quotient graph $\mathcal{Q}=(\mathcal{V}^{\prime},\mathcal{E}^{\prime})$ where each node in $\mathcal{V}^{\prime}$ represents a cluster and each edge in $\mathcal{E}^{\prime}$ an interconnection between two clusters. Our procedure leverages Theorem III.2 and increases the frequency differences between pairs of connected clusters until constraint (7c) is satisfied. The procedure consists of four steps:

1. (i)

If $\lambda_{\text{max}}(\Xi(A,\omega_{\text{av}}))<1$, then $\mu^{*}=[\mu_{1}^{*},\dots,\mu_{n}^{*}]$, with $\mu_{i}^{*}=\omega_{\text{av},i}-\omega_{i}$ is an optimal correction to (7). Otherwise, proceed to the next step. 2. (ii)

Construct a depth-first spanning tree $\mathcal{T}_{\mathcal{Q}}$ of $\mathcal{Q}$ rooted at $r=\arg\min_{k}\omega_{\text{av}}^{(k)}$.444Notice that such a spanning tree always exists, since $\mathcal{Q}$ is connected. 3. (iii)

Assign the frequency $\omega(k,\alpha)=\omega^{(r)}+k\alpha$, $\alpha>0$, to each node of each cluster in $\mathcal{T}_{\mathcal{Q}}$ of depth $k$, $k=1,2,\dots,k_{\text{max}}$, where $k_{\text{max}}$ denotes the height of $\mathcal{T}_{\mathcal{Q}}$.555Given a connected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and a spanning tree $\mathcal{T}$ of $\mathcal{G}$ rooted at $r\in\mathcal{V}$, the depth of a node $v\in\mathcal{V}$ is the length of the path in $\mathcal{T}$ from $r$ to $v$, and the height of $\mathcal{T}$ is the maximum depth among the nodes in $\mathcal{V}$. Let $\omega(\alpha)=[\omega_{1}(\alpha),\dots,\omega_{n}(\alpha)]^{\mathsf{T}}$ denote the resulting vector of modified frequencies. 4. (iv)

Find the smallest $\alpha^{*}$ satisfying $\lambda_{\text{max}}(A,\omega(\alpha^{*}))<1$. Then, $\mu^{*}=[\mu_{1}^{*},\dots,\mu_{n}^{*}]$, with $\mu_{i}^{*}=\omega_{i}(\alpha^{*})-\omega_{i}$, is a (suboptimal) solution to (7).

IV Robustness of the control framework

In this section, we show that the control framework described in Section III, and in fact the stability property of the cluster synchronization manifold $\mathcal{S}_{\mathcal{P}}$, is robust to perturbations of the network parameters. That is, small changes in the oscillators’ natural frequencies and network weights yield a small deviation from cluster-synchronized trajectories. In light of this, the proposed control mechanism lends itself to practical applications, where the network parameters are not known exactly and the neural dynamics is subject to noise.

Consider the dynamics (1) with perturbed parameters:

[TABLE]

where $\tilde{\omega}_{i}=\omega_{i}+\delta\omega_{i}$ and $\tilde{a}_{ij}=a_{ij}+\delta a_{ij}$. Notice that, if $\delta\omega_{i}=0$ and $\delta a_{ij}=0$, the dynamics (8) is equivalent to (1). From (8), the perturbed intra-cluster difference dynamics of nodes $i,j\in\mathcal{P}_{k}$, with $k\in\{1,\dots,m\}$, reads as:

[TABLE]

where $\delta_{ij}=\delta\omega_{j}-\delta\omega_{i}+\sum_{z=1}^{n}[\delta a_{jz}\sin(\theta_{z}-\theta_{j})$ $-\delta a_{iz}\sin(\theta_{z}-\theta_{i})]$. Finally, let $\delta$ be the vector of all $\delta_{ij}$ that affect the nominal intra-cluster dynamics as in (IV).

We are now ready to present the main result of this section, which resorts to the prescriptive stability condition derived in [22] that we recall in the Appendix B for completeness.

Theorem IV.1

(Robustness of cluster synchronization)* Assume that the network weights satisfy Theorem .1, and consider any pair of nodes $i,j\in\mathcal{P}_{k}$, $k\in\{1,\dots,m\}$. Then, for some finite $T>0$ and for all initial conditions such that $|\theta_{j}(0)-\theta_{i}(0)|<\varepsilon$, with $\varepsilon>0$ sufficiently small, the solution to the perturbed dynamics (8) satisfies*

[TABLE]

where $\gamma=\max_{\theta\in[0,2\pi)}\|\delta\|$, and $c$ is a constant that depends only on the network weights.

Proof:

In the first part of the proof, we combine the Lyapunov functions for the isolated clusters $\mathcal{P}_{k}$, $k=1,\dots,m$, into a Lyapunov function for the intra-cluster differences dynamics of the whole network. In the second part of the proof, we show that such Lyapunov function satisfies certain bounds, so that the application of [30, Lemma 9.2] suffices to prove the claimed statement.

We let $x_{ij}=\theta_{j}-\theta_{i}$, and $S$, $x_{\text{intra}}$ and $J_{k}$ be as in the Appendix B. To combine the Lyapunov functions of the isolated clusters, we note that if $S$ is an $M$-matrix, then, along the lines of [22, Proof of Theorem 3.2], the origin of the nominal intra-cluster dynamics of $x_{\text{intra}}$ is locally exponentially stable with Lyapunov function

[TABLE]

where $P_{k}\succ 0$ satisfies $J_{k}^{\mathsf{T}}P_{k}+P_{k}J_{k}=-I$, and $d_{k}>0$ are such that $DS+S^{\mathsf{T}}D\succ 0$, with $D=\mathrm{diag}(d_{1}\dots,d_{m})$ [30].

Consider now the Lyapunov function (11), and notice that: $c_{1}\|x_{\text{intra}}\|^{2}\leq\sum_{k=1}^{m}d_{k}x_{\text{intra}}^{(k)\mathsf{T}}P_{k}x_{\text{intra}}^{(k)}\leq c_{2}\|x_{\text{intra}}\|^{2},$ with $c_{1}=md_{\text{min}}\min_{k}\lambda_{\text{min}}(P_{k})$ and $c_{2}=md_{\text{max}}\max_{k}\lambda_{\text{max}}(P_{k})$. Further, in the ball of radius $r$ of the origin $\mathcal{B}_{r}=\{x_{\text{intra}}\;:\;\|x_{\text{intra}}\|<r,\,\dot{V}(x_{\text{intra}})<0\}$, it holds that $\dot{V}(x_{\text{intra}})\leq-c_{3}\|x_{\text{intra}}\|^{2}$, with $c_{3}=\lambda_{\text{min}}(DS+S^{\mathsf{T}}D)/2$. To see this, consider the derivative of the Lyapunov function $V(x_{\text{intra}})$ along the trajectories of the nominal system. Then, from [30, §9.5] and [22, proof of Theorem 3.2], the following inequality holds in $\mathcal{B}_{r}$:

[TABLE]

and $c_{3}$ follows. Further, since $\|\partial V/\partial x_{\text{intra}}\|=\|2x_{\text{intra}}^{\mathsf{T}}P_{k}\|<2\lambda_{\text{max}}(P_{k})\|x_{\text{intra}}\|$, we have $\|\partial V/\partial x_{\text{intra}}\|\leq c_{4}\|x_{\text{intra}}\|$, with $c_{4}=2c_{2}$. Finally, once the constants $c_{1},c_{2},c_{3}$, and $c_{4}$ are computed, the definition of $x_{\text{intra}}$ and [30, Lemma 9.2] conclude the proof. ∎

Importantly, Theorem IV.1 can be used to provide a quantitative bound on the asymptotic value of $|\theta_{j}-\theta_{i}|$. In fact, we can compute the constant $c$ in (10) by exploiting [30, Lemma 9.2] and $c_{1},c_{2},c_{3},c_{4}$ derived in the above proof.

V Control of functional connectivity in an empirically-reconstructed brain network

We conclude this paper with the application of the control mechanism presented in Section III to the brain network estimated in [1], which is publicly available at http://umcd.humanconnectomeproject.org/umcd. In these data, structural connectivity is proportional to large-scale connection pathways between cortical regions, and the gray matter is subdivided into $n=66$ cortical regions ($33$ per hemisphere). To show the effectiveness of our proposed method in enforcing desired functional connectivity by means of arbitrary synchronization patterns, we partition the structural brain network in $3$ clusters, i.e. $\mathcal{P}=\{\mathcal{P}_{1},\mathcal{P}_{2},\mathcal{P}_{3}\}$, each one comprising $22$ regions that do not belong to any known functionally connected resting-state network. The three clusters are highlighted with different colors in Fig. 3 and Fig. 3. Furthermore, following our goal of providing a method to enhance the synchronization properties of a diseased or damaged brain, we simulate the effects of brain damage, e.g., a stroke, by damping the connectivity of one cluster [31]. That is, we weaken the intra-cluster connections of the first cluster by a scaling factor $10^{-2}$ to echo reduced structural connectivity, and we show that our technique can in fact recover stability of the cluster synchronization manifold associated with the desired network partition.666Specifically, weakening the intra-cluster connections of one cluster is likely to make $\mathcal{S}_{\mathcal{P}}$ unstable [23].

Before presenting our results, we describe the methodology used to simulate human rs-fMRI functional connectivity.

V-A Simulation of functional connectivity

The brain’s neural activity is simulated through a network of coupled Kuramoto oscillators, where we randomly draw the natural frequencies of each oscillator from a uniform distribution in the range $[0,60]$ $[$Hz$]$ so as to include all meaningful neural frequency bands [25]. We set the initial phases in the interval $[0,0.5]$ $[$rad$]$. The Kuramoto phases act as an input to the neurovascular coupling, which is modeled by the Balloon-Windkessel hemodynamic process [33], and whose output is the blood-oxygen-level dependent (BOLD) signal that is measured by rs-fMRI.

The neuronal activity $z_{i}$ of the $i$-th brain region produces an increase in a vasodilatory signal $s_{i}$, which is subject to auto-regulatory feedback. The inflow $f_{i}$ responds in proportion to this signal with concomitant changes in blood volume $\mu_{i}$ and deoxyhemoglobin content $q_{i}$. Mathematically, the dynamics of these quantities reads as:

[TABLE]

The oxygen extraction is a function of the flow $E(f,\rho)=1-(1-\rho^{\/f})$ where $\rho$ denotes the resting oxygen extraction fraction. The biophysical parameters $\kappa,\gamma,\tau,\alpha,$ and $\rho$ are exhaustively treated in [33]. Finally, the BOLD signal is described as a static nonlinear function:

[TABLE]

where $V_{0}=0.02$ denotes the resting blood volume fraction, and $k_{1}=7\rho_{i}$, $k_{2}=2$, $k_{3}=2\rho_{i}-0.2$. Following [9], we choose $z_{i}=\sin(\theta_{i})$. Further, to account for the presence of background noise in the brain, we add white noise to the neural activity $z_{i}$ with variance $10^{-2}$. We simulate $2$ minutes of BOLD signals and process the timeseries as explained below in order to compute functional connectivity estimates that closely resemble that of human rs-fMRI recordings.

To reduce the effect of spurious correlations from small and non-physiological high-frequency components, we filter the synthetic BOLD signals through a low-pass filter. Consequently, to improve the correspondence between resting-state correlations and anatomical connectivity, we process all of the simulated regional BOLD signals by a global signal regression [34] that averages the timeseries of all regions by removing spontaneous oscillations common to the whole brain. Next, we discard the first $40$ seconds of all timeseries to eliminate the effect of initial transients. Finally, we compute the Pearson correlation of the filtered and regressed signals to obtain the synthetic functional connectivity. A pipeline describing the above process is illustrated in Fig. 4.

V-B Application of the clustering control mechanism

In the remainder of this section, we apply the control method proposed in Section III. We first solve the minimization problem (4) to find the optimal correction matrix $\Delta^{*}$ to be applied to $A_{\text{original,SC}}$ such that condition (C$2$) for the invariance of $\mathcal{S}_{\mathcal{P}}$ is satisfied. We choose to constrain the corrective action on a set of edges $\mathcal{E}_{\mathcal{H}}=\mathcal{E}\cup\tilde{\mathcal{E}}$ that includes the original set $\mathcal{E}$ and a minimal set $\tilde{\mathcal{E}}$ of randomly selected edges such that problem (4) is feasible (see Remark 1). Fig. 3 and 3 illustrate the corrective action $\Delta^{*}$ and the network matrix $A_{\text{corrected,SC}}=A_{\text{original,SC}}+\Delta^{*}$, respectively.

We proceed with the frequency tuning technique for invariance and stability of $\mathcal{S}_{\mathcal{P}}$ to $A_{\text{corrected,SC}}$ so that conditions (C$1$) and (C$3$) are satisfied. The first step involves computing the mean natural frequency $\omega_{\text{av}}^{(k)}$ among all oscillators belonging to the same cluster $\mathcal{P}_{k}$: $\omega_{\text{av}}^{(1)}=199.2$, $\omega_{\text{av}}^{(2)}=182.9$ and $\omega_{\text{av}}^{(3)}=115.4$ $[$rad/s$]$. Next, we apply the procedure proposed in Section III-B. We plot in Fig. 5 the spanning tree of the quotient graph $\mathcal{T}_{\mathcal{Q}}$, and in Fig. 5 the optimal $\alpha^{*}$ computed in step (iv) of our procedure. The final natural frequencies are $\omega^{(1)}=131.8$, $\omega^{(2)}=126.4$ and $\omega^{(3)}=115.4$ $[$rad/s$]$. Notice that, although our frequency tuning procedure is sub-optimal, the outcome values remain well within the range of brain activity frequency bands and, based on numerical results, outperform the results of Matlab’s fmincon function.

Finally, by following the pipeline described in the previous subsection, we compute the desired functional connectivity pattern, which we show in Fig. 6. Notably, the functional connectivity of the desired clusters is strong and the correlations between different clusters are negligible. Thus, the proposed method to control synchronization patterns of oscillatory neural activity lends itself to a physiologically plausible framework and shows rather promising results.

VI Conclusion

In this work, we propose a minimally invasive technique to obtain robust synchronization patterns in sparse networks of heterogeneous Kuramoto oscillators. To the best of our knowledge, this is the first attempt at blending mathematically rigorous methods with physiological models of brain activity with the goal of steering whole-brain synchronization dynamics. Specifically, we cast a constrained optimization problem whose solution not only satisfies mathematical conditions for invariance and stability of an arbitrary cluster synchronization manifold, but also meets biological constraints. We decompose the complete optimization problem into two simpler subproblems, and provide efficient methods to solve them. When applying our technique to correct the network parameters of empirically-reconstructed anatomical brain data, we find that our solution, although suboptimal, provides a result that is well within the range of physiologically plausible parameters. Additionally, we show that cluster synchronization is robust to small parameter mismatches and numerical inaccuracies. This result complements previous prescriptive studies on cluster synchronization and enables the use of our framework in practical situations.

A Instrumental result for the proof of Theorem III.1

Lemma .1

Consider a network $\mathcal{G}=(\mathcal{V},\mathcal{E})$, and an arbitrary (nontrivial) partition $\mathcal{P}=\{\mathcal{P}_{1},\dots,\mathcal{P}_{m}\}$ of $\mathcal{V}$. Let $W\in\mathbb{R}^{n\times n}$ and consider the minimization problem

[TABLE]

The minimizer of the problem (12) subject to (12a) is

[TABLE]

Proof:

We prove the result via the method of Lagrange multipliers. The Lagrangian of (12) subject to (12a) is $\mathcal{L}(Z,\Lambda)=\ \|Z-W\|_{\text{F}}^{2}+\mathds{1}^{\mathsf{T}}(\Lambda\circ\bar{V}_{\mathcal{P}}^{\mathsf{T}}ZV_{\mathcal{P}})\mathds{1}=\ \mathrm{tr}((Z-W)^{\mathsf{T}}(Z-W))+\mathrm{tr}(\Lambda^{\mathsf{T}}\bar{V}_{\mathcal{P}}^{\mathsf{T}}ZV_{\mathcal{P}}),$ where $\Lambda\in\mathbb{R}^{(n-m)\times m}$ is a matrix of Lagrange multipliers associated with Constraint (12a), and in the last equation we used that $\mathds{1}^{\mathsf{T}}(A\circ B)\mathds{1}=\mathrm{tr}(A^{\mathsf{T}}B)$. Equating the partial derivatives of $\mathcal{L}$ to zero yelds:

[TABLE]

We next pre- and post-multiply both sides of (14) by $\bar{V}_{\mathcal{P}}^{\mathsf{T}}$ and $V_{\mathcal{P}}$, respectively, and obtain

[TABLE]

where in the second implication we used $V_{\mathcal{P}}^{\mathsf{T}}V_{\mathcal{P}}=I_{n}$, $\bar{V}_{\mathcal{P}}^{\mathsf{T}}\bar{V}_{\mathcal{P}}=I_{n-m}$, and $\bar{V}_{\mathcal{P}}^{\mathsf{T}}V_{\mathcal{P}}=0$, and in the last one we used (15). Finally, (13) follows by substituting (16) into (14). ∎

B $M$-matrix condition for local stability of $\mathcal{S}_{\mathcal{P}}$

We now recall a stability condition established in [22]. We let $x_{ij}=\theta_{j}-\theta_{i}$ denote the phase difference between oscillators $i$ and $j$, and $x_{\text{intra}}=[x_{\text{intra}}^{(1)},\dots,x_{\text{intra}}^{(m)}]^{\mathsf{T}}$ denote a smallest set of intra-cluster differences akin to (see [22]), where $x_{\text{intra}}^{(k)}$ contains only phase differences of oscillators in $\mathcal{P}_{k}$.777The definition in [22] is given for undirected graphs. However, it is straightforward to extend all definitions to digraphs, provided that the quotient graph $\mathcal{Q}$ defined in Section III-B is strongly connected. This is not a restrictive assumption in the context of structural brain networks. Further, let $J_{k}$ be the Hurwitz stable Jacobian matrix of the intra-cluster phase difference dynamics of $x_{\text{intra}}^{(k)}$.

Theorem .1

(Sufficient condition on network weights for the stability of $\mathcal{S}_{\mathcal{P}}$ [22])* Let $\kappa=2\max_{r}|\mathcal{P}_{r}|-2$, and*

[TABLE]

with $k,\ell\in\{1,\dots,m\}$, $i\in\mathcal{P}_{k}$. Define the $m\times m$ matrix $S$:

[TABLE]

where $X_{k}\succ 0$ satisfies $J_{k}^{\mathsf{T}}X_{k}+X_{k}J_{k}=-I$. If $S$ is an $M$-matrix, then the cluster synchronization manifold is locally exponentially stable.

Theorem .1 shows that the cluster synchronization manifold is stable when the intra-cluster dynamics are sufficiently more attractive than the inter-cluster couplings.