All experiments were conducted in accordance with the ethical standards of the Wake Forest University institutional review board and with the Helsinki Declaration of 1975. Functional brain networks of 5 healthy volunteers were constructed according to 
. For each subject, 120 fMRI full-brain volumes were acquired over approximately 5 minutes. Images were corrected for motion, normalized to the MNI (Montreal Neurological Institute) space, and re-sliced to 4×4×5 mm voxel size using SPM99 (Wellcome Trust Centre for Neuroimaging, Longdon, UK). From these volumes, one time series was extracted for each of the 15,996 voxels encompassing all of the gray matter of the cerebrum. Images were corrected for physiological noise by band-pass filtering to eliminate signal outside of the range of 0.009–0.08 Hz 
, and mean time courses from the entire brain, the deep white matter, and the ventricles were regressed from the filtered time series. In the past, the practice of global mean regression has been under scrutiny due to the propensity to produce artificial deactivations, particularly in the white matter and cerebrospinal fluid (CSF) 
. It is important to note, however, that failure to regress the mean signal will prevent detection of true deactivations that are known to occur in the brain. Additionally, the regions that are highly sensitive to these artifacts (white matter and CSF) are not considered in the present work. A full discussion on this topic can be found in 
The time series in each voxel was correlated with every other voxel using the Pearson's correlation coefficient. These correlation values were then represented in a correlation matrix summarizing the functional relationships between every pair of voxels. A threshold was applied to the correlation matrix, above which voxel pairs were said to be connected. This resulted in a binary adjacency matrix where 1 indicated the presence of a link and 0 indicated the absence. The threshold was defined such that the relationship between the number of nodes N
and average number of connections between nodes k
was consistent across subjects. Specifically, the ratio of log(N
) to log(k
) was the same across subjects 
. This threshold resulted in a link density of approximately 0.0015, where density is the ratio of the number of links present in the network to the number of possible links. This density is consistent with the size-density relationship of many self-organized networks described in 
. Moreover, links defined by this threshold represented correlations that are approximately 3 standard deviations above the mean. depicts the process of generating the functional brain networks.
Generating a functional brain network.
Each functional brain network was selectively attacked at the nodes with the highest centrality. In particular, the top 5% highest centrality nodes were removed from the network, along with any links directly connected to those regions. After the removal of the nodes, the respective centrality measure was recalculated and another set of top 5% nodes were identified. This process was repeated until all nodes in the network were removed. Four centrality metrics were utilized, namely, degree centrality, leverage centrality, eigenvector centrality, and betweenness centrality. Degree centrality defines highly central nodes to be those having a high number of links connected to that node. Leverage centrality relates the degree of a node to that of its immediate neighbors. In particular, nodes with higher degrees than their neighbors are considered highly central to their local neighborhood 
. Eigenvector centrality evaluates centrality based on the centrality of immediately connected neighbors, and therefore a node connected to nodes with high degree is highly central by association 
. Betweenness centrality defines the importance of a node by the number of shortest paths between pairs of nodes on which the node lies. In this way, high betweenness nodes facilitate the exchange of information along the most efficient trajectories 
. Formulations for these metrics are provided in 
. In addition to targeted attacks, we also conducted random attacks by iteratively removing 5% of nodes randomly at each step.
After attacking the networks, changes in the network structure were evaluated by assessing three network characteristics: local efficiency (Eloc
), global efficiency (Eglob
), and the size of the giant component (S
). Local and global efficiency are used to infer the efficacy of information exchange through a network by studying its topology 
. Local efficiency quantifies the extent to which nodes communicate with immediate neighbors and can be thought of as an indication of regional specificity. Global efficiency quantifies the extent to which nodes communicate with distant nodes, and indicates the efficacy of information exchange throughout the entire network. As nodes are removed, the network may fragment into isolated subgraphs. The size of the giant component is defined to be the largest connected subgraph, and may be used to indicate the extent of fragmentation.
The impact on dynamics was evaluated using two models. The first is an equation-based spreading activation model described in 
. This model injects signal into a network, and allows the signal to spread through links and decay according to model parameters. The equation governing the spread of activation is given in Equation 1
If N is the number of nodes in the network, St
is an N×1 vector describing the signal at time t, Et
is an N×1 vector containing the external signal injected at time t, γ is the relaxation rate of the signal (0≤γ≤1), α is the relative amount of activity that flows from a node to its neighbors per unit time (α>0), and R
is the N×N connectivity matrix. R
was constructed by eliminating all negative connections in the correlation matrix, setting the diagonal of the matrix to 0, and normalizing the matrix such that each column sums to 1. Therefore R
contained only weighted (normalized) positive connections from the original correlation matrix. External signal, E
, was only present at time t
0, where the 50 seed nodes were set to 1, and all other nodes were 0. The seed nodes for the external signal were randomly selected from the population of nodes that were not deleted. The equation was iterated for 100 time steps. This spreading activation model was tested on the original network and the networks with nodes removed, where 5% through 80% of the nodes were removed in increments of 5%. By examining the total activation in the system over the course of the simulation, we evaluated the impact of removal of highly central nodes on the ability of information to spread through the network. Here, total activation is defined to be the sum of activity values across all nodes in the network at a given time during the simulation. This procedure was performed on 5 subjects. Additionally, the impact of targeting low degree nodes was examined in a single subject in order to further investigate the findings in 
, where the targeted removal of low degree nodes had a greater impact on the dynamics of a network containing coupled oscillators than high degree nodes. For this experiment, we removed nodes that were the top 5% through 30% highest centrality nodes as well as the 5% through 30% lowest centrality nodes, in increments of 5%. Seed nodes were again randomly selected from the pool of remaining nodes in the networks.
Varying the ratio α/γ results in a phase change in the spreading activation model. When α/γ is small, the total activation in the system decays to zero over time (referred to as Phase I), but as α/γ increases, the system enters a regime where the activation builds exponentially in a small component of the system, referred to as Phase II 
. We chose α
1 and tuned gamma until the original networks exhibited Phase II behavior, resulting in α/γ
Changes in dynamics were also evaluated by embedding a coarser form of each network into an agent-based model called the agent-based brain-inspired model (ABBM) described in 
. An agent-based model is a collection of agents that interact with one another by following simple rules. The rules used here were inspired by the work of Stephen Wolfram 
, who has been a major contributor to the study of cellular automata. In this case, agents are represented by the nodes of the functional brain network, and links in the network represent communication pathways between agents. Each agent possesses a state, which can be either on or off, and may update its state based on the states of all connected neighbors by following one of Wolfram's Rules. Due to the computational demand of this model, these networks were constructed by parcellating the brain volume of each subject into 90 anatomical regions using the AAL (automated anatomic labeling) atlas 
. The time series of all voxels belonging to a particular ROI were averaged in order to create 90 ROI time series. These time series were cross-correlated to construct a 90×90 ROI correlation matrix containing positive and negative connection weights. A threshold was applied to these networks to preserve only strong positive or negative connections while preventing fragmentation. Therefore, positive and negative weighted links were present in the ROI networks. The process of creating the ROI networks and the mechanisms underlying the ABBM are described in full in a prior publication 
These ROI networks were selectively attacked by removing 10% of the nodes (9 regions) with the highest centrality, at random, or with the lowest centrality. Slight modifications to the centrality metrics were necessary in order to calculate these metrics in the weighted, signed correlation matrix. Degree was calculated as the sum of the absolute value of the weights of all links belonging to a node. Leverage and eigenvector centrality, which depend only on the degree of the node and its connected neighbors, were calculated using this definition of degree. The weighted form of betweenness was calculated on the absolute value of the correlation matrix using the MATLAB BGL package (http://dgleich.github.com/matlab-bgl/
The impact on dynamics was evaluated by testing the ability of the attacked agent-based model to solve the density classification problem, a problem originally utilized to evaluate whether a one-dimensional cellular automaton (CA) could support computation 
. A CA can be thought of as belonging to a class of agent-based models, where agents are spatially embedded as adjacent cells. The goal of the density classification problem is to find a rule that can determine whether greater than half of the cells in a CA are initially in the on state. If the majority of cells are on (i.e. density >50%), then by the final iteration of the CA, all cells should be in the on state. Otherwise, all cells should be turned off. The system should be able to do this from any random initial configuration of node states. The key is that each node receives input from only a few other nodes in the network. Each node must decide based on this limited information whether to turn on or off in the next time step, resulting in network-wide cooperation without the luxury of network-wide communication. The rule and model parameters that must be used in order to perform this task are identified using a search optimization technique known as genetic algorithms. We have demonstrated that the ABBM is able to perform the density classification task with a high level of accuracy across a range of densities, while null models with randomized connectivity are not successful, indicating that the topology of the brain network is amenable to computation. Here we wished to determine how targeted removal of high centrality nodes would impact performance on this task.
contains a summary of treatments of the functional brain networks used in each procedure for evaluating network structure and dynamics.
Summary of networks used to evaluate network topology and dynamics.