Dynamical Graph Theory Networks Methods for the Analysis of Sparse Functional Connectivity Networks and for Determining Pinning Observability in Brain Networks

Neuroimaging in combination with graph theory has been successful in analyzing the functional connectome. However almost all analysis are performed based on static graph theory. The derived quantitative graph measures can only describe a snap shot of the disease over time. Neurodegenerative disease evolution is poorly understood and treatment strategies are consequently only of limited efficiency. Fusing modern dynamic graph network theory techniques and modeling strategies at different time scales with pinning observability of complex brain networks will lay the foundation for a transformational paradigm in neurodegnerative diseases research regarding disease evolution at the patient level, treatment response evaluation and revealing some central mechanism in a network that drives alterations in these diseases. We model and analyze brain networks as two-time scale sparse dynamic graph networks with hubs (clusters) representing the fast sub-system and the interconnections between hubs the slow sub-system. Alterations in brain function as seen in dementia can be dynamically modeled by determining the clusters in which disturbance inputs have entered and the impact they have on the large-scale dementia dynamic system. Observing a small fraction of specific nodes in dementia networks such that the others can be recovered is accomplished by the novel concept of pinning observability. In addition, how to control this complex network seems to be crucial in understanding the progressive abnormal neural circuits in many neurodegenerative diseases. Detecting the controlling regions in the networks, which serve as key nodes to control the aberrant dynamics of the networks to a desired state and thus influence the progressive abnormal behavior, will have a huge impact in understanding and developing therapeutic solutions and also will provide useful information about the trajectory of the disease. In this paper, we present the theoretical framework and derive the necessary conditions for (1) area aggregation and time-scale modeling in brain networks and for (2) pinning observability of nodes in dynamic graph networks. Simulation examples are given to illustrate the theoretical concepts.

1. Introduction

Novel mathematical paradigms such as graph theoretical techniques can capture the brain connectivity and its topology (Fornito et al., 2012; Giessing and Thiel, 2012; Zeng et al., 2012). New descriptors of complex networks are able to quantify induced changes in topology or network organization or serve as theory-driven biomarkers to be used in disease prediction at the level of the individual subject.

While most graph networks applied to brain research are static graph networks and can not capture the dynamical processes governing the time evolution of neurodegenerative diseases, a new paradigm in brain research—dynamical graph networks—is required to advance this field and overcome the obstacles posed by static graph theory in terms of disease prediction and evolution, and its associated connectivity changes.

Biologically-relevant brain networks represent complex large-scale dynamical systems capturing the interaction of a large number of subsystems. Standard analysis methods and control strategies become very difficult and novel techniques need to be developed to (1) analyze the dynamical behavior of these large-scale systems or (2) to observe their states.

To address the first issue of how to control these networks, we need to employ a model simplification resulting in a model of lower complexity, that is easier to handle, and that will provide a simplified synthesis procedure for design problems at a reduced computational complexity. Balanced truncation is known as a popular method for model reduction since it is relatively simple and the quality of the reduced model is guaranteed. The interpretation of most balancing techniques is based on the concept of past and future energy. The most important contribution was the balancing for stable minimal linear systems (Moore, 1981). It is based on a state—space point of view of employing the well—known observability and controllability Gramians and related to the past input energy (controllability) and future input energy (observability). The idea behind transforming a system into balanced form is to easily detect and remove a state component of the initial system to obtain a reduced–order model. The importance of a component is based on Hankel singular values which determine if the output energy of a certain component is small and thus difficult to observe and if the input energy to reach this state is large. While for linear systems finding a balancing coordinate transformation via solutions of the controllability and observability Lyapunov equations is quite easy, for nonlinear systems these equations are almost impossible to solve and thus balancing becomes in general not a simple task (Lall et al., 2002; Scherpen, 1993). These techniques have been applied to the analysis of gene regulatory networks (Meyer-Bäse, 2008), however they are not quite efficient in terms of model reduction for large-scale networks since they involve the computationally expensive operation of matrix factorization. In addition, this method does not preserve the connection structure between subsystems and erases the neurological connection between the subsystem state variables. Therefore, a network topology-preserving mechanism to provide model reduction is required for large-scale networks. We will present an area aggregation and time-scale modeling for sparse brain networks with densely interconnected hubs and externally sparse interconnections between these hubs (Tahmassebi et al., 2017). In Biyik and Arcak (2008) it was shown that the neurons in the hubs synchronize on the fast time-scale and as aggregated neurons determine the slow dynamics of the neural network. We derive a simplified dynamic two-time scale representation for brain connectivity networks assuming linear connections between the nodes. The eigenvalues of the state matrices of the slow and fast system will provide important information about the dynamic evolution of brain connectivity networks at different stages of neurodegenerative diseases. The structural parameters of these networks expressed by the node and area parameter will unveil changes in the sparsity patterns in the course of the disease associated with the disease.

