The control of brain network dynamics across diverse scales of space and time

TL;DR

This study investigates how different brain regions control neural dynamics across various spatial and temporal scales, revealing that regions controlling easy transitions operate on short scales and slow times, while those controlling difficult transitions operate on long scales and fast times.

Contribution

It provides a novel analysis linking brain control profiles with spatial and temporal dynamics using diffusion tensor imaging and Laplacian eigenspectrum.

Findings

Control of easy transitions linked to short spatial and slow temporal scales.

Control of difficult transitions linked to long spatial and fast temporal scales.

Linear models explain activity propagation and control profiles across brain regions.

Abstract

The human brain is composed of distinct regions that are each associated with particular functions and distinct propensities for the control of neural dynamics. However, the relation between these functions and control profiles is poorly understood, as is the variation in this relation across diverse scales of space and time. Here we probe the relation between control and dynamics in brain networks constructed from diffusion tensor imaging data in a large community based sample of young adults. Specifically, we probe the control properties of each brain region and investigate their relationship with dynamics across various spatial scales using the Laplacian eigenspectrum. In addition, through analysis of regional modal controllability and partitioning of modes, we determine whether the associated dynamics are fast or slow, as well as whether they are alternating or monotone. We find that brain regions that facilitate the control of energetically easy transitions are associated with activity on short length scales and slow time scales. Conversely, brain regions that facilitate control of difficult transitions are associated with activity on long length scales and fast time scales. Built on linear dynamical models, our results offer parsimonious explanations for the activity propagation and network control profiles supported by regions of differing neuroanatomical structure.

Tables

Table 1. Table 7: (Table VII in the manuscript.) We calculate primary network metrics for both the brain network and a Barabasi-Albert graph with the same number of nodes ( N = 234 𝑁 234 N=234 ) and the same normalization (see Section I.B): the metrics of clustering coefficient, characteristic path length, modularity, and coreness, respectively Rubinov and Sporns ( 2010 ) . Uncertainty measurements are standard errors; modularity (a heuristic) is calculated over 1000 trials maximizing a common modularity quality function Newman ( 2004 ) .

	Brain network	Barabasi-Albert graph
Clustering coefficient Watts and Strogatz (1998)	$\left(5.5\times {10}^{-4}\right)\pm \left(1.5\times {10}^{-4}\right)$	$\left(4.3\times {10}^{-3}\right)\pm \left(1.3\times {10}^{-5}\right)$
Characteristic path length	$0.75$	$0.75$
(binary) Watts and Strogatz (1998)
Characteristic path length	$2.1\times {10}^{-3}$	$4.1\times {10}^{-3}$
(weighted) Rubinov and Sporns (2010)
Modularity Newman (2004)	$0.49\pm \left(1.6\times {10}^{-4}\right)$	$0.021\pm \left(2.5\times {10}^{-5}\right)$
Coreness Borgatti and Everett (2000)	$0.37$	$0.22$

Equations

$$\mathbf{x}(t+1)=\mathbf{Ax}(t)+{\mathbf{B}}_{\mathcal{K}}{\mathbf{u}}_{\mathcal{K}}(t),$$

$$B_{\mathcal{K}}=\left[\begin{array}{ccc}{e}_{k_1}& \cdots & e_{k_m}\end{array}\right],$$

$${\mathbf{W}}_{\mathcal{K}}=\sum _{\tau =0}^{\infty }{\mathbf{A}}^{\tau }{\mathbf{B}}_{\mathcal{K}}{\mathbf{B}}_{\mathcal{K}}^{\mathsf{T}}{\mathbf{A}}^{\tau }.$$

$$\text{Trace}({\mathbf{W}}_{\mathcal{K}})$$

$$=\sum _{\tau =0}^{\infty }\text{Trace}\left({\mathbf{A}}^{\tau }{\mathbf{B}}_{\mathcal{K}}{\mathbf{B}}_{\mathcal{K}}^{\mathsf{T}}{\mathbf{A}}^{\tau }\right)$$

$$=\sum _{\tau =0}^{\infty }\text{Trace}\left({\mathbf{A}}^{2\tau }{\mathbf{B}}_{\mathcal{K}}{\mathbf{B}}_{\mathcal{K}}^{\mathsf{T}}\right)$$

$$=\left(\sum _{\tau =0}^{\infty }{\mathbf{A}}^{2\tau }\right)\left({\mathbf{B}}_{\mathcal{K}}{\mathbf{B}}_{\mathcal{K}}^{\mathsf{T}}\right)$$

$$={\left(\mathbf{I}-{\mathbf{A}}^{2\tau }\right)}^{-1}\left({\mathbf{B}}_{\mathcal{K}}{\mathbf{B}}_{\mathcal{K}}^{\mathsf{T}}\right).$$

$${\varphi }_i=\sum _j(1-{\xi }_j^2(A))v_{ij}^2,$$

Full text

