Control energy assessment of spatial interactions among macro-scale brain networks.

Many recent studies have revealed that spatial interactions of functional brain networks derived from fMRI data can well model functional connectomes of the human brain. However, it has been rarely explored what the energy consumption characteristics are for such spatial interactions of macro-scale functional networks, which remains crucial for the understanding of brain organization, behavior, and dynamics. To explore this unanswered question, this article presents a novel framework for quantitative assessment of energy consumptions of macro-scale functional brain network's spatial interactions via two main effective computational methodologies. First, we designed a novel scheme combining dictionary learning and hierarchical clustering to derive macro-scale consistent brain network templates that can be used to define a common reference space for brain network interactions and energy assessments. Second, the control energy consumption for driving the brain networks during their spatial interactions is computed from the viewpoint of the linear network control theory. Especially, the energetically favorable brain networks were identified and their energy characteristics were comprehensively analyzed. Experimental results on the Human Connectome Project (HCP) task-based fMRI (tfMRI) data showed that the proposed methods can reveal meaningful, diverse energy consumption patterns of macro-scale network interactions. In particular, those networks present remarkable differences in energy consumption. The energetically least favorable brain networks are stable and consistent across HCP tasks such as motor, language, social, and working memory tasks. In general, our framework provides a new perspective to characterize human brain functional connectomes by quantitative assessment for the energy consumption of spatial interactions of macro-scale brain networks.

INTRODUCTION

Recently, an increasing number of studies from neuroscience research have revealed that functional brain networks intrinsically interact on multiple spatiotemporal scales (Baajour et al., 2020; Balajoo, Asemani, Khadem, & Soltanian‐Zadeh, 2020; Jiang et al., 2018; Li et al., 2014; Yuan et al., 2018; Zhu et al., 2016). Large‐scale neural interactions allow the brain to dynamically segregate and intergrade to maintain consciousness or respond to external stimuli, which might require substantial amount of energy consumption during activity shifts (Watts, Pocock, & Claudianos, 2018). In fact, human brain is an organ with high energy requirements (Aiello & Wheeler, 1995). Moreover, energetic demands are not uniformly distributed throughout the whole brain but instead increase in localized regions dependent on neural activity. Thus, characterization of macro‐scale functional brain network interactions and assessment of energy consumption of these interactions are of great importance to better understand the brain's function.

After decades of active research using the fMRI technique, functional brain networks and their interactions have been well described with rest and task data (Baajour et al., 2020; Elbich, Molenaar, & Scherf, 2019; Yuan et al., 2018; Zhang et al., 2019, 2020). Thus, it is more likely to achieve a better understanding of the relationship between energy consumption and functional network interactions. In the literature, a variety of approaches have been proposed to analyze the energy characteristics of functional networks, such as computation of the maximum entropy model (Ashourvan, Gu, Mattar, Vettela, & Bassetta, 2017; Ezaki, Watanabe, Ohzeki, & Masuda, 2017; Gu et al., 2018; Kang, Pae, & Park, 2017; Watanabe et al., 2014) and detection of the cerebral energy metabolism (Liang, Zou, He, & Yang, 2013; Lord, Expert, Huckins, & Turkheimer, 2013; Noack, Manjesh, Ruszinko, Siegelmann, & Kozma, 2017). Specifically, the maximum entropy model defines an energy function for the system with a maximum entropy probability distribution. Assuming that the resting‐state network activity can be captured in terms of transitions among locally stable states with local energy minima, that is, attractor states, the maximum entropy model was used to evaluate the transition rates between attractors of the default‐mode network (DMN) and the frontoparietal network (FPN) by calculating the energy barrier between each pair of attractors (Watanabe et al., 2014). An energy landscape in which local energy minima represent attractors characterized by specific patterns of community allegiance was constructed (Ashourvan et al., 2017). Except for the maximum entropy methods, detection of metabolism was also introduced to examine the energy consumption of human brain. More evidences demonstrated that regions with a high degree of functional connectivity have high energy efficiency (Tomasi, Wang, & Volkow, 2013) and require greater blood supply (Liang et al., 2013).

In addition to the fMRI, some other brain imaging techniques, such as MEG (Krzeminski et al., 2020) and PET (Tomasi et al., 2013, 2017), have been utilized to analyze the energy consumption of human brain. The study in (Krzeminski et al., 2020) suggested that energy consumption of attractors in the maximum entropy model becomes higher in Juvenile myoclonic epilepsy (JME) patients. The work in Tomasi et al., (2017) measured the energy demand of resting state fMRI networks using FDG‐PET. Association among the absolute glucose metabolism, temporal metabolic connectivity (t‐MC) and resting‐state fMRI was analyzed. The results implied that mismatch occurred between metabolic and functional connectivity patterns, which reflected differences in the temporal characteristics of the glucose metabolism and brain activation.

In recent years, studies in the neuroscience field have suggested that a promising paradigm to explain a complex brain network system is the linear network control theory (Gu et al., 2015; Honey et al., 2009; Kim et al., 2018). For instance, it was employed to derive the minimum energy required to control the brain network system, which was divided into two kinds of brain regions, called drivers and nondrivers, respectively (Kim et al., 2018). The drivers are directly controlled by external stimuli, while the nondrivers are influenced by the drivers through connectivities and interactions between them. By calculating the energy consumption of controlling nondrivers by drivers, it was suggested that brain networks are organized along energetically favorable principles.

Nevertheless, the existing methods might have some limitations on the following aspects. (1) They mainly focused on the energy cost of some local ROIs (Ezaki et al., 2017; Krzeminski et al., 2020; Watanabe et al., 2014). Only a small amount of methods computed the energy cost for ROIs in the whole brain (Gu et al., 2018; Kang et al., 2017), where spatial interactions between different ROIs have not yet been fully considered. (2) The previous studies built the energy landscape from the viewpoint of either the brain function (Ashourvan et al., 2017; Ezaki et al., 2017; Kang et al., 2017; Watanabe et al., 2014) or the brain structure (Gu et al., 2017; Kim et al., 2018) individually. Correspondence between functional and structural connections has not been well examined from the viewpoint of the energy consumption. (3) Most of the previous studies used resting state fMRI to analyze the energy consumption of brain networks (Ashourvan et al., 2017; Ezaki et al., 2017; Gu et al., 2018; Kang et al., 2017; Krzeminski et al., 2020; Noack et al., 2017; Tomasi et al., 2017; Watanabe et al., 2014). Depicting what the energy consumption characteristics of brain networks are also needs to be addressed in task‐related studies.

On the other hand, modeling of the macro‐scale network connectivity and interaction is another fundamental issue for characterization of energy consumption of functional brain networks. Recent studies suggested that spatial overlaps of functional networks derived from fMRI data can reveal a fundamental organizational principle of the human brain (e.g., Harris & Mrsic‐Flogel, 2013; Xu et al., 2016; Yuan et al., 2018). Our previous study suggested that the dictionary learning can reconstruct robust functional networks (Li et al., 2014; Lv et al., 2015; Zhang et al., 2017; Zhao et al., 2016) and the brain function can be characterized by spatial overlaps/interactions of connectome‐scale brain networks (Jiang et al., 2018; Yuan et al., 2018).

