Topological Principles of Control in Dynamical Network Systems

TL;DR

This paper develops analytical methods linking network topology to control energy in dynamical systems, with applications to brain connectomes across species, revealing insights into energy efficiency and robustness.

Contribution

It introduces closed-form expressions relating network structure to control energy, enabling topological design and targeted interventions in complex networks.

Findings

Connectomes of complex species are wired to minimize control energy.

Removing single edges can alter the brain's control profile.

Humans exhibit efficient and robust network control properties.

Abstract

Networked systems display complex patterns of interactions between a large number of components. In physical networks, these interactions often occur along structural connections that link components in a hard-wired connection topology, supporting a variety of system-wide dynamical behaviors such as synchronization. While descriptions of these behaviors are important, they are only a first step towards understanding the relationship between network topology and system behavior, and harnessing that relationship to optimally control the system's function. Here, we use linear network control theory to analytically relate the topology of a subset of structural connections (those linking driver nodes to non-driver nodes) to the minimum energy required to control networked systems. As opposed to the numerical computations of control energy, our accurate closed-form expressions yield general structural features in networks that require significantly more or less energy to control, providing topological principles for the design and modification of network behavior. To illustrate the utility of the mathematics, we apply this approach to high-resolution connectomes recently reconstructed from drosophila, mouse, and human brains. We use these principles to show that connectomes of increasingly complex species are wired to reduce control energy. We then use the analytical expressions we derive to perform targeted manipulation of the brain's control profile by removing single edges in the network, a manipulation that is accessible to current clinical techniques in patients with neurological disorders. Cross-species comparisons suggest an advantage of the human brain in supporting diverse network dynamics with small energetic costs, while remaining unexpectedly robust to perturbations. Our results ground the expectation of a system's dynamical behavior in its network architecture.

Equations

$$\left[\begin{array}{c}{\dot{\mathbf{x}}}_{\text{d}}\\ {\dot{\mathbf{x}}}_{\text{nd}}\end{array}\right]=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_{\text{d}}\\ {\mathbf{x}}_{\text{nd}}\end{array}\right]+\left[\begin{array}{c}{I}_N\\ 0\end{array}\right]\mathbf{u},$$

$$\text{E}(\mathbf{u})=\sum _{i=1}^N\underset{{\text{E}}_i}{\underbrace{{\int }_0^T{u}_i(t{)}^{2}dt}},$$

$$\text{E}(\mathbf{u})=12({\mathbf{x}}_{\text{nd}}^{\ast }-\frac{1}{2}{A}_{\text{21}}{\mathbf{x}}_{\text{d}}^{\ast }{)}^T(A_{\text{21}}{A}_{\text{21}}^T{)}^{-1}({\mathbf{x}}_{\text{nd}}^{\ast }-\frac{1}{2}{A}_{\text{21}}{\mathbf{x}}_{\text{d}}^{\ast })+{\mathbf{x}}_{\text{d}}^{\ast T}{\mathbf{x}}_{\text{d}}^{\ast },$$

$$\text{E}(\mathbf{u})=12\frac{{\mathbf{v}}_1^T\text{adj}(Q){\mathbf{v}}_1}{\det (Q)}+{\mathbf{v}}_2^T{\mathbf{v}}_2,$$

$$\text{E}(\mathbf{u})=12\left(\frac{{\sum }_{i=1}^M{w}_i{c}_i^2}{{\sum }_{i=1}^M{w}_i}\right)\left(\sum _{k=1}^M\frac{1}{\Vert {\mathbf{a}}_k{\Vert }^2\sin ({\theta }_k{)}^2}\right)+{\mathbf{v}}_2^T{\mathbf{v}}_2,$$

$$\frac{\partial }{\partial A_{\text{21}}}\det (Q)=2\det (Q)(Q^{-1}{A}_{\text{21}}),$$

$$\dot{\mathbf{x}}=A\mathbf{x}.$$

$$\left[\begin{array}{c}{\dot{\mathbf{x}}}_{\text{d}}\\ {\dot{\mathbf{x}}}_{\text{nd}}\end{array}\right]=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_{\text{d}}\\ {\mathbf{x}}_{\text{nd}}\end{array}\right]+\left[\begin{array}{c}{I}_N\\ 0\end{array}\right]\mathbf{u},$$

$$\text{E}(\mathbf{u})=\sum _{i=1}^N\underset{{\text{E}}_i}{\underbrace{{\int }_0^T{u}_i(t{)}^{2}dt}},$$

$$\begin{array}{rl}{\mathbf{x}}_{\text{d}}(t)& ={\int }_0^t\mathbf{u}(\tau )d\tau ,\text{ and }\\ {\mathbf{x}}_{\text{nd}}(t)& =A_{21}{\int }_0^t{\mathbf{x}}_{\text{d}}(\tau )d\tau =A_{21}{\int }_0^t{\int }_0^{\tau }\mathbf{u}({\tau }_1)d{\tau }_{1}d\tau .\end{array}$$

$${\int }_0^t{\mathbf{x}}_{\text{d}}(\tau )d\tau =\mathbf{C},$$

$$\mathcal{C}=\left[\begin{array}{cc}B& AB\end{array}\right]=\left[\begin{array}{cc}{I}_N& 0\\ 0& A_{21}\end{array}\right],$$

$$E_i$$

$$=\underset{x_{\text{di}}}{\min }{\int }_0^T{\dot{x}}_{\text{di}}(t{)}^{2}dt,$$

$${\int }_0^T{x}_{\text{di}}(t)dt$$

$$\underset{x_{\text{di}}}{\min }{\int }_a^{b}F(t,x_{\text{di}},\dot{x}{}_{\text{di}})dt,$$

$${\int }_a^{b}G(t,x_{\text{di}},\dot{x}{}_{\text{di}})dt=0,$$

$$\frac{d}{dt}\frac{\partial }{\partial \dot{x}}(F+\lambda G)$$

$${\ddot{x}}_{\text{di}}^{\ast }(t)$$

$$x_{\text{di}}^{\ast }(t)$$

$${\int }_0^1{\dot{\eta }}^2(F+\lambda G{)}_{\dot{x}\dot{x}}dt=2{\int }_0^1{\dot{\eta }}^{2}dt\ge 0,$$

$$E_i$$

$$={\int }_0^{1}4a_i{t}^2+4a_i{b}_{i}t+b_i^{2}dt$$

$$=\frac{4}{3}{a}_i+2a_i{b}_i+b_i^2.$$

$$x_{\text{di}}(T=1)=a_i+b_i$$

$${\int }_0^{T=1}{x}_{\text{di}}(t)dt=\frac{1}{3}{a}_i+\frac{1}{2}{b}_i$$

$$a_i$$

$$b_i$$

$$E_i=12C_i^2-12C_i{x}_{\text{di}}+4x_{\text{di}}^2$$

$$f(C)=E(\mathbf{u})$$

$$=\sum _{i=1}^{N}12C_i^2-12C_i{x}_{\text{di}}+4x_{\text{di}}^2$$

$$=12{\mathbf{C}}^T\mathbf{C}-12{\mathbf{C}}^T{\mathbf{x}}_{\text{d}}+4{\mathbf{x}}_{\text{d}}^T{\mathbf{x}}_{\text{d}},$$

$$g(\mathbf{C})$$

$$={\mathbf{x}}_{\text{nd}}^{\ast }-A_{\text{21}}\mathbf{C}$$

Full text

I Introduction

Network systems are composed of interconnected units that interact with each other on diverse temporal and spatial scales newman2010networks . The exact patterns of interconnections between these units can take on many different forms, and those forms can dictate how the system functions newman2003structure . Indeed, specific features of network topology – such as small-worldness Watts1998 and modularity simon1962architecture – can give rise to properties like control efficiency Latora2001 ; Vragovic2005 and robustness against component failure Joyce2013 that are highly desirable for both natural and man-made systems. While the relationship between interconnection architecture and dynamics is observed ubiquitously across technological, social, and physical systems, it provides particularly important insights into the functional capabilities of biological systems such as the brain. Here, the topology of interconnection pattterns between neural units is thought to support optimal information processing bassett2006small ; bassett2016small ; sporns2016modular , both at the level of individual cells bettencourt2007functional and at the level of meso-scale brain areas bassett2006adaptive ; bullmore2009complex .

Despite the observation that network topology and system function are related to one another, there still remains very little understanding of the exact mechanisms driving this relationship pasqualetti2014controllability . Gaining such an understanding would have far reaching implications for the analysis, modification, and control of interconnected complex systems newman2006structure . This relationship could be exploited for personalized therapeutics barabasi2011network that would significantly enhance clinical outcomes for patients by mapping the network topology in musculoskeletal murphy2016structure , gene regulatory conaco2012functionalization ; henzler2013staged ; norton2016detecting , and central nervous bassett2016network systems. For example, treatments for drug-resistant epilepsy currently range from surgical resection or laser ablation of epileptic tissue to stimulation designed to stem seizure progression. Yet, these interventions are often complicated by the presence of epileptic tissue in areas of the brain that are essential for motor and language function Nune2015 . An understanding of the specific role of regions and connections in brain networks could inform more targeted and less invasive therapies to make the epileptic state energetically unfavorable to reach, or impossible to maintain ching2012distributed ; khambhati2015dynamic ; khambhati2016virtual .

