A Physarum Centrality Measure of the Human Brain Network

The most important goals of brain network analyses are to (a) detect pivotal regions and connections that contribute to disproportionate communication flow, (b) integrate global information, and (c) increase the brain network efficiency. Most centrality measures assume that information propagates in networks with the shortest connection paths, but this assumption is not true for most real networks given that information in the brain propagates through all possible paths. This study presents a methodological pipeline for identifying influential nodes and edges in human brain networks based on the self-regulating biological concept adopted from the Physarum model, thereby allowing the identification of optimal paths that are independent of the stated assumption. Network hubs and bridges were investigated in structural brain networks using the Physarum model. The optimal paths and fluid flow were used to formulate the Physarum centrality measure. Most network hubs and bridges are overlapped to some extent, but those based on Physarum centrality contain local and global information in the superior frontal, anterior cingulate, middle temporal gyrus, and precuneus regions. This approach also reduced individual variation. Our results suggest that the Physarum centrality presents a trade-off between the degree and betweenness centrality measures.


Results
Spatial distribution of hub nodes and bridge edges. The network hubs of three centrality measures (one standard deviation above the mean), i.e., H(C D,node ), H(C B,node ), and H(C P,node ), were identified according to each centrality map, i.e., C D,node , C B,node , and C P,node . It is noted that the anatomical locations of the obtained hubs were adopted from the predefined template 48 (Table S1). Accordingly, the H(C D,node ) measure appeared in four cortical regions (the precuneus, middle temporal gyrus, superior frontal gyrus [dorsolateral], and postcentral gyrus) in a bilaterally symmetric fashion, and in three other regions (the right precentral gyrus, right calcarine fissure, and left middle occipital gyrus) (Fig. 1A). Furthermore, H(C B,node ) appeared in four cortical regions (the precuneus, superior frontal gyrus [dorsolateral], precentral gyrus, and postcentral gyrus) in a bilaterally symmetric fashion, as well as in three other regions (the left superior frontal gyrus [medial], left anterior cingulate and paracingulate gyrus, and left middle occipital gyrus) (Fig. 1B). Equivalently, H(C P,node ) appeared in nine cortical regions (the precuneus, superior frontal gyrus [dorsolateral], precentral gyrus, postcentral gyrus, left superior frontal gyrus [medial], right calcarine fissure, left middle temporal gyrus, left middle occipital gyrus, and in the left anterior cingulate and paracingulate gyri) (Fig. 1C). Figure 2 shows the cumulative distributions of C B,node and C P,node with a degree distribution in a continuous manner. Accordingly, C B,node and C P,node accounted for 83.47% and 63.67% of the top 50% most connected nodes, respectively. The results indicate that C P,node is homogenously  Cumulative distributions of C B,node and C P,node with degree. Nodes were sorted so that the node with the highest value moved to one, and the node with the lowest centrality value moved to the last index (x-axis). The cumulative distribution of C B,node is shown in blue, the cumulative distribution of C P,node is shown in red, and the degree distribution is shown in green. The network bridge edges are identified when the network edges are greater than one standard deviation (SD) above the mean of each edge centrality measure map. Their Jaccard index is also shown with overlapped bridge edges.
was positively correlated with C B,node (R 2 = 0.939, P = 1.8606e −48 ), and C D,node was positively correlated with C B,node (R 2 = 0.687, P = 2.6774e −21 ). It was noted that the Jaccard index and regression analysis, including C P,node , exhibited a strong tendency to acquire higher values compared to the values of other models. Notably, most of the network hub regions overlapped to some extent. In all three centrality measures, three cortical regions (the precuneus, superior frontal gyrus, and postcentral gyrus) appeared in a bilaterally symmetric fashion, while two regions (the left middle occipital gyrus and right precentral gyrus) appeared in a lateralized manner. We also calculated J(B(C B,edge ), B(C P,edge )) on the edge level to discover common efficient communication paths between two different measures. Its value was estimated to equal 0.791. Notably, the results suggest that the overlapped paths may be core paths, and they may have an important role in the efficient information flow across brain regions.

