Role of Graph Architecture in Controlling Dynamical Networks with Applications to Neural Systems

I. NETWORK TOPOLOGY AND CONTROLLABILITY

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 ∈ ℝ be the weight associated with the edge (i, j) ∈ ℰ, 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 x : ℝ_≥0 → ℝⁿ 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 read as

where x_d and x_nd are the state vectors of the driver and non-driver nodes, A₁₁ ∈ ℝ^N×N, M = n − N, A₁₂ ∈ ℝ^N×M, A₂₁ ∈ ℝ^M×N, A₂₂ ∈ ℝ^M×M, I_N is the N-dimensional identity matrix, and u: ℝ_≥0 → ℝ^N is the control input.

We will use the word controllable to refer to networks that are point-to-point controllable at time T ∈ ℝ_≥0 if, for any pair of states ${\mathit{x}}_{\mathrm{d}}^{\ast }$ and ${\mathit{x}}_{\mathrm{nd}}^{\ast }$, there exists a control input u for the dynamics Eq. (1) such that ${\mathit{x}}_{\mathrm{d}}(T)={\mathit{x}}_{\mathrm{d}}^{\ast }$ and ${\mathit{x}}_{\mathrm{nd}}(T)={\mathit{x}}_{\mathrm{nd}}^{\ast }$. For a detailed discussion and rigorous conditions for the controllability of a system with linear dynamics, see [30]. We define the energy of u as

where u_i is the i-th component of u. The energy of u_i can be thought of as a quadratic cost that penalizes large control inputs.

In the context of the brain, we approximate 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, see [23, 31, 32]). We specifically study the empirical inter-areal meso-scale connectomes of the mouse (112 brain regions, example schematic in Fig. 1g, h) from the Allen Brain Institute, the Drosophila (49 brain regions) [29], and a set of human connectomes (116 brain regions) interconnected by white matter tracts (for empirical details regarding connectivity estimates, see supplementary results XA).

II. PREDICTING CONTROL ENERGY

We seek an accurate, tractable relationship between the energy required to drive a network to a specific state and its connectivity. We begin with the original, non-simplified network (Fig. 2a) involving edges between all nodes, and consider dynamics along the simplified network (Fig. 2b) involving only edges from the driver to the non-driver nodes (for a conceptual schematic of the full and simplified Drosophila connectome, see Fig. 2c–d). We then derive an approximation of the minimum control energy (Lemma X.2 – X.4) by assuming that x_d(0) = 0, x_nd(0) = 0 (Assumption 1), and A₁₁ = 0, A₁₂ = 0, and A₂₂ = 0 (Assumption 2) in Eq. (1), which reads as

We make Assumption 1 because we are interested in the change in brain state through control, and consider initial conditions x_d(0) = 0, x_nd(0) = 0 to be a neutral baseline. Because Eq. (2) only involves edges from driver to non-driver nodes, we call Eq. (2) a first-order approximation to the minimum control energy of the non-simplified network Eq. (1). Importantly, this approximation requires at least as many driver nodes as non-driver nodes for $A_{21}{A}_{21}^T$ to be invertible (i.e. N ≥ M). To assess the accuracy of our expression, we look to classic results in the mathematical theory of systems and control [30], where the spectral properties of the reachability Gramian $W_R(0,T)={\displaystyle {\int }_0^T{e}^{At}B{B}^T{e}^{A^{T}t}dt}$ quantify the minimum amount of energy (Section XI A 2) to control the non-simplified network Eq. (1).

In these brain networks, we observe that the first-order energy approximation is accurate across a range of parameters, which are the magnitude of the adjacency matrix (given by the magnitude of the largest eigenvalue, c = ‖λ_max‖ after multiplying A by a constant scalar), and the fraction d of nodes selected as non-driver nodes (Fig. 2e–g). The error remains below approximately 5% for scaling c < 1.5 and non-driver fraction d < 0.4 (Fig. 2e–g). In this paper, we will use these connectomes scaled such that c = ‖λ_max‖ = 1, and non-driver fraction d ≤ 0.4, to ensure generalizability of our findings to the non-simplified versions of these same networks.

III. DETERMINANT OF THE DRIVER-TO-NON-DRIVER NETWORK

After deriving a closed-form approximation for the minimal energy to control a network, we seek a physical interpretation of the mathematical features that predict the control energy. We let $Q=A_{21}{A}_{21}^T$, and write Eq. (2) as