The brain is an inherently networked system that displays incredibly rich and complex dynamics Bassett and Sporns (2017); Breakspear (2017). Building accurate models of those dynamics remains a key challenge that is fundamental to the field of neuroscience, with the potential to inform personalized medicine by predicting a patient’s disease progression and response to therapy Bansal et al. (2018); Proix et al. (2017); Falcon et al. (2016); Schirner et al. (2015). Efforts to build such models necessarily depend on the development of interdisciplinary approaches informed by dynamical systems theory, statistical physics, and network science as well as substantial knowledge of the intricacies of the underlying biology Bassett et al. (2018). Because of the multiscale nature of the system Breakspear and Stam (2005); Betzel and Bassett (2017), there are some regimes of function or modalities of measurement whose dynamics are well-fit by one mathematical model and others whose dynamics are well-fit by another mathematical model Friston (2005, 2008); Deco et al. (2008); Breakspear (2017). On the one hand, large-scale measurements of brain activity can display a wave-like nature and other characteristics consistent with linear, diffusive dynamics instantiated on a physically embedded system Honey et al. (2009); Fernández Galán (2008); Raj et al. (2012); Aquino et al. (2012); Viventi et al. (2011). On the other hand, fine time scale measurements of brain activity can display oscillatory characteristics consistent with nonlinear, synchronization dynamics thought to support distributed communication and computation Kopell et al. (2014); Fries (2015); Palmigiano et al. (2017); Kirst et al. (2016).

Evidence for the first type of dynamics comes, for example, from studies of both large- and small-scale neuronal processes whose intrinsic or endogeneous activity can be partially explained by simple linear models Honey et al. (2009); Fernández Galán (2008). Moreover, wave-like activity has been observed using functional magentic resonance imaging (fMRI) in the visual cortex, where traveling wave responses to visual stimuli are observed to propagate out from the fovea Aquino et al. (2012). Notably, such traveling waves with various timescales have been observed across different regions of the brain Muller et al. (2018), and it has been posited that the direction and wavelength of the wave could encode information transmitted between regions Heitmann et al. (2013). In addition to the outward propagation of activity, spiral waves also occur frequently in the neocortex in vivo, both during pharmacologically induced oscillations and during sleep-like states Huang et al. (2010), while seizures may manifest as recurrent spiral waves that propagate in the neocortex Viventi et al. (2011). Such multiplicity of phenomena motivates the use of dynamical models that emphasize the spatially embedded nature of brain networks Stiso and Bassett (2018), and the extended versus transient responses that embedding can support Chaudhuri et al. (2015). An example is the use of finite element modeling to simulate the propagation of brain stimulation Cancelli et al. (2016); Lopes et al. (2012), where effects are distributed on a localized part of brain tissue Metwally et al. (2012). Another example is that of diffusive models, where network diffusion over long time scales has been used to model disease progression in the brain for dementia Raj et al. (2012), or to predict longitudinal patterns of atrophy and metabolism in Alzheimer’s disease Raj et al. (2015).

Evidence for the second type of dynamics comes from studies of both large- and small-scale ensembles that produce rhythmic or oscillatory activity Kopell et al. (2014); Fries (2015); Jones (2016). For example, scalp electrodes in humans can be used to measure rhythms of certain frequencies, which in turn have been associated with different cognitive processes and behavioral responses Palmigiano et al. (2017); Kirst et al. (2016). For example, attentional control via the inhibition of behaviorally irrelevant stimuli and motor responses has been associated with the synchronization of $\alpha$ and $\beta$ rhythms between right inferior frontal and primary sensory neocortex Sacchet et al. (2015). Similarly, higher frequency $\gamma$ rhythms in local synchronization have been observed during visual responses, while lower frequency $\beta$ rhythms reflecting coherence over parietal and temporal cortices have been observed during activities that required more multi-modal sensory integration von Stein et al. (1999). Initial efforts sought to explain such phenomena using a simple neuronal model with conduction delays, where excitatory units describe pyramidal cells and inhibitory units describe interneurons Kopell et al. (2000). More recent work suggests that the $\beta$ rhythms specifically can emerge from the integration of nearly synchronous bursts of excitatory synaptic drive targeting proximal and distal dendrites of pyramidal neurons, where the defining feature of a $\beta$ event was a strong distal drive to supragranular and infragranular layers of cortex that lasted one $\beta$ period Sherman et al. (2016). Notably, the laminar architecture of the neocortical network has also proven critical for explaining $\gamma$ rhythms Lee and Jones (2013).

Although empirical studies find pervasive evidence of these two distinct types of dynamics, there is little theoretical framework to model how the spatiotemporal characteristics of these dynamics can inform the control properties of specific brain regions, which are interconnected by a fixed structural network. To address this question, recent efforts have begun to study how such distinct dynamics can be guided or controlled by local energy input Tang and Bassett (2018). For the control of those features of the dynamics that are well-approximated by linear systems theory, efforts have stipulated a linear model of the dynamics dependent upon the structural network implicit in white matter connectivity Gu et al. (2015); Taylor et al. (2015). This framework of network control has been useful to understand brain anatomy and regional function across multiple species including the nematode C. elegans Yan et al. (2017a), the fly Drosophila Kim et al. (2018), the mouse Kim et al. (2018), the macaque Gu et al. (2015), and the human Tang et al. (2017). Network control theory has also been used to explore the control profiles of brain regions based on their network connectivity Pasqualetti et al. (2014a), and hence to map control propensities on to canonically understood cognitive functions Gu et al. (2015); Tang et al. (2017); Cornblath et al. (2018) and their alterations in psychiatric disorders Jeganathan et al. (2018). For the control of those features of the dynamics that require nonlinear systems theory, network control is more difficult Motter (2015); Muldoon et al. (2016), but some initial work has begun to assess the utility of control principles for seizure abatement Taylor et al. (2015) and to determine routes along which information can be propagated by transient synchrony Palmigiano et al. (2017); Kirst et al. (2016). In extending these ongoing efforts, a key challenge remains to build models that can predict not only the generic control properties of a system, but also the length scales over which control-induced dynamics will propagate, and to link these length scales to the time scales of dynamical transitions or control of activity enacted by regional drivers.