The second issue of observability of neural states can be achieved over synchronization. Synchronization plays an important role in the analysis of neural networks in neurodegenerative diseases. For many neurodegenerative diseases it is very important to obtain some information about some neural states in order to recover the others. This new concept of “pinning observability,” first proposed in Yu et al. (2014), refers to observing a small number of neurons such that the states of the other neurons can be recovered. Differently from the concept of pinning controllability (Chen et al., 2007; Liu et al., 2011; Song and Cao, 2010; Tang et al., 2012), the dynamics of the neurons can be heterogeneous. This property is extremely appealing to brain connectivity networks where we aim to obtain insight into the dynamics of several regions by observing only a few of them. We will derive for these networks a general criterion for synchronization and then some decoupled conditions for pinning observability taking into account their heterogeneous architecture.

In the present paper, we will present a reduced-model approximation over time for large-scale brain networks and derive pinning observability conditions for competitive neural networks and illustrate in an example the theoretical analysis.

2. Reduced-Model Approximation Over Time

Many brain connectivity networks exhibit a heterogeneous structure with densely linked nodes in an area but with sparse connections between these areas. The network is viewed as an interconnected graph with links between the areas which are viewed as nodes in the graph. Thus, the architecture can be described concisely by two main parameters (Chow and Kokotovic, 1985): the node parameter d and the area parameter δ. The node parameter is given as

$\begin{array}{ll}d=\frac{c^E}{c^I}\ll 1& (1)\end{array}$

where c^E is the number of external links of the node with the largest number of external links to nodes outside its area and c^I is the number of internal links of the node with the smallest number of internal links to nodes inside its area. d needs to be a small number. The area parameter is given as

$\begin{array}{ll}\delta =\frac{{\gamma }^E}{m{c}^I}\ll 1& (2)\end{array}$

where γ^E is the number of external links of the area with the largest number of external links and m is the minimal number of nodes found in an area.

While the node and area parameters describe the static behavior and the sparsity pattern of the network, a more interesting approach is to analyze the dynamical behavior of such a network. In order to obtain a reduced-model approximation we view the brain connectivity graph as a structured representation with dense areas (clusters) and sparse interconnections between these areas (Biyik and Arcak, 2008). The time behavior of networks with densely linked neurons in an area but with sparsely connected areas can be described dynamically by a two-time scale system where the neurons within the same area synchronize on the fast time-scale because the dense internal links allow the neural potential in the same area to quickly reach an equilibrium. During the fast time scale the exchange with the other areas is slow because of the sparsity of the links and this becomes significant only over a longer period thus leading to a slow time-scale. This coupled dynamics leads to a reduced-order model describing the long-term behavior of the overall network.

The dynamical formulation will be that the activation states within a group synchronize with each other and achieve a reference group velocity υ(t):

$\begin{array}{l}\underset{t\to \infty }{\lim }|x_i(t)-x_j(t)|=0,i,j=1,\mathrm{\dots },N\\ \underset{t\to \infty }{\lim }|{\dot{x}}_i-\upsilon (t)|=0,i=1,\mathrm{\dots },N& (3)\end{array}$

The synaptic connections d_ik are defined as