Differences of hub nodes and bridge edges.
A post-hoc analysis was performed after the analysis of variance (ANOVA) test on the network hub regions based on three Z-transformed centrality measures (C D,node vs. C B,node vs. C P,node ). Table 2 shows the details of the ANOVA and post-hoc analyses. In the five cortical hub regions (the left precentral gyrus, right superior frontal gyrus [dorsolateral part], left anterior cingulate, and left and right postcentral gyri), C D,node was significantly lower than C B,node a C P,node , but C B,node and C P,node were not significantly different. In the three cortical hub regions (the right precentral gyrus and left superior frontal gyrus [dorsolateral and medial part]), only C P,node was significantly higher than C D,node . In the four cortical hub regions (the calcarine fissure, right precuneus, and left and right middle temporal gyri), C D,node was significantly higher than C B,node and C P,node , and C P,node was significantly higher than C B,node (C D,node > C P,node > C B,node ). A similar tendency was observed in the left precuneus. Table 3 shows the differences of the bridge edges based on the C B,edge and C P,edge values. B(C P,edge ) contained 46 additional bridge edges which B(C B,edge ) did not have. The additional B(C P,edge ) mainly connected with hub nodes, such as the calcarine fissure and the middle temporal gyrus, rather than B(C B,edge ). However, B(C B,edge ) had only 18 additional bridge edges which B(C P,edge ) did not have (Table 4).
Individual variability of hub nodes. The coefficient of variation (CV) was calculated in network hub regions, including the left and right precentral gyri, left and right superior frontal gyri (dorsolateral part), left superior frontal gyrus (medial part), left anterior cingulate, left calcarine fissure, left middle occipital gyrus, left and right postcentral gyri, left and right precuneus, and left and right middle temporal gyri (Table 5). A two-tailed t-test was performed to determine the statistical significance of the differences in CV values among the three centrality measures (C D,node , C B,node and C P,node ). Accordingly, it was found that CV (C D,node ) was significantly lower than CV (C B,node ) (P = 2.5610e −19 , t-test) and CV(C P,node ) (P = 1.2673e −12 , t-test). Additionally, CV(C P,node ) was also significantly lower than CV(C B,node ) (P = 8.2899e −18 , t-test). Furthermore, the values of CV(C B,node ), CV(C D,node ), and CV(C P,node ) were found to lie in the ranges of 0.5777-0.6776, 0.1478-0.1778, and 0.2378-0.3138, respectively.