In this article, in order to overcome the above‐mentioned limitations and fully explore energy consumption characteristics of the macro‐scale network connectivity and interaction, we present a novel framework for characterizing and assessing the energy consumption of functional networks during their spatial interactions throughout the whole brain. Specifically, a method integrating dictionary learning and hierarchical clustering is designed to obtain macro‐scale consistent brain network templates from tfMRI data. Then, spatial interactions between functional networks in the whole brain are modeled by a networked system and the minimum control energy is computed based on the driver‐to‐nondriver optimal control scheme (Kim et al., 2018). Afterwards, energy consumption patterns of brain networks are comprehensively analyzed across different tasks and individuals. Moreover, the relationship between the energy consumption of functional networks and their structural connection patterns is examined. Extensive experimental results on four HCP tfMRI datasets reveal inherent energy consumption characteristics of macro‐scale network interactions.

Compared with the previous studies, the main contributions of the article are as follows. (1) Energy consumption of the functional network interaction in the whole brain is considered, which would bring a more comprehensive understanding for energy characteristics of human brain. (2) Energy consumption is analyzed from viewpoints of both the functional interaction and structural connection, and correspondence between them is established. (3) Task‐related fMRI data are used to characterize the energy consumption of brain networks. Consistency of energy consumption patterns across different tasks is examined to further reveal regularity of the human brain function.

MATERIALS AND METHODS

Overview

As shown in Figure 1, the computational pipeline of the proposed framework consists of four main steps. In the first step, the whole‐brain tfMRI data of each subject in HCP Q1 and Q3 releases are decomposed by sparse coding (Lv et al., 2015), to extract temporal dynamics and the corresponding spatial profiles of functional network components, represented by the dictionary and sparse weighting coefficients, respectively. In the second step, a hierarchical clustering method is designed to compute the common group‐wise functional networks (CGFNs) from the network components learned in the first step. In the third step, functional connectivities and interactions among brain networks of individuals are modeled by a weighted graph. At last, the control energy of each node in the graph is calculated by the linear optimal control theory (Kim et al., 2018). In addition, energetically favorable nodes of each individual are identified and their energy consumption distributions are further assessed on the group‐wise level to characterize the energy consumption of CGFNs statistically.

FIGURE 1

Overall computational pipeline of the proposed framework. (a) The temporal dynamics and the corresponding spatial maps of brain network components are computed from the tfMRI data of each subject by the dictionary learning and sparse representation. (b) Common group‐wise functional networks (CGFNs) are obtained by a hierarchical clustering algorithm. (c) Functional connectivities and interactions among CGFNs are represented by a weighted graph model, where the nodes and edges denote the CGFNs and their spatial interactions, respectively. (d) The control energy consumption of each CGFN is calculated by the linear optimal control theory and the energy consumption characteristics are comprehensively analyzed and assessed

For the detailed description of the proposed method in this article, the following definitions and notations are used.

input matrix of tfMRI of th subject.

dictionary.

) atom in the dictionary.

coefficient matrix.

the th common group‐wise functional network.

the th functional network component.

functional connectivity and interaction graph with a node set and an edge set .

the edge between nodes and in .

the weight of the edge .

the state vector of nodes in .

control input of the dynamics of .

control energy of the dynamics of .

the state vector of the driver nodes.

the state vector of the nondriver nodes.

‐1 norm.

norm.

In addition, correspondences between abbreviations used in this article and their full names are listed.

CGFNs: common group‐wise functional networks.

EMF/ELF: energetically most/least favorable.

EMFN/ELFN: EMF/ELF nodes.

EMFN‐I/ELFN‐I set: EMF/ELF node set of individuals.

J‐EMFN‐I/J‐ELFN‐I set: joint EMF/ELF node set of individuals.

EMFN‐G/EMLN‐G set: group‐wise EMF/ELF node set.

J‐EMFN‐G/J‐ELFN‐G set: group‐wise joint EMF/ELF node set.

LFD‐I/HFD‐I set: individual node set with a low/high fiber density.

LFD‐G/HFD‐G set: group‐wise node set with a low/high fiber density.

RBO: rank biased overlap.

Data acquisition and preprocessing

We use the publicly released high‐quality tfMRI and DTI data in Human Connectome Project (HCP) (Q1 and Q3 releases) to develop and evaluate the proposed method. Four tfMRI datasets including motor, language, social, and working memory tfMRI data are collected from 61 subjects in the HCP Q1 release as the test bed. Note that, although the HCP Q1 release contains 68 subjects, we select 61 subjects of them since data of the other seven subjects are incomplete. Correspondingly, the DTI data of 61 subjects are acquired to examine the relationship between energy consumption of functional networks and their structural connections. Besides, another 61 subjects are randomly chosen from the HCP Q3 release for reproducibility verification. The main acquisition parameters of the tfMRI data are as follows: 90 × 104 × 72 dimension, 220 mm FOV, 72 slices, TR = 0.72 s, TE = 33.1 ms, flip angle = 52°, BW = 2290 Hz/Px, in‐plane FOV = 208 × 180 mm, and 2.0 mm isotropic voxels. The data have been preprocessed by the HCP pipeline (Glasser et al., 2013) and we perform a further refinement for these data by the FSL tools including additional spatial smoothing with FWMH = 4 mm and high‐pass filtering (>0.01 Hz) to remove global drift. (Barch et al., 2013; Smith et al., 2013). For the DTI data, the main acquisition parameters are 168 × 144 × 111 dimension, 111 slices, TR = 5,520 s, TE = 89.5 ms, flip angle = 78°, BW = 1488 Hz/Px, in‐plane FOV = 210 × 180 mm, b values = 1000, 2000, and 3000 s/mm2, 1.25 mm isotropic voxels. The DTI data are preprocessed by the HCP diffusion pipeline (Glasser et al., 2013), including the normalization of b0 image intensity across runs, correction for EPI susceptibility and eddy‐current‐induced distortions, gradient nonlinearities, subject motion, and application of a brain mask. In addition, FA maps are obtained by fitting a diffusion tensor model using the FSL‐FDT toolkit (Behrens et al., 2003).

Sparse representation of tfMRI data

Assume that the tfMRI time series (of dimension , where is the number of timepoints and is the number of voxels) of the th subject in each task constitute the input signal matrix . We decompose by an ‐1 regularized online dictionary learning algorithm (Lv et al., 2015; Mairal, Bach, Ponce, & Sapiro, 2010) to obtain a sparse representation of . Specifically, dictionary learning is to learn an over‐complete dictionary , composed of “atoms” ), and a weighting coefficient matrix , such that each input signal vector is represented by a linear combination of atoms, that is, or . The dictionary learning can be solved by a minimization problem. The cost function is composed of two terms, that is, the reconstruction error and the sparsity measure . The weight regularizes the trade‐off between two terms. denotes the constraint on the dictionary . By imposing a sparsity requirement on , we can obtain the sparse representation of . In this article, the online dictionary learning algorithm (Mairal et al., 2010) is adopted to solve the minimization problem in Equation (1) subject to the constraint in Equation (2). The obtained dictionary and weighting coefficient characterize the underlying temporal dynamics and spatial maps of functional network components of , respectively (Lv et al., 2015). In this way, network components for the total number of subjects across four tasks in both the HCP Q1 and Q3 releases are obtained.

Computation of common group‐wise functional networks