Here, we perform initial numerical experiments to address this challenge by relating new developments in the linear control of regional activity in brain networks and their predicted time scales, with the length scales implicit in nonlinear dynamical models of inter-regional synchronization. We use the framework of network control Liu and Barabási (2016) to probe the relation between control and dynamics in brain networks constructed from diffusion tensor imaging data in a young adult subset of the Philadelphia Neurodevelopmental Cohort, a large community based sample of youth Satterthwaite et al. (2016); Calkins et al. (2015). Use of this large sample allows us to probe relations between control and dynamics in an ensemble of brain networks in which control propensities vary extensively Tang et al. (2017). The paper is organized as follows. Following a description of the formalism in Sec. I, we turn in Sec. II to a review of network controllability metrics and their predictions regarding the control profiles of certain brain regions (see the Appendix for details regarding brain network construction from neuroimaging data). In Sec. III, we consider metrics for the characterization of spatial scales of synchronization in network dynamics Atasoy et al. (2016) and for the control of such synchronization He et al. (2014). In Sec. IV, we consider metrics for the characterization of temporal scales of control in network dynamics, specifically assessing the role of modal controllers in driving transient versus extended temporal modes of brain activity Pasqualetti et al. (2014a). In Sec. V, we compare the relation between control of time scales and control of length scales that we observe in brain networks to those that we observe in networks constructed via the rules of preferential attachment, gaining insight into the dependence of spatiotemporal control on the underlying topology. In Sec. VI, we place our findings within the broader context of empirical results from neuroimaging experiments. We also discuss the implications of our findings for our understanding of how various controllers are associated with distinct types of dynamics, and in turn how distinct dynamics are associated with different regional controllers. We close with a brief discussion of potential future directions.

I Mathematical Framework

In this section, we describe a formalism for building network models of brain architecture and states, and for studying network controllability of brain dynamics.

I.1 Network models of brain architecture and states

Here, we study 190 brain networks describing the white matter connectivity of youth in the Philadelphia Neurodevelopmental Cohort. For each network composed of $N=234$ cortical and subcortical brain regions, $\mathbf{A}\in\mathbb{R}^{N\times N}$ is a symmetric and weighted adjacency matrix whose elements indicate the number of white matter streamlines connecting two different brain regions indexed with $(i,j)$. For further details regarding the construction of brain networks from neuroimaging data see the Appendix, and for prior work in this data set see Tang et al. (2017); Baum et al. (2017, 2018); Cornblath et al. (2018).

More formally, a networked system such as the brain can be represented by the graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ and $\mathcal{E}$ are the vertex and edge sets, respectively. Let $a_{ij}$ be the weight associated with the edge $(i,j)\in\mathcal{E}$, and define the weighted adjacency matrix of $\mathcal{G}$ as $A=[a_{ij}]$, where $a_{ij}=0$ whenever $(i,j)\not\in\mathcal{E}$. We associate a real value (state) with each node, collect the node states into a vector (network state), and define the map $x:\mathbb{N}_{\geq 0}\rightarrow\mathbb{R}^{n}$ to describe the evolution (network dynamics) of the network state over time.

Next we seek to define an appropriate equation of state. Based on prior work demonstrating the utility of simple linear models in predicting intrinsic brain dynamics across spatial scales Honey et al. (2009); Fernández Galán (2008), we employ a simplified noise-free linear discrete-time and time-invariant network model Gu et al. (2015):

[TABLE]

where $\mathbf{x}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{N}$ describes the state (i.e., a measure of the electrical charge, oxygen level, or firing rate) of brain regions over time, and $\mathbf{A}\in\mathbb{R}^{N\times N}$ is the structural connectivity described in the previous section. Hence the size of the vector $\mathbf{x}$ is given by $N$, and the value of $\mathbf{x}$ describes the brain activity of that region.

I.2 Network controllability of brain dynamics

To assess controllability of this networked system, we must first ensure its stability and define input points for the injection of control energy. We note first that the diagonal elements of the matrix $\mathbf{A}$ satisfy $A_{ii}=0$. Next we note that to assure Schur stability, we divide the matrix by $1+\xi_{0}(\mathbf{A})$, where $\xi_{0}(\mathbf{A})$ is the largest eigenvalue of $\mathbf{A}$. The input matrix $\mathbf{B}_{\mathcal{K}}$ identifies the control points $\mathcal{K}$ in the brain, where $\mathcal{K}=\{k_{1},\dots,k_{m}\}$ and

[TABLE]

and $e_{i}$ denotes the $i$-th canonical vector of dimension $N$. The input $\mathbf{u}_{\mathcal{K}}:\mathbb{N}\rightarrow\mathbb{R}^{|\mathcal{K}|}$ denotes the control strategy, i.e. the control signal injected into the network via the nodes $\mathcal{K}$. More precisely, it is a map from the set of natural numbers (representing the time instant), and the set of real vectors of length $|\mathcal{K}|$ (representing the actual control vector).