Existing paradigms for exploring the mechanisms by which a complex network topology drives observable dynamics come from diverse intellectual fields and are built on varying assumptions. One paradigm stems from the field of nonlinear dynamics, and deals with attractors (states a system naturally tends to) and basins of attraction (regions of initial states that naturally fall into an attractor) Sprott2015 . This is an effective tool for understanding meso-scale networks that display nonlinear dynamics, and for defining perturbations of state trajectories to force a system to transition from one basin to another Cornelius2013 . A principal challenge of this approach is that analytical solutions explaining mechanisms of network control remain sparse. An alternative paradigm involves calculating statistical correlations between observed function and structure using graph theoretical metrics such as network communicability estrada2008communicability , modularity newman2006modularity , or search information trusina2005communication ; rosvall2005searchability ; sneppen2005hide . This approach has been used in the context of brain networks to predict – in a statistical sense – patterns of activity from patterns of wiring goni2014resting , the impact of focal lesions on distant brain areas following stroke crofts2011network , and the role of modules in facilitating adaptive functions such as learning Bassett2011 ; mantzaris2013dynamic ; bassett2013task ; bassett2014cross ; bassett2015learning . The principal challenge of this approach is that graph statistics themselves do not constitute mechanisms in a philosophical craver2005beyond or mathematical sense bassett2016mattar .

A promising paradigm that meets these challenges is linear network control theory kalman1963mathematical ; kalman1963controllability ; lin1974structural , which assumes that the dynamics of a system are state-dependent, and that they are linear with respect to the state. In this framework, the state of a system at a given time is a function of the previous state, the structural network linking units of that system, and any input to the system provided in the form of control energy (see next section for explicit mathematical definitions). From this paradigm arises the possibility of identifying driver nodes in the system liu2011controllability ; campbell2015topological ; ruths2014control : units that have the potential to influence the system along diverse control strategies. Also, within this framework we can compute optimal inputs that move the system from one state to another with minimal cost. This formulation has been particularly useful in understanding a variety of networked systems, including the human brain where control points facilitate diverse cognitive strategies Gu2015 ; gu2016optimal , facilitate efficient activation and deactivation patterns during endogeneous activity Betzel2016 , and inform optimal targets for brain stimulation to alter neural activity muldoon2016stimulation .

While the identification of control points and optimal trajectories is computationally tractable, basic intuitions about the network properties that enhance control – either locally around a given node or globally throughout the whole system – have remained elusive. In this study, we sought to identify the key topological features that determine network controllability, and use these principles to understand and modify the complex connectivity of meso-scale brain networks in a dynamically meaningful way. To reach these goals, we formulate a linear control problem that examines a subset of a network’s edges: the connections that link driver nodes to non-driver nodes. Using this bipartite subgraph of the entire network, we explore the mathematical theory of fundamental driver$\to$non-driver dynamics to the extent of explicitly understanding the role of every node and edge in a networked system, and modifying those objects to alter control energy in an exactly predictable manner. We show that the intuitions and solutions for the control of the bipartite subgraph provide excellent estimates of the control of the full network. Our results include analytical derivations of expressions relating a network’s minimum control energy to its connection topology, along with an intuitive geometric representation to visualize this relationship. While our mathematical contributions are applicable to any complex network system whose dynamics can be approximated by a linear model, we illustrate the utility of the formulation in the context of brain networks estimated from the mouse brain (made publically available by the Allen Brain Institute) oh2014mesoscale ; rubinov2015wiring , the drosophila brain shih2015connectomics , and the human brain (Fig. 1d–f). Specifically, we use the analytical expressions to (i) understand key patterns and principles of connectivity that determine a network’s control profile, (ii) describe the implications of the connectivity of brain networks on their control profiles, and (ii) explicitly modify the control properties of the brain by performing energetically favorable edge deletions, thereby informing potential clinical interventions. Together, these results offer fundamental insights into key patterns of connections between brain regions that directly impact their minimum control energy, providing a link between the structure and function of neural systems.

II Network Topology and Controllability

We are particularly interested in understanding how a network’s specific topological features (edge connections, edge weights, weight distributions) affect the energy required to control the network. We are also interested in how that same topology facilitates or inhibits certain states from being reached. To this end, we develop our analysis tools based on a simplified network model, which effectively reveals hidden relations between network topology and control energy. Finally, we validate the predictive power of our results on network models representing the mouse, drosophila, and human brains.

We consider a network represented by the directed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{1,\dots,n\}$ and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ are the sets of network vertices and edges, respectively. Let $a_{ij}\in$ be the weight associated with the edge $(i,j)\in\mathcal{E}$, and let $A=[a_{ij}]$ be the weighted adjacency matrix of $\mathcal{G}$. We associate a real value (state) with each node, collect the nodes’ states into a vector (network state), and define the map $\bm{x}:_{\geq 0}\rightarrow^{n}$ to describe the evolution (dynamics) of the network state over time (Fig. 1a–c). We assume that a subset of $N$ nodes, called drivers, is independently manipulated by external controls and, without loss of generality, we reorder the network nodes such that the $N$ drivers come first. Thus, the network dynamics with controlled drivers read as

[TABLE]

where ${\bm{x}}_{\textup{d}}$ and ${\bm{x}}_{\textup{nd}}$ are the state vectors of the driver and non-driver nodes, $A_{11}\in^{N\times N}$, $M=n-N$, $A_{12}\in^{N\times M}$, $A_{21}\in^{M\times N}$, $A_{22}\in^{M\times M}$, $I_{N}$ is the $N$-dimensional identity matrix, and $\bm{u}:_{\geq 0}\rightarrow^{N}$ is the control input.

We refer to the network as controllable at time $T\in_{\geq 0}$ if, for any pair of states ${\bm{x}}_{\textup{d}}^{*}$ and ${\bm{x}}_{\textup{nd}}^{*}$, there exists a control input $u$ for the dynamics Eq. (1) such that ${\bm{x}}_{\textup{d}}(T)={\bm{x}}_{\textup{d}}^{*}$ and ${\bm{x}}_{\textup{nd}}(T)={\bm{x}}_{\textup{nd}}^{*}$. For a detailed discussion and rigorous conditions for the controllability of a system with linear dynamics, see TK:80 . Finally, we define the energy of $\bm{u}$ as

[TABLE]

where $u_{i}$ is the $i$-th component of $\bm{u}$. The energy of $u_{i}$ can be thought of as a quadratic cost that penalizes large control inputs to drive the system. In what follows, we characterize how a network’s topology determines the control energy needed for a given control task.

We approach this problem by approximating the interactions between brain regions as linear, time invariant dynamics, where a stronger structural connection between two regions represents a stronger dynamic interaction (for empirical motivation for this approximation, see Gu2015 ; galan2008network ; honey2009predicting ). From these dynamics, we identify physically meaningful interpretations of the underlying mathematical features, which determine the controllability of the simplified network only containing connections from drivers $\to$ non-drivers (Fig. 1h). To justify this approximation, we show that the simplified network well approximates the control dynamics along the full network (Fig. 1g) for a wide range of model parameters. From this simplification, we derive a closed-form expression of the minimum control energy that shows that the similarity in these driver $\to$ non-driver connections scales the energy required to control the network. Given this scaling, we show that dissimilarly connected regions are easiest to control, while similarly connected regions are most difficult to control. Finally, we use these principles to gain insight into the relationship between the structural connectivity and function of brain networks, and make a few strategic topological modifications to affect a profound reduction in control energy. We conclude by discussing the utility of these insights in informing interventions to modulate brain dynamics.

III Results

III.1 Control energy is well predicted by direct connections between driver and non-driver nodes

We seek an accurate, tractable relationship between the topology of a network and the energy required to drive the network from one state to another. To find this relationship, we look to classic results in the mathematical theory of systems and control TK:80 , where the spectral properties of the reachability Gramian $W_{R}(0,T)=\int_{0}^{T}e^{At}BB^{T}e^{A^{T}t}$ quantify the minimum amount of energy (Section V.3) to control the network Eq. (1). Explicit formulas and bounds for the eigenvalues of the Gramian are intricate and hard to derive, and approximate formulations are often preferred.