To create a common reference space of group‐wise functional network templates for all subjects across four tasks and across the HCP Q1 and Q3 releases, we use a hierarchical clustering method to cluster all the network components. It has been demonstrated in our previous study (Yuan et al., 2018) that the hierarchical clustering can generate meaningful intermediate results with much less computational cost. In this article, the clustering method contains two layers, as shown in Figure 1. In the first layer, the affinity propagation algorithm (Frey & Dueck, 2007) is utilized to cluster the number of learned network components in one task of the HCP Q1 (or Q3) release across all number of subjects to obtain task‐related functional networks. In the second layer, the results of the first layer, that is, cluster centers of all four tasks in both the HCP Q1 and Q3 releases are pooled together, and further clustered via the spectral clustering algorithm (Zhao et al., 2016) to yield common group‐wise functional networks . In both layers, the similarity between two networks (or network components) is measured by the spatial overlap rate as the summation of the minimum activation value in each voxel between two networks (or network components) over the summation of the averaged activation value in each voxel (Zhao et al., 2016). where is the activation value of the kth voxel in the ith functional network (or network component), and denotes the voxel number. More specifically, the activation value in the affinity propagation clustering equals to the weighting coefficient of the sparse representation, while in the spectral clustering is computed by the average weighting coefficient of the kth voxel in all network components of the ith cluster obtained by the affinity propagation clustering. It is worthwhile to point out that, two kinds of different clustering techniques are employed in two layers, respectively. The reason is that, without a priori information, it's unlikely to pre‐specify a well‐turned cluster number for a large number of network components (24,400 in our study) in the first layer. The affinity propagation clustering algorithm does not require the number of clusters to be pre‐specified. Thus, the affinity propagation clustering is used in the first layer. In contrast, for the spectral clustering algorithm in the second layer, a priori information about the cluster number is available as 199 for group‐wise consistent resting‐state networks (Zhao et al., 2016). Therefore, the cluster number is assigned to be 200 in the second layer of our method for computing functional network templates using tfMRI data. The main reasons contain two aspects. (1) Our previous studies on generation of fine‐granularity functional networks (Lv et al., 2015; Zhao et al., 2016; Zhao, Ge, & Liu, 2018) have revealed that most functional networks occur in both resting states and task performances, and only a small number of networks are strictly task‐evoked. (2) The work in (Zhao et al., 2016) created 199 initial clusters, and 55 of them were irreproducible. At last, 144 networks remained as templates. Likewise, 135 network templates were obtained using convolutional neural networks (Zhao et al., 2018). Based on these results, we infer that 200 would be suitable as the initial cluster number in the second layer, since if we choose a greater cluster number, many of the obtained networks are irreproducible or irrational, which may result in an excessive computational cost.

Modeling of functional connectivity and interaction of CGFNs

In this article, we represent the functional connectivities and interactions between CGFNs by a weighted graph , where denotes the set of nodes and is the set of edges, as shown in Figure 1c. The edge between a pair of nodes and represents the functional connectivity and interaction between them. For each subject, the node set contains functional network components , which have the highest similarity with the network templates . For two network components and , we compute the spatial overlap rate using Equation (3) as the weight of the edge between them. The weight measures to what extent the spatial interaction occurs between each pair of network components of individuals and reflects the strength of functional connectivity and interaction between them.

Note that, more than one functional network component would be similar to a CGFN in some cases. But we only choose the network component that has the highest similarity with the CGFN as its correspondence. That is, a CGFN corresponds to only one network component for each task of a subject. Likewise, one network component can correspond to at most one CGFN. Therefore, the node in the is unique and no repetitive edge occurs for a same node.

Minimum control energy of CGFNs in functional connectivity and interaction

It has been suggested in (Kim et al., 2018) that some brain regions, called drivers, are independently manipulated by external controls and can guide the brain states along diverse trajectories. In contrast, the other brain regions are driven by the drivers and cannot be directly influenced by external inputs. The drivers are the functional networks that first respond to the external task stimuli, and then activate and regulate nondriver functional networks to perform the tasks together. The experimental results in this article (Section 3.3) have indicated that limbic networks are more inclined to act as the drivers to influence other networks, while visual networks are more likely to be the nondrivers to make the whole brain use the energy more efficiently.

From the perspective of the linear network control theory, the driver nodes are directly controlled by external stimuli, while the nondriver nodes are influenced by the driver nodes through functional connectivities and interactions between brain networks. Thus, how much energy the brain consumes when the brain networks are driven is of particular importance for understanding the functional principles of the human brain.

To explicitly characterize and assess the control energy characteristics of brain networks, we identify the energetically favorable CGFNs in the . Following the idea of (Kim et al., 2018), we approximate functional connectivities and interactions between brain networks as a linear, time‐invariant system. As a result, the describing spatial overlaps among brain networks can be modeled by: where denotes the state vector of nodes in the . The state represents the level of activity in each brain network (Kim et al., 2018). Assuming that nodes in the are drivers and the other nodes are nondrivers, denotes the control input vector, which is related to the external stimuli. and are the state transition matrix and input transform matrix, respectively. is composed of the weights and measures the spatial overlap rates between all pairs of brain network components. Equation (4) can be rewritten by: where and are the state vectors of the driver and nondriver nodes, respectively. , , , and . is a M‐dimensional identity matrix, and . Note that, we assume in this article that spatial overlap rates of brain networks do not change over time. As a result, the system denoted by (4) and (5) is linear time invariant. If the spatial overlap rates change, the will be modeled by a linear time‐varying system, where the state transition matrix is time‐varying. The energy consumption analysis on this kind of system will be a future study.

It has been proven that, if the system represented by (5) is controllable, for any desired state , there exists a control input , such that at time (Hespanha, 2009). The control energy of is Following the assumption in (Kim et al., 2018), the original graph can be simplified by a new graph , which only involves edges from the driver nodes to the nondriver nodes. That means , , and . Accordingly, the minimum control energy of the dynamics along the simplified graph is It needs to be pointed out that, (7) is a first‐order approximation to (6). This approximation requires at least as many driver nodes as nondriver nodes for to be invertible (i.e., ). Based on this assumption, the study in (Kim et al., 2018) further considered a range of nondriver fractions from 0.1 to 0.4. The analysis results showed that the linear network control model holds for this range and the first‐order energy approximation is sufficiently accurate (the error remains below approximately 5%), and thus the fraction of nondriver nodes was chosen as 0.2. In our study, we follow the conclusion in (Kim et al., 2018) and also set the fraction of nondriver nodes to be 0.2.

The minimum control energy in (7) can be rewritten by a more concise form where , , . is the adjugate matrix of , and is the determinant of , which acts as a scaling factor for the total energy.

Like the work in (Kim et al., 2018), in this article, we set the initial states to be a neutral baseline. In order to assess the energy consumption characteristics of brain networks as fully as possible, we compute the control energy required to drive brain networks from initial states to 10,000 random desired states, which are set to be and , i.e., each element in and is a random number between and .

Identification of energetically favorable functional networks

Utilizing the expression in (8), we identify nondrivers that minimize and maximize this energy, which are defined as the energetically most favorable (EMF) nodes and energetically least favorable (ELF) nodes, respectively, from the perspective of the linear optimal control. We denote the energetically most favorable nodes and energetically least favorable nodes by the abbreviations EMFN and ELFN, respectively. In what follows, we discuss how to compute the individual and group‐wise EMFN and ELFN, respectively.

First, for each subject in each task, the node set is composed of network components. Every time we select one node as the nondriver, and the other nodes are used as the drivers to jointly control the nondriver to reach the given desired states. For each desired state, the energy consumption is calculated using (8). We compute the total energy consumption for all 10,000 desired states. In this way, we go through all nodes and obtain number of total energy consumption values. We sort the nodes in ascending order with respect to these energy consumption values, and select top 20% of nodes as the EMF node set of individuals (EMFN‐I set). The last 20% of nodes are chosen as the ELF node set of individuals (ELFN‐I set).