We can now study the controllability of this dynamical system, which refers to the possibility of driving the state of the system to a specific target state by means of an external control input Kalman et al. (1963). Classic results in control theory ensure that controllability of the network (1) from the set of nodes $\mathcal{K}$ is equivalent to the controllability Gramian $\mathbf{W}_{\mathcal{K}}$ being invertible, where

[TABLE]

Consistent with Gu et al. (2015); Tang et al. (2017); Cornblath et al. (2018); Jeganathan et al. (2018), we use this framework to choose control nodes one at a time, and thus the input matrix $\mathbf{B}_{\mathcal{K}}$ in fact reduces to a one-dimensional vector, e.g., $\mathbf{B}_{\mathcal{K}}=\begin{pmatrix}1&0&0&...\end{pmatrix}^{T}$ when the first brain region is the control node. In this case, $\mathcal{K}$ simply describes this control node, i.e. the controllability Gramian can be indexed by the $i$-th control node that it describes: $\mathbf{W}_{i}$.

II Controllability metrics for brain networks

Within the network controllability framework, we study two different control strategies that describe the ability to move the network into different states defined as patterns of regional activity. While the brain certainly displays non-linear activity, we study a linear model of this more complex nonlinear process which is valid in a neighborhood of the linearization point (see Fig. 1a). Modeling of brain activity in large-scale regional networks shows that the linear approximation provides fair explanatory power of resting state fMRI BOLD data Honey et al. (2009). Further, studies of this controllability framework using non-linear oscillators connected with coupling constants estimated from white matter tracts shows a good overlap with the linear approximation Muldoon et al. (2016).

The energy landscape of the linear model is quadratic, with at most one energy minimum (at the origin), corresponding to the only state that can be maintained without cost. Note that this origin of the linear system corresponds to the linearization point, i.e. the equilibrium configuration of the nonlinear model from which the linear model is obtained. Using the linear model described in the previous section, we can introduce two controllability metrics. Intuitively, average controllability describes the ease of transition to many states nearby on an energy landscape, while modal controllability describes the ease of transition to a state distant on this landscape (see Fig. 1a).

Average controllability of a network equals the average input energy from a set of control nodes and over all possible target states. As a known result, average input energy is proportional to $\text{Trace}(\mathbf{W}_{\mathcal{K}}^{-1})$, the trace of the inverse of the controllability Gramian. Instead and consistent with Gu et al. (2015); Tang et al. (2017); Cornblath et al. (2018), we adopt $\text{Trace}(\mathbf{W}_{\mathcal{K}})$ as a measure of average controllability for two main reasons: first, $\text{Trace}(\mathbf{W}_{\mathcal{K}}^{-1})$ and $\text{Trace}(\mathbf{W}_{\mathcal{K}})$ satisfy a relation of inverse proportionality, so that the information obtained from the two metrics are correlated with one another and, second, $\mathbf{W}_{\mathcal{K}}$ is typically very ill-conditioned even for coarse network resolutions, so that $\text{Trace}(\mathbf{W}_{\mathcal{K}}^{-1})$ cannot be accurately computed even for small brain networks. It should be noted that $\text{Trace}(\mathbf{W}_{\mathcal{K}})$ encodes a well-defined control metric, namely the energy of the network impulse response or, equivalently, the network $H_{2}$ norm Kailath (1980). As discussed above, when a brain region $i$ serves as a control node, the resulting Gramian can be indexed as $\mathbf{W}_{i}$, in order to compute the regional average controllability.

To calculate the trace of the Gramian, we avoid the numerical difficulties associated with solving the Lyapunov equation Hewer and Kenney (1988); Gahinet et al. (1990); Hammarling (1982); Sorensen and Zhou (2002). Given the stability of our matrix $\mathbf{A}$, we can use the following:

[TABLE]

Here we used some properties of our system: that ${\mathcal{K}}$ contains only one node and $\mathbf{A}$ is symmetric.

Modal controllability refers to the ability of a node to control each time-evolving mode of a dynamical network, and can be used to identify states that are difficult to control from a set of control nodes. Modal controllability is computed from the eigenvector matrix $V=[v_{ij}]$ of the network adjacency matrix $\mathbf{A}$. By extension from the PBH test Kailath (1980), if the entry $v_{ij}$ is small, then the $j$-th mode is poorly controllable from node $i$. Following Pasqualetti et al. (2014b), we define

[TABLE]

as a scaled measure of the controllability of all $N$ modes $\xi_{0}(\mathbf{A}),\dots,\xi_{N-1}(\mathbf{A})$ from the brain region $i$ – allowing the computation of regional modal controllability. Regions with high modal controllability are able to control all of the dynamic modes of the network, and hence to drive the dynamics towards hard-to-reach configurations.

We study these metrics in the 190 structural brain networks derived from diffusion tensor imaging data (see Fig. 1b). On a group-representative brain network constructed by averaging all 190 subject-specific networks, we calculate these controllability metrics and show their distribution across the brain (see Fig. 1c.i). Consistent with prior work in both this and other data sets Gu et al. (2015); Tang et al. (2017), we find that regions with high average controllability tend to be located in network hubs associated with the default-mode system, a set of regions that tend to be active during intrinsic processing. Interestingly and again consistent with prior work in this and other data sets Gu et al. (2015); Tang et al. (2017); Jeganathan et al. (2018); Wu-Yan et al. (2018), regions with high modal controllability tend to display low average controllability (see Fig. 1c.ii), and are predominantly found in regions of the brain that become active during tasks that demand high levels of executive function or cognitive control.