$\begin{array}{ll}{d}_{ik}=\{\begin{array}{lll}\hfill 1,& :\hfill & \text{if the i-th node is the positive end of the k-th link}\hfill \\ \hfill -1,& :\hfill & \text{if the i-th node is the negative end of the k-th link}\hfill \\ \hfill 0& :\hfill & \text{otherwise}\hfill \end{array}& (4)\end{array}$

We assume having a neural network of N nodes and M total links yielding thus a N × M incidence matrix D describing this neural network. The formulated objectives in Equation (3) are achieved by the following network architecture

$\begin{array}{ll}{\dot{x}}_i=-{\displaystyle \sum _{k=1}^M{d}_{ik}}{f}_k({\theta }_k)+\upsilon (t),i=1,\cdots ,N& (5)\end{array}$

where the difference variable θ_k is the difference between the positive and negative ends of k-th link, i.e.,

$\begin{array}{ll}{\theta }_k:={\displaystyle \sum _{l=1}^N{d}_{lk}}{x}_l=\{\begin{array}{lll}{x}_i-x_j\hfill & :\hfill & \text{if i is the positive end,}\hfill \\ x_j-x_i\hfill & :\hfill & \text{if j is the positive end},\hfill \end{array}& (6)\end{array}$

and f_k(θ_k) being uniformly monotone and fulfilling the sector condition f_k(θ_k)θ_k > 0 for all θ_k ∈ R.

In matrix representation this can be expressed as

$\begin{array}{ll}\dot{x}=-Df\left(D^{T}x\right)+1_N\upsilon (t)& (7)\end{array}$

with $f(\theta )={\left[f_1({\theta }_1),\cdots ,f_M({\theta }_M)\right]}^T$ and 1_N being the N-vector of ones.

We reorder the incidence matrix D as

$\begin{array}{ll}D=\left[D_I|D_E\right]& (8)\end{array}$

such that $D_I=\text{diag}(D_1^I,\dots ,D_r^I)$ correspond to the r areas.

We assume that we have an N-node network with r internally dense regions but sparsely connected. Area α has m_α neurons with α = 1, 2, ⋯, r and the vector $x^{\alpha }={\left[x_1^{\alpha }\cdots x_{m_{\alpha }}^{\alpha }\right]}^T$ contains all neural activities in area α. Then we define the slow variable (Biyik and Arcak, 2008; Chow and Kokotovic, 1985) representing an aggregate region as

$\begin{array}{ll}{y}_{\alpha }:={\displaystyle \sum _{i=1}^{m_{\alpha }}\frac{x_i^{\alpha }}{m_{\alpha }}}=\frac{1}{m_{\alpha }}{u}_{\alpha }^T{x}^{\alpha },& (9)\end{array}$

i.e., y_α is the average of the components of x^α, and in matrix representation, we obtain

$\begin{array}{ll}y=M_a^{-1}{U}^{T}x& (10)\end{array}$

with M_a = diag(m₁, m₂, ⋯, m_r) and U = diag(u₁, u₂, ⋯, u_r) where each u_α = 1_mα is an m_α-vector of 1's.

The fast variable z_α is given as the transformation of the differences between the activation of the neurons in the same region (Biyik and Arcak, 2008; Chow and Kokotovic, 1985)

$\begin{array}{ll}{z}_{\alpha }=Q_{\alpha }{x}^{\alpha }& (11)\end{array}$

where the (m_α − 1) × m_α matrix Q_α is a difference matrix. In Chow and Kokotovic (1985) Q_α is given by

$\begin{array}{ll}{Q}_{\alpha }=\left(\begin{array}{ccccc}-1& 1& 0& \cdots & 0\\ -1& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1& 0& 0& \cdots & 1\end{array}\right)& (12)\end{array}$

while in Biyik and Arcak (2008) it is chosen to have orthonormal rows with the vector 1_N as a null vector:

$\begin{array}{ll}{Q}_{\alpha }=\left(\begin{array}{ccccc}-1+(n-1)v& 1-\upsilon & -\upsilon & \cdots & -\upsilon \\ -1+(n-1)v& -\upsilon & 1-\upsilon & \cdots & -\upsilon \\ ⋮& ⋮& \cdots & \ddots & ⋮\\ -1+(n-1)v& -\upsilon & \cdots & -\upsilon & 1-\upsilon \end{array}\right)& (13)\end{array}$

where n = m_α and $\upsilon =(n-\sqrt{n})/(n(n-1))$. In this work we will use the latter version of Q_α given in (Equation 13). This has the advantage that the pseudoinverse of Q_α is simply $Q_{\alpha }^T$.

In matrix form we have the fast z (N − r)-vector with $z={\left[z_1^T{z}_2^T\cdots z_r^T\right]}^T$

$\begin{array}{ll}z=Qx& (14)\end{array}$

with Q = diag(Q₁, Q₂, ⋯, Q_r) being an (N − r) × N block diagonal matrix.

We have thus defined a new transformation of the original neural activation X into a slow and fast variable (Biyik and Arcak, 2008)

$\begin{array}{ll}\left(\begin{array}{c}y\\ z\end{array}\right)=\left(\begin{array}{c}{M}_a^{-1}{U}^T\\ Q\end{array}\right)x& (15)\end{array}$

Similarly to the matrix Q describing the connections' differences between the nodes in one area, we introduce the matrix R, an (r − 1) × r matrix also with orthonormal rows and with null vector 1_r, that models the interarea dynamics

$\begin{array}{ll}\tilde{y}=Ry& (16)\end{array}$

where $\tilde{y}$ is the area difference variable. For further computations, we define a new matrix C as $C=R{M}_a^{-1}{R}^T$.

By splitting the nonlinear vector f(θ) from Equation (7) into two components and U^TD and QD we obtain the block diagonal matrices representations

$\begin{array}{ll}{U}^{T}D=\left(\begin{array}{c}0|\Sigma \end{array}\right),f=\left(\begin{array}{c}{f}_I\\ f_E\end{array}\right),QD=\left(\begin{array}{c}A|B\end{array}\right)& (17)\end{array}$

where $\Sigma =U^T{D}_E$, A is the matrix of the internal links over all areas and B corresponds to the external links. This yields to a new singular perturbed model

$\begin{array}{l}{\dot{\tilde{y}}}_i=-CR\Sigma f_E({\Sigma }^T{R}^T\tilde{y}+B^{T}z)\\ \dot{z}=-A{f}_I(A^{T}z)-B{f}_E({\Sigma }^T{R}^T\tilde{y}+B^{T}z)& (18)\end{array}$

The above model (Equation 18) applies to large-scale brain networks for sufficiently small network parameters δ and d. The theoretical results from Biyik and Arcak (2008) state that the neurons in the same region synchronize in fast time-scale leading to a substitute aggregate neuron in the slow time-scale.

Since most brain connectivity graph networks are modeled as linear systems, we can derive a reduced model approximation over time. We assume linear connections for the brain network and replace the nonlinearity f by the identity function and set υ(t) = 0. Then $\dot{x}$ = Kx = (K^I + K^E)x with K^E being the external and K^I being the internal connection matrix.

The new linear singular perturbed system is given as

$\begin{array}{ll}\left(\begin{array}{c}\dot{y}\\ \dot{z}\end{array}\right)=\left(\begin{array}{cc}\widetilde{A_{11}}& \widetilde{A_{12}}\\ \widetilde{A_{21}}& \widetilde{A_{22}}\end{array}\right)·\left(\begin{array}{c}y\\ z\end{array}\right)& (19)\end{array}$

where

$\begin{array}{ll}\left(\begin{array}{cc}\widetilde{A_{11}}& \widetilde{A_{12}}\\ \widetilde{A_{21}}& \widetilde{A_{22}}\end{array}\right)=\left(\begin{array}{cc}G{K}^{E}U& G{K}^E{Q}^T\\ Q{K}^{E}U& Q(K^I+K^E)Q^T\end{array}\right)& (20)\end{array}$

and $G=M_a^{-1}{U}^T$. Note that in Equation (20) we chose Q using the Q_α matrices given in Equation (13). If we had chosen to use (Equation 12) instead, the rows of Q would not have been orthonormal and the Q^T in $\widetilde{A_{12}}$ and $\widetilde{A_{21}}$ would be replaced with Q⁺ = Q^T(QQT)⁻¹.

We can determine the time-scale model by defining the fast and slow time-scales

$\begin{array}{ll}{t}_f=c^{I}t\text{ and }{t}_s=\delta t_f& (21)\end{array}$

and rescaling the matrices A_ij as

$\begin{array}{l}{A}_{11}=\frac{\widetilde{A_{11}}}{c^I\delta }{A}_{12}=\frac{\widetilde{A_{12}}}{c^I\delta }\\ A_{21}=\frac{\widetilde{A_{21}}}{c^{I}d}{A}_{22}=\frac{\widetilde{A_{22}}}{c^I}& (22)\end{array}$

leading to a new system

$\begin{array}{l}\frac{dy}{d{t}_s}=A_{11}y+A_{12}z\\ \frac{\delta dz}{d{t}_s}=d{A}_{21}y+A_{22}z& (23)\end{array}$

The above results are summarized in the theorem proven in Chow and Kokotovic (1985).

Theorem 1. Chow and Kokotovic (1985) There are δ* and d* such that for all 0 < δ ≤ δ*, 0 < d ≤ d* the system in Equation (23) has r slow eigenvalues and n − r fast eigenvalues. The fast and slow subsystems are given as

$\begin{array}{l}\frac{d{y}_s}{d{t}_s}=(A_{11}-d{A}_{12}{A}_{22}^{-1}{A}_{21})y_s=A_0{y}_s,y_s(0)=y(0)\\ \frac{d{z}_f}{d{t}_f}=A_{22}{z}_f,z_f(0)=z(0)+d{A}_{22}^{-1}{A}_{21}y(0).& (24)\end{array}$

In Chow and Kokotovic (1985) a simplified formulation of the slow system was presented, the so-called aggregate system, given as

$\begin{array}{ll}{M}_a\frac{d{y}_s}{dt}=K_a{y}_s& (25)\end{array}$

with $K_a=U^T{K}_{E}U$.

We illustrate the time-scale separation properties and the aggregation procedure first on a general example and then on a structural and functional connectivity network from healthy and dementia patients. The reduced-order model will deliver important dynamic parameters derived from the state matrices of the slow and fast subsystem that are different for healthy and dementia subjects. In addition, it determines sparsity parameters based on the area and node parameters which are different at each specific stage of the disease.

2.1. Examples

Example 1. Consider an 8-node brain network that is partitioned into two areas as shown in Figure 1.

The connection matrix for Figure 1 are

$\begin{array}{ll}{K}_1^I=\left(\begin{array}{cccc}-3& 1& 1& 1\\ 1& -2& 1& 0\\ 1& 1& -3& 1\\ 1& 0& 1& -2\end{array}\right),K_2^I=\left(\begin{array}{cccc}-3& 1& 1& 1\\ 1& -3& 1& 1\\ 1& 1& -3& 1\\ 1& 1& 1& -3\end{array}\right),& (26)\end{array}$

$\begin{array}{ll}{K}_{12}^E=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right)={\left(K_{21}^E\right)}^T,& (27)\end{array}$

$\begin{array}{ll}{K}_{11}^E=\text{diag}\left(0,0,-1,-1\right),\text{ and }{K}_{22}^E=\text{diag}(-1,0,0,-1).& (28)\end{array}$

We have c^I = 2, c^E = 1, d = 0.5, γ^E = 2, δ = 0.25. Further we have for the slow subsystem according to the Equation (24)

$\begin{array}{ll}\frac{d{y}_s}{d{t}_s}=A_0{y}_s& (29)\end{array}$

with $A_0=\left(\begin{array}{cc}-0.7634& 0.7634\\ 0.7634& -0.7634\end{array}\right)$ and having the following two eigenvalues 0 and −1.5267. These two eigenvalues are close to those two of 0 and −2 of the matrix A₁₁. Thus we could show that the system can be in the long-term correctly approximated by the slow subsystem.

FIGURE 1

Figure 1. An 8-node, 2-area network.

Example 2. We apply the theoretical results on functional (FDG-PET) and structural (MRI) connectivity graphs (Ortiz et al., 2015) for control (CN), mild cognitive impairment (MCI) and Alzheimer's disease (AD) subjects. For the structural data, the connections in the graph show the inter-regional covariation of gray matter volumes in different areas while in the case of functional data, the connections do not show the correlation in activity but in the glucose uptake between different regions. (Ortiz et al., 2015) considered only 42 out of the 116 from the AAL in the frontal, parietal, occipital and temporal lobes. The nodes in the graphs represent the regions while the links show if a connection is existing between these regions or not.

Figure 2 shows the clusters found on the functional data for (A) controls, (B) MCI and (3) AD (Ortiz et al., 2015). We perform a time-scale modeling and area aggregation with two main areas on the three functional networks from Figure 2. For CN and MCI, we can apply Theorem 1, however for AD we are not able to obtain an area aggregation since the conditions in Theorem 1 are not satisfied.

FIGURE 2

Figure 2. Areas of the connectivity graph for functional data for (A) controls, (B) MCI, and (C) AD.

The results of the in-depth dynamical analysis are shown in Table 1. The controls show smaller node and area parameters than the MCIs. But most importantly, both the exact as well as the rigid aggregate model as shown in Equation (25) show smaller eigenvalues for the controls than the MCIs.

TABLE 1

Table 1. Area aggregation parameters and time-scale modeling for functional connectivity graphs.

Similarly, we perform a time-scale modeling and area aggregation with two main areas on the three structural networks from Figure 3. For CN and MCI, we can apply Theorem 1, however for AD we are not able to obtain an area aggregation again because the conditions are not met.

FIGURE 3

Figure 3. Areas of the connectivity graph for structural data for (A) controls, (B) MCI, and (C) AD.

Figure 3 shows the clusters found on the structural data for (A) controls, (B) MCI and (3) AD (Ortiz et al., 2015). The results of the in-depth dynamical analysis are shown in Table 2. Similarly, the controls show smaller area parameters than the MCIs and smaller eigenvalues. However, the node parameter is larger in the case of controls.

TABLE 2

Table 2. Area aggregation parameters and time-scale modeling for structural connectivity graphs.

While the results obtained through static graph analysis revealed the loss of strong connections in AD and MCI compared to CN, the dynamic graph analysis reveals different slow modes between MCI and CN. The CN have smaller eigenvalues than the MCI for both functional and structural data and those eigenvalues remain operative. The contribution of the larger eigenvalues over time decreases quickly. The range of the eigenvalues for each subject represents an important biomarker for disease prediction. It is worth noting that an area aggregation is not possible for the ADs in both structural and functional networks since the conditions of Theorem 1 are not fulfilled. By providing an area and node parameter, we are able to add to the dynamic biomarkers additional static graph descriptors. The reduced-order time scale modeling provides a change in the sparsity pattern for MCI and AD patients compared to the CN, and shows higher values for the node and area parameters.

It is worth mentioning that the slow or aggregate variable represents in Markov chains the probability for a Markov process to be in a group of states.

3. Pinning Observability

Brain connectivity networks have lots of nodes and connections between them, and it is effectively impossible to observe the states of all nodes such that the network can be recovered by reaching synchronization between the original network and a reconstructed network. Pinning observability is introduced as a new technique to reduce the number of observable nodes, and at the same time to be able to recover the states of other nodes. We give a new criterion for synchronization via pinning observability for nonlinear brain networks and derive decoupled conditions for pinning observability for brain connectivity networks.

We consider in the following the general neural network equation describing the temporal evolution of the neural activation states for the ith neuron of an N–neuron network:

$\begin{array}{ll}{\dot{x}}_i=A_i{x}_i+{\displaystyle \sum _{j=1}^n{\tilde{d}}_{ij}}\overline{f}(x_j)i=1,\mathrm{\dots },N& (30)\end{array}$

where $x_i={(x_{i1}(t),\cdots ,x_{in}(t))}^T\in R^n$ is the state vector of node i, $A_i\in R^{n\times n}$ is a matrix, and $\overline{f}(x_i)$ is the neuron's output. ${\tilde{d}}_{ij}$ represents a synaptic connection parameter between the ith neuron and the jth neuron and is defined as the matrix $\tilde{D}=({\tilde{d}}_{ij})$.

Pinning observability is applied only to a small number l of nodes and thus we obtain a pinning observable network supposing these first l nodes are selected:

$\begin{array}{ll}{\dot{\tilde{x}}}_i=A_i{\tilde{x}}_i+{\displaystyle \sum _{j=1}^n{\tilde{d}}_{ij}}\overline{f}({\tilde{x}}_j)+u_{i}i=1,\mathrm{\dots },N& (31)\end{array}$

where

$\begin{array}{ll}{g}_i(x_j)={\displaystyle \sum _{j=1}^N{\tilde{d}}_{ij}}\overline{f}(x_j)& (32)\end{array}$

and

$\begin{array}{ll}{u}_i=-d_i({\tilde{x}}_i(t)-x_i(t))1,\cdots ,l& (33)\end{array}$

are n-dimensional linear feedback controllers with the control gains d_i > 0, i = 1, ⋯, l and d_i = 0 for i = l + 1, ⋯, N.

By subtracting the pinned system Equation (31) and the original system Equation (30), we obtain the error dynamical system

$\begin{array}{ll}{\dot{e}}_i(t)=A_i{e}_i(t)+g_i(\widetilde{x_1}(t),\cdots ,\widetilde{x_N}(t))-g_i(x_1(t),\cdots ,x_N(t))& (34)\end{array}$

and the error signal

$\begin{array}{ll}{e}_i(t)=\widetilde{x_i}(t)-x_i(t).& (35)\end{array}$

The following assumptions and lemmas are needed to obtain the main result of the paper.

Assumption A Assume a scale-free graph network and let $\tilde{D}={({\tilde{d}}_{ij})}_{N\times N}$ be the coupling configuration matrix. Then there exists a symmetric matrix Γ such that

$\begin{array}{l}{\displaystyle \sum _{i=1}^N{(\widehat{x_i}-\widehat{w_i})}^T}(g_i(\widehat{x_1}(t),\cdots ,\widehat{x_N}(T))-g_i(\widehat{w_1}(t),\cdots ,\widehat{w_N}(t)))\\ \le {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{l}_{ij}{(\widehat{x_i}-\widehat{w_i})}^T}}\Gamma \left(\widehat{x_j}-\widehat{w_j}\right).& (36)\end{array}$