Discussion
In this study, we proposed a novel methodological framework for defining the importance of network nodes and edges using the Physarum model. Other centrality measures, such as C D and C B , assume that information flows in a network only through the paths that are associated with the shortest connections, but C P considers all possible information flows between brain regions.
Many previous studies have detected brain network hubs and bridges using various measures, such as C D and C B 8,9,[21][22][23][24][25][26][27]51 . These measures of centrality helped the interpretation of the meaning of nodes and edges in the network 3,5 . Accordingly, C D,node , usually defined as the number of connections of the target node, quantified the local properties without global information flow. Furthermore, C B,node identified the node that played an important role with the use of information based on the global flow patterns and on the shortest path concept, while C B,edge used a similar approach to that used by C B,node at the edge level. Equivalently, C P and C B used the global information flow. However, C P used all the possible paths from the Physarum model instead of the shortest path concept and was shown to be affected by local characteristics, such as C D .
As shown in Fig. 2, the C P,node was homogenously distributed among all network nodes compared to C B,node . However, the C P,node considered the optimal paths from the Physarum model independent of the assumption used by other centralities according to which the information flow of a network only spread through the shortest connecting paths. Previous studies have shown the existence of a communication scheme that contradicted the assumption that only the shortest connections are used 37,52 . These studies have shown that C P can be uniformly distributed compared to C B .
The Jaccard index was used to examine the overlap ratio between the sets of each centrality measure. As shown in Table 1, the Jaccard indices between hub sets based on C P,node and C B,node had higher values than those associated with other combinations. As shown in Fig. 4, similar Jaccard index patterns were observed between C P,node and other measures in a continuous manner. The Jaccard index value estimated between bridge sets based on C P,edge and C B,edge also yielded higher values. Thus, the network hub regions determined based on C P,node possessed  Table 1. Jaccard indices between network hubs from three centrality measures. The Jaccard index of the hub regions is the ratio of the number of overlapping hub nodes to the total number of hub nodes based on any two centrality measures. The value of the Jaccard index varies from zero (no overlap) to one (perfect overlap).
www.nature.com/scientificreports www.nature.com/scientificreports/ local and global network properties. Based on global information, some regions of the precentral gyrus and superior frontal and anterior cingulate gyri were defined as network hubs. However, these regions were not considered as network hubs based on local information, such as C D,node . In previous studies, these regions were classified as multimodal and functional hubs that are parts of cognitive resting-state networks, such as the default mode 5,25 . In addition, these regions were also defined as network hubs in other species, like in macaques and cats 14,22,53 . Some studies have found that the high FA values in the superior frontal gyrus were associated with post-traumatic stress disorder 54 , and exhibited decreased blood oxygen level-dependent activation of the superior frontal gyrus during a working memory task in individuals with schizotypal personality disorders 55 . FA plays an important role in the detection of network hub regions in global communication processes in the superior frontal and anterior cingulate regions 21 . As shown in Table 2, C D,node has lower values than C B,node and C P,node in the superior frontal and anterior cingulate regions.
Although similar patterns were observed in network hubs (Table 1 and Fig. 1) and bridges ( Fig. 3) based on C B and C P because they both used global information, and because their core paths exhibited similar patterns, some regions, such as the calcarine fissure and the middle temporal gyrus, could not be detected based on C B,node , which measures only the shortest path between network regions. Notably, the middle temporal gyrus is a meaningful network hub 25,56 . The association of the middle temporal gyrus is reduced on voxel-based DTI measures 57 , and network efficiency and centrality in the middle temporal gyrus have been shown to be disrupted in individuals with Alzheimer's disease 58,59 . The grey matter volume is reduced in the middle temporal gyrus in individuals with schizotypal personality disorders 60,61 . The bridges based on C P,edge (Table 3) also yielded more connections with the calcarine fissure and the middle temporal gyrus compared to the bridges based on C B,edge ( Table 4). The precuneus plays an important role in the brain network, thus suggesting that it has mutual connections with other areas 56,62 . Specifically, the precuneus was connected with parietal regions that were related to visuo-spatial information processing 63 . Both C B,node and C P,node could detect the precuneus as a network hub (Fig. 1B,C). However,   www.nature.com/scientificreports www.nature.com/scientificreports/ network hubs based on C P,node reflected the important network properties of the precuneus in a better manner compared to C B,node (Table 2). Additionally, network bridges based on C P,edge also included more connections with precuneus than C B,edge ( Table 3).

Hub regions F-value P-value
The process of competition to find the optimal paths-instead of the shortest paths-from the Physarum model required increased information flow. Accordingly, it would be helpful to enhance the flow information efficiently across different brain regions. Centrality measures identified based on the shortest path assumption  www.nature.com/scientificreports www.nature.com/scientificreports/ have been used in many brain network analyses, such as computer viruses, news, rumors, or infections 32-35 , but they have not been used in most real networks. In the brain network, there are some considerations against the shortest path assumption because it is difficult to elucidate the mechanism of an action potential that encodes the route and its destinations 52 . The shortest path assumption can also lead to nonresilient communication or loss of information 52,64 . Accordingly, the communication model with multiple connected paths is likely to be more appropriate, and can produce various alternative paths, which increase the efficiency, robustness, and resilience of the brain network 50,65 . The Physarum model has been suggested to combine the flux of tubular networks and competing edges through many possible paths. Therefore, we conclude that the Physarum model can improve the efficiency, robustness, and resilience of the brain network.
It is important to investigate the common features and variability of network centralities across subjects, and it is also critical to minimize the intrasubject variability 66,67 . The coefficient of variation (CV) is computed to describe the variation across subjects. As shown in Table 5, the CV of C P,node was significantly lower than that of C B,node . This indicated that the Physarum network hubs were more consistent throughout the dataset. Notably, network hubs based on C D,node yielded the lowest CV values compared to those based on centrality measures. This   www.nature.com/scientificreports www.nature.com/scientificreports/ is because C D,node is a relatively simple method, and local information is less variable than global information. Therefore, C P,node may capture the characteristics of local information to detect network hubs.
Many previous studies on brain network analyses have used various predefined atlases to define network nodes. The choice of the atlas defining the network nodes affects the network measures 66 . The spatial location of highly connected brain regions can be different depending on which atlas is used 66,68 . While the automated anatomical labeling (AAL) atlas was used in this study to compare and interpret the existing network hub and bridge results obtained in the previous studies 8,25,69 , some rigorous experiments with different atlases will be needed in future studies to compare the effects of atlas selection.
In this study, we illustrated a novel methodological framework for the identification of influential nodes and connections of the human brain network based on C P . This model has not been previously employed in a brain network. Comparison of the validation results between C P and other network centrality measurements indicated that C P contained local and global information. Additionally, this measure was not based on the assumption that the information flow of a network spread only through the shortest connections. Accordingly, C P could reduce the within-individual variation and detect some regions and connections that are related to post-traumatic stress disorder, schizotypal personality disorder, and Alzheimer's disease. Therefore, it would be helpful to apply this measure to individuals with neurological disorders that could provide biologically meaningful network results.