III Controllability of spatial scales: Network dynamics in harmonic waves

In contrast to the time-varying or transient dynamics that can be described using models mentioned in the previous section, recent observations of brain dynamics also includes sustained patterns of collective activity across relatively large areas of the brain Muller et al. (2018); Roberts et al. (2019). Such wave patterns instead call for a different description, one better suited by a framework developed for networked oscillators and synchronized activity.

Specifically, the previous linear model (Eq. 1) is not relevant to wave-like dynamics, even as it is supported by other empirical observations Honey et al. (2009); Fernández Galán (2008); Yan et al. (2017a). Therefore, to complement our previous efforts to study the control of brain state transitions via a simple linear model of network dynamics, here we turn to the question of how to control the synchrony of harmonic waves Atasoy et al. (2016) which naturally arise in the same underlying white matter architecture. Notably, extensive prior work in the dynamical systems literature has considered ways in which to measure the synchronizability of a networked system. To be clear, here we define the synchronizability as a measure of the ability of a network to persist in a single synchronous state $\mathbf{s}(t)$, i.e. $\mathbf{x}_{1}(t)=...=\mathbf{x}_{n}(t+1)=\mathbf{s}(t)$. The master stability function (MSF) allows analysis of the stability of this synchronous state without detailed specification of the properties of the dynamical units (see Fig. 1d; Pecora and Carroll (1998)). Within this framework, linear stability depends upon the positive eigenvalues $\{\lambda_{i}\},i=1,...,N-1$ of the Laplacian matrix $L$ defined by $L_{ij}=\delta_{ij}\sum_{k}A_{ik}-A_{ij}$. The synchronous state is stable when the master stability function is negative for all positive eigenvalues of the graph Laplacian (see Fig. 1d).

For the purposes of control, it is useful to note that the structure of the graph Laplacian provides information regarding the most likely modes of standing waves: the vibrational states of a dynamical system in which the frequency of vibration is the same for all elements. These modes are given by the eigenvectors of the graph Laplacian. For example, the eigenvector $\phi_{1}$ of the smallest positive Laplacian eigenvalue $\lambda_{1}$ is an odd mode associated with large-scale waves Atasoy et al. (2016); Pecora and Carroll (1998) (Fig. 2a). In contrast, the eigenvector $\phi_{N-1}$ of the largest Laplacian eigenvalue $\lambda_{N-1}$ is an even mode associated with small-scale waves (Fig. 2b). Recent developments in push-pull control capitalize on these relations to control the transition into and out of synchrony by targeting network nodes that have a large weight for a given extremal eigenvector of the Laplacian He et al. (2014).

To investigate the control of synchrony in the harmonic waves of brain networks, we calculate the eigenvectors of the Laplacian for each subject, and then we average the eigenvectors across all of the 190 subjects. We observe that the strength of control for harmonic waves of a network can be understood in terms of the network’s predisposition to average and modal control. First, we observe that the regional strength of the large-scale waves $|\phi_{1}^{j}|$ is positively correlated with regional modal controllability (Spearman correlation coefficient $\rho=0.27$, $df=233$, $p<1\times 10^{-4}$; Fig. 2c). This relation suggests that regions that enable synchronous behavior over long distances are also predicted to be effective in moving the brain to energetically distant states (see a list of the highest and lowest such regions in Table 1 of the Appendix). Second, we observe that the regional strength of the small-scale waves $|\phi_{N-1}^{j}|$ is positively correlated with regional average controllability (Spearman correlation coefficient $\rho=0.95$, $df=233$, $p<1\times 10^{-5}$; Fig. 2d). This relation suggests that regions that enable synchronous behavior over short distances are also predicted to be effective in moving the brain to energetically nearby states (see a list of the highest and lowest such regions in Table 2 of the Appendix). To confirm that our results do not depend on the stage at which we average over subjects, we also calculate these correlations for the Laplacian spectra obtained from a group-averaged network (procedure described in Section II). We find that the regional strength of the large-scale waves in the group network is still positively correlated with regional modal controllability (Spearman correlation coefficient $\rho=0.28$, $df=233$, $p<1\times 10^{-4}$), while the regional strength of the small-scale waves in the group network is still positively correlated with regional average controllability (Spearman correlation coefficient $\rho=0.86$, $df=233$, $p<1\times 10^{-5}$). These results demonstrate that our findings do not depend unduly on the stage at which we average over subjects.

IV Controllability of temporal scales: Network dynamics in fast or slow, alternating or monotone modes