The following two lemmas will be useful.

Lemma 4. For any vectors x, y ∈ R^N and positive definite matrix G ∈ R^N×N, the following matrix inequality holds

$\begin{array}{ll}2x^{T}y\le x^{T}Gx+y^T{G}^{-1}y.& (37)\end{array}$

Lemma 5. Schur complement (Boyd et al., 1994). The following linear matrix inequality (LMI)

$\begin{array}{ll}\left(\begin{array}{cc}Q(x)& S(x)\\ S{(x)}^T& R(x)\end{array}\right)>0& (38)\end{array}$

with Q(x) = Q(x)^T, R(x) = R(x)^T is equivalent to either of the following equivalent conditions

$\begin{array}{l}(1)Q(x)>0,R(x)-S{\left(x\right)}^{T}Q{\left(x\right)}^{-1}S(x)>0\\ (2)R(x)>0,Q(x)-S(x)R{\left(x\right)}^{-1}S{\left(x\right)}^T>0& (39)\end{array}$

We can now state a result concerning the two networks Equations (30) and (31).

Theorem 6. Under Assumption A the two networks Equations (30) and (31) are globally synchronized if the following condition is satisfied

$\begin{array}{ll}(A+A^T)/2-D\otimes I_n+(L\otimes \Gamma +L^T\otimes \Gamma )/2<0& (40)\end{array}$