where ${\mathit{v}}_1={\mathit{x}}_{\mathrm{nd}}^{\ast }-\frac{1}{2}{A}_{21}{\mathit{x}}_{\mathrm{d}}^{\ast }$ and ${\mathit{v}}_2={\mathit{x}}_{\mathrm{d}}^{\ast }$, and adj(Q) is the adjugate matrix of Q. We notice that the determinant of Q acts as a scaling factor for the total energy. This insight is useful because of the geometric interpretation of a Gram matrix determinant. Specifically, let a_i ∈ ℝ^1×N be the i-th row of A₂₁ (which we will call the weight vector), representing weights 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 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, and their corresponding geometric parallelotopes in Fig. 3e–g with weight-vector a₁ in gray and a₂ in tan. We also compute the distribution of control energy required to drive each network from initial states x_d = 0, x_nd = 0 to 10,000 random final states ${\mathit{x}}_{\mathrm{nd}}^{\ast }\in {(-1,1)}^M$, ${\mathit{x}}_{\mathrm{d}}^{\ast }\in {(-1,1)}^N$ in Fig. 3h. As the non-drivers x_nd1, x_nd2 become more similarly connected, the total area of the parallelotope (and corresponding Gram determinant) decreases (Fig. 3e–g), and the control energy increases (Fig. 3h). We note that this determinant relationship persists for any number of nodes where N > M. We conclude that the similarity between weight-vectors generally scales the control energy through det(Q), allowing us to analyze and modify the connectivity of a network with respect to its control energy.

IV. IDENTIFYING ENERGETICALLY FAVORABLE CONTROL NODES

Here, we further explore the idea of “similarity” between connections a_i, to quantify the impact of each individual non-driver on the control energy.

V. COMPLEX BRAIN NETWORKS ARE ENERGETICALLY FAVORABLE

Given the relationship between a network’s connectivity and minimum control energy in Eq. (4), we seek to understand if brain networks are organized along energetically favorable principles. Fundamentally, we ask how well a network’s specific set of connectivity features ║a_k║ and sin(θ_k) combine to minimize the topology-dependent energy term $\left({\displaystyle {\sum }_{k=1}^M\frac{1}{{\Vert {\mathit{a}}_k\Vert }^2\sin {({\theta }_k)}^2}}\right)$. In networks that are not designed along these energetic principles, we expect to see no particular relationship between ║a_k║ and sin(θ_k). In networks that minimize the topology dependent energy term, we expect a compensatory effect, where non-drivers with small angles have large magnitudes, and vice versa.

To explore the relationship between ║a_k║ and sin(θ_k) in brain networks, we selected 10,000 random permutations of non-drivers in each of the Drosophila, mouse, and 10 human connectomes, at non-driver fraction d. For each permutation, we calculated ║a_k║ and sin(θ_k) for every non-driver. Then, we averaged ║a_k║ and sin(θ_k) for each non-driver across all permutations, giving us an averaged magnitude ║a_k║ and sin(θ_k) for each brain region in each network. Finally, we plotted the averaged sin(θ_k) versus ka_kk for all brain regions in each network for d = 0.2 (Fig. 5a–c). We find little relationship between the averaged ║a_k║ and sin(θ_k) in the Drosophila (Spearman ρ = −0.25, p = 0.0748), a moderate negative relationship in the mouse (ρ = −0.36, p = 0.000125), and a strong negative relationship in the human (ρ = −0.73, p ≈ 0). This ordering holds for a wide range of non-driver fractions (Fig. 5d). We graphically demonstrate how this negative sin(θ_k) versus ║a_k║ relation might arise in networks, using a simple 5-node network with two communities of 3 and 2 strongly interconnected sets of nodes (Fig. 5d-f), which has a strong negative relationship (Fig. 5h).

VI. NETWORK MANIPULATION TO FACILITATE CONTROL

Here, we consider 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 [33], metabolic [34], and gene regulatory [35] networks. Specifically, we quantify the effect of modifying each edge weight on the determinant in Lemma X.5 as

and compute the decrease in control energy as a result of deleting edges that maximally increase the determinant.