To complement the findings reported in the previous section explicating the spatial scales of control, we next turn to a consideration of the time scales of control. We note that in a discrete-time system, the eigenvalues of the network adjacency matrix $\mathbf{A}$ give the rate of decay of the various supported modes of activity. These modes of a linear time-invariant system are the elementary components of the dynamic response of the system to stimuli. Each mode gives rise to an elementary response (or trajectory), and the combination of the modes generates more complex responses and trajectories. Distinct processes and time scales of activity are observed in the brain, where rhythms of certain frequencies have been associated with different cognitive processes and behavioral responses Palmigiano et al. (2017); Kirst et al. (2016). For instance, higher frequency $\gamma$ rhythms in local synchronization emerge during visual responses, while lower frequency $\beta$ rhythms reflecting coherence over parietal and temporal cortices have been observed during activities that required more multi-modal sensory integration von Stein et al. (1999); Sherman et al. (2016). In our discrete time linear model, the eigenvalues with large magnitude decay slowly while the eigenvalues with small magnitude decay quickly. The relationship between eigenvalues and the rate of decay of associated modes of activity suggests that the control of such modes can be analyzed separately for different time scales. While our prior expression for modal controllability (Eq. 4) summed over all $N$ modes of the adjacency matrix $\mathbf{A}$, we can instead partition the sum into different groups based on the eigenvalues $\xi_{j}$ of the modes. For instance, a region $i$ can be evaluated for its controllability of fast modes with $\phi_{i}^{\textrm{fast}}=\sum_{j}^{\textrm{small }|\xi_{j}|}v_{ij}^{2}$ and for its controllability of slow modes with $\phi_{i}^{\textrm{slow}}=\sum_{j}^{\textrm{large }|\xi_{j}|}v_{ij}^{2}$. Note that in this case, all the eigenvalues have magnitude smaller than 1. Further, in the above expression we have also removed the scaling factor in Eq. 4 of $\left(1-\xi_{j}^{2}(A)\right)$ as this factor is unnecessary for the calculation of control time-scales; in supporting calculations, we also verified that the inclusion or exclusion of this factor does not alter our subsequent findings.

To investigate regional propensities for the control of fast versus slow modes, we considered the group-representative network constructed by averaging all 190 subject-specific brain networks. We examined the spectrum of this representative network (see Fig. 3), and note that there are three well-separated groups of eigenvalues. In the center, there is a densely populated cluster of magnitude between 0 and 0.2 (fast), which have at least double the modes in each bin of width 0.1 as compared to all the bins outside this cluster. On the two ends, there are sparse modes with magnitude above 0.6 (slow) that have a clear gap from all the other modes. In between these above mentioned groups, there is a lightly populated cluster of eigenvalue magnitude between 0.2 and 0.6, which are distinct from the other groups due to the reasons mentioned. Given this visually trimodal distribution, it is hence natural to partition the eigenvalues into three groups that differ in their relevant time scales: $|\xi_{j}|<0.2$ (fast), $0.2<|\xi_{j}|<0.6$ (medium), and $|\xi_{j}|>0.6$ (slow).

Now, it is also important to note that eigenvalues carry additional information: the sign of the eigenvalue determines the nature of the system’s response, where a mode with a positive eigenvalue decays monotonically with each time step, while a mode with a negative eigenvalue has a decaying response that alternates between positive and negative values at each time step. To allow sensitivity to these monotone and alternating modes, we can restrict the range in the sum of $\phi_{i}^{\textrm{fast}}$ and $\phi_{i}^{\textrm{slow}}$ to the positive or to the negative eigenvalues. By doing so, we now have four groups of interest: $0<\xi_{j}<0.2$ (fast monotone), $\xi_{j}>0.6$ (slow monotone), $-0.2<\xi_{j}<0$ (fast alternating), and $\xi_{j}<-0.6$ (slow alternating). By separately identifying these diverse modes, we can determine whether and how brain regions’ differing controllability of particular time scales might relate to their differing controllability of various length scale modes.

Beginning with the control of only monotone dynamics, we find that support for slow monotone dynamics (eigenvalues $\xi_{j}>0.6$) differs across the brain with greatest support located in the prefrontal cortex and temporoparietal junction (Fig. 4a; see a list of the highest and lowest such regions in Table 3 of the Appendix). Moreover, we find that the degree to which a regional controller supports slow monotone dynamics is strongly and positively correlated with regional average controllability (Spearman correlation coefficient $\rho=0.99$, $df=233$, $p<1\times 10^{-4}$; Fig. 4c), providing support for energetically easy state transitions. Similarly, we also find that support for fast monotone dynamics (eigenvalues $0<\xi_{j}<0.2$) differs across the brain with greatest support located in subcortical areas (Fig. 4b; see a list of the highest and lowest such regions in Table 4 of the Appendix). In addition, we find that the degree to which a regional controller supports fast monotone dynamics is strongly and positively correlated with regional modal controllability ($\rho=0.59$, $df=233$, $p<1\times 10^{-4}$; Fig. 4d), providing support for energetically more difficult state transitions.

When studying the control of alternating dynamics, we find that support for slow alternating dynamics (eigenvalues $\xi_{j}<-0.6$) differs across the brain with greatest support located in the temporoparietal junction, inferior frontal gyrus, and precuneus (Fig. 5a; see a list of the highest and lowest such regions in Table 5 of the Appendix). Moreover, we find that the degree to which a regional controller supports slow alternating dynamics is strongly and positively correlated with regional average controllability ($\rho=0.83$, $df=233$, $p<1\times 10^{-4}$; Fig. 5c), providing support for energetically easy state transitions. Lastly, we also find that support for fast alternating dynamics (eigenvalues $-0.2<\xi_{j}<0$) differs across the brain with greatest support located in several midline structures (Fig. 5b; see a list of the highest and lowest such regions in Table 6 of the Appendix). In this case, we find that the degree to which a regional controller supports fast alternating dynamics is strongly and positively correlated with regional modal controllability ($\rho=0.24$, $df=233$, $p<2\times 10^{-4}$; Fig. 5d), providing support for energetically more difficult state transitions.