where A = diag(A₁, …, A_N), D = diag(d₁, …, d_l, 0, …, 0), $\tilde{D}={({\tilde{d}}_{ij})}_{N\times N}$ and I_n being an n-dimensional identity matrix.

Proof. Consider the Lyapunov function candidate

$\begin{array}{ll}V(t)=\frac{1}{2}{\displaystyle \sum _{i=1}^N{e}_i^T}(t)·e_i(t).& (41)\end{array}$

We take the derivative of V(t) along the trajectories of Equation (34) and obtain

$\begin{array}{l}\dot{V}(t)={\displaystyle \sum _{i=1}^N{e}_i^T}(t)·\stackrel{.}{e_i}(t)\\ ={\displaystyle \sum _{i=1}^N{e}_i^T}(t)(A_i{e}_i(t)+g_i(\widetilde{x_1}(t),\cdots ,\widetilde{x_N}(t))\\ -g_i(x_1(t),\cdots ,x_N(t))).& (42)\end{array}$

From Assumption A we have

$\begin{array}{l}{\displaystyle \sum _{i=1}^N{e}_i^T}(t)\left(g_i(\widetilde{x_1},\cdots \widetilde{x_N},t)-g_i(x_1,\cdots x_N,t)\right)\\ \le {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{\tilde{d}}_{ij}}}{e}_i^T(t)\Gamma e_j(t).& (43)\end{array}$