To make progress on this problem, we begin by calling the network involving edges between all nodes a non-simplified network (Fig. 2a), and the network involving only the edges from the driver to the non-driver nodes a simplified, first-order network (Fig. 2b). We then derive an accurate approximation of the minimum control energy (Lemma V.2 - V.4) by assuming that ${\bm{x}}_{\textup{d}}(0)=0$, ${\bm{x}}_{\textup{nd}}(0)=0$ (Assumption 1), and $A_{11}=0$, $A_{12}=0$, and $A_{22}=0$ (Assumption 2) in Eq. (1), which reads as

[TABLE]

where ${\bm{x}}_{\textup{d}}^{*}$ and ${\bm{x}}_{\textup{nd}}^{*}$ are the desired final states of the driver and non-driver nodes, respectively. We make Assumption 1 because we are interested in the change in brain state through control, and consider initial conditions ${\bm{x}}_{\textup{d}}(0)=0$, ${\bm{x}}_{\textup{nd}}(0)=0$ to be a neutral baseline.

Because the expression Eq. (2) involves only the edges in the simplified network from the driver to the non-driver nodes, that is, the matrix $A_{21}$, we say that Eq. (2) is a first-order approximation to the minimum control energy of the non-simplified network Eq. (1). One topological feature that impacts the accuracy of the first-order energy approximation is the fraction of nodes that are selected as non-drivers: the drivers of a simplified first-order network (Fig. 2b) can only control the non-drivers through the direct driver $\to$ non-driver connections, while a non-simplified network (Fig. 2a) can additionally use indirect non-driver $\to$ non-driver connections. In a network with few non-drivers, there are relatively few indirect connections, so most of the control dynamics rely upon the direct connections; in a network with many non-drivers, there are more indirect connections for the drivers to utilize, making the first-order control energy Eq. (2) a worse approximation. A second topological feature that affects the accuracy of the first-order approximation is the relative scaling of the adjacency matrix $A$, given by a constant $c$ times every element in $A$. For matrix $A$ multiplied by some scaling coefficient $c$, the natural dynamics along the network without control obey $\dot{\bm{x}}=cA\bm{x}$, which implies $\ddot{\bm{x}}=cA\dot{\bm{x}}=c^{2}A^{2}\bm{x}$, and $\bm{x}^{[k]}=c^{k}A^{k}\bm{x}$. In simplified networks, $A^{k}=0$ for $k>1$. Hence, for small $c$, the natural dynamics along the simplified and non-simplified networks are similar. The larger the scaling coefficient $c$, the more the non-simplified network dynamics deviate from the simplified network dynamics.

We analyzed the accuracy of this first-order approximation, and its dependence on the scaling coefficient $c$ and fraction of non-driver nodes, in brain networks of several species: the empirical meso-scale mouse connectome consisting of 112 interconnected brain regions from the Allen Brain Institute, a drosophila connectome consisting of 49 interconnected brain regions shih2015connectomics , and a set of human connectomes consisting of 116 brain regions interconnected by white matter tracts estimated using a 705-direction diffusion imaging scan over 55 minutes (for empirical details regarding connectivity estimates, see Methods; for a conceptual schematic of the full and simplified Drosophila connectome, see Fig. 2c–d). We normalized each connectome by the magnitude of its largest eigenvalue Betzel2016 , and selected a range of scaling coefficients $c$ and fractions of non-driver nodes $d$. For each combination of scaling coefficient $c_{i}$ and fraction of non-driver nodes $d_{j}$, we selected 1000 random permutations of drivers and non-drivers, and computed the minimum energy required to drive the simplified and non-simplified networks from initial states ${\bm{x}}_{\textup{d}}=0$, ${\bm{x}}_{\textup{nd}}=0$ to random final states ${\bm{x}}_{\textup{nd}}^{*}\in(-1,1)^{M},{\bm{x}}_{\textup{d}}^{*}\in(-1,1)^{N}$. We then calculated the median magnitude of the percent error for the control energies between the simplified and non-simplified networks (Fig. 2e–g for drosophila, mouse, and human, respectively).

We observed that the percent error magnitude between the control energies along the non-simplified versus simplified networks was similar across species. In general, the error remained below approximately 5% for scaling coefficients $c<1.5$, and fraction of non-driver nodes $d<0.4$, and below 10% for $c<3.3$, and $d<0.39$ (Fig. 2e–g), confirming that the first-order energy approximation is accurate within a range of scaling coefficients and non-driver fractions for these empirical connectomes. More generally, the accurate, closed-form expression relating a network’s topology and control energy can be used to extract the underlying topological features that determine controllability. For the remainder of this paper, we will use the same connectomes (drosophila, mouse, and human) at a scaling coefficient of $c=1$, and non-driver fraction $d\leq 0.4$, to ensure generalizability of our findings to the non-simplified versions of these same networks. For any examples requiring a specific fraction of non-drivers, we will use a fraction of 0.2, corresponding to $M=10$ in drosophila, $M=22$ in mouse, and $M=23$ in human connectomes.

III.2 Determinant of the driver-to-non-driver connection matrix scales the control energy

After deriving a closed-form approximation for the minimal total energy required to drive a network from one state to another, we next sought to provide a physical interpretation of the mathematical features that drive the control energy. First, we let $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}$, and notice that Eq. (2) can be rewritten as

[TABLE]

where $\bm{v}_{1}={\bm{x}}_{\textup{nd}}^{*}-\frac{1}{2}{A}_{\textup{21}}{\bm{x}}_{\textup{d}}^{*}$ and $\bm{v}_{2}={\bm{x}}_{\textup{d}}^{*}$, and $\text{adj}(Q)$ is the adjoint matrix of $Q$. We notice that, independently of the vectors $v_{1}$ and $v_{2}$, the determinant of $Q$ acts as a scaling factor for the total energy. This insight is useful because of the intricate geometric interpretation of a Gram matrix determinant. Specifically, let $\bm{a}_{i}\in^{1\times N}$ be the $i$-th row of ${A}_{\textup{21}}$ (which we will call the weight vector), representing the connections from all $N$ drivers to the $i$-th non-driver node (Fig. 3a). Then, the determinant of the gram matrix $Q$ is equal to the squared volume of the parallelotope formed by all $\bm{a}_{i}$.

To gain an intuition for these results, we show a simple system with 3 drivers and 2 non-drivers with varying network topologies in Fig. 3b–d, their corresponding geometric parallelotopes in Fig. 3e–g, with $\bm{a}_{1}$ as the vector of gray-colored weighted connections into ${\bm{x}}_{\textup{nd1}}$, and $\bm{a}_{2}$ as the tan-colored connections into ${\bm{x}}_{\textup{nd2}}$. The control task has initial states ${\bm{x}}_{\textup{d}}(0)=0$, ${\bm{x}}_{\textup{nd}}(0)=0$, and final states ${\bm{x}}_{\textup{d}}^{*}=0$, ${x}_{\textup{nd1}}^{*}=1,{x}_{\textup{nd2}}^{*}=-1$. These final states hard-code the empirical observation that the most common dynamics of intrinsic activity patterns in large-scale human brain networks are anti-correlated activation states chen2016 ; Medaglia2015a ; Reddy2017 , often referred to task-positive and task-negative fox2005 . We see that as the total area of the parallelogram shrinks from Fig. 3e–g, the total control energy to move the non-drivers increases in Fig. 3h–j.

Intuitively, if two non-drivers ${x}_{\textup{nd1}},{x}_{\textup{nd2}}$ are very similarly connected to the drivers, it is difficult to drive one of them independently from the other. Geometrically we show in Fig. 3 how similarity between two driver $\to$ non-driver connections $\bm{a}_{1},\bm{a}_{2}$ decreases the volume of the corresponding parallelotope, thereby decreasing the determinant in Eq. (3) and inversely scaling the control energy. We see that this relationship between the determinant and similarity of weight vectors $\bm{a}_{i}$ persists for any number of drivers and non-drivers where the first-order network is a good approximation. We conclude that the similarity in distribution of weights directed into non-driver nodes scales the control energy through the determinant of $Q$, where more similarly connected non-drivers require more energy to control differentially. This relationship is significant because we can now analyze and modify the connectivity of a network, knowing the topological features that determine its control.

III.3 Using connection topology to identify energetically favorable control nodes

In the previous section, we derived an approximate relationship between a network’s topology and its minimum control energy, and showed that the similarity in driver to non-drivers connections $\bm{a}_{i}$ changed the determinant of the gram matrix $Q$, thereby scaling the control energy. Here, we further explore this idea of “similarity” between connections $\bm{a}_{i}$, in order to quantify the impact of each individual non-driver on the control energy. With this knowledge, we can begin asking questions about the most or least controllable regions in a network, and how to modify specific connections in a network to improve controllability.

C.1. Derivation of the Main Topological Contributors to Control Energy. Our analysis is rooted in the intuition that the edge weights $\bm{a}_{i}$ that maximize the parallelotope volume, thereby facilitating network control, are large in magnitude and orthogonal to each other. Let $\lambda_{i}$ and $\bm{e}_{i}$ be the eigenvalues and eigenvectors of the matrix $Q$ in Eq. (3). We derive in Lemma V.6 the equivalent, alternative control energy expression written as