V Dependence Upon Graph Architecture

The results described thus far provide striking relations between the temporal and spatial scales of network dynamics that a brain region can control. An interesting and important question to ask is whether and to what degree these relations hold true in graph models with vastly different topologies. In choosing a particularly dissimilar comparator, we consider the fact that human brain networks are known to display high clustering and short path length Bassett and Bullmore (2016, 2006), community structure Baum et al. (2017), and core-periphery structure Betzel et al. (2018). A notable graph model that stands in stark contrast to this architecture is the scale free model Albert and Barabasi (2002), also sometimes referred to as the preferential attachment model Price (1976) or the Barabasi-Albert model Barabasi and Albert (1999), which displays low clustering and longer path length, less community structure, and less core-periphery structure. For ease of comparison, we provide a table contrasting these primary network metrics in both networks, in the Appendix (see Table VII).

We therefore study a synthetically constructed Barabasi-Albert network with the same number of nodes ($N=234$) and the same average strength as the group-representative brain network in our empirical data set. As before, we calculate controllability metrics and the Laplacian eigenvectors to investigate the spatial harmonics. We find that the relationships between controllers of different length scales with average and modal controllers remain similar across both brain and Barabasi-Albert graphs. Specifically, the nodes in the Barabasi-Albert graph most relevant for the large-scale mode $|\phi_{1}^{j}|$ tend to be the nodes with high modal controllability: $\rho=0.64$, $df=233$, $p<1\times 10^{-4}$. Meanwhile, the nodes in the Barabasi-Albert graph most relevant for the small-scale mode $|\phi_{N-1}^{j}|$ tend to be the nodes with high average controllability: $\rho=0.77$, $df=233$, $p<1\times 10^{-5}$. Such positive correlations were also found in the brain network (Section III).

To study the control of modes with different time-scales, we employ the same partition that we used in the brain network, i.e. groups of $|\xi_{j}|<0.2$, $0.2<|\xi_{j}|<0.6$, and $|\xi_{j}|>0.6$ for both positive and negative eigenvalues, respectively (see Fig. 6). Interestingly, we find that the distribution of eigenvalues differs significantly in the Barabasi-Albert network in comparison to the group-representative brain network, with the striking absence of slow alternating modes. When we consider only the fast alternating modes, we observe that the relationship between the controllers of these fast modes and the modal controllers is inverted in the Barabasi-Albert network in comparison to the group-representative brain network (see Fig. 7). In the brain, the degree to which a regional controller supports fast alternating dynamics tracks the degree to which that same regional controller supports energetically more difficult transitions. In the Barabasi-Albert network, the degree to which a regional controller supports fast alternating dynamics is anti-correlated with the degree to which that same regional controller supports energetically more difficult transitions (Spearman correlation coefficient $\rho=-0.91$, $df=233$, $p<1\times 10^{-4}$).

VI Discussion

The human brain displays a variety of different types of dynamics. Here we use two different models to compare the control profiles of various brain regions: diffusion dynamics and network control, respectively. The two dynamical models we use capture different properties of the same brain network, hence it is useful to compare and contrast these different dynamics produced by the same system. The model for synchronizability uses the Master Stability Function (MSF) formalism based on a model of connected oscillators, whereas our controllability analysis relies on a linear discrete time model. Note that while the MSF approach originally studied continuous time systems, there have been extensions to discrete time systems Porfiri (2011). In the latter case, a modified graph Laplacian is used that has the same eigenvalues as a regular Laplacian operator (and the eigenvectors are similar up to a constant shift) Sun et al. (2009). Hence, our two models are related by the structure of the network, which appears in the interconnections among the oscillators in the first model, and as a weighted adjacency matrix in the second model. Through these analyses we find that brain regions with high average controllability are also predicted to effectively control dynamics on short length scales, while regions of high modal controllability are predicted to do the converse. In addition, regions of high average controllability are also predicted to effectively control dynamics on slow scales, while regions of high modal controllability are predicted to do the converse. Collectively, our results expand our understanding of the relations between brain network dynamics and control, and the relative roles that different brain regions can play in controlling the diverse length and time scales of neural dynamics.

Pertinent Theoretical Considerations. There are two important theoretical considerations relevant to our work: the relation between linear and nonlinear network control, and the nature of the Master Stability Function. The first consideration is crucial, particularly given recent results regarding the difficulty of applying linear control models to nonlinear systems Jiang and Lai (2019). While linear models are certainly valid only within the domain of linearization, we note that the modelling of brain activity in large-scale regional networks shows that the linear approximation provides fair explanatory power of resting state fMRI BOLD data Honey et al. (2009). Further, studies of this controllability framework using non-linear oscillators connected with coupling constants estimated from large-scale white matter structural connections shows a good overlap with the linear approximation Muldoon et al. (2016). While the model we employ is a discrete-time system, this controllability Gramian is statistically similar to that obtained in a continuous-time system Gu et al. (2015), through the comparison of simulations run using MATLAB’s lyap function. Regarding the second, we note that we have plotted a typical example of a Master Stability Function (MSF) for a network of oscillators schematically in Fig. 1d; however, specific details will depend on the dynamics of individual nodes and the connectivity between them. The shape of the MSF for various families of dynamical systems is typically convex for generic oscillator systems, including chaotic oscillators that have stable limit cycles Huang et al. (2009).