First, for each species and each of a range of non-driver fractions, we randomly selected 2,000 permutations of non-drivers. For each permutation, we extracted the block matrix A₂₁, calculated 2 det(Q)(Q⁻¹)A₂₁, and found the element a_ij ≠ 0 yielding the largest increase in det(Q) based on Eq. (6). We then simulated an edge deletion by setting a_ij = 0, and repeated the process to obtain networks of 1, 2, 3, and 4 deleted edges. Finally, we computed the percent change in control energy required to bring the non-simplified network from initial states x_nd(0) = 0, x_d(0) = 0, to final states ${\mathit{x}}_{\mathrm{nd}}^{\ast }\in {(-1,1)}^M$, ${\mathit{x}}_{\mathrm{d}}^{\ast }\in {(-1,1)}^N$ before and after edge deletion (Fig. 6a–d).

As can be seen in Fig. 6a, the removal of one edge can sometimes lead to more than a 10% average reduction in control energy, while the removal of four edges (Fig. 6d) can sometimes lead to more than a 30% reduction. Across most non-driver fractions, 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 these increasingly complex species are already energetically favorably wired, they may not experience as much improvement after modification.

VII. CONTRIBUTION AND FUTURE DIRECTIONS

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 better their outcomes [36]. While some theoretical predictions exist in nonlinear network systems [15], the majority of recent advances have been made in the context of linear control [21, 22]. Nevertheless, basic intuitions regarding how edge weights impact control have remained elusive. Although spectral analysis of a network’s controllability Gramian [30] yields theoretically useful information about the overall behavior of the network under control [37], 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, where nodes and edges often have known functions [38] that may modulate or influence one other.

A distinct advantage of our approach is the focus on a physically meaningful topological understanding of the principles governing network control. We map control behavior to network topology through a simplified network only involving connections from driver to non-driver nodes. This simplification hard-codes the fact that energy can be transmitted directly from drivers to non-drivers along walks of length unity, and is motivated by recent work demonstrating that relatively sparse network representations of complex biological systems [39, 40] can contain much of the information needed to understand the system’s structure and dynamics [41, 42]. Our results inform our understanding of how much first-order connections contribute to the overall dynamics of our network control systems. Moreover, they inform the development of analytical constraints on the accessible state space of a networked system, 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 [43, 44]. While many initial studies have examined unconstrained state spaces [23, 25, 26], understanding viable states and state trajectories is critical for the translation of these ideas into the clinic [45]. Further, by formally quantifying the contribution of the network connectivity to the control energy, we 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 [46].

Finally, we make strategic, task-agnostic edge deletions that maximally increase the determinant and observe that, even in an overdetermined, unsimplified system (N > M), a single edge deletion could produce 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 [47, 48]. Moreover, these results also suggest that minor, targeted structural changes through concussive injury can lead to drastic changes in overall brain function [49, 50], via altering the controllability landscape of the brain [24]. We further observed 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 natural to perform a similar analysis that takes into account the specific tasks v₁, v₂ by taking the derivative of the full energy term E_total with respect to A₂₁, which would optimize the network topology for a specific task, as studied in more detail in [25].

To achieve the most meaningful comparison between species, we only analyzed weighted meso-scale whole brain networks. As such, we did not include binary neuronal connectomes (e.g., C. elegans), and binary or partial connectomes (e.g., macaque). As more connectomes become available, we hope to further explore the role of species complexity on network controllability. Until then, we consider the comparison of energetically favorable connectivity between species to be a preliminary excursion into a nuanced evolutionary phenomena. As demonstrated in the significant percent change in energy after edge deletion, we emphasize that uncertainty in network connectivity has the potential to yield substantial changes in average control energy. Finally, we note that while methodological limitations prevent us from resolving excitatory versus inhibitory connectivity, all results are directly applicable to networks with signed elements. Further important theoretical considerations and methodological limitations pertinent to our approach, linear model of dynamics, optimality of control trajectories, and empirical data sets are discussed in the SI.

In closing, we note that the natural direction in which to take this work will be to use higher-order approximations of this framework found in the supplement to gain intuition for the role of complex network topologies (e.g. self-loops, cycles) in controlling networks. Moreover, it would be interesting to apply this reduced framework to random graphs and other well-known benchmarks – both from a mathematical perspective [51] and also in the context of neural systems [52, 53] – 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 [54], material science [55], and other contexts where optimal control of physical systems is of paramount importance.

Supplementary Material