Substituting Equation (43) into Equation (42), we have

$\begin{array}{l}\dot{V}(t)\le {\displaystyle \sum _{i=1}^N{e}_i^T}(t)·(A_i-d_i{I}_n)·e_i(t)\\ +{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{l}_{ij}}}{e}_i^T(t)\Gamma e_j(t)\\ =e^T(t)\left((A+A^T)/2-D\otimes I_n+\tilde{D}\otimes \Gamma \right)e(t).& (44)\end{array}$

As shown in Yu et al. (2014), it is possible to identify the number of nodes that can be observed. The following corollary gives a condition to check for each node whether it can be controlled or not without involving the other nodes.

Synchronized if there exists a constant c > 0 such that the following condition is satisfied

Corollary 7. Under the Assumption A, the two networks Equations (30) and (31) are globally synchronized if there exists a constant c > 0 such that the following condition is satisfied

$\begin{array}{l}(A_i+A_i^T)/2-d_i{I}_n+{\tilde{d}}_{ii}\Gamma +\frac{1}{2}{c}_i{\displaystyle \sum _{j=1,j\ne i}^N|}{\tilde{d}}_{ij}|\Gamma {\Gamma }^T\\ +\frac{1}{2}{\displaystyle \sum _{j=1,i\ne j}^N\frac{1}{c_j}}|{\tilde{d}}_{ji}|I_n<0& (45)\end{array}$