As a complex biological system, human brain would work in a concurrent way, where a number of functional networks are simultaneously driven. Thus, we consider the joint EMFN‐I (J‐EMFN‐I) set and joint ELFN‐I (J‐ELFN‐I) set, where multiple nondriver nodes are driven, to better capture the concurrent characteristics of functional networks. Here, we consider the J‐EMFN‐I set as an example and the J‐ELFN‐I can be obtained in the same way. Computation of the J‐EMFN‐I set for each individual contains three steps.In the same way, we can obtain the J‐ELFN‐I set including nodes with maximum energy consumption values by the greedy search, whose number is also 20% of .

Initialization. We select the node with the minimum energy consumption value from the EMFN‐I set and add it into the nondriver set.

Greedy search for a new EMF node. Assuming the current nondriver and driver sets consist of and numbers of nondriver and driver nodes, respectively, every time we add one node from the driver set into the nondriver set temporarily, and then compute the energy consumption corresponding to the case that number of driver nodes drive number of nondriver nodes. We go through all the nodes in the driver set and obtain energy consumption values. Then, we select the node with the minimum energy consumption value and really add it into the nondriver set. At this time, the new nondriver and driver sets contain and nodes, respectively.

Repeat step (b) until the number of nodes in the nondriver set reaches 20% of . Then, all nodes in the nondriver set constitute the J‐EMFN‐I set.

After identifying the individual EMFN/ELFN sets, we continue to define the group‐wise EMFN/ELFN sets. To this end, we first examine the EMFN‐I set for all individuals in each task. The group‐wise EMFN (EMFN‐G) set is composed of the nodes that occur in the EMFN‐I set of at least 20 subjects (i.e., 33% of all 61 subjects). In other words, each node in the EMFN‐G set is energetically most favorable for at least 20 subjects. Likewise, the group‐wise ELFN (ELFN‐G) set contains nodes that occur in the ELFN‐I set of at least 20 subjects. Correspondingly, we construct the group‐wise joint EMFN/ELFN (J‐EMFN‐G/J‐ELFN‐G) sets with the nodes that occur in the J‐EMFN‐I/J‐ELFN‐I sets of at least 20 subjects.

RESULTS AND DISCUSSIONS

Functional brain networks and their connectivity

To infer functional network components (Section 2.3), the dictionary size is set to 400 and is set to 1.5 according to our prior experiences (Lv et al., 2015). In fact, the spatial maps and temporal patterns of brain networks are not very sensitive to the dictionary size from 100 to 600 (Lv et al., 2015). In addition, our extensive observations show that as long as the sparsity level is within a reasonable range (in this study, from 1.2 to 2.0), the sparse representation results are largely the same. It is worthwhile to point out that, in this article, we use the tfMRI data of the same 68 subjects in HCP Q1 as those in (Lv et al., 2015). The only difference is that we only use the data from four tasks, instead of all seven tasks. Therefore, we adopt the same parameters as (Lv et al., 2015) for dictionary learning.

After learning the individual network components of all subjects across four tasks in the HCP Q1 and Q3 releases, we use the proposed hierarchical clustering method to cluster these components. In the first layer, 2,089, 2,032, 2,058, and 2,240 clusters are obtained in the motor, language, social, and working memory tasks of the HCP Q1 release, respectively. Likewise, 2,137, 2,097, 2,057, and 2,235 clusters are generated in the first layer for the HCP Q3 release, respectively. In the second layer, 200 initial clusters are obtained with the spectral clustering. We evaluate them as follows. For each initial cluster, a template is computed by averaging all the spatial maps of individual network components within this cluster. Then we find the correspondence between the template and the network components of each subject. The network component with the maximum spatial overlap rate to the template is taken as its correspondence. If the spatial overlap rate between a template and its correspondence is <0.2 in more than 50% subjects, the template is taken as irreproducible and thus is eliminated. Finally, 110 templates in total remain and form the CGFNs, that is, .

To examine correspondence between the resultant CGFNs and the network components of subjects, we give a box plot of the number of subjects related to each CGFN in four tasks of HCP Q1, as shown in Figure 2, where the horizontal axis corresponds to the task and the vertical axis represents the number of subjects related to each CGFN. As can be seen from Figure 2, most CGFNs correspond to a large number of subjects in each task. In total, 66, 60, 63, and 64 CGFNs are generated by network components of more than 90% subjects in the motor, language, social, and working memory tasks, respectively.

FIGURE 2

The number of subjects related to CGFNs in four tasks

We would like to point out that, a large‐scale functional network overlap is a general property of brain functional organization (Lv et al., 2015; Xu et al., 2016). To delicately characterize spatial overlaps among brain networks needs a sufficient number of them. If a too small number of networks are involved, spatial coverage areas of some networks would be too large, which may result in that different networks are merged by mistake and spatial interactions between them are ignored. Thus, a reasonably large number of brain networks are of benefit to precisely capture spatial interaction patterns of the functional organizations. In this article, 110 fine‐grained functional network templates are obtained. Figure 3 shows spatial maps of some CGFNs, which are consistent with 10 well‐defined RSNs (Smith et al., 2009). Spatial maps of all 110 CGFNs are shown in Figure S1. The results in Figure 3 suggest that the sparse representation and hierarchical clustering can generate meaningful and high‐quality functional networks.

FIGURE 3

Correspondence between the obtained CGFNs and 10 well‐studied RSNs

In addition to the 10 CGFNs that correspond to the well‐defined RSN templates, for the left 100 CGFNs, some of them are likely to be the fine‐grained subdivisions of the RSN templates. For example, #12, #20, and #71 CGFNs in the supplementary materials might be the fine‐grained subdivisions of #6 RSN template. On the other hand, some CGFNs might be new brain networks, for example, #91 and #78 CGFNs, which need to be further studied.

The functional connectivity matrices of the 110 CGFNs for four tasks are shown in Figure 4. It is obtained by computing the average connectivity strength (spatial overlap rate) between all pairs of brain network components over all subjects. It can be found that the functional connectivity of brain networks presents different strength levels across four tasks. Especially, brain networks have stronger functional connectivity in motor and social tasks.

FIGURE 4

Functional connectivity of brain networks in four tasks

Assessment of minimum control energy on individual level

We first examine the difference on MCE values when one nondriver is driven by drivers. The violin plots in Figure 5 show distributions of the maximum and minimum MCE values of all subjects in four task of the HCP Q1 release. Ratios of the maximum to minimum MCE values are 5.35, 4.42, 5.57, and 5.22 in the motor, language, social, and working memory tasks, respectively. We perform the t‐test for examining the difference between the maximum and minimum MCE values. The p values are 2.4206e‐33, 5.9993e‐45, 1.1191e‐23, and 1.4447e‐39 for motor, language, social, and working memory tasks, respectively. It can be seen that the energy consumption is significantly distinct for driving different brain networks in all four tasks.

FIGURE 5

Violin plots of maximum and minimum MCE values of all subjects in four tasks

Then, we compute the MCE of each network component when it is driven separately for each subject in four tasks of the HCP Q1 release and obtain four EMFN‐I sets, each of which contains 22 nodes. We count the nodes that occur in all these four EMFN‐I sets for each subject. The average number is 2.7 and its standard deviation is 1.5. Likewise, the average number of common ELF nodes occurring in all four tasks is 5.8 and its standard deviation is 1.6. The average number of common ELF nodes is greater than that of common EMF nodes, which means that the ELF nodes are more stable across different tasks and subjects.