[TABLE]

where $w_{i}=\prod_{j\neq i}^{M}\lambda_{j}$, $c_{i}=\bm{e}_{i}^{T}\bm{v}_{1}$, and $\theta_{k}$ is the angle formed between $\bm{a}_{k}$ and the parallelotope formed by $\bm{a}_{j\neq k}$. For $N$ drivers and $M$ non-drivers, we can visualize the $M$ weight vectors $\bm{a}_{k}$ as forming a parallelotope in an $N$-dimensional space. The variable $\theta_{k}$ then represents the angle formed between $\bm{a}_{k}$ and the paralellotope formed by the remaining $M-1$ vectors $\bm{a}_{j\neq k}$. An example with $N=3,M=2$ is shown in Fig. 3e–g, where $\theta_{1}=\theta_{2}$ is the angle between the tan and gray vectors.

Here, we have segregated the control energy into a task-based ($\frac{\sum_{i=1}^{M}w_{i}c_{i}^{2}}{\sum_{i=1}^{M}w_{i}}$) and topology-based ($\sum_{k=1}^{M}{\frac{1}{\|\bm{a}_{k}\|^{2}\sin(\theta_{k})^{2}}}$) term. The task-based term is a weighted average of eigenvalue products $w_{i}$, weighted by the eigenvector composition $c_{i}$ of a specific task $\bm{v}_{1}$. The topology-based term is a sum of elements, where each element is an inverse squared function of the magnitude and angle of each weight vector. This segregation allows us to analyze the topology separate from the specific control task, and shows that each non-driver additively contributes to the total control energy. We note that for a non-driver $i$, this contribution is smallest when $\|\bm{a}_{i}\|$ and $\theta_{i}$ are large.

C.2. Most and least energetically favorable driver-nondriver sets in brain connectomes. To support this discussion, we used the expression Eq. (4) that segregated the control task from the topology in tandem with a greedy algorithm to find the sets of $M$ non-drivers that minimized and maximized this topology term. We first calculated the topology term for all permutations of 4 non-driver regions, and found the sets of 4 regions that minimized and maximized the topology-based term. Then, we iteratively appended individual regions that minimized and maximized the term until we reached $M$ non-drivers. We defined the most and least energetically favorable networks to be the selections of $N$ driver and $M$ non-driver nodes that minimize and maximize this topology term, respectively.

As an example using a non-driver fraction of 0.2, we show the distribution of magnitudes $\|\bm{a}_{i}\|$ of each driver to non-driver weight vector in drosophila, mouse, and human for both the energetically most and least favorable networks (Fig. 4a–c, respectively). We observe that the energetically least favorable networks have significantly weaker driver to non-driver connections than the energetically most favorable networks via a two-sample $t$-test between the most and least favorable networks in the drosophila ($t(8)=7.19$, $p=1.08\times 10^{-6}$), mouse ($t(20)=7.20$, $p=7.39\times 10^{-9}$), and human ($t(21)=8.22$, $p=1.93\times 10^{-10}$). We also show the corresponding distributions of angles between each driver to non-driver weight vector $\bm{a}_{i}$ and the parallelotope formed by the remaining $M-1$ non-drivers given by $\theta_{i}$ in Eq. (4) in the drosophila, mouse, and human connectomes (Fig. 4d–f, respectively). We observe that the angles $\theta_{i}$ in the energetically least favorable networks are significantly smaller than the angles in the energetically most favorable networks via a two-sample $t$-test between the most and least favorable networks in the drosophila ($t(8)=3.67$, $p=1.80\times 10^{-3}$), mouse ($t(20)=6.04$, $p=3.43\times 10^{-7}$), and human ($t(21)=3.71$, $p=5.72\times 10^{-4}$).

Next, we demonstrate the utility and robustness of these topological features for network control by computing the minimum control energy along the non-simplified network using the driver and non-driver designations from the simplified network in Eq. (4) for a range of non-driver fractions. For each $M$, we calculated the $\log_{10}$ energy from 2000 random control tasks with initial states ${\bm{x}}_{\textup{nd}}(0)=0$, ${\bm{x}}_{\textup{d}}(0)=0$, and final states ${\bm{x}}_{\textup{nd}}^{*}\in(-1,1)^{M},{\bm{x}}_{\textup{d}}^{*}\in(-1,1)^{N}$ on the most energetically favorable network, the least energetically favorable network, and a set of randomly chosen $N$ drivers and $M$ non-drivers. We show the control energies across the tasks for the drosophila, mouse, and human connectomes (see Fig. 4g–i, respectively). We also show the means and standard deviations of the control energy across the random tasks for each $M$ (see Fig. 4j–l). As can be seen across all three species, the most energetically favorable networks require around 0.5–1 order of magnitude less control energy than the random networks, and 2.5–4 orders of magnitude less control energy than the least energetically favorable networks. This difference indicates an energetic advantage for some configurations of drivers and non-drivers over others.

III.4 Brain Networks of Increasingly Complex Species have More Energetically Favorable Topological Relationships

Given the relationship between a network’s topology and minimum control energy in Eq. (4), we seek to understand if brain networks are organized along energetically favorable principles. Previously, we showed that the control energy was inversely proportional to the squared product of the magnitudes $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ of the driver $\to$ non-driver connections. Here, we show that brain networks of increasingly complex species more effectively balance these two topological features to yield robust and energetically favorable connectivities.

Fundamentally, we are asking how well a network’s specific set of topological components $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ combine to minimize the topology dependent energy term $\sum_{k=1}^{M}{\frac{1}{\|\bm{a}_{k}\|^{2}\sin(\theta_{k})^{2}}}$. In networks that are not designed along these energetic principles, we expect to see no particular relationship between $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$. In networks that minimize the topology dependent energy term, we expect to see a compensatory effect, where non-drivers with small angles have large magnitudes, and non-drivers with small magnitudes have large angles. While it may seem intuitive that pairing large magnitudes with large angles would yield lower energies, this strategy also requires smaller magnitudes to pair with smaller angles, yielding disproportionately large energy contributions.

To explore the relationship between $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ in brain networks, we selected 10,000 random permutations of non-drivers in the drosophila, mouse, and 10 human connectomes, at a non-driver fraction of 0.2. For each permutation, we calculated $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ for every non-driver. Then, we averaged $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ for each non-driver across all permutations, giving us an averaged magnitude $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ for each brain region in the drosophila, mouse, and each of 10 humans. Finally, we plotted the averaged $\sin(\theta_{k})$ vs $\|\bm{a}_{k}\|$ for all brain regions in the drosophila, mouse, and the average across all 10 human subjects (Fig. 5a–c). We find little relationship between the averaged $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ in the drosophila (Spearman $\rho=-0.25$, $p=0.0748$), a moderate negative relationship in the mouse ($\rho=-0.36$, $p=0.000125$), and a strong negative relationship in the human ($\rho=-0.73$, $p\approx 0$).

To graphically demonstrate how this negative $\sin(\theta_{k})$ vs. $\|\bm{a}_{k}\|$ relation might arise in networks, we show a simple 5 node network with two communities of 3 and 2 strongly interconnected sets of nodes (Fig. 5d). We first look at the average $\|\bm{a}_{k}\|$ and $\sin(\theta_{k})$ between a specific node in the 3 strongly interconnected set (Non-Driver 1, colored light blue in Fig. 5e) and all four permutations of one other non-driver (designated Non-Driver 2) as shown in Fig. 5e. We see that, because Non-Driver 1 is a member of 3 strongly interconnected nodes, it has on average 1.5 strong driver $\to$ non-driver edges across the 4 permutations, and 2/4 configurations where it and Non-Driver 2 are similarly connected (small angle). In contrast, we look at a specific node in the 2 strongly interconnected set (Non-Driver 1, colored light blue in Fig. 5f), and all four permutations of one other non-driver (Non-Driver 2, Fig. 5f). Here, we find that Non-Driver 1 only has an average of 0.75 strong driver $\to$ non-driver connections, and only 1/4 similarly connected configurations. Hence, on average, a non-driver among 3 strongly interconnected nodes will have stronger driver $\to$ non-driver connections (larger $\|\bm{a}_{k}\|$) and greater number of similarly connected configurations (smaller $\sin(\theta_{k})$) than a non-driver among 2 strongly interconnected nodes.

III.5 Network manipulation to facilitate control

In the previous section, we segregated the energy contribution of network topology into two key components (magnitude and angle), and found that brain networks of more complex species were organized in more energetically favorable ways. Here, we look to extend this concept to network modifications that lead to lower control energies. We focus on the effects of edge deletion since it is often useful in the study of biological systems such as brain Alstott2009 , metabolic Aristidou1994 ; Patnaik1994 , and gene regulatory Sander2014 networks. Perhaps counter-intuitively, we show that the deletion of certain edges improves the general controllability of the network.