The proof is based on using the same Lyapunov function as in the above Theorem.

Thus we are able to give a much simpler condition with this fixed constant c: if d_i = 0 and the condition given in the below equation is satisfied for a node i, then the node may not be controlled. Otherwise, if this condition is not satisfied for node i, then the node can be controlled.

$\begin{array}{l}{\Xi }_i(c)=\left(\begin{array}{cc}(A_i+A_i^T)/2-d_i{I}_n+l_{ii}{I}_n+\frac{c}{2}{\sum }_{j=1,j\ne i}^N{l}_{ij}{I}_n& \sqrt{{\sum }_{j=1,i\ne j}^N|l_{ij}|}{I}_n\\ \sqrt{{\sum }_{j=1,i\ne j}^N|l_{ij}|}{I}_n& -2c{I}_n\end{array}\right)<0,\\ i=1,2,\cdots ,N.\end{array}$

The above equation shows a simple modality to find the few nodes that have to be controlled such that synchronization is reached.

The following examples are given to show the application of the theory to a simple neural network. We illustrate the concept of pinning observability first on a general example describing a network with nonlinear coupling and then on a network processing olphactoric stimuli with linear coupling. Both examples employ a simplified condition Equation 46 for pinning observability developed in the above theoretical framework.

3.1. Examples

Example 3. In the following we give a numerical example to elucidate the theoretical results in pinning observability.

The simplest example consists of three neurons. Let N = 3, a₁ = 2, a₂ = 1.6, a₃ = 0.5, c = 1 $\tilde{D}=\left(\begin{array}{ccc}0.7& -0.2& -0.1\\ -0.2& 0.5& -0.3\\ -0.1& -0.3& 0.8\end{array}\right)$, and let the nonlinearity f be a sigmoid function with K = 1.

The condition to be tested for each neuron is given by Theorem 6.

In case the condition is satisfied for a fixed parameter c and d_i = 0, then the neuron i can not be controlled. Node 3 is the one that needs to be controlled and d₁ = 1.5.

Example 4. The processing of olphactoric stimuli was shown to involve a network of cortical and subcortical regions (Nigri et al., 2013). The resulting connectivity graph in fMRI studies shows a small number of hubs in this network suggesting that they have a predominant role in information gathering. Applying the theory of pinning observability on this networks reveals that there are few nodes in the graphs that can be controlled, and that these nodes are surprisingly not hubs but have sparse connections. Figure 4 shows the distribution of these controlling nodes in the olphactoric network.

FIGURE 4

Figure 4. Nodes that can be controlled in the connectivity graph during olphactoric stimulation marked as a donut (Nigri et al., 2013).

4. Discussion

In this paper, we introduced some novel dynamic graph theory techniques to the analysis of the dynamical behavior of brain connectivity networks. We applied the new concepts of time-scale modeling for sparse networks and pinning observability on brain networks of heterogeneous architecture. We considered graphs with densely linked nodes in an area but with sparse connections between these areas. We have shown that the nodes in a dense area synchronize on the fast time-scale while the dense areas become aggregate nodes on the slow time scale. For the time scale modeling, we have derived new local models that describe the fast and slow dynamics of brain connectivity networks assuming linear connections between the nodes. This new paradigms provides us with important disease descriptors showing changes over the disease trajectory such as the modes of the dynamic system and the sparsity patterns.

Observing a small number of nodes in brain connectivity network and recovering the states of the others is of major interest in a large-scale network. This was achieved through pinning observability. We formulated a new criterion for synchronization for nonlinear brain networks and derived decoupled simplified conditions for determining the small number of observable states in the network. Examples are given to elucidate the theoretical results for both new concepts.

While static graph theory shows the changes in graph measures at certain points in time and the differences between disease and normal control, the derived results may have important implication for understanding and controlling the evolution of neurodegenerative diseases that may further lead to better therapeutic interventions. Thus by describing the dynamic of the aggregate areas and the resulting time-scale modeling and determining the observable nodes in a brain network, a new research avenue is opened that allows more detailed study of the differences between disease and healthy groups and which provides a wider variety of characteristic parameters over time for those groups.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication. AM: manuscript writing, research and revision, RR, IAI, UM, ML, AS, and KP: research and revision.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We wish to thank Elena Gonzalez Avidad for her design contribution in this work. We acknowledge funding through a Fulbright Award and Marie Sklodowska-Curie actions (MSCA-IF-GF-656886).