At last, we count the nodes that are energetically most favorable in some tasks but energetically least favorable in the other tasks. To this end, we compute for each subject the union of four EMFV‐I sets as well as the union of four ELFV‐I sets, and then compute the intersection of the two resultant union sets. The average number of vertices that are energetically most favorable in some tasks but energetically most favorable in the other tasks is 5.4 and its standard deviation is 2.6. For instance, Figure 6 shows the radar chart of the variation of energy consumption roles of #8 subject across four tasks. The three layers from the inside to outside of the radar chart represent the EMF, NEMF‐NELF (neither energetically most favorable nor energetically least favorable) and ELF roles of energy consumption, respectively. That is to say, the nondriver set contains three categories of nodes, EMF, NEMF‐NELF, and ELF nodes, from the viewpoint of roles of energy consumption. The 11 vertices of the radar chart represent 11 CGFNs of #8 subject. Four colored curves correspond to the four tasks, respectively. For example, the red curve of the motor task intersects with the axis of #107 CGFN at the EMF layer, indicating that #107 CGFN is an EMF node in the motor task. The four curves show how the energy consumption roles of 11 CGFNs of #8 subject shift with the change of task designs. It can be seen that brain networks can exchange their energy consumption roles between different tasks.

FIGURE 6

Variation of energy consumption roles of #8 subject across four tasks

To further quantitatively examine the consistency of energy consumption roles between different tasks, for each task, we arrange three energy consumption roles of the CGFNs to yield a vector. Then, we compute the Spearman's rank correlation coefficient between two vectors corresponding to any two tasks, as shown in Table 1. The Spearman's rank correlation coefficient measures the strength and direction of the association between two ranked variables and reflects their consistency of ranking. It varies from to . The value of 1 means a perfect association between two ranks. The value of 0 means no association and the value of means a perfect negative association. Because all the p values in Table 1 are larger than 0.05, the consistency of energy consumption roles is low between different tasks. That means, distributions of the energy consumption roles have significant differences across tasks on the individual level.

TABLE 1

Spearman's rank correlation coefficients of energy consumption roles between different tasks

	Motor	Language	Social	WM
Motor	/	0.18	0.57	0.21
Language	0.18	/	0.39	0.20
Social	0.57	0.39	/	0.17
WM	0.21	0.20	0.17	/

Assessment of minimum control energy on group‐wise level

We first consider the case that N − 1 driver nodes drive one nondriver node. Table 2 shows the EMFN‐G nodes that are EMF in at least 60% subjects in the motor task of the HCP Q1 release, where the frequency denotes the percentage of subjects for which an EMFN‐G node occurs in their EMFN‐I sets. The spatial maps of these nodes are shown in Figure 7. These networks are mostly responsible for visual processing and attention. It is reasonable that the least energy is consumed to drive these networks as they usually need to be evoked promptly to perceive and react in a cognitive task. On the other hand, the other networks have strong interactions with these networks to make the whole brain function more efficient.

TABLE 2

EMFN‐G nodes that are EMF in at least 60% subjects in the motor task

CGFN No.	Frequency (%)
#18	70.5
#21	82.0
#24	86.9
#27	65.6
#29	68.9
#44	60.7
#72	67.2
#106	75.4
#108	68.9

FIGURE 7

Spatial maps of the EMFN‐G nodes that are EMF in at least 60% subjects