Ultimately, we seek to understand if, and to what degree, edge lesioning may be used therapeutically to improve a network’s control profile. To that end, we quantified the effect of modifying each edge weight on the determinant, deleted edges that maximally increased the determinant, and demonstrated the corresponding impact on the total control energy. First, we let $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}$ be as in Eq. (3) and derived (Lemma V.5) that

[TABLE]

which characterized the rate of change of the determinant of $Q$ with respect to a change in any particular edge weight in ${A}_{\textup{21}}$. This metric was of particular utility in assessing the sensitivity of the control energy with respect to changes of the network’s edge weights.

Then, we used Eq. (5) to make informed modifications to the network topology to increase or decrease the network control energy. First, we randomly selected 10,000 permutations of non-drivers at a non-driver fraction of 0.2, and designated the remaining regions as drivers. For each permutation, we extracted the block matrix ${A}_{\textup{21}}$, calculated $2\det(Q)(Q^{-1}){A}_{\textup{21}}$, and found the element $a_{ij}\neq 0$ yielding the largest change based on Eq. (5). We then simulated an edge deletion by setting $a_{ij}=0$, and we repeated the process to obtain networks of 1, 2, 3, and 4 deleted edges. Finally, for each permutation, we performed a random control task on the non-simplified network with initial states ${\bm{x}}_{\textup{nd}}(0)=0$, ${\bm{x}}_{\textup{d}}(0)=0$, and final states ${\bm{x}}_{\textup{nd}}^{*}\in(-1,1)^{M},{\bm{x}}_{\textup{d}}^{*}\in(-1,1)^{N}$, and calculated the percent change in minimum energy between the pre- and post-modified networks for drosophila, mouse, and human connectomes Fig. 6a–d. As can be seen in Fig. 6a, the removal of a single weight can sometimes lead to more than a 50% reduction in control energy, while the removal of four edges (Fig. 6d) can sometimes lead to more than an 80% reduction in control energy. We also note that the human connectome, which we showed was already energetically favorable wired, also had the smallest percent decrease in energy from edge-deletion.

Here we use the scaling property of the determinant to make small, strategic topological changes that drastically reduce the control energy. We show that the deletion of one edge can yield as much as a 50% reduction in energy. We also show that, for the same topological modification, the drosophila experienced greater energy reduction than the mouse, which also experienced greater energy reduction than the human. This corresponds to the previous finding where, because brain networks of increasingly complex species are already energetically favorably wired, they may not experience as much improvement after modification. We emphasize that this was a purely topologically-motivated modification that did not cater to the specific end-state of the control task. This analysis may also yield a measure of how fragile or robust a specific network is to topological disruption, and introduces the perspective of targeted improvements in network controllability through topological changes.

IV Discussion

The control of networked systems is a critical frontier in science, mathematics, and engineering, as it requires a fundamental understanding of the mechanisms that drive network dynamics and subsequently offers the knowledge necessary to intervene in real-world systems to enhance or better their outcomes motter2015network . While some theoretical predictions are beginning to be made in nonlinear network systems Cornelius2013 , the overwhelming majority of recent advances have been made in the context of linear control liu2011controllability ; campbell2015topological ; ruths2014control . Yet, despite the significant efforts, some very basic intuitions regarding how edge weights impact control – either locally around a certain node or globally throughout the whole system – have remained elusive. Here, we sought to address this gap by segregating network nodes into either drivers or non-drivers, and examining the dynamic interactions between them. We show that for a wide range of parameters, the minimum energy required to control a network is mostly a function of these directed bipartite connections, offering a fundamental theory of driver$\to$non-driver interactions. We apply this framework to the inter-areal connectomes of the mouse oh2014mesoscale ; rubinov2015wiring , drosophila shih2015connectomics , and human to demonstrate that the predictions of the theory derived from the bipartite subgraph hold for the fully weighted, directed network. The work thus offers important insights into network analysis and design by delineating key topological principles in under-actuated network systems.

More specifically, we presented an equivalent first-order energy expression segregating the control energy into task-based and topology-based terms. We built upon our finding that similarly connected non-driver nodes decrease the gram determinant and increase the control energy by quantifying the connection similarity between non-driver $i$ and non-drivers $j\neq i$ in this topology term. Then we selected non-drivers with the topological goal of minimizing or maximizing this similarity, and showed that on average, networks with non-drivers that minimized this similarity required 2.5 – 4 orders of magnitude less energy than those that maximized this similarity. We also showed that the connectomes of increasingly complex species were topologically organized to be energetically favorable. Ultimately, we have shown that there is an inherent topological contribution to the total control energy, which can be used to inform and analyze the selection of drivers and non-drivers. We concluded by using these principles to show that the deletion of one edge in many thousands of potential edges can yield as much as a 50% reduction in energy, providing insights into how to target improvements in network controllability through fine-scale changes in topology.

IV.1 A novel conceptualization of network control

A distinct advantage of our approach is the focus on a physically meaningful topological understanding of the principles governing network control. Although spectral analysis of a network’s controllability Gramian TK:80 yields theoretically useful information about the overall behavior of the network under control pasqualetti2014controllability ; yan2015spectrum , it is not obvious how specific patterns of connectivity or selections of driver and non-driver nodes contribute to this behavior. Understanding this relationship is crucial when analyzing empirical biological networks such as the brain compared to purely mathematical networks, because the nodes and edges of a brain network often have known functions and perform known computations Lanteaume2007 ; Burgess2002 , and we are interested in understanding how these functions and computations modulate or influence one other.

We address this gap in understanding between the control behavior and topology of networks through a simplified network only involving connections from driver to non-driver nodes. This simplification is quite generally motivated by recent work demonstrating that relatively sparse network representations of complex biological systems olhede2014network ; park2015sparse ; liu2016sparse can contain much of the information needed to understand the system’s structure and dynamics clauset2008hierarchical ; yang2016predicting ; pan2016predicting ; zhu2015information ; lu2015toward ; navlakha2012network . More specifically, the simplification hard-codes the fact that energy can be transmitted directly from drivers to non-drivers along walks of length unity. By simplifying the complex network structure into driver $\to$ non-driver interactions, we reached a powerful closed-form bilinear mathematical approximation relating a network’s topology to the total control energy. We found that for a range of matrix scales and fractions of non-driver nodes, this simplification well approximates the minimum control energy of the unsimplified network. This implies that control dynamics along the first-order connections from driver $\to$ non-driver nodes dominate the dynamics along other connections within the viable range of parameters. These results inform our understanding of how much first-order connections contribute to the overall dynamics of generic network control systems goni2014resting .

We used this approach to demonstrate that the similarity and strength of connections from driver to non-driver nodes are key topological principles that govern controllability. To reach this conclusion, we posited a geometric perspective of network control by showing that the control energy was inversely proportional to the determinant of the gram matrix $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}$, and that this determinant was equal to the squared volume of the parallelotope formed by the vector of connections $\bm{a}_{k}$ into the $k$-th non-driver node. This principle allows us to make hypotheses about the function of a network given its structure. Brain regions that have very similar connection distributions form a parallelotope with a very small volume, implying that most state changes will be prohibitively energetically costly. Hence, the subspace of energetically viable state changes will be of fairly low dimension. However, very differentially connected regions will form a parallelotope with much larger volume, allowing the subspace of energetically viable state changes to be of higher dimension. These results inform the development of analytical constraints on the accessible state space of a networked system Cornelius2013 , particularly informing the set of states within which one might seek to push the brain using stimulation paradigms common in the treatment of neurological disorders and psychiatric disease chen2014harnessing ; chrysikou2011noninvasive . While many initial studies have examined unconstrained state spaces Gu2015 ; Betzel2016 ; muldoon2016stimulation , understanding viable states and state trajectories is critical for the translation of these ideas into the clinic bassett2017emerging .

Indeed, to provide further insights into the potential utility of these approaches in informing interventions in brain systems, we formally quantified the contribution of the network topology to the control energy as a sum of contributions from each non-driver, where each contribution was a function of the magnitude and angle of $\bm{a}_{k}$. This formulation allowed us to identify brain regions that were inherently costly to control, and therefore potentially should be chosen as stimulation targets with careful consideration. Moreover, we were able to use a greedy algorithm to find the set of $M$ non-drivers that maximized or minimized this topology-based energy component. Through this topology-based energy contribution, we obtained approximate answers to questions such as ”which sets of brain regions are easiest versus most difficult to control,” and ”which non-drivers should be designated drivers to most improve controllability.” Importantly, these insights lay the groundwork for the optimization of stimulation sites in neural systems, a problem that has received very little theoretical treatment, and is considered one of the current critical challenges in neuroengineering johnson2013neuromodulation .

