The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure\DisableLigatures
[f]encoding = *, family = *
The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure
Shi Gu1,2, Matthew Cieslak3, Benjamin Baird4, Sarah F. Muldoon5, Scott T. Grafton3, Fabio Pasqualetti6, Danielle S. Bassett2,7,*
1 Applied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, PA, 19104 USA
2 Department of Bioengineering, University of Pennsylvania, Philadelphia, PA, 19104 USA
3 Department of Psychological and Brain Sciences, University of California, Santa Barbara, CA, 93106 USA
4 Center for Sleep and Consciousness, University of Wisconsin - Madison, Madison, WI 53706
5 Department of Mathematics and CDSE Program, University at Buffalo, SUNY, Buffalo, NY 14260
6 Department of Mechanical Engineering, University of California, Riverside, CA, 92521 USA
7 Department of Electrical & Systems Engineering, University of Pennsylvania, Philadelphia, PA, 19104 USA
* To whom correspondence should be addressed: firstname.lastname@example.org.
A critical mystery in neuroscience lies in determining how anatomical structure impacts the complex functional dynamics of human thought. How does large-scale brain circuitry constrain states of neuronal activity and transitions between those states? We address these questions using a maximum entropy model of brain dynamics informed by white matter tractography. We demonstrate that the most probable brain states – characterized by minimal energy – display common activation profiles across brain areas: local spatially-contiguous sets of brain regions reminiscent of cognitive systems are co-activated frequently. The predicted activation rate of these systems is highly correlated with the observed activation rate measured in a separate resting state fMRI data set, validating the utility of the maximum entropy model in describing neurophysiologial dynamics. This approach also offers a formal notion of the energy of activity within a system, and the energy of activity shared between systems. We observe that within- and between-system energies cleanly separate cognitive systems into distinct categories, optimized for differential contributions to integrated versus segregated function. These results support the notion that energetic and structural constraints circumscribe brain dynamics, offering novel insights into the roles that cognitive systems play in driving whole-brain activation patterns.
How does the complex interconnection structure between large-scale brain regions impact how we think? Does this structure inform how we move between different thoughts or actions? We address these questions using a simple mathematical model of brain dynamics that is informed by empirical estimates of anatomical interconnection structure between brain regions. Our results suggest that while interconnection structure plays an important role in constraining and predicting brain dynamics, an additional layer of explanatory power is offered by considering energetic constraints on those dynamics. These data offer a richer landscape of mechanisms that enhance our understanding of how we may move from one thought to the next.
A human’s adaptability to rapidly changing environments depends critically on the brain’s ability to carefully control the time within (and transitions among) different states. Here, we use the term state to refer to a pattern of activity across neurons or brain regions . The recent era of brain mapping has beautifully demonstrated that the pattern of activity across the brain or portions thereof  differs in different cognitive states . These variable patterns of activity have enabled the study of cognitive function via the manipulation of distinct task elements , the combination of task elements [4, 5], or the temporal interleaving of task elements [6, 7]. Such methods for studying cognitive function are built on the traditional view of mental chronectomy , which suggests that brain states are additive and therefore separable in both space and time (although see  for a discussion of potential caveats).
Philosophically, the supposed separability and additivity of brain states suggests the presence of strong constraints on the patterns of activations that can be elicited by the human’s environment. The two most common types of constraints studied in the literature are energetic constraints and structural constraints . Energetic constraints refer to fundamental limits on the evolution  or usage of neural systems , which inform the costs of establishing and maintaining functional connections between anatomically distributed neurons . Such constraints can be collectively studied within the broad theory of brain function posited by the free energy principal – a notion drawn from statistical physics and information theory – which states that the brain changes its state to minimize the free energy in neural activity [14, 15]. The posited preference for low energy states motivates an examination of the time within and transitions among local minimums of a predicted energy landscape of brain activity [16, 17].
While energetic costs likely form critical constraints on functional brain dynamics, an arguably equally important influence is the underlying structure and anatomy linking brain areas. Intuitively, quickly changing the activity of two brain regions that are not directly connected to one another by a structural pathway may be more challenging than changing the activity of two regions that are directly connected [18, 13]. Indeed, the role of structural connectivity in constraining and shaping brain dynamics has been the topic of intense neuroscientific inquiry in recent years [19, 20, 21, 22, 23]. Evidence suggests that the pattern of connections between brain regions directly informs not only the ease with which the brain may control state transitions , but also the ease with which one can externally elicit a state transition using non-invasive neurostimulation .
While energy and anatomy both form critical constraints on brain dynamics, they have largely been studied in isolation, hampering an understanding of their collective influence. Here, we propose a novel framework that combines energetic and structural constraints on brain state dynamics in a free energy model explicitly informed by structural connectivity. Using this framework, we map out the predicted energy landscape of brain states, identify local minima in the energy landscape, and study the profile of activation patterns present in these minima. Our approach offers fundamental insights into the distinct role that brain regions and larger cognitive systems play in distributing energy to enable cognitive function. Further, the results lay important groundwork for the study of energy landscapes in psychiatric disease and neurological disorders, where brain state transitions are known to be critically altered but mechanisms driving these alterations remain far from understood [26, 27].
Materials and Methods
Human DSI Data Acquisition and Preprocessing
Diffusion spectrum images (DSI) were acquired for a total of subjects (mean age 22.65.1 years, 24 female, 2 left handed) along with a weighted anatomical scan at each scanning session . Of these subjects, were scanned ones, was scanned twice, and were scanned three times, for a total of 61 scans.
DSI scans sampled directions using a half shell acquisition scheme with a maximum value of and an isotropic voxel size of mm. We utilized an axial acquisition with the following parameters: s, ms, slices, FoV (,, mm). All participants volunteered with informed written consent in accordance with the Institutional Review Board/Human Subjects Committee, University of California, Santa Barbara.
DSI data were reconstructed in DSI Studio (www.dsi-studio.labsolver.org) using -space diffeomorphic reconstruction (QSDR) . QSDR first reconstructs diffusion weighted images in native space and computes the quantitative anisotropy (QA) in each voxel. These QA values are used to warp the brain to a template QA volume in MNI space using the SPM nonlinear registration algorithm. Once in MNI space, spin density functions were again reconstructed with a mean diffusion distance of mm using three fiber orientations per voxel. Fiber tracking was performed in DSI Studio with an angular cutoff of , step size of mm, minimum length of mm, spin density function smoothing of , maximum length of mm and a QA threshold determined by DWI signal in the CSF. Deterministic fiber tracking using a modified FACT algorithm was performed until streamlines were reconstructed for each individual.
Structural Network Construction
Anatomical scans were segmented using FreeSurfer  and parcellated according to the Lausanne 2008 atlas included in the connectome mapping toolkit . A parcellation scheme including regions was registered to the B0 volume from each subject’s DSI data. The B0 to MNI voxel mapping produced via QSDR was used to map region labels from native space to MNI coordinates. To extend region labels through the gray/white matter interface, the atlas was dilated by mm. Dilation was accomplished by filling non-labeled voxels with the statistical mode of their neighbors’ labels. In the event of a tie, one of the modes was arbitrarily selected. Each streamline was labeled according to its terminal region pair.
From these data, we built structural brain networks from each of the diffusion spectrum imaging scans. Consistent with previous work [13, 32, 33, 34, 35, 24, 36, 25, 37], we defined these structural brain networks from the streamlines linking large-scale cortical and subcortical regions extracted from the Lausanne atlas . We summarize these estimates in a weighted adjacency matrix whose entries reflect the structural connectivity between region and region (Fig. 1A).
Following , here we use an edge weight definition based on the quantitative anisotropy (QA). QA is described by Yeh et. al (2010) as a measurement of the signal strength for a specific fiber population in an ODF [38, 39]. QA is given by the difference between and the isotropic component of the spin density function (SDF, ) scaled by the SDF’s scaling constant. Along-streamline QA was calculated based on the angles actually used when tracking each streamline. Although along-streamline QA is more specific to the anatomical structure being tracked, QA is more sensitive to MRI artifacts such as B1 inhomogeneity. QA is calculated for each streamline. We then averaged values over all streamlines connecting a pair of regions, and used this value to weight the edge between the regions.
Resting state fMRI data
As an interesting comparison to the computational model, we used resting state fMRI data collected from an independent cohort composed of 25 healthy right-handed adult subjects (15 female), with a mean age of 19.6 years (STD 2.06 year). All subjects gave informed consent in writing, in accordance with the Institutional Review Board of the University of California, Santa Barbara. Resting-state fMRI scans were collected on a 3.0-T Siemens Tim Trio scanner equipped with high performance gradients at the University of California, Santa Barbara Brain Imaging Center. A T2*-weighted echo-planar imaging (EPI) sequence was used (TR=2000 ms; TE=30 ms; flip angle=90\unichar176; acquisition matrix=6464; FOV=192 mm; acquisition voxel size = 333.5 mm; 37 interleaved slices; acquisition length=410s).
We preprocessed the resting state fMRI data using an in-house script adapted from the workflows described in detail elsewhere [40, 41]. The first four volumes of each sequence were dropped to control for instability effects of the scanner. Slice timing and motion correction were performed in AFNIusing 3dvolreg and FreeSurfer\unichar8217s BBRegister was used to co-register functional and anatomical spaces. Brain, CSF, and WM masks were extracted, the time series were masked with the brain mask, and grand-mean scaling was applied. The temporal derivative of the original 6 displacement and rotation motion parameters was obtained and the quadratic term was calculated for each of these 12 motion parameters, resulting in a total of 24 motion parameters which were regressed from the signal. The principal components of physiological noise were estimated using CompCor (aCompCor and tCompCor) and these components were additionally regressed from the signal. The global signal was not regressed. Finally, signals were low passed below 0.1 Hz and high passed above 0.01 Hz in AFNI. To extract regional brain signals from the voxel-level time series, a mask for each brain region in the Lausanne2008 atlas was obtained and FreeSurfer was used to individually map regions to the subject space. A winner-takes-all algorithm was used to combine mapped regions into a single mask. The resulting signal for each region was then extracted in FreeSurfer using mrisegstats.
Following data preprocessing and time series extraction, we next turned to extracting observed brain states. Importantly, physiological changes relating to neural computations take place on a time scale much smaller than the time scale of BOLD image acquisition. Thus, we treat each TR as representing a distinct brain state. To maximize consistency between the model-based and data-based approaches, we transformed the continuous BOLD magnitude values into a binary state vector by thresholding regional BOLD signals at 0. From the set of binary state vectors across all TRs, we defined activation rates in a manner identical to that described for the maximum entropy model data.
Defining an Energy Landscape
We begin by defining a brain state both intuitively and in mathematical terms. A brain state is a macroscopic pattern of BOLD activity across regions of the brain (Fig. 1A). For simplicity, here we study the case in which each brain region can be either active () or inactive (). Then, the binary vector respresents a brain state configuration.
Next, we wish to define the energy of a brain state. We build on prior work demonstrating the neurophysiological relevance of maximum entropy models in estimating the energy of brain states in rest and task conditions [42, 43]. For a given state , we write its energy in the second order expansion:
where represents an interaction matrix whose elements indicate the strength of the interaction between region and region . If , this edge decreases the energy of state , while if , this edge increases the energy of state . The column sum of the structural brain network, , is the strength of region . In the thermodynamic equilibrium of the system associated with the energy defined in Eqn , the entropy of the system is maximized and the probability of the configuration is .
The choice of the interaction matrix is an important one, particularly as it tunes the relative contribution of edges to the system energy. In this study, we seek to study structural interactions in light of an appropriate null model. We therefore define the interaction matrix to be equal to the modularity matrix  of the structural brain network:
for , where is the adjacency matrix, , and . This choice ensures that any element measures the difference between the strength of the edge in the brain and the expected strength of that edge in an appropriate null model (here given as the Newman-Girvan null model ). If the edge is stronger than expected, it will decrease the energy of the whole system when activated, while if the edge is weaker than expected, it will increase the energy of the whole system when activated.
Discovering Local Minima
The model described above provides an explicit correspondence between a brain’s state and the energy of that state, in essence formalizing a multidimensional landscape on which brain dynamics may occur. We now turn to identifying and characterizing the local minima of that energy landscape (Fig. 1C). We begin by defining a local minimum: a binary state is a local minimum if for all vectors satisfying , which means that the state realizes the lowest energy among its neighboring states within the closed unit sphere. We wish to collect all local minima in a matrix with
where is the number of local minima and is the number of nodes in the structural brain network (or equivalently the cardinality of the adjacency matrix ).
Now that we have defined a local minimum of the energy landscape, we wish to discover these local minima given the pattern of white matter connections represented in structural brain networks. To discover local minima of , we first note that the total number of states is , which – when – prohibits an exhaustive analysis of all possibilities. Moreover, the problem of finding the ground state is an NP-complete problem , and thus it is unrealistic to expect to identify all local minima of a structural brain network. We therefore choose to employ a clever heuristic to identify local minima. Specifically, we combine the Metropolis Sampling method  and a steep search algorithm using gradient descent methods. We identify a starting position by choosing a state uniformly at random from the set of all possible states. Then, we step through the energy landscape via a random walk driven by the Metropolis Sampling criteria (see Algorithm Discovering Local Minima). At each point on the walk, we use the steep search algorithm to identify the closest local minimum.
[H] Let be the vector obtained by changing the value of the -th entry of \Fort = 1 : N Randomly select an index Set with probability and otherwise \While is not a local minimum Set Set Set
Here, are the sampled states, are the sampled local minima, and is the temperature parameter which can be absorbed in . In the context of any sampling procedure, it is important to determine the number of samples necessary to adequately cover the space. Theoretically, we wish to identify a number of samples following which the distribution of energies of the local minima remains stable. Practically speaking, we choose 4 million samples in this study, and demonstrate the stability of the energy distribution in the Supplement. A second important consideration is to determine the initial state, that is the state from which the random walk begins. Here we choose this state uniformly at random from the set of all possible states. However, this dependence on a uniform probability distribution may not be consistent with the actual probability of states in the energy landscape. We therefore must ensure that our results are not dependent on our choice of initial state. To ensure independence from the initial state, we dismiss the first local minima identified, and we demonstrate in the Supplement that this procedure ensures our results are not dependent on the choice of initial state.
Characterizing Local Minima
Following collection of local minima, we wished to characterize their nature as well as their relationships to one another. First, we estimated the radius of each local minimum as the Hamming distance from the minimum to the closest sampled point on energy landscape. Next, calculated the Hamming distance from each local minimum to the first sampled local minimum, a second quantification of the diversity of the energy lanscape that we traverse in our sampling. Finally, we quantify how diverse the observed local minima are by calculating the pairwise normalized mutual information of each pair of local minima.
Next, we wished to understand the role of different regions and cognitive systems in the minimal energy states. Cognitive systems are sets of brain regions that show coordinated activity profiles in resting state or task-based BOLD fMRI data . They include the visual, somatosensory, motor, auditory, default mode, salience, fronto-parietal, cingulo-opercular, dorsal and ventral attention systems, as well as subcortical areas. Here, the specific association of regions of interest to cognitive systems are exactly as listed in  and based originally on results described in . We characterize the roles of these systems in the local minima by assessing their activation rates, as well as the utilization energies required for communication within and between systems.
Intuitively, we define the activation rate of a node as the average activation state of that node over all the local minima. Formally, the activation rate for region is defined as
where indexes over states, and recall is the number of local minima. The computed activation rate offers a prediction of which regions are more versus less active across the local minima (that is, the brain’s locally “stable ” states), and can be directly compared with the resting state activation rate estimated from empirical fMRI data.
To complement the information provided by the activation rates, we also defined the energetic costs associated with utilizing within- and between-systems interactions. We note that each cognitive system is a subnetwork of the whole brain network. We use the index set to indicate the set of nodes associated with the cognitive system, and thus gives the number of nodes in the system. Then, for a given state , the within-system energy measures the cost associated with the set of interactions constituting the subnetwork. The between-system energy measures the cost associated with the set of interactions between the subnetwork and all other nodes in the whole network. Formally, we define
where measures the within-system energy, measures the between-system energy, and the normalization coefficients , are chosen by considering the number of the corresponding interactions.
Permutation Tests for State Association
For a given local minimum configuration , we associate it with system ,
where is the configuration pattern of system such that the corresponding regions for system are activated and the others not and where “NMI” refers to the Normalized Mutual Information, which is used to measure the similarities between the two states. To obtain the null distribution, for each local minimum configuration in the collection , we permute the configuration at each position of to achieve a random configuration with the same activation rate, and we then compute the associated percentage in each system. Then we repeat this procedure to generate samples and construct the null distribution of the probability of being configured as each system pattern. Considering the large size of the state collection, the variance of the samples in the null distribution will be small. We pick here. See Fig. 4 D for the results.
Local Minima in the Brain’s Energy Landscape
By sampling the energy landscape of each structural connectivity matrix, we identified an average of approximately 450 local minima or low energy brain states: binary vectors indicating the pattern of active and inactive brain regions (see Methods). On average across the 61 scans, 144 brain regions were active in a given local minimum, representing 61.70% of the total (standard deviation: 6.04%). This large percentage suggests that the brain may utilize a broad set of rich and diverse activations to perform cognitive functions , rather than activation patterns localized to small geographic areas.
To quantify this diversity, we examined the location of minima on the energy landscape, the size of the basins surrounding the minima, and the mutual information between minima. First, we estimated the distance from the first local minima identified to all subsequent minima (see Methods; Fig. 2A). We observe an order of magnitude change in the distance between the first and second local minima, and the first and last local minima, suggesting that local minima span a broad geographic domain in the energy landscape. Interestingly, these minima differ not only in their location on the energy landscape, but also in the size of the basins surrounding them. We estimate basin size by calculating the radius of each local minimum (see Methods) and show that the distribution of radii follows a power-law, with the majority of minima displaying a small radius, and only a few minima displaying a large radius (Fig. 2B). More specifically, we fit the function – where is a constant – to the data using a statistically principled approach [51, 52]. We identified an for , and the -value for the goodness of fit was indicating that the power law was an accurate fit to the data. As a final quantification of minima diversity, we estimated the normalized mutual information between every pair of local minima, as an intuitive measurement of the similarity between anatomical compositions of the minima. We observe that the probability distribution of normalized mutual information between minima pairs is heavy tailed, indicating that most minima pairs display very dissimilar anatomical compositions, and only a few minima pairs display similar anatomical compositions (Fig. 2C).
From a neurophysiological perspective, it is also important to note that these local minima displayed significant local structure. Specifically, we found that regions within known cognitive systems tended to be active together. The probability that regions were co-active is 48.22%, which was significantly greater than that expected in a null distribution (associated was ; see Methods). These results indicate that the structural connectivity between brain regions, and the assumption of energy conservation, predict that regions that belong to the same cognitive system will tend to be co active with one another during diverse cognitive functions. Indeed, these predictions are consistent with previous studies of functional neuroimaging data demonstrating that groups co-active regions tended to align well with known cognitive systems [53, 54].
Activation Rates of Cognitive Systems
Given the alignment of activation patterns with cognitive systems, we next asked whether certain cognitive systems were activated more frequently than others. To address this question, we studied the activation rate of each cognitive system, which measures how frequently the regions in the cognitive system participated in the set of states identified as local minima. Intuitively, if the activation rate is high, the system is more likely to be active in diverse brain states. We observed that systems indeed showed signficantly different activation rates (Fig. 2D). Sensorimotor systems (auditory, visual, somatosensory) tended to display the lowest activation rates, followed by higher order cognitive systems (salience, attention, fronto-parietal, and cingulo-opercular), and subcortical structures. The system with the largest activate rate was the default mode system, suggesting that activation of this system is particularly explicable from structural connectivity and the assumption of energy conservation. The unique role of the default mode system is consistent with predictions from network control theory that highlight the optimal placement of default mode regions within the network to maximize potential to move the brain into many easily reachable states with minimal energetic costs .
It is important to determine whether this activation rate is driven by simple properties of the structural connectivity matrix that do not depend on assumptions of energy conservation. To address this question, we next assessed the relationship between a simple summary statistic of the structural connectivity matrix – the strength, or weighted degree, of a brain region – and the predicted activation rate drawn from the maximum entropy model. We observed that the activation rate was not well predicted by the weighted degree on average over brain regions (see Supplement). These data suggest that the additional assumption of energy conservation produces a set of brain states that cannot be predicted from simple statistics of structural connectivity alone.
Relating Predicted Activation Rates to Rates Observed in Functional Neuroimaging Data
Before exercising the model further to explore how energy is utilized in the brain, we wished to quantify the relationships between the theoretically predicted activation rates, and activation rates observed empirically in functional neuroimaging data. Specifically, we studied fMRI data acquired in a separate set of healthy adult human subjects.
Next, we wished to directly quantify similarities between the predicted activation rates and those observed empirically in resting state fMRI. In the resting state data, we observed that highly active regions were located in broad swaths of frontal and parietal cortex, as well as medial prefrontal, precuneus, and cingulate (Fig. 3A). This pattern of high activation is consistent with the so-called “default-mode” of resting state brain function . In our maximum entropy model, we observed that the areas predicted to have high activation rates show a broad similarity to those observed empirically in the resting state (Fig. 3B). Indeed, we observed that the empirical resting activation rate of brain regions is significantly correlated with the activation rate predicted from the maximum entropy model (Fig. 3C; Pearson correlation coefficient , ). These results suggest that the modeling framework we use here has significant similarities to observable features of resting state brain dynamics. However, it is important to mention that there are also noticeable differences between the two maps: the predicted activation rates are strong along the medial wall, while the resting state activation rates extend to larger sections of lateral cortices.
Utilization Energies of Cognitive Systems
We next turned to exercising our model to further understand the potential constraints on brain state dynamics. Specifically, we asked how cognitive systems utilized the minimal energy presumably available to them. Intuitively, this question encompasses both how energy is utilized by within-system interactions, and how energy is utilized by between-system interactions. We therefore defined the within-system energy, which measures the cost associated with the set of interactions constituting the cognitive system, and the between-system energy, which measures the cost associated with the set of interactions between cognitive systems. We observed a fairly strong dissociation between these two variables: cognitive systems that display a large within-system energy are not necessarily those that display a large between-system energy (see Fig. 4A and B). Indeed, within- and between-system energies are not significantly correlated across systems (Pearson’s correlation coefficient and ), suggesting that these two variables offer markers of distinct constraints.
Moreover, we observed that the 2-dimensional plane mapped out by the within- and between-system energies of all brain regions revealed the presence of 4 surprisingly distinct clusters (Fig. 4C) that are not explicable by simple statistics such as network degree (see Supplement). Each cluster represents a unique strategy in energy utilization that is directly reflected in its activation pattern; in other words, each cluster offers a distinct balance between the energetic costs of within-system interactions and the energetic costs of between-system interactions. The central cluster, displaying high within-system energies but low between-system energies, is composed of subcortical and fronto-parietal systems. A high between-system energy cone emanating from this central cluster is composed of predominantly primary and secondary sensorimotor cortices in somatosensory, visual, and auditory systems. A second cone emanating from the central cluster with a slightly lower between-system energy is composed predominantly of regions in the default mode system. The final cone emanating from the central cluster with between-system energies near zero is composed predominantly of regions in the dorsal and ventral attention systems. These results suggest that sensorimotor, default mode, attention, and cognitive control circuits display differential preferences for energy utilization: regions in attentional systems share less energy with other networks than regions in sensorimotor systems, while the default mode maintains an intermediate balance.
The clear differences in the energies utilized by different cognitive systems and by between-system interactions begs the question of whether the brain cares about these energies. Does the brain prefer smaller within-system energies, smaller between-system energies, or some balance between the two? To address this question, we studied the ensemble of local minima, and asked which systems were commonly expressed. Specifically, for each local minimum, we determined which system was most activated, assigned the minimum to that system, and performed this assignment for all minima. We observed that 3 systems were represented at higher percentages than expected in a permutation-based null model (see Methods): the default mode system, the visual system, and the somatosensory system (Fig. 4D). Importantly, these three systems represent the systems with the highest between system energies (Fig. 4C), but are indistinguishable from other cognitive systems in terms of within-system energy. These results suggest that the brain may prefer high integration between systems over low integration, and that the constraint of between-system energies is more fundamental to brain function than the constraint of within-system energies.
In this paper, we address the question of how large-scale brain circuitry and distinct energetic constraints produce whole-brain patterns of activity. We build our approach on a maximum entropy model of brain dynamics that is explicitly informed by estimates of white matter microstructure derived from deterministic tractography algorithms. The model allows us to study minimal energy states, which we observe to be composed of co-activity in local spatially-contiguous sets of brain regions reminiscent of cognitive systems. These systems are differentially active, and activity patterns are significantly correlated with the observed activation rate measured in a separate resting state fMRI data set. Finally, we exercise this model to ask how cognitive systems utilize the minimal energy presumably available to them. We find that the energy utilized within and between cognitive systems distinguishes 4 classes of energy utilization dynamics, corresponding to sensorimotor, default mode, attention, and cognitive control functions. These results suggest that diverse cognitive systems are optimized for differential contributions to integrated versus segregated function via distinct patterns of energy utilization. More generally, the results highlight the importance of considering energetic constraints in linking structural connectivity to observed dynamics of neural activity.
The Role of Activation vs. Connectivity in Understanding Brain Dynamics
As the interest in understanding structural brain connectomics has blossomed in the last several years [56, 57], it has not been accompanied by an equally vivid interest in linking subsequent insights to the more traditional notions of brain activation profiles . Indeed, the fields of systems and cognitive neuroscience have instead experienced a pervasive divide between the relatively newer notions of connectome mapping and the relatively traditional yet highly effective notions of brain mapping , which have led to powerful insights into neural function in the last quarter century . This divide is at least in part due to the fact that graph theory and network-based methods on which the field of connectomics is based have few tools available to link node properties (activity) with edge properties (connectivity) . While technically explicable, however, the conceptual divide between these fields can only lead to their detriment, and synergistic efforts are necessary to develop a language in which both activity and connectivity can be examined in concert . This study offers one explicit mathematical modeling framework in which to study the relationships between activation profiles across the whole brain and underlying structural connectivity linking brain regions. Complementary approaches include the model-based techniques formalized in the publicly available resource The Virtual Brain [62, 63].
In this study, we observed that brain regions affiliated with known cognitive systems – including somatosensory, visual, auditory, default mode, dorsal and ventral attention, fronto-parietal, and cingulo-opercular – also tend to be active together with one another in low energy brain states. Indeed, these theoretical results are consistent with previous studies of functional neuroimaging data demonstrating that groups of co-active regions tended to align well with known cognitive systems [53, 54]. This correspondence is particularly interesting when one considers how these cognitive systems were initially defined: and that is, based on strong and dense functional connectivity . Thus, our results point to a fundamental mapping between activity and connectivity: regions that are active together tend to be functionally connected together. The presence of such a map has been empirically observed in the resting state (though not in task ), in both healthy controls and in people with schizophrenia where the map appears to be fundamentally altered in its nature [64, 65, 66]. Here we offer a mechanism for this mapping based on a combined consideration of energy- and connectivity-based constraints.
Critical Importance of Energy Constraints
The quest to understand and predict brain dynamics from the architecture of underlying structural connectivity is certainly not a new one. In fact, there have been concerted efforts over the last decade and more to identify structural predictors of the resting state BOLD signal. Seminal contributions have included the observations of statistically significant correlations between structural connectivity estimated from diffusion imaging data and functional connectivity estimated from fMRI , as well as extensions of these correlations that account for long distance paths along white matter tracts  and spectral properties of structural matrices . The question of how brain structure constrains a wide range of brain states (beyond simply the resting state) is a very open area of inquiry. Moreover, this question is particularly challenging to address with empirical data because it is difficult to obtain data from humans in more than a handful of task states . For this reason, computational models play a very important role in offering testbeds for the development of theories linking structure to ensembles of brain states, which can in turn offer testable predictions. Our results suggest that an understanding of the relationship between brain structure and function is perhaps ill-constrained when examining connectivity alone. The additional assumption of energy conservation produces a set of brain states that cannot be simply predicted from statistics of structural connectivity, perhaps offering a mechanism for the large amount of unexplained variance in prior predictions [67, 22].
Our results are built on the formalism of the maximum entropy model, which is predicated on pairwise statistics . However, emerging evidence suggests that some neurophysiological phenomenon are better studied in the framework of simplicial complexes rather than dyads . For example, integrate and fire neurons exposed to common fluctuating input display strong beyond-pairwise correlations that cannot be captured by maximum entropy models . Similar arguments can also be made for co-activation patterns in BOLD fMRI [53, 54, 73]. It will be interesting in future to determine the role energetic and structural constraints on observed higher-order functional interactions during human cognitive function.
A second important consideration is that the maximum entropy model is appropriate for systems at equilibrium. Therefore, the local minima identified may not accurately represent the full class of states expected to be elicited by daily activity. Instead, the local minima identified here are expected to more accurately represent the set of states expected to appear as a person rests in the so-called default mode of brain function, which is thought to lie near a stable equilibrium . Such an interpretation is consistent with our findings that the activation rates predicted by the maximum entropy model are strongly correlated with the activation rates observed in resting state fMRI data.
The analyses presented in this study produce information regarding an underlying energy landscape through which the brain is predicted to move. The existence of such a landscape motivates the very interesting question of how the brain transitions between states. In sampling this landscape, we have used a simple random walk in an effort to extract a large ensemble of possible brain states, measured by local minima. However, the question of which walks a healthy (or diseased) brain might take through this landscape remains open. Such walks or dynamical trajectories may be determined by energetic inputs to certain regions of the brain , either by external stimuli or by neuromodulation . In this context, network control theory may offer explicit predictions regarding the optimal dynamic trajectories that the brain may take through a set of states to move from an initial state to a target state with little energetic resources [24, 74]. In addition to inputs to single regions, changes in a cognitive task – for example elicited by task-switching paradigms – may also drive a specific trajectory of brain states. Indeed, it is intuitively plausible that the asymmetric switch costs observed between cognitively effortful and less effortful tasks [75, 76] may be explained by characteristics of the energy landscape defined by structural connections between task-activated brain regions.
Stability of the Energy Distribution with respect to the Number of Samples. We plot the probability distribution of the energy for the first 2000, 4000, 6000 and 8000 samples. We observe that the shapes of the probability distributions are qualitatively consistent. We confirm this qualitative observation with Kolmogorov-Smirnov tests (see Supplemental text).
Activation Rate is Poorly Predicted by Regional Degree and Energy. (A) Scatterplot of weighted regional degree versus activation rate. (B) Scatterplot of regional energy versus activation rate. We observe that the activation rate is not well predicted by either regional energy or weighted degree.
Relationship Between Utilization Energies and Degree. (A) Scatterplot of within-system connectivity versus between-system connectivity, for individual brain regions that are color-coded by cognitive system. (B) The same data presented in panel (A) except only for the default mode, attention, and task control systems, demonstrating the indistinguishability of default mode and attention systems. (C) Scatterplot of the between- and within- system energy. (D) The same data presented in panel (C) except only for the detault mode, attention, and task control systems. We observe that cognitive systems are more clearly separated in the 2-dimensional space of the within- and between-system energies than in the 2-dimensional space of the within- and between-system connectivity. Across all four panels, data points indicate brain regions, and color of data points indicates which cognitive system that region is affiliated with (see legend for map from color to cognitive system).
Simulated Activation Rate is Significantly Correlated with the Rates Observed in Resting State Functional Neuroimaging Data. We observe that the activation rates estimated from resting state fMRI data are significantly positively correlated with the activation rates estimated from the local minima of the maximum entropy model (Pearson’s correlation coefficient , ). Each data point represents a brain region, with either observed or predicted activation rates averaged over subjects.
S1 Appendix. Supplementary Information for “The Energy Landscape of Neurophysiological Activity Implicit in Brain Network Structure”.
We thank Ankit Khambhati, Lia Papadopoulus, and Evelyn Tang for helpful comments on an earlier fersion of this manuscript. D.S.B., S.G., and R.F.B. acknowledge support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the Army Research Laboratory and the Army Research Office through contract numbers W911NF-10-2-0022 and W911NF-14-1-0679,the National Institute of Mental Health (2-R01-DC-009209-11), the National Institute of Child Health and Human Development (1R01HD086888-01), the Office of Naval Research, and the National Science Foundation (BCS-1441502, BCS-1430087, and PHY-1554488). FP acknowledges support from BCS-1430280. B.B. was supported by the National Institutes of Health under Ruth L. Kirschstein National Research Service Award F32NS089348 from the NINDS.
- Tang YY, Rothbart MK, Posner MI. Neural correlates of establishing, maintaining, and switching brain states. Trends in cognitive sciences. 2012;16(6):330–337.
- Mahmoudi A, Takerkart S, Regragui F, Boussaoud D, Brovelli A. Multivoxel pattern analysis for fMRI data: a review. Comput Math Methods Med. 2012;2012:961257.
- Gazzaniga MS, editor. The cognitive neurosciences. MIT Press; 2013.
- Szameitat AJ, Schubert T, Muller HJ. How to test for dual-task-specific effects in brain imaging studies: an evaluation of potential analysis methods. Neuroimage. 2011;54(3):1765–1773.
- Alavash M, Hilgetag CC, Thiel CM, Giessing C. Persistency and flexibility of complex brain networks underlie dual-task interference. Hum Brain Mapp. 2015;36(9):3542–3562.
- Ruge H, Jamadar S, Zimmermann U, Karayanidis F. The many faces of preparatory control in task switching: reviewing a decade of fMRI research. Hum Brain Mapp. 2013;34(1):12–35.
- Muhle-Karbe PS, De Baene W, Brass M. Do tasks matter in task switching? Dissociating domain-general from context-specific brain activity. Neuroimage. 2014;99:332–341.
- Donders FC. On the speed of mental processes. Acta Psychol. 1969;30:412–431.
- Mattar MG, Cole MW, Thompson-Schill SL, Bassett DS. PLoS Comput Biol. 2015;11(12):e1004533.
- Bullmore E, Sporns O. The economy of brain network organization. Nature Reviews Neuroscience. 2012;13(5):336–349.
- Niven JE, Laughlin SB. Energy limitation as a selective pressure on the evolution of sensory systems. J Exp Biol. 2008;211:1792–1804.
- Attwell D, Laughlin SB. An energy budget for signalling in the grey matter of the brain. J Cereb Blood Flow and Metab. 2001;21:1133–1145.
- Bassett DS, Greenfield D, Meyer-Lindenberg A, Weinberger DR, Moore SW, Bullmore ET. Efficient physical embedding of topologically complex information processing networks in brains and computer circuits. PLoS Comput Biol. 2010;6(4):e1000748.
- Friston K, Kilner J, Harrison L. A free energy principle for the brain. Journal of Physiology-Paris. 2006;100(1):70–87.
- Friston K. The free-energy principle: a unified brain theory? Nature Reviews Neuroscience. 2010;11(2):127–138.
- Moreno-Bote R, Rinzel J, Rubin N. Noise-induced alternations in an attractor network model of perceptual bistability. Journal of neurophysiology. 2007;98(3):1125–1139.
- Tsodyks M, Pawelzik K, Markram H. Neural networks with dynamic synapses. Neural computation. 1998;10(4):821–835.
- Achard S, Bullmore E. Efficiency and cost of economical brain functional networks. PLoS computational biology. 2007;3(2):e17.
- Honey C, Sporns O, Cammoun L, Gigandet X, Thiran JP, Meuli R, et al. Predicting human resting-state functional connectivity from structural connectivity. Proceedings of the National Academy of Sciences. 2009;106(6):2035–2040.
- Honey CJ, Thivierge JP, Sporns O. Can structure predict function in the human brain? Neuroimage. 2010;52(3):766–776.
- Deco G, Jirsa VK. Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J Neurosci. 2012;32(10):3366–3375.
- Goni J, van den Heuvel MP, Avena-Koenigsberger A, Velez de Mendizabal N, Betzel RF, Griffa A, et al. Resting-brain functional connectivity predicted by analytic measures of network communication. Proc Natl Acad Sci U S A. 2014;111(2):833–838.
- Becker CO, Pequito S, Pappas GJ, Miller MB, Grafton ST, Bassett DS, et al. Accurately Predicting Functional Connectivity from Diffusion Imaging. arXiv;1512:02602.
- Gu S, Pasqualetti F, Cieslak M, Telesford QK, Yu AB, Kahn AE, et al. Controllability of structural brain networks. Nat Commun. 2015;6:8414.
- Muldoon SF, Pasqualetti F, Gu S, Cieslak M, Grafton ST, Vettel JM, et al. Stimulation-based control of dynamic brain networks. arXiv. 2016;In Prep.
- Ravizza SM, Moua KC, Long D, Carter CS. The impact of context processing deficits on task-switching performance in schizophrenia. Schizophr Res. 2010;116(2–3):274–279.
- Wylie GR, Clark EA, Butler PD, Javitt DC. Schizophrenia patients show task switching deficits consistent with N-methyl-d-aspartate system dysfunction but not global executive deficits: implications for pathophysiology of executive dysfunction in schizophrenia. Schizophr Bull. 2010;36(3):585–594.
- Cieslak M, Grafton ST. Local termination pattern analysis: a tool for comparing white matter morphology. Brain Imaging Behav. 2014;8(2):292–299.
- Yeh FC, Tseng W. NTU-90: a high angular resolution brain atlas constructed by -space diffeomorphic reconstruction. Neuroimage. 2011;58(1):91–99.
- Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis. I. Segmentation and surface reconstruction. NeuroImage. 1999;9(2):179–194.
- Hagmann P, Cammoun L, Gigandet X, Meuli R, Honey CJ, Wedeen VJ, et al. Mapping the structural core of human cerebral cortex. PLoS Biology. 2008;6(7):e159.
- Bassett DS, Brown JA, Deshpande V, Carlson JM, Grafton ST. Conserved and variable architecture of human white matter connectivity. Neuroimage. 2011;54(2):1262–1279.
- Hermundstad AM, Bassett DS, Brown KS, Aminoff EM, Clewett D, Freeman S, et al. Structural foundations of resting-state and task-based functional connectivity in the human brain. Proc Natl Acad Sci U S A. 2013;110(15):6169–6174.
- Hermundstad AM, Brown KS, Bassett DS, Aminoff EM, Frithsen A, Johnson A, et al. Structurally-constrained relationships between cognitive states in the human brain. PLoS Comput Biol. 2014;10(5):e1003591.
- Biol PC. Resolving structural variability in network models and the brain. 10;(3):e1003491.
- Small-World Propensity in Weighted, Real-World Networks. arXiv. 2015;1505:02194.
- Sizemore A, Giusti C, Bassett DS. Classification of weighted networks through mesoscale homological features. arXiv. 2015;1512:06457.
- Yeh FC, Wedeen VJ, Tseng WY. Generalized Q-Sampling Imaging. Medical Imaging, IEEE Transactions on. 2010;29(9):1626–1635.
- Tuch DS. Q-ball imaging. Magnetic Resonance in Medicine. 2004;52(6):1358–1372.
- Baird B, Smallwood J, Gorgolewski KJ, Margulies DS. Medial and lateral networks in anterior prefrontal cortex support metacognitive ability for memory and perception. The Journal of Neuroscience. 2013;33(42):16657–16665.
- Satterthwaite TD, Elliott MA, Gerraty RT, Ruparel K, Loughead J, Calkins ME, et al. An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. Neuroimage. 2013;64:240–256.
- Watanabe T, Hirose S, Wada H, Imai Y, Machida T, Shirouzu I, et al. A pairwise maximum entropy model accurately describes resting-state human brain networks. Nat Commun. 2013;4:1370.
- Watanabe T, Hirose S, Wada H, Imai Y, Machida T, Shirouzu I, et al. Energy landscapes of resting-state brain networks. Front Neuroinform. 2014;8:12.
- Newman MEJ. Modularity and community structure in networks. Proceedings of the National Academy of Sciences of the United States of America. 2006;103(23):8577–8696.
- Clauset A, Newman ME, Moore C. Finding community structure in very large networks. Physical review E. 2004;70(6):066111.
- Cipra BA. The Ising model is NP-complete. SIAM News. 2000;33(6):1–3.
- Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. The journal of chemical physics. 1953;21(6):1087–1092.
- Manning CD, Raghavan P, Schütze H, et al. Introduction to information retrieval. vol. 1. Cambridge university press Cambridge; 2008.
- Power JD, Cohen AL, Nelson SM, Wig GS, Barnes KA, Church JA, et al. Functional network organization of the human brain. Neuron. 2011;72(4):665–678.
- Evolution of brain network dynamics in neurodevelopment. arXiv. 2016;.
- Clauset A, Shalizi CR, Newman MEJ. Power-law distributions in empirical data. SIAM Review. 2009;51(4):661–703.
- Y V, Clauset A. Power-law distributions in binned empirical data. Annals of Applied Statistics. 2014;8(1):89 –119.
- Crossley NA, Mechelli A, Vertes PE, Winton-Brown TT, Patel AX, Ginestet CE, et al. Cognitive relevance of the community structure of the human brain functional coactivation network. Proc Natl Acad Sci U S A. 2013;110(28):11583–11588.
- Crossley NA, Mechelli A, Scott J, Carletti F, Fox PT, McGuire P, et al. The hubs of the human connectome are generally implicated in the anatomy of brain disorders. Brain. 2014;137(Pt 8):2382–2395.
- Raichle ME, MacLeod AM, Snyder AZ, Powers WJ, Gusnard DA, Shulman GL. A default mode of brain function. Proceedings of the National Academy of Sciences. 2001;98(2):676–682.
- Sporns O, Tononi G, Kötter R. The human connectome: A structural description of the human brain. PLoS Comput Biol. 2005;1(4):e42.
- Sporns O. The human connectome: a complex network. Ann N Y Acad Sci. 2011;1224:109–125.
- Bassett DS, Yang M, Wymbs NF, Grafton ST. Learning-induced autonomy of sensorimotor systems. Nat Neurosci. 2015;18(5):744–751.
- Sporns O. Cerebral cartography and connectomics. Philos Trans R Soc Lond B Biol Sci. 2015;370(1668):20140173.
- Zeki S. Introduction: cerebral cartography 1905–2005. Phil Trans R Soc B. 2005;360:651\unichar8211–652.
- Newman MEJ, Clauset A. Structure and inference in annotated networks. arXiv. 2015;1507:04001.
- Schirner M, Rothmeier S, Jirsa VK, McIntosh AR, Ritter P. An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data. Neuroimage. 2015;117:343–357.
- Ritter P, Schirner M, McIntosh AR, Jirsa VK. The virtual brain integrates computational modeling and multimodal neuroimaging. Brain Connect. 2013;3(2):121–145.
- Bassett DS, Nelson BG, Mueller BA, Camchong J, Lim KO. Altered resting state complexity in schizophrenia. Neuroimage. 2012;59(3):2196–2207.
- Zalesky A, Fornito A, Egan GF, Pantelis C, Bullmore ET. The relationship between regional and inter-regional functional connectivity deficits in schizophrenia. Hum Brain Mapp. 2012;33(11):2535–2549.
- Yu Q, Sui J, Liu J, Plis SM, Kiehl KA, Pearlson G, et al. Disrupted correlation between low frequency power and connectivity strength of resting state brain networks in schizophrenia. Schizophr Res. 2013;143(1):165–171.
- Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, Meuli R, et al. Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A. 2009;106(6):2035–2040.
- Becker C, Pequito S, Pappas GJ, Miller MB, Grafton ST, Preciado V. Accurately Predicting Functional Connectivity from Diffusion Imaging. arXiv. 2016;1512:02602.
- Cole MW, Bassett DS, Power JD, Braver TS, Petersen SE. Intrinsic and task-evoked network architectures of the human brain. Neuron. 2014;83(1):238–251.
- Bialek W, Cavagna A, Giardina I, Mora T, Silvestri E, Viale M, et al. Statistical mechanics for natural flocks of birds. Proc Natl Acad Sci U S A. 2012;109(13):4786–4791.
- Giusti C, Ghrist R, Bassett DS. Two’s company, three (or more) is a simplex: Algebraic-topological tools for understanding higher-order structure in neural data. Trends in Cognitive Sciences. 2016;Under Consideration.
- Leen DA, Shea-Brown E. A Simple Mechanism for Beyond-Pairwise Correlations in Integrate-and-Fire Neurons. J Math Neurosci. 2015;5(1):30.
- Fox PT, Lancaster JL, Laird AR, Eickhoff SB. Meta-analysis in human neuroimaging: computational modeling of large-scale databases. Annu Rev Neurosci. 2014;37:409–434.
- Pasqualetti F, Zampieri S, Bullo F. Controllability Metrics, Limitations and Algorithms for Complex Networks. IEEE Transactions on Control of Network Systems. 2014;1(1).
- Wu S, Hitchman G, Tan J, Zhao Y, Tang D, Wang L, et al. The neural dynamic mechanisms of asymmetric switch costs in a combined Stroop-task-switching paradigm. Sci Rep. 2015;5:10240.
- Davidson MC, Amso D, Anderson LC, Diamond A. Development of cognitive control and executive functions from 4 to 13 years: evidence from manipulations of memory, inhibition, and task switching. Neuropsychologia. 2006;44(11):2037–2078.