References

Biyik, E., and Arcak, M. (2008). Area aggregation and time-scale modeling for sparse nonlinear networks. Sys. Cont. Lett. 57, 142–149. doi: 10.1016/j.sysconle.2007.08.003

CrossRef Full Text | Google Scholar

Boyd, S., El Ghaoui, L., Feron, E., and Balakrishnan, V. (1994). Linear Matrix Inequalities in System and Control Theory. SIAM.

Google Scholar

Chow, J. H. and Kokotovic. P. (1985). Time scale modeling of sparse dynamic networks. IEEE Trans. Automat. Cont. 30, 714–722. doi: 10.1109/TAC.1985.1104055

CrossRef Full Text | Google Scholar

Chen, T., Liu, X., and Lu, W. (2007). Pinning complex networks by a single controller. IEEE Trans. Circ. Syst. I 54, 1317–1326. doi: 10.1109/TCSI.2007.895383

CrossRef Full Text | Google Scholar

Fornito, A., Zalesky, A., Pantelis, C., and Bullmore, E. (2012). Schizophrenia, neuroimaging and connectomics. Neuroimage 62, 2296–2314. doi: 10.1016/j.neuroimage.2011.12.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Giessing, C., and Thiel, C. (2012). Pro-cognitive drug effects modulate functional brain network organization. Front. Behav. Neurosci. 6, 1–9. doi: 10.3389/fnbeh.2012.00053

PubMed Abstract | CrossRef Full Text | Google Scholar

Lall, S., Marsden, J., and Glavaski, S. (2002). A subspace approach to balanced truncation for model reduction of nonlinear control systems. Int. J. Robust Nonlin. Control. 12, 519–535. doi: 10.1002/rnc.657

CrossRef Full Text | Google Scholar

Liu, Y. Y., Slotine, J. J., and Barabasi, A. (2011). Controllability of complex networks. Nature 473, 167–173. doi: 10.1038/nature10011

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer-Bäse, A. (2008). Gene regulatory networks simplified by nonlinear balanced truncation. SPIE Proc. Ser. 6979:69790C. doi: 10.1117/12.777292

CrossRef Full Text | Google Scholar

Moore, B. (1981). Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Trans. Auto. Cont. 26, 17–32. doi: 10.1109/TAC.1981.1102568

CrossRef Full Text | Google Scholar

Nigri, A., Ferraro, S., D'incerti, L., Critchley, H., Bruzzone, M., and Minati, L. (2013). Connectivity of the amygdala, piriform, and orbitofrontal cortex during olfactory stimulation: a functional MRI study. Neuroreport 24, 171–175. doi: 10.1097/WNR.0b013e32835d5d2b

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortiz, A., Munilla, J., Alvarez-Illan, I., Gorriz, M., and Ramirez, J. (2015). Exploratory graphical models of functional and structural connectivity patterns for Alzheimer's disease diagnosis. Front. Comput. Neurosci. 9:132. doi: 10.3389/fncom.2015.00132

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Q., and Cao, J. (2010). On pinning synchronization of directed and undirected complex dynamical networks. IEEE Trans. Circ. Sys. I 57, 672–680. doi: 10.1109/TCSI.2009.2024971

CrossRef Full Text | Google Scholar

Scherpen, J. (1993). Balancing for nonlinear systems. Sys. Cont. Lett. 21, 143–153. doi: 10.1016/0167-6911(93)90117-O

CrossRef Full Text | Google Scholar

Tahmassebi, A., Pinker-Domenig, K., Wengert, G., Lobbes, M., Stadlbauer, A., Romero, F. J., et al. (2017). “Dynamical graph theory networks techniques for the analysis of sparse connectivity networks in dementia,” in Proceeding of SPIE 10216, Smart Biomedical and Physiological Sensor Technology XIV (Anaheim, CA). doi: 10.1117/12.2263555

CrossRef Full Text | Google Scholar

Tang, Y., Wang, Z., Gao, H., Swift, S., and Kurths, J. (2012). A constrained evolutionary computation method for detecting controlling regions of cortical networks. IEEE/ACM Trans. Comput. Biol. Bioinformat. 9, 1569–1581. doi: 10.1109/TCBB.2012.124

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, W., Wen, G., Lü, L., and Yu, X. (2014). Pinning observability in complex networks. IET Cont. Theory Appl. 8, 2136–2144. doi: 10.1049/iet-cta.2013.0969

CrossRef Full Text | Google Scholar

Zeng, L., Shen, H., Liu, L., Wang, L., Li, B., Fang, P., et al. (2012). Identifying major depression using whole-brain functional connectivity: a multivariate pattern analysis. Brain 1498, 1498–1507. doi: 10.1093/brain/aws059

CrossRef Full Text | Google Scholar