Finally, we used the scaling principle of the gram determinant to make strategic, task-agnostic edge deletions that maximally increased the determinant. We saw that even the deletion of 1 edge occasionally produced a 50% reduction in total control energy, while the deletion of 4 edges occasionally produced an 80% reduction in total control energy. The first insight was that, even in an overdetermined, unsimplified system ($N>M$), a single edge deletion could produce such a profound improvement in the general controllability of a network. This sensitivity suggests that dynamical networks such as the brain can produce fairly drastic changes in dynamical behavior given minute changes in physiological topology, consistent with observations of critical dynamics in human and animal neurophysiology bassett2006adaptive ; rubinov2011neurobiologically ; shew2015adaptation ; deco2012ongoing . Moreover, these results also suggest that minor, targeted structural changes through concussive injury can lead to drastic changes in overall brain function hart2016graph ; caegenberghs2016mapping ; horn2016altered , via altering the controllability landscape of the brain gu2016optimal . The second interesting insight was that these topological modifications were task-agnostic edge deletions, signifying that even in a linear regime, the presence of an unfavorable edge can have a profoundly negative impact on the controllability of a network. We note that it is trivial to perform a similar analysis that takes into account the specific tasks $\bm{v}_{1},\bm{v}_{2}$ by taking the derivative of the full energy term $E_{total}$ with respect to ${A}_{\textup{21}}$, which would optimize the network topology for a specific task, as studied in more detail in Betzel2016 .

IV.2 Cross-species comparison of controllability in structural brain networks

Emerging neurotechnologies are uncovering the richness of brain connectivity with unprecedented detail bassett2016network . The fine-scale maps produced by these efforts make quantiative cross-species comparison tractable in a way that was not possible in previous years heuvel2016comparative ; reid2016cross . Here, we take advantage of these advances to address the question of whether and how species differ in brain network controllability. Importantly, the more general question of whether brain network architecture in different species harbors similar or distinct topological attributes is not a new one bassett2006small ; bassett2016small ; bassett2010efficient ; kaiser2006nonoptimal . The majority of work has focused on reporting cross-species similarities, while few have addressed the question of evolutionary drivers explaining differences across species capable of more or less complex function herculano2012remarkable . Here we find a monotonic gradation in the mapping between edge strength and topological angle across the 3 species, suggesting that increasingly complex species are more energetically favorably wired. Interestingly, the human, in addition to being most energetically favorably wired, also had the smallest percent decrease in energy following edge-deletion. These results point to an advantage of the human brain in supporting diverse network dynamics with small energetic costs, while remaining unexpectedly robust to perturbations. It will be interesting in the future to expand this analysis to the connectomes of other species as they become available, and also to examine brain network robustness in individuals sustaining traumatic brain injury gu2016optimal .

IV.3 Utility in informing control of brain networks

In this paper, we apply our theoretical framework to real-world data collected via tract-tracing in drosophila and mouse, and via diffusion imaging in human. Such brain networks represent particularly important contexts in which to understand control bassett2016network ; bassett2017emerging . Even at the microscale of individual neurons, neural control engineering schiff2011neural seeks to identify minimal energy control nabi2013minimum , that can potentially be used to explain homeostatic mechanisms controlling abonormal bursts of activity wiles2016autaptic or to develop exogeneous control strategies to terminate bursts wilson2014hamilton . Open questions surround the role of symmetries whalen2012observability ; whalen2015observability or synchronizability tang2016structural within the network that may constrain these control properties. At the macroscale of centimeter-sized brain regions, network control offers a novel perspective on how the brain controls itself, a characteristic known to cognitive neuroscientists as cognitive control Gu2015 , as well as which mental states might be preferred Betzel2016 , and how both features may be altered following traumatic brain injury gu2016optimal . Indeed, in clinical populations the need for realistic, low-energy control systems is particularly pressing, for example to restore neural function in Parkinson’s disease santaniello2015therapeautic , to supress bursting activity in coma ching2012neurophysiological ; ching2013real , and to control the distributed propagation of seizures ching2012distributed in epilepsy, which has become known as an inherently network-based disorder burns2014network ; khambhati2015dynamic ; khambhati2016virtual . Our results offer a novel framework in which to address these questions and challenges, as well as tools to design interventions (the strengthening and weakening of neural connections) to facilitate optimal control and enhance therapeutic benefit bassett2017emerging .

IV.4 Theoretical considerations and methodological limitations

The work is built on several important assumptions that bring with them significant theoretical considerations. First, we only consider first-order connections from driver nodes to non-driver nodes. In other words, we study paths of length $1$ between drivers and non-drivers; we do not study the propagation of control energy through longer paths or walks in the network. This is an inherent limitation of the work, as it is clear from prior studies that the strength of long walks through the network has a non-trivial impact on the control energy Betzel2016 ; gu2016optimal . To address this limitation, we demonstrate that this simplification enables us to better understand the dependence of control energy on network topology. Moreover, we show that these first-order control dynamics can offer reasonable approximations for non-simplified networks constructed from real-world neuroimaging data in drosophila, mouse, and human.

Second, we assume that these networks adhere to linear dynamics liu2011controllability ; muller2011few ; yan2015spectrum , which may limit the applicability of the results to (i) linear systems, or (ii) nonlinear systems for which control is sought in short time intervals, enabling a linearization of the dynamics around the operating point luenberger1979introduction . Moreover, in the context of our application to understand neural architecture, this choice is consistent with prior brain network control studies, such as Gu2015 ; Betzel2016 ; muldoon2016stimulation , as well as prior mathematical modeling studies on human neuroimaging data galan2008network ; honey2009predicting .

Third, we assume that the input functions are chosen in such a way as to minimize both the control energy and the distance of the current state from the target state Betzel2016 ; gu2016optimal . Importantly, these assumptions represent the best case scenario, offering a lower-bound on the relationships between topology and control energy. Additional variables of interest to examine in future studies include the time scales of control and the tortuosity of the trajectory.

Finally, as with all empirical data, the connectomes of the drosophila, mouse, and human are fundamentally estimates of the true large-scale brain connectivity, and emerging neurotechnologies will likely offer increasingly accurate estimates. A specific limitation of the human connectome that is important to mention is that the diffusion imaging data must be submitted to sophisticated tractography algorithms to construct region-to-region estimates of structural connectivity. While evolving at a swift pace, these algorithms may still report spurious tracts or fail to report existing tracts (thomas2014anatomical, ; reveley2015superficial, ; pestilli2014evaluation, ). However, these issues are somewhat mitigated by the fact that we use an exceptionally high-resolution scan, capitalizing on a multiband sequence taking place over 55 minutes, and estimating diffusion over 705 directions, thereby increasing the resolution of the data by an order of magnitude over most existing studies (which estimate diffusion over 30-64 diffusion directions).

IV.5 Conclusion and future directions

In closing, we note that the natural direction in which to take this work will be to include higher order interactions in the bipartite framework, and further to expand the bipartite framework to include driver$\to$driver and non-driver$\to$non-driver interactions. Moreover, it would be interesting to apply this reduced framework to random graphs and other well-known benchmarks – both from a mathematical perspective bollobas1985random and also in the context of neural systems klimm2014resolving ; sizemore2016classification – to better understand the phenotypes present in those graph ensembles. Third and finally, informing the design of new networks with these tools may be particularly useful in neuromorphic computing pedroni2016mapping ; pfeil2013six , material science papadopoulos2016evolution ; giusti2016topological , and other contexts where optimal control of physical systems is of paramount importance.

IV.6 Acknowledgements

JK acknowledges support from NIH T32-EB020087. JMS and DSB acknowledge support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the U.S. Army Research Laboratory and the U.S. Army Research Office through contract numbers W911NF-10-2-0022 and W911NF-14-1-0679, the National Institute of Health (2-R01-DC-009209-11, 1R01HD086888-01, R01-MH107235, R01-MH107703, R01MH109520, 1R01NS099348 R21-M MH-106799, and T32-EB020087), the Office of Naval Research, and the National Science Foundation (BCS-1441502, CAREER PHY-1554488, BCS-1631550, and CNS-1626008). AEK and JMV acknowledge support from the U.S. Army Research Laboratory contract number W911NF-10-2-0022. FP acknowledges support from the National Science Foundation (BCS-1430280 and BCS 1631112). The content is solely the responsibility of the authors and does not necessarily represent the official views of any of the funding agencies.

V Supplement

V.1 Connectome Data