Likewise, the ELFN‐G nodes that are ELF in at least 60% subjects in the motor task are given in Table 3. Their spatial maps are shown in Figure 8. Most of them are limbic networks (e.g., #1, #3, #15, #35, #63, and #79), which are directly related to emotion, motivation, memory, and valuation. It is reasonable that high energy consumption is required to drive these networks as they are responsible for high‐level mental functions and neural activities. Viewed from another perspective, high energy consumption means that these networks are not easy to be driven by other brain networks. Since brain networks are organized along energetically favorable principles (Kim et al., 2018), it can be inferred that the networks in Figure 8 are more inclined to be drivers to make the whole brain use the energy more efficiently.

TABLE 3

ELFN‐G nodes that are ELF in at least 60% subjects in the motor task

CGFN No.	Frequency (%)
#1	67.2
#3	83.6
#15	86.9
#23	68.8
#26	60.7
#33	67.2
#35	68.9
#55	91.8
#58	81.7
#61	63.9
#63	90.2
#79	75.4

FIGURE 8

Spatial maps of the ELFN‐G nodes that are ELF at least 60% subjects

Then, we compute the minimum control energy of J‐EMFN‐G and J‐ELFN‐G nodes, respectively, when multiple nondriver nodes are driven. Tables 4 and 5 show the J‐EMFN‐G/J‐ELFN‐G sets and the frequency that a node occurs in the J‐EMFN‐I/J‐ELFN‐I sets.

TABLE 4

J‐EMFN‐G set of the motor task in HCP Q1 release

CGFN No.	Frequency (%)
#7	36.1
#12	34.4
#20	32.8
#27	32.8
#35	41.0
#70	32.8
#71	36.1
#77	41.0
#83	37.7
#84	32.8
#88	32.8
#101	37.7
#102	34.4
#107	41.0

TABLE 5

J‐ELFN‐G set of the motor task in HCP Q1 release

CGFN No.	Frequency (%)
#1	42.6
#3	63.9
#15	60.6
#23	45.9
#33	44.3
#35	63.9
#45	41.0
#52	34.4
#55	82.0
#58	60.6
#59	32.7
#61	45.9
#63	68.9
#69	54.1
#70	36.1
#73	44.3
#79	45.9
#87	34.4
#109	34.4

It can be observed from Tables 2 to 5 that the EMFN‐G and J‐EMFN‐G sets, which contain 26 and 14 nodes, respectively, share only four same nodes, while the ELFN‐G and J‐ELFN‐G sets, which contain 21 and 19 nodes, respectively, share 17 same nodes. That means, almost all J‐ELFN‐G nodes occur in the ELFN‐G set. Therefore, the ELF nodes present higher stability on the group‐wise level than the EMF nodes.

The statistical results of the minimum control energy of EMF and ELF nodes in the other three tasks are similar to those in the motor task on both the individual and group‐wise levels. To further examine the consistency of the group‐wise minimum control energy across four tasks, we compute elements in the intersections of EMFN/ELFN sets across four tasks, as shown in Table 6. The first column gives the intersection of EMFN‐G sets of four tasks and the second column presents that of ELFN‐G sets of four tasks. Their spatial maps are shown in Figures 9 and 10, respectively. It can be found that 17 and 18 nodes occur as the EMF and ELF nodes in all four tasks, respectively, while the other nodes have different characteristics of the minimum control energy in four tasks. For instance, the ELFN‐G set in the motor task contains 21 nodes, where 18 nodes (about 85.7%) behave as ELFN‐G nodes in all four tasks, and only one node (#46 CGFN) is NEMF‐NELF in the other three tasks. The other two nodes (#52 and #70 CGFNs) are ELF in the motor, language, and working memory tasks, and are NEMF‐NELF in the social task. Likewise, the EMFN‐G set in the motor task is composed of 26 nodes, where 17 nodes (about 65.4%) are EMF nodes in all four tasks, and three nodes (#80, #83, and #84 CGFNs) are NEMF‐NELF in the other three tasks. Two nodes (#51 and #102 CGFNs) are EMF in the motor and language tasks, and are NEMF‐NELF in the working memory and social tasks. One node (#95 CGFN) is EMF in the motor, social, and working memory tasks, and is NEMF‐NELF in the language task. The other three nodes (#10, #47, and #53 CGFNs) are EMF in the motor, language, and working memory tasks, and are NEMF‐NELF in the social task.

TABLE 6

Intersections of group‐wise EMFN/ ELFN sets across four tasks in HCP Q1 release

Intersection of EMFN‐G	Intersection of ELFN‐G	Intersection of J‐EMFN‐G	Intersection of J‐ELFN‐G
#16	#1	#84	#1
#18	#3	#102	#3
#21	#15	/	#15
#24	#23	/	#23
#27	#26	/	#33
#29	#33	/	#35
#30	#35	/	#55
#38	#45	/	#58
#44	#55	/	#61
#49	#58	/	#63
#72	#61	/	#69
#75	#63	/	#73
#86	#68	/	/
#103	#69	/	/
#106	#73	/	/
#108	#77	/	/
#110	#79	/	/
/	#87	/	/

FIGURE 9

Spatial maps of 17 common EMF nodes in all four tasks

FIGURE 10

Spatial maps of 18 common ELF nodes in all four tasks

The 17 common EMF CGFNs in the first column of Table 6 and Figure 9 contain all 9 CGFNs in Table 2 and Figure 7, which are mainly visual related cortices. Except these nine CGFNs, most of the other eight ones are also visual and attention networks. Visual networks can be easily driven by other networks (e.g., limbic networks shown in Figure 8), since they require lower energy consumptions. Thus, it can be inferred that visual networks are more likely to be the nondrivers, because brain networks are organized along energetically favorable principles. In contrast, the 18 common ELF CGFNs in the second column of Table 6 and Figure 10 contain all 12 CGFNs in Table 3 and Figure 8, most of which are limbic networks. Except these 12 CGFNs, the other six CGFNs' spatial maps are also among the limbic networks and frontal‐temporal networks. The limbic networks are primarily responsible for controlling emotional drives and memory formation, while the frontal‐temporal networks are related to high‐level brain activities. Thus, they consume more energy and are difficult to be driven. We can infer that these networks might act as the drivers in the brain functional organization.

The third column in Table 6 shows the intersection of the J‐EMFN‐G sets in four tasks and the fourth column gives the intersection of the J‐ELFN‐G sets in four tasks. It can be seen from Table 6 that more ELF CGFNs occur in all four tasks (the second and fourth columns), which implies the ELF CGFNs are more stable and consistent than the EMF ones. From the visualization results of EMF and ELF CGFNs in Figures 7 and 8, we can find that most EMF CGFNs are distributed around the cerebral cortex, while the ELF CGFNs mainly lie deep inside the brain. Although activation of the cerebral cortex regions is more likely to be influenced by different task stimuli, information is always transmitted through internal cerebral nerve fibers. Therefore, the ELF CGFNs have more stable energy consumption across different tasks.

Furthermore, we compute the intersection sets between the EMFN‐G and ELFN‐G sets across four tasks, as shown in Table 7. We can find that all these intersection sets are empty. In other words, on the group‐wise level, the nodes that are ELF in some tasks are very likely to be either ELF or NEMF‐NELF in the other tasks, rather than being EMF. Likewise, the nodes that are EMF in some tasks are not likely to be ELF in the other tasks. Therefore, the energy consumption roles of each CGFN remain relatively fixed across different tasks on the group‐wise level, although they exhibit large differences on the individual level.

TABLE 7

Intersection sets between the EMFN‐G and ELFN‐G sets across four tasks

	EMFN‐G in motor task	EMFN‐G in language task	EMFN‐G in social task	EMFN‐G in WM task
ELFN‐G in motor task	∅	∅	∅	∅
ELFN‐G in language task	∅	∅	∅	∅
ELFN‐G in social task	∅	∅	∅	∅
ELFN‐G in WM task	∅	∅	∅	∅

Validation for correspondence between energy consumption and structural fiber density

We conduct experiments for examining the relationship between energy consumption characteristics of the CGFNs and the structural connectomes of the corresponding brain regions. We first compute the fiber density of each functional network component for each subject by the ratio of the number of fibers spatially covered by the network component to the number of nonzero voxels in it. Then, we define four sets involving the fiber density on both individual and group‐wise levels. (1) We sort the 110 nodes of each subject in ascending order with respect to the fiber density values in each task, and select top 20% of nodes as the individual node set with a low fiber density (LFD‐I). (2) Likewise, the last 20% of nodes constitute the individual node set with a high fiber density (HFD‐I). (3) The group‐wise node set with a low fiber density (LFD‐G) is composed of the nodes that occur in the LFD‐I sets of at least 20 individuals (33% of all 61 subjects). (4) The group‐wise node set with a high fiber density (HFD‐G) contains the nodes that occur in the HFD‐I sets of at least 20 individuals.

Figure 11 shows four group‐level MCE curves of the motor task in the HCP Q1 release, where the horizontal and vertical axes denote the number of nondriver nodes and the average MCE values over all subjects, respectively. The green and red curves represent variations of the maximum and minimum MCE values as the number of nondriver nodes gradually increases when computing the J‐EMFN‐I and J‐ELFN‐I sets, respectively. It can be seen that, the MCE value corresponding to 22 ELF nondriver nodes is about one order of magnitude higher than that corresponding to 22 EMF nondriver nodes. This implies that the energy consumption of brain networks presents significant differences. The black and blue curves denote the group‐level MCE evolutions as the HFD‐I/LFD‐I nodes are gradually added into the nondriver node set in descending/ascending order with respect to the density values, respectively. From Figure 11, we know that the MCE of HFD‐I nodes is larger than that of LFD‐I ones.

FIGURE 11

Group‐level MCE evolutions versus the number of nondriver nodes

We count the nodes that simultaneously occur in the EMFN‐I and LFD‐I sets (both of them contain 22 nodes) for each subject, as shown in Figure 12a. Figure 12b shows the number of common nodes that occur in both the ELFN‐I and HFD‐I sets (they also contain 22 nodes) for each subject. Likewise, Figure 12c,d shows the number of common nodes shared by the J‐EMFN‐I and LFD‐I sets, and that shared by the J‐ELFN‐I and HFD‐I sets, respectively. Comparing Figure 12a–c, it can be seen that the EMFN‐I set shares more common nodes with the LFD‐I set than the J‐EMFN‐I set. Likewise, the ELFN‐I set shares more common nodes with the HFD‐I set than the J‐ELFN‐I set.

FIGURE 12

Number of common nodes in the sets with different energy consumption values and fiber densities. (a) Number of common nodes in the EMFN‐I and LFD‐I sets. (b) Number of common nodes in the ELFN‐I and HFD‐I sets. (c) Number of common nodes in the J‐EMFN‐I and LFD‐I sets. (d) Number of common nodes in the J‐ELFN‐I and HFD‐I sets

Afterwards, we examine the correspondence between the MCE and fiber density on the group‐wise level. Table 8 gives the numbers of elements in intersections between the EMFN‐G/ELFN‐G sets and LFD‐G/HFD‐G sets for the motor task when one node is driven. The intersection between the EMFN‐G and LFD‐G sets and that between the ELFN‐G and HFD‐G sets contain 15 and 13 nodes, respectively. In comparison, the EMFN‐G and HFD‐G sets share no node, so do the ELFN‐G and LFD‐G sets. These results further suggest that the brain networks with a higher fiber density are less energetically favorable, while the brain networks with a lower fiber density are more energetically favorable. Table 9 gives the numbers of elements in intersections between the J‐EMFN‐G/J‐ELFN‐G sets and LFD‐G/HFD‐G sets for the motor task when multiple nondriver nodes are driven. From Table 9, the J‐EMFN‐G set shares only three common nodes with the LFD‐G set. In contrast, the J‐ELFN‐G set shares 11 common nodes with the HFD‐G set. This observation is basically consistent with the analysis results of one nondriver case, and also demonstrates that the ELF nodes have higher stability than the EMF ones.

TABLE 8

Correspondence between the MCE and fiber density in the motor task when one node is driven

	EMFN‐G	ELFN‐G
LFD‐G	15	0
HFD‐G	0	13

TABLE 9

Correspondence between the MCE and fiber density in the motor task when multiple nodes are driven

	J‐EMFN‐G	J‐ELFN‐G
LFD‐G	3	0
HFD‐G	2	11

Validation for relationship among functional connectivity, energetic consumption and fMRI amplitude

In the first and second columns of Table 6, we have obtained two intersections of group‐wise EMFN/ELFN sets across four tasks in the HCP Q1 release. Here, we call the elements in these two intersections as common EMF/ELF CGFNs of four tasks. We first examine the average fMRI amplitudes of the common EMF and ELF CGFNs over all subjects in the HCP Q1 release and give the box plots in Figure 13. The t‐test is performed to verify whether the fMRI amplitude has a significant difference between the EMF and ELF brain networks. The p values are 9.5105e‐8, 1.5146e‐7, 1.8873e‐9, and 1.533e‐8 for motor, language, social, and working memory tasks, respectively. These results demonstrate that the EMF brain networks have significantly higher fMRI amplitudes than the ELF ones.

FIGURE 13

The average fMRI amplitudes of common EMF and ELF CGFNs in four tasks. (a) Motor task. (b) Language task. (c) Social task. (d) Working memory task

Furthermore, we use the Spearman's rank correlation coefficient to measure the correlation between the energy consumption and fMRI amplitude. The Spearman's rank correlation coefficients are −0.85, −0.86, −0.87, and −0.88 for four tasks, respectively. Correspondingly, the p values are 9.1829e‐11, 2.35e‐11, 1.2686e‐11, and 5.216e‐12, respectively. The results show a very strong negative correlation between the energy consumption and fMRI amplitude in all four tasks.

Then, we calculate the average functional connection degree, that is, the strength of spatial interactions, of the CGFNs according to the functional connectivity matrix given in Figure 4. The distribution of the MCE is shown in Figure 14, where three axes denote the CGFN index, the average functional connection degree of each CGFN in the motor task, and the corresponding group‐wise MCE value, respectively. It can be found that different brain networks have diverse functional connection degrees as well as energy consumption values.

FIGURE 14

The distribution of the MCE versus the functional connection and CGFN index

Figure 15 shows the box plots of the average functional connection degree of the common ELF and EMF CGFNs in four tasks. The t‐test is also performed to verify whether the average functional connection degree has a significant difference between the EMF and ELF brain networks. The p values are 5.9669e‐14, 1.9443e‐13, 2.1014e‐10, and 2.8507e‐11 for motor, language, social, and working memory tasks, respectively, which show that the EMF brain networks have significantly larger functional connection degrees than the ELF ones.

FIGURE 15

The average functional connection degree of common EMF and ELF CGFNs in four tasks. (a) Motor task. (b) Language task. (c) Social task. (d) Working memory task

In addition, we use the Spearman's rank correlation coefficient to measure the correlation between the energy consumption and average functional connection degree. The Spearman's rank correlation coefficients are −0.94, −0.89, −0.96, and −0.88 for four tasks, respectively. Correspondingly, the p values are 1.0047e‐16, 8.6416e‐13, 1.0468e‐19, and 3.3957e‐12, respectively. The results show a very strong negative correlation between the energy consumption and average functional connection degree in all four tasks.

From Figures 13 and 15, the EMF brain networks have both higher fMRI amplitudes and functional connection strength, which is consistent with the result that brain networks with a high functional connectivity are energy efficient in glucose metabolism (Tomasi et al., 2013). These observations suggest that the human brain is organized in a complicated network topology, which would be optimized and functionally re‐organized for a specific task to facilitate network control. To support complex even variable functional connectivities requires efficient use of energy (Kim et al., 2018). Thus, energetically favorable organization of brain networks may dominate the functional connectivity architecture of human brain.

To further explore the functional and structural connectivity characteristics of the EMF/ELF CGFNs, we define the structural connection degree of a brain network for ith subject by: where denotes the number of fibers whose end points lie in the #j and #k brain network components of the ith subject, is the number of voxels activated in the #j network component of the ith subject. For each task, the average structural connection degree of each brain network over all subjects is For different tasks of each subject, the functional network components corresponding to the same CGFN have slight differences. Therefore, for a given CGFN, structural connections of the functional network components are not completely the same in four different tasks. Figure 16 shows the box plots of the average structural connection degree of common EMF and ELF CGFNs over all subjects for four tasks. The t‐test for the difference of the average structural connection degree between EMF and ELF CGFNs is performed, and the p values are 3.0e‐5, 1.0e‐3, 5.2e‐6, and 2.3e‐5 for motor, language, social, and working memory tasks, respectively. It can be seen that values of the ELF CGFNs are significantly larger than those of the EMF ones. Likewise, we use the Spearman's rank correlation coefficients to measure the correlation between the energy consumption and average structural connection degree. They are 0.80, 0.67, 0.82, and 0.70 for four tasks, respectively. Correspondingly, the p values are 7.6088e‐9, 1.2e‐5, 1.8951e‐9, and 3.0e‐6, respectively. The results show a strong correlation between the energy consumption and average structural connection degree.

FIGURE 16

The average structural connection degree of common EMF and ELF CGFNs in four tasks. (a) Motor task. (b) Language task. (c) Social task. (d) Working memory task

Comparing Figures 15 and 16, we can know that the ELF brain networks have stronger structural connectivity but weaker functional connectivity, and the finding is opposite for the EMF brain networks.

Reproducibility verification for energy consumption characteristics of brain networks

In this experiment for reproducibility verification, we use the HCP Q3 release to compute and assess the energy consumption of CGFNs on both the individual and group‐wise levels. We can obtain similar results to the HCP Q1 release. For instance, Table 10 lists the elements in the intersections of EMFN/ELFN sets across four tasks in the HCP Q3 release. Comparing Table 10 with Table 6, we find that intersections of EMFN‐G sets in the HCP Q1 and Q3 releases share six common nodes and intersections of ELFN‐G sets in two releases share 10 common nodes. In addition, intersections of J‐EMFN‐G sets in two releases share no node and intersections of J‐ELFN‐G sets in two releases share seven common nodes. These results also verify that the energy consumption of the ELF nodes is more stable and reproducible than that of the EMF nodes across four tasks and across two releases.

TABLE 10

Intersections of group‐wise EMFN/ELFN sets across four tasks in HCP Q3 release

Intersection of EMFN‐G	Intersection of ELFN‐G	Intersection of J‐EMFN‐G	Intersection of J‐ELFN‐G
#6	#1	#5	#3
#29	#3	#77	#15
#72	#15	/	#23
#86	#23	/	#26
#92	#25	/	#35
#95	#26	/	#40
#101	#32	/	#50
#103	#33	/	#55
#106	#35	/	#58
#108	#40	/	#63
/	#55	/	#68
/	#58	/	#77
/	#62	/	/
/	#63	/	/
/	#68	/	/
/	#71	/	/
/	#76	/	/
/	#77	/	/
/	#97	/	/

To quantitatively verify the reproducibility of the analysis results on energy consumption of brain networks, we compare three groups of experimental results between the HCP Q1 and Q3 releases.

Comparison of the MCE of CGFNs. Specifically, we calculate for each CGFN in each task of the HCP Q1 release the average MCE value of all 61 subjects when the CGFN acts as a nondriver node separately. We go through all CGFNs and obtain N (N = 110) number of average MCE values, which constitute an N‐dimensional vector . In the same way, we can also obtain an N‐dimensional vector for each task in the HCP Q3 release. Then, we calculate the Spearman's rank correlation coefficient between and , which is 0.78, 0.75, 0.77, and 0.79 for four tasks, respectively. Correspondingly, the p values are 1.2903e‐23, 7.9825e‐21, 1.5323e‐22, and 1.127e‐24, respectively. The results verify consistency of the average energy consumption between two HCP releases.

Comparison of ELFN‐G sets. We examine the consistency of both the ELFN‐G and HFD‐G nodes between the HCP Q1 and Q3 releases. Specifically, we sort the ELFN‐G nodes in the motor task of the HCP Q1 and Q3 releases in descending order with respect to the MCE values, respectively, as shown in the first and third columns of Table 11. Then, we sort the HFD‐G nodes in the motor task of two HCP releases in descending order with respect to the fiber density values, respectively, and the results are given in the second and fourth columns of Table 11. The top three nodes in the ELFN‐G set of the HCP Q1 release are #15, #58, and #35 CGFNs, respectively. They also have high rankings in the HFD‐G set of the HCP Q1 release, as well as in the ELFN‐G and HFD‐G sets of the HCP Q3 release. Their spatial maps are shown in Figure 17. To give a quantitative evaluation, we choose the rank biased overlap (RBO) similarity to compare the results in two releases. The RBO measures the similarity between two ranked lists which may or may not contain the same items. It ranges between 0 and 1. Two ranked lists with a high RBO value are very similar, whereas a low RBO value indicates dissimilar lists. We compute the RBO similarity between the first and third columns of Table 11. The RBO value is 0.66 for the motor task. Likewise, the RBO similarities are 0.63, 0.71, and 0.64 for the language, social, and working memory tasks, respectively. The results further demonstrate the consistency and reproducibility of our analysis on the energy consumption of brain networks.

TABLE 11

Rankings of ELFN‐G and HFD‐G sets in the motor task of HCP Q1 and Q3 releases

Motor task in HCP Q1		Motor task in HCP Q3
ELFN‐G	HFD‐G	ELFN‐G	HFD‐G
#15	#15	#35	#58
#58	#35	#26	#15
#35	#71	#68	#1
#26	#58	#15	#71
#1	#1	#1	#68
#61	#32	#23	#26
#79	#68	#40	#18
#23	#26	#76	#35
#87	#33	#8	#21
#69	#62	#33	#24
#33	#23	#58	#97
#68	#25	#71	#110
#52	#34	#77	#54
#45	#61	#62	#44
#73	#17	#12	#77
#70	#87	#25	#81
#77	/	#32	#51
#46	/	#97	#32
/	/	#46	#90
/	/	#84	#84
/	/	/	#83

FIGURE 17

Spatial maps of three CGFNs with high rankings in the ELFN‐G and HFD‐G sets of two releases

3) Comparison of EMFN‐G sets. We also sort the EMFN‐G and LFD‐G nodes in the motor task of the HCP Q1 and Q3 releases in ascending order with respect to the MCE and fiber density values, respectively, as given in Table 12. The RBO similarities between two releases are 0.56, 0.33, 0.53, and 0.68 for four tasks, respectively. These results further suggest that the EMF CGFNs have lower regularity and stability than the ELF CGFNs.

TABLE 12

Rankings of EMFN‐G and LFD‐G sets in the motor task of HCP Q1 and Q3 releases

Motor task in HCP Q1		Motor task in HCP Q3
EMFN‐G	LFD‐G	EMFN‐G	LFD‐G
#21	#103	#21	#50
#24	#38	#6	#53
#106	#37	#106	#5
#29	#47	#92	#6
#72	#60	#99	#48
#18	#5	#72	#78
#108	#92	#86	#103
#27	#53	#108	#94
#44	#11	#29	#47
#47	#89	#24	#61
#110	#98	#41	#60
#75	#102	#44	#80
#38	#80	#95	#52
#86	#85	#53	#85
#103	#24	#103	#75
#30	#21	#94	#41
#16	#108	#80	#76
#10	#96	#85	#99
#53	#72	#27	#66
#51	#29	#39	#34
#49	#27	#43	#98
#95	#64	#75	#92
#80	#51	#47	#2
#83	#95	#101	#57
#102	#94	#20	#37
#84	#84	#64	#108
/	/	#82	#95
/	/	#50	#96

CONCLUSION AND DISCUSSION

In this article, a novel method for analysis and assessment of energy consumption of functional brain networks was proposed. From a technical perspective, a novel computational framework was designed to create macro‐scale common group‐wise functional networks and model their functional interactions through a weighted graph. Then, the energy consumption for driving different brain networks was calculated from the viewpoint of the linear optimal control theory, and assessed on both individual and group‐wise levels. The framework has been applied and evaluated on four HCP task fMRI datasets and interesting results have been achieved.

From a neuroscientific perspective, our experimental results on both the HCP Q1 and Q3 releases revealed the following interesting observations. First, the obtained macro‐scale CGFNs are in agreement with the existing literature studies, which might provide a common brain network space for various network composition patterns of different cognitive states. Moreover, some fine‐grained subdivisions of well‐defined RSN templates and some new brain networks have been found in the obtained macro‐scale CGFNs, which may facilitate to understand new brain function in task performance. Second, utilization of spatial interactions to model the functional connectivity of brain networks is reasonable and can yield meaningful and consistent energy assessment results. Both energetically most favorable and energetically least favorable functional networks coexist in human brain. Their energy consumption characteristics have significant differences. In particular, the energetically least favorable brain networks require higher energy consumption to be driven and are more stable and consistent across different tasks. Moreover, brain networks potentially take on completely different roles during their spatial interaction process, in the sense of linear network control. For instance, limbic networks require high energy consumption to be driven during spatial interactions of brain networks. They might act as drivers in the brain functional organization. In contrast, visual networks can be easily driven by others and are more likely to be the nondrivers. Third, significant differences exist in energy consumption of brain networks between different subjects. However, group‐wise energy consumption is basically stable and consistent across four tasks. Especially, the energetically least favorable brain networks present remarkable regularity and reproducibility across different tasks as well as different HCP releases. Fourth, energetically most favorable brain networks have significantly higher functional connection degrees, which implies that brain networks with strong functional connectivity are energy efficient to facilitate human brain control in a specific task. Fifth, correspondences between the energy consumption and fiber density of brain networks were found. It is inspiring that correspondences between the least favorable energy consumption and the high fiber density were well established and interpreted according to current neuroscientific knowledge in the field, which suggests that brain networks with complex fiber structures are more stable and more unlikely to be driven during functional brain interactions, and thus might take on more complicated cognitive functions and more important roles involved in information transmission among human brains. The good interpretation and reproducibility suggest the validity of the quantitative assessments and analyses on energy consumptions of brain networks, which provides new understanding of the human brain function.

FIGURE S1. Spatial maps of 110 CGFNs

Click here for additional data file.