Separately, it is worth noting that there are several important and interesting questions at the intersection of control theory and network neuroscience. First, early work focused on the question of whether a brain network is controllable from a single region, providing initial evidence in the affirmative Gu et al. (2015). A subsequent numerical study cast doubt on that evidence Tu et al. (2018), while two even later studies provided additional theoretical evidence in the affirmative RE: Warnings and caveats in brain controllability (2019); Menara et al. (2017). Notably, the structural controllability of networks depends upon their topology, being quite different for diverse graph models including the Erdős-Rényi, ring lattice, Watts-Strogatz, modular, random geometric, and Barabasi-Albert model RE: Warnings and caveats in brain controllability (2019); Wu-Yan et al. (2018), the latter of which we also study here. Second, several studies address the question of whether a set of control points (network nodes) can be identified that can effectively and reliably induce a particular change in system dynamics Yan et al. (2017b); Zañudo et al. (2017). Third, another body of work focuses on the question of how to guide the brain from one state of activity to another state of activity given the underlying structural connections along which activity can propagate Gu et al. (2017); Betzel et al. (2016). Work in this vein has offered validated predictions of the effects of brain stimulation Stiso et al. (2019), and explanations of individual differences in cognitive control Cui et al. (2019). This list of questions is far from exhaustive, and we expand upon it in this work by probing the time scales and spatial scales relevant for control in brain networks.

Future Directions. Our work provides an initial foray into the combined study of brain network dynamics and control, and how the two depend upon the length scales and time scales of brain states and state transitions. From a mathematical perspective, it would be interesting to study these relations in other graph models with different topologies, and also to work to find analytical forms for the relations we observe in numerical experiments. From a neurophysics perspective, future efforts could extend our work into other age groups, or into clinical groups of patients with psychiatric disease or neurological disorders where such relations might be altered. Efforts could also determine whether our findings in humans are recapitulated in other animals, and across different spatial resolutions of imaging data. Another avenue for further investigation could be the use of directed brain networks, such as those that are available from tract-tracing studies of macaques, which would display complex eigenvalues thereby providing additional oscillatory time-scales for more extensive study. Lastly, this trimodal nature of the data (see Fig. 3) appears to arise when studying group-representative data, which led to our choice of time-scale bins. Future work could extend our study to the single-subject level, where a distinct and subject-specific binning procedure might prove helpful.

VII Acknowledgements

D.S.B. and E.T. acknowledge support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the ISI Foundation, the Paul Allen Foundation, the Army Research Laboratory (W911NF-10-2-0022), the Army Research Office (Bassett-W911NF-14-1-0679, Grafton-W911NF-16-1-0474, DCIST- W911NF-17-2-0181), the Office of Naval Research, the National Institute of Mental Health (2-R01-DC-009209-11, R01 – MH112847, R01-MH107235, R21-M MH-106799), the National Institute of Child Health and Human Development (1R01HD086888-01), National Institute of Neurological Disorders and Stroke (R01 NS099348), and the National Science Foundation (BCS-1441502, BCS-1430087, NSF PHY-1554488 and BCS-1631550). The content is solely the responsibility of the authors and does not necessarily represent the official views of any of the funding agencies.

VIII Appendix

All data were acquired from the Philadelphia Neurodevelopmental Cohort (PNC), a large community-based study of brain development. This resource is publicly available through the Database of Genotypes and Phenotypes. Each subject provided their informed consent according to the Institutional Review Board of the University of Pennsylvania who approved all study protocols. These subjects had no gross radiological abnormalities that distorted brain anatomy, no history of inpatient psychiatric hospitalization, no use of psychotropic medications at the time of scanning, and no medical disorders that could impact brain function. Each included subject also passed both manual and automated quality-assessment protocols for DTI Roalf et al. (2016) and T1-weighted structural imaging Vandekar et al. (2015), and had low in-scanner head motion (less than 2mm mean relative displacement between b=0 volumes). Here we study a subset of the full group previously reported in Tang et al. (2017); specifically, we consider the 190 young adults aged 18 to 22 years.

Structural connectivity was estimated using 64-direction DTI data. The diffusion tensor was estimated and deterministic whole-brain fiber tracking was implemented in DSI Studio using a modified FACT algorithm, with exactly 1,000,000 streamlines initiated per subject after removing all streamlines with length less than 10mm Gu et al. (2015). A 234-region parcellation Daducci et al. (2012) was constructed from the T1 image using FreeSurfer. Parcels were dilated by 4mm to extend regions into white matter, and registered to the first non-weighted (b=0) volume using an affine transform. Edge weights $A_{ij}$ in the adjacency matrix were defined by the number of streamlines connecting each pair of nodes end-to-end. All analyses were replicated using an alternative edge weight definition, where weights are equal to the number of streamlines connecting each node pair divided by the total volume of the node pair, as well as using probabilistic fiber tracking methods. A schematic for structural connectome construction is depicted in Fig. 1b.