Drosophila Connectome. The full reconstruction of the Drosophila connectome can be found in the FlyCircuit 1.1 database chiang2011three ; shih2015connectomics . This database contains images of 12,995 neurons, as well as their projections, that are characteristic of the Drosophila female. In this database, each neuron was labeled using green fluorescent protein (GFP) and its location was estimated from 3-dimensional images that were co-registered to a template using a rigid linear transform. To obtain a mesoscale representation of this fine-scale data, neurons were assigned to one of 49 local populations, based on their morphology and known functions. Following the original work from Shih and colleagues, we treated each of these 49 populations as the nodes of the network, and we treated the directed, weighted edges between populations as network edges shih2015connectomics .

Mouse Inter-Areal Connectome Data. In addition to drosophila connectome, we also analyzed the inter-areal connectome of the mouse. In particular, we use the exact network studied in rubinov2015wiring , which was reconstructed from original tract-tracing data recently released by the Allen Brain Institute oh2014mesoscale . The entire brain was separated into 112 regions, which we treat as network nodes. Each pair of regions was then linked by directed edges that encoded the presence or absence of inter-regional projections. The weight of each edge was defined by the number of projections normalized by the volumes of the two regions being connected.

Human Diffusion Imaging Data. Ten healthy adult human subjects (m) were imaged as part of an ongoing data collection effort at the University of Pennsylvania; the subjects provided informed consent in writing, in accordance with the Institutional Review Board of the University of Pennsylvania. All scans were acquired on a Siemens Magnetom Prisma 3 Tesla scanner with a 64-channel head/neck array at the University of Pennsylvania. Each data acquisition session included both a diffusion spectrum imaging (DSI) scan as well as a high-resolution T1-weighted anatomical scan. The diffusion scan was 730-directional with a maximum $b$-value of 5010s/mm2 and TE/TR = 102/4300 ms, which included 21 $b=0$ images. Matrix size was 144$\times$144 with a slice number of 87. Field of view was 260$\times$260mm2 and slice thickness was 1.80mm. Acquisition time per DTI scan was 53:24min, using a multi-band acceleration factor of 3. The anatomical scan was a high-resolution three-dimensional T1-weighted sagittal whole-brain image using a magnetization prepared rapid acquisition gradient-echo (MPRAGE) sequence. It was acquired with TR = 2500 ms; TE=2.18 ms; flip angle = 7 degrees; 208 slices; 0.9mm thickness.

DWI is highly sensitive to subject movement(Yendiki:2013ez, ), which can cause significant distortions in the reconstructed ODFs if not corrected. Motion correction is typically applied by determining an affine or non-linear transform to align each DWI volume to a reference derived from the high-signal $b=0$ images. The high b-values used in DSI present a problem for this approach, as the low signal in many of the volumes leads to poor registration. To address this, we interspersed $b=0$ volumes in the scan sequence, one for every 35 volumes. An initial average template was produced by averaging the $b=0$ images together and then improved by registering the $b=0$ images to the initial template and re-averaging. Each $b=0$ was finally re-registered to the improved template, and then each volume in the DSI scan was then motion corrected by applying the transformation calculated for the closest $b=0$ volume. Motion correction also impacts the effective $b$-matrix directions since the rotated images are no longer aligned with the scanner; therefore the transforms applied to motion correct each volume were also used to rotate the corresponding $b$-vectors.(Leemans:2009jn, ) The processing pipeline was implemented using Nipype(Gorgolewski:2011ex, ) with registration performed using the Advanced Normalization Tools (ANTs)(Avants:2011kk, ).

Using DSI-Studio (http://dsi-studio.labsolver.org), orientation density functions (ODFs) within each voxel were reconstructed from the corrected scans using GQI (Yeh:2010jq, ). We then used the reconstructed ODFs to perform a whole-brain deterministic tractography using the derived QA values in DSI-Studio (Yeh:2013fa, ). We generated 1,000,000 streamlines per subject, with a maximum turning angle of 35 degrees(Bassett2011structure, ) and a maximum length of 500mm(Cieslak2014, ). By holding the number of streamlines between participants constant, we use the number of streamlines that connect brain region pairs as an estimate of the strength of the connection and examine individual variability in structural connectivity(Griffa:2013gn, ).

To examine the relationship between structural connectivity and individual differences in learning rate, we constructed networks for each subject where nodes are atlas regions and edges are the measured connection strength between region pairs (Hagmann:2008gd, ). The nodes of the network were derived from spatially-defined regions of a brain atlas. We chose the anatomically-defined AAL atlas, originally developed in Statistical Parametric Mapping (SPM) (TzourioMazoyer:2002bi, ), which divides each brain hemisphere into 45 regions. We used a version in MNI-space that was then warped into subject-specific space using ANTs. Edges of the network were derived from streamlines that started and ended between the region pair and excluded streamlines that passed through one or both of the regions.

V.2 Mathematical Framework

Here we reiterate our mathematical notation and assumptions in greater detail, and provide lemmas for the main results. Consider a network represented by the directed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{1,\dots,n\}$ and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ are the sets of network vertices and edges, respectively. Let $a_{ij}\in$ be the weight associated with the edge $(i,j)\in\mathcal{E}$, and let $A=[a_{ij}]$ be the weighted adjacency matrix of $\mathcal{G}$. We associate a real value (state) with each node, collect the nodes’ states into a vector (network state), and define the map $\bm{x}:_{\geq 0}\rightarrow^{n}$ to describe the evolution (dynamics) of the network state over time. We let the network dynamics be linear and time invariant, as described by the equation

[TABLE]

We are particularly interested in characterizing how the network structure $\mathcal{G}$ influences the control properties of the dynamical system (6). To this aim, we assume that a subset of $N$ nodes, called drivers, is independently manipulated by external controls and, without loss of generality, we reorder the network nodes such that the $N$ drivers come first. Thus, the network dynamics with controlled drivers read as

[TABLE]

where ${\bm{x}}_{\textup{d}}$ and ${\bm{x}}_{\textup{nd}}$ are the state vectors of the driver and non-driver nodes, $A_{11}\in^{N\times N}$, $M=n-N$, $A_{12}\in^{N\times M}$, $A_{21}\in^{M\times N}$, $A_{22}\in^{M\times M}$, $I_{N}$ is the $N$-dimensional identity matrix, and $\bm{u}:_{\geq 0}\rightarrow^{N}$ is the control input.

We say that the network is controllable at time $T\in_{\geq 0}$ if, for any pair of states ${\bm{x}}_{\textup{d}}^{*}$ and ${\bm{x}}_{\textup{nd}}^{*}$, there exists a control input $u$ for the dynamics (1) such that ${\bm{x}}_{\textup{d}}(T)={\bm{x}}_{\textup{d}}^{*}$ and ${\bm{x}}_{\textup{nd}}(T)={\bm{x}}_{\textup{nd}}^{*}$. We refer the interested reader to TK:80 for a detailed discussion and rigorous conditions for the controllability of a system with linear dynamics. Finally, we define the energy of $u$ as

[TABLE]

where $u_{i}$ is the $i$-th component of $\bm{u}$. In what follows we characterize how the network topology and weights determine the control energy needed for a given control task. We restrict our analysis to a class of bipartite networks, as specified in the following assumptions. We remark that our assumptions, although restrictive, allow us to thoroughly predict how driver $\to$ non-driver connections facilitate or inhibit network control even in more complex network models (see Section V.3).

Assumption 1**.**

The network initial state satisfies ${\bm{x}}_{\textup{d}}(0)=0$ and ${\bm{x}}_{\textup{nd}}(0)=0$. $\square$

Assumption 2**.**

The network $\mathcal{G}$ contains only edges from the drivers to the non-drivers, that is, $A_{11}=0$, $A_{22}=0$, and $A_{12}=0$. Thus, the dynamics (1) simplify to ${\dot{\bm{x}}}_{\textup{d}}=\bm{u}(t)$, and ${\dot{\bm{x}}}_{\textup{nd}}=A_{21}{\bm{x}}_{\textup{d}}$. $\square$

An example of the original network (Fig. 7a) and the network satisfying our assumption (Fig. 7b) are shown. From assumptions 1 and 2 we readily observe that

[TABLE]

We see that ${\bm{x}}_{\textup{nd}}$ provides an integral constraint to ${\bm{x}}_{\textup{d}}$, and represent the specific values of the constraints as

[TABLE]

where $\bm{C}\in^{N\times 1}$ is an $N$-dimensional vector of real-valued constants. Furthermore, the set of controllable states can be characterized as follows. For a matrix $M$, let $\operatorname{Im}(M)$ and $\text{Rank}(M)$ denote the image and rank of $M$, respectively CDM:01 .

Lemma V.1**.**

(Controllability)* The network (1) is controllable if and only if $\text{Rank}(A_{21})=M$. Furthermore, the set of controllable states is $\operatorname{Im}\left(\begin{bmatrix}I_{N}&0\\ 0&A_{21}\end{bmatrix}\right)$.*

Proof.

Notice that the controllability matrix of (1) is

[TABLE]

and recall that a state is controllable if and only if it belongs to the range space of the controllability matrix. ∎

Lemma V.2**.**

(Minimum Energy Control Input)* The ${i}^{\textup{th}}$ driver trajectory ${x}_{\textup{d{i}}}(t)$ that minimizes the control energy takes the form ${x}_{\textup{d{i}}}(t)=a_{i}t^{2}+b_{i}t$*

Proof.

Recall from (1) and assumption 2 that ${\dot{x}}_{\textup{d{i}}}(t)=u_{i}(t)$. We minimize the energy

[TABLE]

where ${x}_{\textup{d{i}}}(0)=0,{x}_{\textup{d{i}}}(T)={x}_{\textup{d{i}}}^{*}$. From (9), ${x}_{\textup{d{i}}}$ is also subject to some integral constraint

[TABLE]

where $C_{i}$ is the $i$th element of $\bm{C}$. We see this naturally takes the form of the isoperimetric problem in the calculus of variations, which finds

[TABLE]

where $F(t,{x}_{\textup{d{i}}},\dot{x}{}_{\textup{d{i}}})=\dot{x}{}_{\textup{d{i}}}(t)^{2}$, $a=0$, and $b=T=1$, constrained by

[TABLE]

where $G(t,{x}_{\textup{d{i}}},\dot{x}{}_{\textup{d{i}}})={x}_{\textup{d{i}}}(t)-C_{i}$. The trajectory ${x}_{\textup{d{i}}}^{*}(t)$ which locally minimizes the cost function must satisfy the necessary (Euler-Lagrange) and sufficient (Jacobi) conditions. The Euler-Lagrange equation reads

[TABLE]

which, after substituting $F$ and $G$, yields

[TABLE]

to give the only extremal solution satisfying assumption 1

[TABLE]

where $\lambda$ is the lagrange multiplier. Because $(F+\lambda G)_{xx}=(F+\lambda G)_{x\dot{x}}=0$, and $(F+\lambda G)_{\dot{x}\dot{x}}=2$, the Jacobi condition becomes

[TABLE]

which holds true for any arbitrary smooth function $\eta$, where $\eta(0)=\eta(1)=0$. As ${x}_{\textup{d{i}}}^{*}(t)=at^{2}+bt$ is the only extremal function, and is also minimum, ${x}_{\textup{d{i}}}$ is the global minimum of the constrained control energy.

∎

Lemma V.3**.**

(Minimum Control Energy)* The required control energy for the $i^{th}$ driver is $E_{i}=12C_{i}^{2}-12C_{i}{x}_{\textup{d{i}}}+4{x}_{\textup{d{i}}}$.*

Proof.

We recall that the energy required to drive ${x}_{\textup{d{i}}}$ is

[TABLE]

We solve for $a_{i}$ and $b_{i}$ via the final state and integral constraint to yields equations

[TABLE]

from which we get

[TABLE]

substituting $a_{i}$ and $b_{i}$ into the equation for $E_{i}$, we get

[TABLE]

∎

Lemma V.4**.**

(Total Control Energy)* The total control energy is $E_{total}=12\bm{v}_{1}^{T}Q^{-1}\bm{v}_{1}+\bm{v}_{2}^{T}\bm{v}_{2}$, where $\bm{v}_{1}={\bm{x}}_{\textup{nd}}^{*}-\frac{1}{2}{A}_{\textup{21}}{\bm{x}}_{\textup{d}}^{*}$, $\bm{v}_{2}={\bm{x}}_{\textup{d}}^{*T}{\bm{x}}_{\textup{d}}^{*}$, and $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}$*