Methods
Subjects and data acquisition. This study used the Human Connectome Project (HCP, https://www. humanconnectome.org/) dataset and included 339 healthy participants (age: 28.2 ± 3.9 years, female: 159, male: 180). Their scans and data were released after they passed the HCP quality control and assurance standards 70 . Table 6 shows the details of these datasets. Data preprocessing. An automated processing-pipeline (CIVET) was used to process T 1 -weighted magnetic resonance (MR) images (http://mcin-cnim.ca/neuroimagingtechnologies/civet/) 71 . The T 1 -weighted MR images were first registered to ICBM152 T1 template in the Montreal Neurological Institute (MNI) space using an affine linear transformation 72 , and were then corrected for intensity nonuniformities owing to magnetic field inhomogeneities using an N3 algorithm 73 . After the removal of tissues unrelated to the brain matter, registered and corrected images were segmented into the white matter, grey matter, cerebrospinal fluid, and background, using an advanced neural-net classifier 71 .
Diffusion Tensor Imaging (DTI) datasets were managed using the FMRIB's software library (http://www.fmrib. ox.ac.uk/fsl). Motion artifacts and eddy current distortions were corrected by normalizing diffusion-weighted images to the baseline image using the affine registration in the FMRIB's linear image registration tool (FLIRT). A diffusion tensor matrix from the corrected diffusion-weighted images was generated based on a simple linear fitting algorithm, and the FA of each voxel was then calculated. DTI tractography was performed in the diffusion MR space using the FACT algorithm 74 , and was implemented using the Diffusion Toolkit (http://trackvis.org/) for the extraction of approximately 100,000 fibers from each subject. An angle of <45° between each fiber tracking step and a minimum/maximum path length of 20/200 mm were set as the terminating conditions. The classified white matter map masked the tractography results to eliminate false positives.

Construction of structural connectivity matrices.
It is important to define the basic elements of networks as edges. Because definitions and processes of constructing network nodes and edges have been described in detail previously, we explained them briefly as follows 69,75,76 (Fig. 5A).
Node definition. We used the AAL atlas 48 with the exception of the cerebellum and subcortical regions to segment the cortical regions into 78 areas, which represent the nodes of the network. Individual T 1 -weighted images were nonlinearly transformed to the ICBM 152 template, and the AAL atlas in the MNI space was then transformed to the T 1 native space using the inverse transformation parameters. Therefore, the individual AAL atlas was defined in the T 1 native space.
Edge definition. Tractography results were used to quantify edges between different AAL regions for individual networks. Individual T 1 -weighted images were coregistered to the baseline image using the affine registration in the FMRIB's FLIRT. Tractography results were transformed into the T 1 native space using the inverse transformation parameters. Fiber tractography results and the AAL atlas thus represented the same individual T 1 native space. Two nodes were considered to be structurally connected when at least three fiber tracts were present between these two nodes 69,75,77,78 . Accordingly, the edge was defined as the mean FA value along the fiber tracts 69,75 . Structural connectivity matrices were then constructed for each individual.  www.nature.com/scientificreports www.nature.com/scientificreports/ Physarum centrality. C P was calculated in two steps using an in-house software implemented in MATLAB (Version R2012b, Mathworks, Natick, MA, USA) (Fig. 5B). Based on the Physarum model, the optimal path was obtained within all pairs of network nodes, and C P was then calculated at each node 43,45 .
Physarum model for path finding. The basic concept underlying the Physarum model is that long and narrow tubes tend to weaken, and short and wide tubes strengthen with the positive feedback of flux in tubes during the competition process in the effort expended to identify the optimal paths. This concept assumes that short and wide tubes are the most effective for fluid transmission of information.
In a Physarum tubular network, each tube segment is denoted as the edge e ij , and its two ends are linked nodes i and j. If the flow along the tube is a Hagen−Poiseuille flow, the flux Q of each edge e ij can be defined as  where s and t are starting and ending nodes, respectively. Thus, the total flux in the brain from the starting to the ending nodes is a fixed constant I 0 during the path-finding process. Therefore, the pressure of each node and flux are calculated using Eqs (1 and 2), respectively. The flux can be updated according to the calculated pressure at all the nodes. The conductivities of the tubes are strengthened by large fluxes based on the positive feedback in the Physarum model, or are weakened by small fluxes when the lengths of the tubes are maintained fixed. The conductivity D ij is thus changed over time and is expressed as, ij ij ij γ = − Figure 5. Flowchart of measurement of Physarum centrality. The process for Physarum centrality (C P ) measurement was assessed in two steps. In step 1, the optimal path using the Physarum model was iteratively calculated within all pairs of network nodes, respectively. In step 2, C P was extracted in each node or edge based on the optimal path within all pairs of network nodes. www.nature.com/scientificreports www.nature.com/scientificreports/ where γ is a decay rate of the tube and f(Q) is usually a simply increasing function with f(0) = 0 43 . The tubes without flux are removed, and the pressure at each node is updated during iterations. This process is repeated until the optimal path is found, thus indicating that as the brain information flow through the path between two nodes increases, the importance of that route increases. Finally, unused paths are removed in the Physarum model. There is a negative correlation between the length of the path and the amount of flux through the path.
Centrality measure. The criticality (c) of each edge e ij is defined as, ij k ij k where Q ij k is the kth final flux through edge e ij , and k are the different path indices between all different pairs of nodes, and c ij indicates the sum of the flux through the edge e ij between all pairs of nodes i and j. Correspondingly, C P,node of node i is defined as, where c ij is the criticality (c) of each edge e ij . In addition, C P,node is defined as the sum of the criticality of each edge e ij attached to i. C P,edge was also calculated by the sum of flux (c) that passed through a given edge from the Physarum model.

Other centrality measures.
In this study, the values of C D,node and C B,node were compared with C P,node .
Equivalently, the value of C D,node 28 of node i is defined as, B node j i k jk jk , where ρ jk is the number of the shortest paths from node j to k, and ρ jk (i) is the number of the shortest paths between nodes j and k that pass through node i. Accordingly, C B,node (i) captures the influence of a node on the information flow between other nodes in the network. A node with a high degree of C B,node indicates increased interconnectivity with other regions in the network. Thus, the value of C B,edge was calculated as the number of the shortest paths that passes through a given edge instead of a node 33 . These measures were calculated using the Brain Connectivity Toolbox (http://www.brain-connectivity-toolbox.net).
Statistical analyses. The Jaccard index (J) values of each set of network hubs or bridges from different centrality measures Were calculated to show how they overlapped quantitatively 79 . J was defined as the ratio of the number of overlapping network hubs or bridges to the total number of these hubs (or bridges) based on any two centrality measures, where A and B are the sets of the hubs or bridges from each centrality measure, ∩ A B is the number of overlapping hubs or bridges, and ∪ A B is the total number of A and B hub or bridge sets. The value of J varies from zero (no overlap) to one (perfect overlap). Linear regression analyses was performed at all 78 network nodes to assess the relationship of three different pairs (C D,node vs. C B,node , C P,node vs. C D,node , and C P,node vs. C B,node ) in a continuous manner.
The Z-transform was applied on the centrality measures to ensure a fair comparison, and an ANOVA test was then performed to determine significant differences among centrality measures (C D,node , C B,node , and C P,node ) from the different concepts. All 78 network nodes were analyzed separately, and P values such that P < 0.05 were considered statistically significant with Bonferroni post-hoc correction. Intersubject variability, which assesses whether centrality measures could be reliably reproduced across all subjects, was quantified based on the coefficient of variation (CV), n n where X n is the normalized centrality value of a network node at the nth subject, N is the total number of datasets, and σ and μ denote the standard deviation and mean, respectively. Equivalently, CV quantifies the central tendency and variability of the samples. Therefore, a lower CV indicates lower intersubject variability and higher consistency across subjects in the group.