Proof.

Here, we use the method of Lagrange multipliers to minimize the total energy $f(\bm{C})$ as a function of $\bm{C}$ given in (9), constrained by $g(\bm{C})$ given by (8). We can write the total energy as

[TABLE]

with $M$ constraining equations, the set of which are given by

[TABLE]

where the $k^{th}$ constraint $g_{k}(\bm{C})$ is given by the $k^{th}$ row of $g(\bm{C})$. The method of Lagrange multipliers defines the Lagrangian given by

[TABLE]

where $\lambda_{k}$ is the $k^{th}$ Lagrange multiplier to compose $\bm{\lambda}\in^{M\times 1}$, and sets the gradient of the Lagrangian to 0

[TABLE]

which allows us to solve for $\bm{C}$ with respect to $\bm{\lambda}$

[TABLE]

By substituting $\bm{C}$ into the total energy equation and grouping terms, we get a preliminary formulation of $E(\bm{u})$ with respect to $\bm{\lambda}$

[TABLE]

To solve for $\bm{\lambda}$, we substitute the expression for $\bm{C}$ into our constraint equations $g(\bm{C})$ to yield

[TABLE]

and solve for $\bm{\lambda}$

[TABLE]

Substituting $\bm{\lambda}$ into $E(\bm{u})$, we get

[TABLE]

∎

Lemma V.5**.**

(Derivative of Gram Matrix)* The determinant of the gram matrix $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}$ with respect to the elements of ${A}_{\textup{21}}$ is $\frac{\partial}{\partial{A}_{\textup{21}}}\det(Q)=2\det(Q)(Q^{-1}{A}_{\textup{21}})$*

Proof.

For ${A}_{\textup{21}}$, we note the matrix determinant derivative identity

[TABLE]

If we set $B=I$, this simplifies to

[TABLE]

We note that

[TABLE]

which ultimately yields

[TABLE]

∎

Lemma V.6**.**

(Gram Vector Decomposition)* For system matrix ${A}_{\textup{21}}$ with linearly independent rows $\bm{a}_{j}$, and symmetric positive definite gram matrix $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}=PDP^{T}$ with eigenvalues $\lambda_{i}$ and eigenvectors $\bm{e}_{i}$, the total control energy can be represented by $E(\bm{u})=12\left(\frac{\sum_{i=1}^{M}w_{i}c_{i}^{2}}{\sum_{i=1}^{M}w_{i}}\right)\sum_{k=1}^{M}{\frac{1}{\|\bm{a}_{k}\|^{2}\sin(\theta_{k})^{2}}}+\bm{v}_{2}^{T}\bm{v}_{2}$, where $w_{i}=\prod_{j\neq i}^{M}\lambda_{j}$, $c_{i}=\bm{e}_{i}^{T}\bm{v}_{1}$, and $\theta_{i}$ is the angle formed by $\bm{a}_{k}$ and the sub-parallelotope formed by the remaining $\bm{a}_{j\neq k}$. .*

Proof.

We recall that the total control energy is given by

[TABLE]

We multiply each $k$ term to find a common denominator to yield

[TABLE]

We note that the left term is just a weighted average, with weights $w_{i}=\prod_{j\neq i}^{M}\lambda_{j}$. We also note that $\sum_{k=1}^{M}\prod_{l\neq i}^{M}\lambda_{l}$ is the $M^{th}$ term of the characteristic polynomial of $Q$, which is equivalent to $\sum_{k=1}^{M}\det(Q_{kk})$, where $Q_{kk}$ represents the $(k,k)$ minor of $Q$. Hence, we write

[TABLE]

We make use of the geometric fact that the determinant of $Q={A}_{\textup{21}}{A}_{\textup{21}}^{T}$ is equal to the squared volume of the parallelotope formed by the rows of ${A}_{\textup{21}}$. We also note that minor $Q_{kk}$ is the gram matrix of ${A}_{\textup{21}}$ after removing $\bm{a}_{k}$, represented by ${A}_{\textup{21}}^{k*}$. Therefore the ratio of the determinants of $Q_{kk}$ and $Q$ becomes the squared ratio of parallelotope volumes with and without $\bm{a}_{k}$.

[TABLE]

Finally, we realize that the contribution of $\bm{a}_{k}$ to the parallelotope volume is by a multiple of $\|\bm{a}_{k}\|\sin(\theta_{k})$, where $\theta_{k}$ is the angle formed by $\bm{a}_{k}$ and the sub-parallelotope, given by $vol({A}_{\textup{21}})=vol({A}_{\textup{21}}^{k*})\|\bm{a}_{k}\|\sin(\theta_{k})$ to yield

[TABLE]

∎

V.3 Validity of the First-Order Approximation

Until now, we have derived several useful closed-form expressions from the first-order minimum energy approximation by making the simplifying assumptions 1, 2. We explore how well this energy approximation holds when we relax the topological assumption 2. As a generalized, analytic closed-form energy solution for non-simplified networks is typically intractable or not informative, we compare the first-order energy approximation to a numerical computation of the unsimplified control energy. From linear control theory, we know that for an LTI system $\dot{\bm{x}}=A\bm{x}+B\bm{u}$ obeying the dynamics in (1), we can define the Reachability Gramian:

[TABLE]

The minimum control energy takes the form

[TABLE]

from which we calculate the percent error between the minimum energy of the simplified vs. unsimplified network. The results of this numerical analysis is shown in Fig. 2b.