Functional control of electrophysiological network architecture using direct neurostimulation in humans

Novel neurotechnologies capable of perturbing the physiological state of neural systems are rapidly gaining popularity for their potential to treat neurological disease and psychiatric disorders (Stacey & Litt, 2008). Chronically implantable devices that stimulate the human brain are clinically approved to treat Parkinson’s disease, essential tremor, dystonia, epilepsy, and obsessive-compulsive disorder and have been investigated for major depressive disorder and Tourette syndrome (Lozano & Lipsman, 2013). Recent human studies have investigated the ability for direct stimulation of cortical and subcortical structures to modulate biomarkers of memory (Ezzyat et al., 2017; Inman et al., 2017), visual perception (Rangarajan et al., 2014; Winawer & Parvizi, 2016), language production (Chang, Kurteff, & Wilson, 2017), somatosensory perception (Muller et al., 2018), sensorimotor function (W. Wang et al., 2013), and subjective experience (Foster & Parvizi, 2017). While neurostimulation is a promising interventional approach to modulate brain state, current practices of calibrating where, when, and how to stimulate the brain are “open-loop” and limited in efficacy—relying on manual and periodic tuning of device parameters to optimize therapy (Morrell, 2011). Automated, “closed-loop” approaches would augment the capability of current stimulation devices to dynamically adjust parameters based on the physiological state of the brain network, monitored in real time (Stanslaski et al., 2012). Undoubtedly, the translational prospect of neurostimulation to manipulate brain networks that generate abnormal rhythms, dysrhythmias, or bursts of activity associated with dysfunction is promising. However, critical gaps in knowledge hinder the development of a robust control policy for next-generation implantable devices.

How does the architecture of the neural system mediate the effect of neurostimulation on neurophysiology and behavior? Network control theory (Pasqualetti, Zampieri, & Bullo, 2014) provides a mathematical framework for mapping the influence of a control signal on the dynamics of an interconnected system. When combined with graph modeling tools from network neuroscience (Bassett & Sporns, 2017), where nodes represent discrete brain regions and edges represent the structural connections between brain regions, control theoretic approaches can elucidate how the brain’s structural architecture of white matter fiber pathways shapes its ability to navigate through a repertoire of dynamical states (Gu et al., 2015). Theoretical rules of controllability prescribe the trajectories through state space elicited by a given control signal (Betzel, Gu, Medaglia, Pasqualetti, & Bassett, 2016; Gu et al., 2017, 2015), and begin to explain why one brain network may be more or less influential on brain dynamics than another (Kim et al., 2018). Recent efforts to test control theoretic predictions of the relationship between controllability and brain activity have relied on in silico models in which neuronal ensembles are interlinked by structural connections measured by human neuroimaging (Muldoon et al., 2016). Despite the promising convergence between theory and model simulation, empirical stimulation data bridging network control and neurophysiology are lacking.

Network control theory accounts for the structural connections that convey modulated brain activity to downstream regions in the network; however, it does not account for the functional rules that govern whether communication between brain regions can occur at a specific point in time. At the millimeter scale, synchronous oscillations in the local field potential are thought to actively gate the transfer of information across the network (Bonnefond, Kastner, & Jensen, 2017; Buzsáki et al., 2012; Canolty & Knight, 2010; Fries, 2015; Schalk, 2015) and are commonly observed during higher order cognitive processing (Buzsáki, 2006). A functional (rather than structural) network representation of the coherence between different ensembles of neurons may capture dynamical states of communication (Deco & Kringelbach, 2016; Hutchison et al., 2013). The neurophysiologic interpretation of these states can depend on the measured frequency range of the functional network (Bassett, Meyer-Lindenberg, Achard, Duke, & Bullmore, 2006; Solomon et al., 2017), which in turn implicates certain types of cells interacting over specific spatial scales (Kopell, Ermentrout, Whittington, & Traub, 2000). Prior studies have examined how these functional networks may reconfigure during higher order cognitive functions such as learning new skills (Bassett, Wymbs et al., 2011; Bassett et al., 2013; Bassett, Yang, Wymbs, & Grafton, 2015; Mattar, Cole, Thompson-Schill, & Bassett, 2015), forming memories (Braun et al., 2015), attending to the environment (Shine, Koyejo, & Poldrack, 2016), and processing language (Chai, Mattar, Blank, Fedorenko, & Bassett, 2016). Complimentary work also posits that reconfiguration of functional networks may underlie neurophysiological abnormalities in patients with epilepsy (Khambhati et al., 2017, 2015), schizophrenia (Bassett, Nelson, Mueller, Camchong, & Lim, 2012; Braun et al., 2016), Parkinson’s disease (Olde Dubbelink et al., 2014; Sang et al., 2015), and stroke (Grefkes & Fink, 2011; L. Wang et al., 2010). While these studies explain changes in functional network reconfiguration when the brain is perturbed en masse, a rigorously quantified map of functional network reconfiguration due to controlled, focal perturbation has not been attained.

Here we seek to elucidate the network control principles by which neurostimulation can alter function and behavior based on constraints prescribed by structural connectivity and spontaneous functional interactions. We measure the electrocorticogram (ECoG) in 94 drug-resistant epilepsy patients undergoing neurostimulation (Figure 1A–B), and we construct structural networks from diffusion imaging data acquired in the same individuals. We also construct functional networks before and after individual stimulation trials using multitaper coherence between sensors (Prieto, Parker, & Vernon, 2009) in distinct frequency bands (Khambhati, Davis, Lucas, Litt, & Bassett, 2016; Kramer et al., 2011) (Figure 1C), and we define brain state before and after stimulation using a previously validated biomarker of memory (Ezzyat et al., 2017). We test four hypotheses. First, we hypothesize that the strength and location of stimulation can differentially drive two separate modes of global versus local control over functional architecture (Muldoon et al., 2016). Intuitively, stimulation to functional hubs—nodes that tend to interact strongly with the rest of the network—may have swiftly attenuated effects due to signal dispersion across many downstream regions, while stimulation to nonhubs may have more localized and targeted effects. Second, we hypothesize that regions with strong baseline functional interaction with the stimulation site are more likely to exhibit altered hub properties following stimulation than brain regions with weak functional interaction with the stimulation site, indicating a functional conduit of stimulation. Third, based on prior data (Betzel, Medaglia, Kahn, et al., 2017), we hypothesize that these functional interactions—particularly in high-frequency bands—colocalize with structural white matter networks (Figure 1D). Fourth, we hypothesize that neurostimulation directed towards modal control points (Gu et al., 2015; Pasqualetti et al., 2014), which tend to be structural nonhubs of a patient’s white matter network thereby minimizing signal dispersion, facilitate a stronger shift in dynamical state associated with memory encoding, a function that is altered in patients with epilepsy (Aarts, Binnie, Smit, & Wilkins, 1984; Holmes & Lenck-Santini, 2006; Uhlhaas & Singer, 2006). Collectively, these analyses will supply a roadmap of the impact of neurostimulation on network physiology, mediated by network structure, and provide fundamental mechanistic insight into the influence of neurostimulation on behavioral state.

Ninety-four patients undergoing intracranial EEG monitoring as part of clinical treatment for drug-resistant epilepsy were included in this study. Data were collected as part of a multicenter project designed to assess the effects of electrical stimulation on memory-related brain function. Data analyzed in this study were collected at the following centers: Thomas Jefferson University Hospital (N = 23), University of Texas Southwestern (N = 23), Mayo Clinic (N = 17), National Institutes of Health (N = 11), Dartmouth-Hitchcock Medical Center (N = 9), Hospital of the University of Pennsylvania (N = 6), Columbia University Medical Center (N = 4), and Emory University Hospital (N = 1). The research protocol was approved by the institutional review board at each hospital and informed consent was obtained from each participant.

Patients undergoing surgical treatment for medically refractory epilepsy believed to be of neocortical origin underwent implantation of intracranial electrodes to localize the seizure onset zone. These procedures were applied after presurgical evaluation with scalp EEG recording of ictal epochs, MRI, PET, and neuropsychological testing suggested that focal cortical resection may be a therapeutic option. Patients were then deemed candidates for implantation of intracranial electrodes to better define epileptic networks. Electrode configurations spanned the surface of the cortex (linear and two-dimensional arrays, each sensor is 2.3 mm diameter spaced 10 mm apart) and subcortical depth (each sensor is 1.5–10 mm apart). All electrode configurations were planned by a multidisciplinary team of neurologists and neurosurgeons at each of the eight medical centers.

Electrodes were anatomically localized using separate processing pipelines for surface and depth electrodes. To localize depth electrodes we first labeled hippocampal subfields and medial temporal lobe cortices in a pre-implant, 2 mm thick, coronal T2-weighted MRI using the automatic segmentation of hippocampal subfields (ASHS) multiatlas segmentation method (Yushkevich et al., 2015). We additionally used whole-brain segmentation to localize depth electrodes not in medial temporal lobe cortices. We next coregistered a postimplant CT with the preimplant MRI using advanced normalization tools (ANTs; Avants, Epstein, Grossman, & Gee, 2009). Electrodes visible in the CT were then localized within subregions of the medial temporal lobe by a pair of neuroradiologists with expertise in medial temporal lobe anatomy. The neuroradiologists performed quality checks on the output of the ASHS/ANTs pipeline. To localize subdural electrodes, we first extracted the cortical surface from a pre-implant, volumetric, T1-weighted MRI using FreeSurfer (Fischl et al., 2004). We next coregistered and localized subdural electrodes to cortical regions using an energy minimization algorithm (Dykstra et al., 2012). For patient imaging in which automatic localization failed, the neuroradiologists performed manual localization of the electrodes.

The electrocorticogram (ECoG) was recorded and digitized at 500 Hz, 512 Hz, 1,000 Hz, 1,024 Hz, or 2,000 Hz depending on clinical considerations at each medical center. Signals were recorded using a referential montage with the reference electrode, chosen by the clinical team, distant to the site of seizure onset.

To study the response of the electrocorticogram to neurostimulation, we used a mapping procedure in which stimulation was delivered to cortical and subcortical brain regions. Patients were not instructed to engage in any other task before or during stimulation. Prior to the start of each mapping session, we selected a pair of adjacent electrodes for stimulation by prioritizing electrodes in brain regions thought to be associated with memory function. For each mapping session, we selected a new stimulation site and patients underwent one or several mapping sessions depending on their availability for testing and the monitoring needs of the clinicians. Prior to the start of a mapping session, we recorded 30 s of ECoG activity as a baseline epoch. During a mapping session, we performed several stimulation trials in which a single trial consisted of the following epochs: (a) a half-second pre-stimulation epoch, (b) a stimulation epoch with variable duration, two consecutive and nonoverlapping half-second post-stimulation epochs, and an inter-stimulation epoch with variable duration. During each stimulation trial, we delivered stimulation using charge-balanced, biphasic, rectangular pulses with a pulse width of 600 μs and combinations of the following parameters: pulse frequency (10, 25, 50, 100, 200 Hz), pulse amplitude (maximum safe amplitude minus 0, 0.5, 1 mA; range of 0.125–3.0 mA across subjects), stimulation duration (250, 500, 1,000 ms), and inter-stimulation interval (2,750–3,250 ms). These stimulation parameter ranges were chosen to be well below the accepted safety limits for charge density (Shannon, 1992) and ECoG was continuously monitored for afterdischarges by a trained neurologist.

We will now describe a typical stimulation experimental session in further detail. In a single session, stimulation location was kept constant across trials, and in each trial stimulation parameters (pulse frequency, pulse amplitude, and duration) were each separately drawn uniformly at random from the aforementioned list of parameters. This procedure produced 45 possible parameter combinations. The sampling distribution of parameter combinations was consistent across patients; however, the number of trials in a given session differed based on a number of factors in the hospital setting of the Epilepsy Monitoring Unit, including patient fatigue, availability, and willingness to participate. In general, across 94 participants we conducted an average of 3 ± 2 experimental sessions (unique stimulation locations), and per session we conducted an average of 1,655 ± 697 trials—more than 36 times the number of possible parameter combinations. We confirmed that the full parameter space was sampled within a session.

To eliminate confounding effects of stimulation on signal quality and saturation, we disregarded ECoG data collected during the stimulation epoch and the 100 ms following stimulation offset. We also employ a conservative electrode screening procedure, in which we discard nonstimulated channels that exhibit evidence of stimulation-related artifact. Specifically, before re-referencing to a common average reference, we use a paired t test to compare the distribution of mean signal amplitude during the pre-stimulation epoch to the distribution of mean signal amplitude during the post-stimulation epoch, for each electrode across stimulation trials. Using a Bonferroni uncorrected p-value threshold of 0.05, we discard electrodes that exhibit significantly elevated raw, mean signal amplitude during each stimulation session.

We analyzed ECoG data collected during the baseline, pre-stimulation, and post-stimulation epochs. The post-stimulation epoch following 100 ms of a buffer period was split into two consecutive and nonoverlapping segments, 0.5 s in duration to assess delayed effects of stimulation—we refer to the first segment as the 100-ms response and the second segment as the 600-ms response. Because of the time constraints on experimentation in the hospital setting of the Epilepsy Monitoring Unit, we utilized the 30 s of ECoG activity recorded before the start of a stimulation session as an indicator of the patient’s baseline state. To quantify the difference between the stimulation-induced effect on network topology in contrast to the spontaneous effects over the passage of time, we constructed a surrogate distribution of stimulation “trials” from the baseline period by sampling-with-replacement time windows equal in duration to the true stimulation trials. In other words, the surrogate distribution consisted of the same number of trials as the number of trials conducted during the stimulation session, and each surrogate baseline trial was associated with a true stimulation trial and was constructed from the same duration of baseline data as its corresponding stimulation trial. This sampling scheme was used to mitigate the limited baseline data that were available for analysis.

ECoG signals were divided into 0.5s, nonoverlapping, time windows—the pre-stimulation epoch consisted of one time window per stimulation trial, and the post-stimulation epoch consisted of two time windows per stimulation trial (100–600 ms post-stimulation and 600–1,100 ms post-stimulation). We applied a common average reference to the artifact-free ECoG signal before constructing functional networks (Burns, Santaniello, Yaffe, Jouny, & Crone, 2014; Khambhati et al., 2016; Kramer et al., 2010, 2011; Towle, Carder, Khorasani, & Lindberg, 1999).

To measure functional interactions between ECoG signals in each time window, we computed spectral coherence, which is a measure of correlation between the power spectra of two signals within a frequency range. Prior studies have shown that coherence is largely independent of the shape of the power spectrum in ECoG signals (Bullock et al., 1995a, 1995b; Towle et al., 1999), and underlies different forms of synchronous interactions between neural populations (Kopell et al., 2000). We constructed functional networks in each time window using multitaper coherence estimation, which defines a graph edge between electrode pairs (graph nodes) as the power spectral similarity of signal activity over a specific frequency band. We applied the mtspec Python implementation (Prieto et al., 2009) of multitaper coherence estimation with time-bandwidth product of five and eight tapers in accord with related studies (Kramer et al., 2011). This procedure resulted in a symmetric adjacency matrix A(t, f) with size N × N, where N is the number of network nodes, or electrode sensors, t is the time window, and f is the frequency band. In this study, we examined network activity in the following four frequency bands: α/θ (5–15 Hz), β (15–25 Hz), low γ (30–40 Hz), and high γ (95–105 Hz). These frequency ranges cover traditional oscillatory classes and have been previously examined for their network topology (Khambhati et al., 2016; Kramer et al., 2011).

An alternate representation of the symmetric, square adjacency matrix A(t, f) is a configuration vector $Â$(t, f), which tabulates all N × N pairwise interactions. Because of symmetry of the adjacency matrix, we unravel the upper triangle of A, resulting in the weights of E = $N(N−1)2$ functional interactions. Thus, $Â$(t, f) is a vector of size E.

In this study, we investigated the effect of neurostimulation on functional network architecture at the scale of network nodes and at the scale of network edges. At the node scale, we first quantified the change in the node strength—a measure of functional “hubness”—of individual network nodes. Specifically, we computed the node strength as k_i(t,f) = $1N−1$ ∑_j∈NA_ij(t,f), where k is the strength of node i and A_ij is the edge weight between nodes i and j. Based on the time-dependent set of node strengths in the network, we computed the change in the mean of node strengths between time windows and the change in the variance of node strengths between time windows. To assess the magnitude of change in node strength for a node between time windows t_n and t_m, we calculated Δk_i(t_n,m, f) = abs(k_i(t_m, f) − k_i(t_n, f)).

At the edge scale, we quantified the amount of change in the configurational pattern of the network edges, or coherences, as described previously in Khambhati et al. (2015). Specifically, we computed the configuration similarity between configuration vectors $Â$(t_n, f) and $Â$(t_m, f), where t_n and t_m are two different time windows, using the Pearson correlation test statistic. Two vectors with a Pearson correlation value closer to 0 are more dissimilar in their configurational pattern of network edges than two vectors with a Pearson correlation value closer to 1.

For the stimulation epoch, we computed global and local network metrics between the pre-stimulation time window and the post-stimulation time window of a stimulation trial. For the baseline epoch, we computed global and local network metrics between time windows separated by an equal length of time as the duration of stimulations in the associated stimulation session.

We collected diffusion-weighted imaging data for a subset of patients from Thomas Jefferson University Hospital (N = 11) and Hospital of the University of Pennsylvania (N = 3) and validated our analysis of the functional network response to neurostimulation.

All scans at Thomas Jefferson University Hospital were acquired with a 3T Philips Achieva with an 8-channel head coil using an echo-planar diffusion-weighted technique. The diffusion scan was 62-directional with a b-value of 3,000 s/mm² and TE/TR = 98/7,251 ms. The matrix size was 96 × 96 and the slice number was 52. The field of view was 230 × 230 mm² and the slice thickness 2.5 mm. Acquisition time was 496 s per DTI scan.

All scans at the Hospital of the University of Pennsylvania were acquired with a 3T Siemens Tim Trio with a 32-channel head coil using an echo-planar diffusion-weighted technique. The diffusion scan was 116-directional with a b-value of 2,000 s/mm² and TE/TR = 117/4,180 ms. The matrix size was 96 × 96 and the slice number was 92. The field of view was 210 × 210 mm² and the slice thickness 1.5 mm. Acquisition time was 506 s per DTI scan.

Based on recent evidence that diffusion imaging is highly sensitive to subject movement (Yendiki, Koldewyn, Kakunoori, Kanwisher, & Fischl, 2014) and to directional eddy currents (Jezzard, Barnett, & Pierpaoli, 1998), we processed data using the FMRIB Software Library (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012). We first created individual masks of the patient brain using BET (Smith, 2002). We next simultaneously corrected for motion effects and eddy current distortions by applying the EDDY correction tool (Andersson & Sotiropoulos, 2016) to the diffusion scans and a b = 0 image collected at the beginning of the scan.

We next reconstructed orientation density functions (ODFs) of the diffusion imaging in each voxel. Specifically, we used DSI Studio (http://www.dsi-studio.labsolver.org) and generalized q-sampling imaging (GQI; Yeh, Wedeen, & Tseng, 2010) to compute the quantitative anisotropy (QA; Yeh, Verstynen, Wang, Fernández-Miranda, & Tseng, 2013) in each voxel. To conduct fiber tractography on the reconstructed diffusion images, we used DSI Studio to generate 1,000,000 streamlines with a maximum turning angle of 35° (Bassett, Brown, Deshpande, Carlson, & Grafton, 2011) and a maximum length of 500 mm (Cieslak & Grafton, 2014). We next used the streamlines to define the structural brain network at the two following spatial resolutions of the Lausanne atlas included in the Connectome Mapping Toolkit (Cammoun et al., 2012) consistent with previous work (Bassett, Brown, et al., 2011; Bassett et al., 2010; Gu et al., 2017, 2015; Hermundstad et al., 2013, 2014; Muldoon et al., 2016): (a) at the fine scale we divided the brain into N = 1,015 cortical and subcortical regions of interest (ROIs) (average ROI volume 1.56 ± 0.06 cm³), and (b) at the coarse scale we divided the brain into N = 234 cortical and subcortical ROIs (average ROI volume 6.78 ± 0.02 cm³). We summarized these measurements in a symmetric and weighted structural adjacency matrix S whose entries S_ij reflect the structural connectivity (quantitative anisotropy) between region i and region j, separately for the fine-scale and coarse-scale parcellation.

We localized electrodes in native subject T1-weighted MRI space to the Lausanne anatomical space by using ANTs (Avants et al., 2009) to register the subject’s T1 image to the subject’s diffusion 0 image via affine transformation and also to register the subject’s T1 image to MNI space (also native Lausanne space) using a nonlinear warp.

To study the architectural constraints of the structural brain network with the functional network response to neurostimulation, we adopted a control theoretic approach known as network controllability. Briefly, the controllability of a networked system refers to its ability to be driven to specific dynamical states upon external input (Kalman, 1963). Recent research efforts have made substantial progress in the development of quantitative heuristics to characterize different strategies for control (Pasqualetti et al., 2014; Pequito et al., 2016). These approaches are now being applied to brain imaging data to understand how structural brain network topology constrains function and behavior (Betzel et al., 2016; Gu et al., 2017, 2015; Kim et al., 2018; Tang et al., 2017).

One control strategy that we investigate in this study is modal controllability—the ability of a network region to feasibly control all the dynamical modes of a system (Pasqualetti et al., 2014). To calculate the modal controllability of an anatomical brain region, we first computed the eigenvector matrix V = [v_ij] of the structural network adjacency matrix S—intuitively, v_ij encodes the ability to control the j-th dynamical mode from region i (Kailath, 1980). Based on our previous work, we defined ϕ_i = ∑_j∈N(1 − $λj2$(S))$vij2$ as a scaled measure of the controllability of all N dynamical modes λ₁(S), …, λ_N(S) from brain region i (Gu et al., 2015; Muldoon et al., 2016; Pasqualetti et al., 2014; Tang et al., 2017). Brain regions with high modal controllability are versatile in their ability to control all dynamical modes of the network, and brain regions with low modal controllability are specific in their ability to control a subset of dynamical modes of the network.

To provide additional insight into the topological properties of structural control points, we evaluated the structural “hubness” of each brain region by computing the structural node strength as k_i = $1N−1$ ∑_j∈NS_ij—similar to the calculation for functional node strength specified earlier.

To relate structural controllability to functional network topology of the stimulated electrodes, we first computed metrics of the structural network topology for 234 or 1015 brain regions defined by the Lausanne anatomical parcellation. The advantage of computing these measures using the anatomical parcellation is the ability to account for whole-brain structural connectivity, including areas that are not directly sampled by the intracranial electrodes. We next assigned intracranial electrodes to the Lausanne brain regions based on a nearest voxel approach. Specifically, we identified the voxel closest to the electrode and assigned the electrode to the brain region containing that voxel. Based on this assignment, we associated values of each structural network metric to the intracranial electrodes.

We examined stimulation-driven changes in dynamical brain state using a classifier of neural activity associated with memory encoding processes that was previously validated on data collected during behavioral experimentation with the same patients recruited in this study (Ezzyat et al., 2017; Kragel et al., 2017). Briefly, in these prior studies a logistic regression classifier was trained to discriminate memory encoding-related changes in spectral power in eight logarithmically spaced frequency bands across intracranial electrodes that are predictive of whether a word was later remembered or forgotten during a free-recall task (Ezzyat et al., 2017; Kragel et al., 2017). In this study, we evaluated the trained memory encoding state classifier on task-free stimulation data of the same patients by measuring spectral power during the pre-stimulation epoch and the post-stimulation epoch, and by computing the change in probability of good memory encoding state for each stimulation trial. We next calculated the average change in probability of good memory encoding state across all stimulation trials of each stimulation mapping session of each patient. These data allowed us to assess the putative effects of different dimensions of the stimulation parameter space on previously validated, neural biomarkers of cognitive function, specifically on memory encoding.

We first ask the question, “How does neurostimulation alter the architecture of functional brain networks?” Based on recent theoretical insights on the costs of forming and breaking connections in structural and functional brain networks (Achard & Bullmore, 2007; Betzel, Medaglia, Papadopoulos, et al., 2017; Bullmore & Sporns, 2012), we expect stimulation to heterogeneously affect existing coherent interactions, strengthening some and weakening others. To test these expectations, we study the average change in the following three measures of network topology across stimulation trials per patient: two at the topological scale of nodes and one at the topological scale of edges (Figure 2A). At the node scale, we first compute the strength, or average coherence, for each network node during the pre-stim epoch and post-stim epoch for each of the four coherence frequency bands. We next examine the change in the mean of node strengths and the change in the variance of node strengths between the pre-stim epoch and the post-stim epoch. Intuitively, a change in the mean of node strengths quantifies the likelihood that nodes exhibit greater frequency-specific functional interaction following stimulation, and a change in the variance of node strengths quantifies the likelihood that nodes exhibit greater heterogeneity in their degree of functional interaction with other nodes in the network. At the edge scale, we compute the configuration similarity (Khambhati et al., 2015): a Pearson correlation between the vector of coherence weights during the pre-stim epoch and the vector of coherence weights during the post-stim epoch. Similarity values near 0 imply a greater change in the configuration of network coherences, and values near 1 imply a lesser change in the configuration of network coherences. We confirmed that the topological measurements at the node scale capture different reconfiguration phenomenon than the topological measurements at the edge scale by observing weak relationships between changes in the mean and variance of node strengths to configuration similarity (Supporting Information: Figure 1-1).

Next, we test our expectation that stimulation heterogeneously affects existing coherent interactions, strengthening some and weakening others, by comparing the change in network topology measures during stimulation trials to passive changes in network topology that occur during the baseline period before the stimulation session begins. We first study changes in the mean and variance of node strengths between the pre-stim epoch and post-stim epoch for each of the four coherence frequency bands (Figure 2B, C). Using a Wilcoxon rank-sum test, Bonferroni correction for multiple comparisons, and effect size based on rank correlation, we examine whether node-level changes in the network 100 ms after stimulation offset are any greater than passive changes observed over an equal duration of spontaneous activity at baseline, before any stimulation, across stimulation sessions over subjects. We find that stimulation leads to a significantly greater change in the mean of node strengths than expected at baseline in the alpha/theta band (Z(247) = 11,448, p = 0.001, rank correlation = 0.63) and to a nonsignificant change in the beta band (Z(247) = 14,006, p = 0.68), in the low gamma band (Z(247) = 13,131, p = 0.13), and in the high gamma band (Z(247) = 14,883, p = 1.0). We also find that stimulation leads to a significantly greater change in the variance of node strengths than expected at baseline in the alpha/theta band (Z(247) = 10,660, p = 0.0001, rank correlation = 0.65) and to a non-significant change in the beta band (Z(247) = 13,422, p = 0.24), in the low gamma band (Z(247) = 14,137, p = 0.84), and in the high gamma band (Z(247) = 14,086, p = 0.76). We find that these effects indeed persist and possibly strengthen in the beta band and low gamma band at least 600 ms after stimulation offset (Supporting Information: Figure 2-1b, c). These results demonstrate that stimulation amenably alters functional network organization in lower alpha/theta band frequencies (5–15 Hz) at the node scale. Specifically, we observe that nodes generally exhibit an increase in low-frequency interaction following neurostimulation. However, changes in node strengths are also heterogeneously distributed across nodes in the network.

We next ask whether stimulation may still alter functional network topology at the edge scale. Using a Wilcoxon rank-sum test and Bonferroni correction for multiple comparisons, we examine whether configurational changes in the network edges 100 ms after stimulation offset are any greater than the passive change observed over an equal duration of spontaneous activity at baseline, before any stimulation, across stimulation sessions over subjects. We find that stimulation leads to a significantly lower configuration similarity (greater reconfiguration) than expected at baseline in the high gamma band (Z(247) = 9,252, p = 2.3 × 10⁻⁶, rank correlation = 0.70) and to a nonsignificant change in the alpha/theta band (Z(247) = 13,543, p = 1.0), in the beta band (Z(247) = 12,820, p = 1.0), and in the low gamma band (Z(247) = 13,502, p = 1.0). We find that these effects indeed persist at least 600 ms after stimulation offset (Supporting Information: Figure 2-1d). These results demonstrate that stimulation amenably alters functional network organization in high gamma band frequencies (95–105 Hz) at the edge scale. Specifically, we observe that functional interactions undergo a change in their configurational pattern in high frequencies following neurostimulation.

Building on our observations of a complex, frequency-dependent network response to stimulation, we next ask, “Do properties of the stimulation signal, such as amplitude, pulse frequency, and duration, influence the extent of functional network reconfiguration?” There is a fundamental gap in knowledge of how different parameters of direct brain stimulation influence brain networks—delaying the therapeutic benefits of stimulation in the treatment of neurological and neuropsychiatric disorders by several months or years (Crowell, Garlow, Riva-Posse, & Mayberg, 2015; Fisher, 2011; Nair & Morrell, 2017). To understand how stimulation parameters influence functional network reconfiguration, we draw stimulation parameters from a predefined list, uniformly at random, for each consecutive trial (Figure 3A), and we compute the stimulation intensity as the product between the three parameters (Figure 3B). Based on prior observations of a relationship between stimulation intensity and volume of tissue activated (Butson & McIntyre, 2008), we hypothesize that stronger stimulation input into the functional network will lead to more widespread change in functional architecture than weaker stimulation input, presumably by penetrating the network along short axonal fibers in the gray matter and long myelinated fibers in the white matter.

To test this hypothesis, we first compute a within-session Spearman’s ρ correlation between stimulation intensity and the three measures of functional network reconfiguration (change in mean of node strengths, change in variance of node strengths, and configuration similarity) for the four coherence frequency bands (Figure 3C–E). Using a one-sample t test, Bonferroni correction for multiple comparisons, and effect size based on Cohen’s d we test whether increasing stimulation intensity drives greater node-level changes in the network 100 ms after stimulation offset (Figure 3C, D). We find that greater stimulation intensity leads to a significant decrease in the mean of node strengths in the high gamma band (t(247) = −6.5, p = 7.0 × 10⁻⁹, Cohen’s d = 0.41), and to a nonsignificant change in the alpha/theta band (t(247) = 1.6, p = 1.0), beta band (t(247) = 0.2, p = 1.0), and low gamma band (t(247) = −2.0, p = 0.53). We also find that greater stimulation intensity leads to a significant decrease in the variance of node strengths in the high gamma band (t(247) = −5.9, p = 7.5 × 10⁻⁵, Cohen’s d = 0.40), and to a nonsignificant change in the alpha/theta band (t(247) = 1.2, p = 1.0), in the beta band (t(247) = −0.7, p = 1.0), and in the low gamma band (t(247) = −1.6, p = 1.0). We find that these effects indeed persist in the high gamma band at least 600 ms after stimulation offset (Supporting Information: Figure 3-1c, d). Our results indicate a robust dependence of high-frequency functional reorganization at the scale of network nodes on stimulation strength. Specifically, greater stimulation intensity disrupts and decreases cohesive node-level interactions in high-frequency bands. We did not observe a similar disruption in node-level architecture in the lower frequency bands.

Logically, we next ask whether stimulation intensity similarly alters the edge-level architecture of the network. Using a one-sample t test and Bonferroni correction for multiple comparisons, we test whether increasing stimulation intensity drives greater configurational change in the network edges 100 ms after stimulation offset (Figure 3E). We find that greater stimulation intensity leads to a significant decrease in the configuration similarity (greater reconfiguration) in the high gamma band (t(247) = −2.8, p = 0.04, Cohen’s d = 0.18), and to a nonsignificant change in the alpha/theta band (t(247) = −2.8, p = 0.06), in the beta band (t(247) = −1.6, p = 1.0), and in the low gamma band (t(247) = −1.2, p = 1.0). We find that these effects dissipate 600 ms after stimulation offset (Supporting Information: Figure 3-1e). Our results indicate that greater stimulation intensity drives greater reconfiguration of the functional topology in high-frequency bands. Additionally, stimulation strength only explains immediate edge-level reconfiguration of network topology and does not exhibit a relationship with later stage edge-level reconfiguration.

Lastly, we asked whether the observed changes in network topology were primarily driven by any single dimension of the stimulation parameter space (trial duration, pulse amplitude, pulse frequency). In contrast to stimulation duration and stimulation amplitude, we observed that changes in the frequency of stimulation significantly drives altered network topology at the node level and at the edge level (Supporting Information: Figure 3-2g–i). Specifically, faster stimulation frequencies may disrupt high-frequency coherence between network nodes, by presumably redistributing coherent edges across the network and reducing the variance in node strengths. Conversely, slower stimulation frequencies may increase node strengths in high-frequency networks by driving less topological reconfiguration of the network edges and simply reinforcing existing functional interactions.

We next build upon our analysis of the influence of stimulation parameters on functional network topology by similarly investigating the role of 248 unique stimulation locations over 94 subjects in the functional brain network (83 depth locations and 165 surface locations, Figure 4A; see Supporting Information, Figure 4-1, for regional distribution of stimulation location). We specifically ask, “Do functional hubs drive more widespread reconfiguration of the functional network than functionally isolated brain areas?” To answer this question, we measure the node strength as the mean coherence of the stimulation node to all other nodes at baseline. We hypothesize that stimulation of a stronger functional hub will lead to greater dispersion of input intensity throughout the network, driving a homogenous network response; stimulation of a weaker functional hub will lead to more targeted dispersion of input intensity to a subset of network nodes, driving a heterogenous network response (Figure 4B).

To test this hypothesis, we first compute a Spearman’s ρ correlation between baseline node strength of the stimulation region and the average of each of the three measures of functional network reconfiguration (change in the mean of node strengths, change in the variance of node strengths, and the configuration similarity) for the four coherence frequency bands (Figure 4C–E). Using a Bonferroni correction for multiple comparisons, we test whether greater node strength of the stimulation region drives greater node-level changes in the network 100 ms after stimulation offset, over stimulation sessions across subjects (Figure 4C, D). We find that stimulation node strength does not significantly influence the mean of node strengths in the alpha/theta band (ρ(246) = −0.02, p = 1.0), in the beta band (ρ(246) = −0.02, p = 1.0), in the low gamma band (ρ(246) = −0.06, p = 1.0), or in the high gamma band (ρ(246) = −0.03, p = 1.0). We also find that the stimulation node does not significantly influence the variance of node strengths in the alpha/theta band (ρ(246) = −0.05, p = 1.0), in the beta band (ρ(246) = −0.03, p = 1.0), in the low gamma band (ρ(246) = −0.14, p = 0.09), or in the high gamma band (ρ(246 = −0.11, p = 0.3). We also do not observe these effects after 600 ms following stimulation offset (Supporting Information: Figure 4-2c, d). Our results indicate that baseline node strength does not play an influential role in altering large-scale organization of network nodes.

We next ask whether the strength of the stimulation node can differentially drive reconfiguration of edge-level architecture of the network. Using a Bonferroni correction for multiple comparisons, we test whether greater node strength of the stimulation region drives greater configurational change in the network edges 100 ms after stimulation offset, over stimulation sessions across subjects (Figure 4E). We find that greater stimulation node strength leads to a significantly greater configuration similarity (lower reconfiguration) in the low gamma band (ρ(246) = 0.17, p = 0.03) and in the high gamma band (ρ(246) = 0.58, p = 2.2 × 10⁻²²), and to a nonsignificant change in the alpha/theta band (ρ(246) = 0.01, p = 1.0) and in the beta band (ρ(246) = 0.13, p = 0.2). We find that these effects persist in the low gamma band and in the high gamma band, and that they strengthen in the beta band at least 600 ms after stimulation offset (Supporting Information: Figure 4-1e). These results suggest that the functional topology of the stimulation region significantly impacts the pattern of coherent interactions in low and high gamma coherence frequency bands. Specifically, stimulation of weaker functional hubs tends to drive a greater change in the pattern of coherent interactions in low gamma networks and in high gamma networks. We find that a location-based rule for using stimulation to control the distributed reconfiguration of functional interactions is most robust for high gamma networks thought to reflect activity associated with synaptic input and short-range interactions.

Combined with our earlier findings on the negative relationship between stimulation intensity and edge reconfiguration, our findings suggest that stimulation of stronger functional hubs may lead to greater attenuation of the stimulation intensity and drive less edge-level reconfiguration than stimulation of weaker functional hubs. Another possible explanation for our findings is that stronger coherent interactions between stimulated hub nodes and the remaining nodes in the network mechanistically constrain the network response to stimulation—which we assess next.

Our findings in the previous section point to an important role of the baseline functional network in constraining the network response to stimulation. Logically, a stimulation node that exhibits stronger coherence with one set of brain regions may be more likely to convey the stimulation input to these brain regions than to another set of brain regions with which it exhibits weaker coherence. We therefore next test the hypothesis that the baseline coherence between the stimulation node and a downstream node predicts the probability that the downstream node will be evoked during a stimulation trial (Figure 5A). In other words, a stronger baseline coherence between the stimulation node and the downstream node may facilitate a greater magnitude change in the node strength of the downstream node.

To address this hypothesis, we first compute the within-session mean and variance of the magnitude change in node strength. For each session, we next compute the Spearman’s ρ correlation between the set of baseline coherence values between the stimulation node and the downstream nodes and the mean of the magnitude changes in node strength of each downstream node (Figure 5B). Intuitively, positive correlation values imply that stronger baseline coherence between the stimulation node and the downstream nodes predicts greater magnitude change in node strength across stimulation epochs. Using a Wilcoxon rank-sum test, Bonferroni correction for multiple comparisons, and effect size based on rank correlation, we find a significant positive correlation values in the alpha/theta band (Z(247) = 4,986, p = 0.003, rank correlation = 0.84), in the beta band (Z(247) = 2,092, p = 5.4 × 10⁻¹⁵, rank correlation = 0.93), in the low gamma band (Z(247) = 2,355, p = 1.4 × 10⁻¹³, rank correlation = 0.92), and we find a nonsignificant positive trend in the high gamma band (Z(247) = 5,928, p = 0.19). We also assess the Spearman’s ρ correlation between the set of baseline coherence values between the stimulation node and the downstream nodes and the variance of the magnitude changes in node strength of each downstream node (Figure 5B). We find a significantly positive correlation in the beta band (Z(247) = 2,721, p = 9.9 × 0⁻¹¹, rank correlation = 0.91) and in the low gamma band (Z(247) = 2,997, p = 2.0 × 10⁻¹⁰, rank correlation = 0.90), and we find a nonsignificant positive trend in the alpha/theta band (Z(247) = 5,993, p = 0.24) and in the high gamma band (Z(247) = 6,340, p = 0.72). We find that these effects persist in the alpha/theta band, in the beta band, and in the low gamma band and strengthen in the high gamma band at least 600 ms after stimulation offset (Supporting Information: Figure 5-1b,c).

Our findings are consistent with the hypothesis that the baseline functional network topology involving the stimulation node may be used to predict the downstream network regions that are most influenced by stimulation. We also find that this predictive capacity is dependent on the frequency band of the coherent measurement: The likelihood of modulating the coherent interactions of a downstream node is less predictable for higher frequencies. This finding suggests that the direct coherence between a stimulation node and a downstream node may be more influential in conveying stimulation input to the downstream node for lower coherence frequency bands. Our analysis highlights a putative mechanism of node-level flexibility, or the ability for a network region to dynamically alter its level of interaction with other regions in the network. Specifically, stronger baseline coherence between the stimulation node and a downstream node tends to predict greater variability with which the downstream node changes its level of interaction with the rest of the network within a stimulation session. Such a rule can guide more principled targeting of network structures to amenably drive flexible reconfiguration of the functional network.

We lastly seek to integrate our observations on stimulation-driven reconfiguration of functional brain networks with first principles theory. Network control theory provides a mathematical framework to model changes in the state of a complex system under a set of constraints prescribed by the structure of that system (Pasqualetti et al., 2014; Yan et al., 2017). For brain networks, network control theory offers an opportunity to model the logical progression of a stimulus input into an anatomically defined structural brain network, the traversal of that input through the network, the resulting change in interregional communication, and an accompanying shift in the dynamical state of the brain that accommodates a change in behavior (Figure 6A; Gu et al., 2017, 2015; Kim et al., 2018; Muldoon et al., 2016; Pequito et al., 2016). The structural topology of the network may confer important control properties to a complex system such as modal controllability, which enables a system to move from its current dynamical state to more difficult-to-reach dynamical states through an efficient expenditure of energy resources (Ashourvan, Gu, Mattar, Vettel, & Bassett, 2017; Gu et al., 2015; Pasqualetti et al., 2014). Recent theoretical inquiry into the relationship between brain structure and function during stimulation posited that stimulation of the structural brain network’s modal controllers may drive a heterogenous change in functional architecture (Muldoon et al., 2016). However, experimental evidence linking the network control theoretic model to brain stimulation and its influence on functional architecture and dynamical brain state via the structural brain network is lacking.

In this study, we have so far shown that stimulation drives a rich and complex functional network response that is dependent on the stimulation input intensity and the stimulation input location. While the input intensity drives the magnitude of functional reconfiguration, the baseline coherence of the input location constrains the spatial specificity of the functional reconfiguration. Yet, how does baseline network topology in the vicinity of the stimulation region relate to the modal control strategy put forth by structural control theory? To answer this question, we construct structural brain networks by applying deterministic tractography to diffusion-weighted imaging data that are parcellated into 1,015 anatomically defined, cortical and subcortical regions of interest (Cammoun et al., 2012) in a subset of 14 participants (see Methods; the other participants did not have diffusion-weighted imaging data). We next assign the intracranial ECoG sensors for a subject to their nearest anatomical ROI based on the shortest spatial distance between a sensor and the ROI centroids. We then measure the structural connectivity between anatomical ROIs and the modal controllability of each anatomical ROI in the structural network, and we assign these values to the intracranial sensors based on their proximity to the nearest ROI.

To establish a clear relation between structural connectivity and baseline functional connectivity (measured by coherence), we calculate the Pearson correlation between the pattern of structural connections and the pattern of baseline coherent interactions between intracranial sensors, for each coherence frequency band and each subject (Figure 6B). This procedure is conducted separately for connections that span the same hemisphere of the brain (intra-hemispheric connections) and for connections that span opposing hemispheres of the brain (inter-hemispheric connections), and correlation values are subsequently z-scored using the Fisher transformation. Using a one-sample t test, we test the null hypothesis that the z-transformed correlation between structural connectivity and functional connectivity is equal to 0. For intra-hemispheric connections, we observe significant negative correlation between structural connectivity and baseline coherence in the low gamma band (t(13) = −3.0, p = 0.04, corrected) and in the high gamma band (t(13) = −3.6, p = 0.01, corrected), and non-significant negative trends in the alpha/theta band (t(13) = −1.5, p = 0.64, corrected) and in the beta band (t(13) = −1.8, p = 0.34, corrected). For inter-hemispheric connections, we observe nonsignificant positive trends between structural connectivity and baseline coherence in the alpha/theta band (t(13) = −0.7, p = 0.48, corrected), in the beta band (t(13) = 0.9, p = 0.37, corrected), in the low gamma band (t(13) = 0.9, p = 0.39, corrected), and in the high gamma band (t(13) = 0.8, p = 0.44, corrected). Generally, inter-hemispheric connections exhibited positive trends between structural connectivity and ECoG functional connectivity, while intra-hemispheric connections exhibited negative trends between structural connectivity and ECoG functional connectivity. This observed relationship was stronger at the finer scale parcellation of the brain than at the coarser scale parcellation (Supporting Information: Figure 6-1). Our results imply that weaker structural connections between brain regions within the same hemisphere are generally associated with stronger coherent interactions in higher frequency bands, which typically reflect local, bottom-up processing associated with synaptic input. This divergent relationship can be explained by an expectedly strong coherence between brain regions that are in close proximity to one another, which are more likely to be weakly linked by structural U-shaped white matter fibers (Schmahmann et al., 2007).

Based on the result that structural connectivity weakly constrains baseline coherent interactions between brain regions, we ask whether a control strategy motivated by modal controllability of the structural brain network predicts stimulation-driven reconfiguration of the functional network topology. To answer this question, we first calculate the modal controllability for the brain region in each stimulation session across subjects. We next compute the Spearman’s ρ correlation between modal controllability of the stimulated brain region and the average, within-session configuration similarity. Using Bonferroni correction for multiple comparisons, we find a significant positive correlation between modal controllability and configuration similarity in the beta band (ρ(24) = 0.61, p = 0.003), in the low gamma band (ρ(24) = 0.65, p = 0.001), and in the high gamma band (ρ(24) = 0.59, p = 0.005), and we find a nonsignificant positive trend in the alpha/theta band (ρ(24) = 0.47, p = 0.05). These effects persist at least 600 ms after stimulation offset (Supporting Information: Figure 6-2a). These results imply that stimulation of stronger modal controllers drives lower reconfiguration of the coherent interactions in the functional network. To contextualize these findings, we draw upon the well-established positive relationship between modal controllability and structurally isolated brain regions (Gu et al., 2015). By targetting strong modal controllers in the structural brain network, stimulation modulates structurally isolated brain areas and drives lower functional reconfiguration than stimulation of weak modal controllers in structurally connected brain areas.

Notably, our results demonstrate a complex and frequency-specific link between structural topology and functional topology—weak structural connections tend to span between brain areas with stronger, high-frequency functional interactions, and stimulation of modal controllers in structurally isolated brain regions tends to limit the extent of functional network reconfiguration. Combined with earlier findings, we put forth a putative sequence of physiological events associated with modal control in which stimulation of strong modal controllers activates strong, local functional hubs that drive less functional reconfiguration of network-wide edges (Figure 4E) and greater change in the node strengths of functionally connected, downstream brain regions (Figure 5B). In contrast, stimulation of weak modal controllers activates weaker, functionally isolated regions that drive distributed functional reconfiguration of network-wide edges (Figure 4E). While these findings establish a link between network control and functional reconfiguration, they do not establish a link between network control and changes in dynamical brain state.

To experimentally examine the relationship between the modal controllability of the stimulation region and the shift in dynamical brain state following stimulation, we leverage a previously documented and validated binary classifier of neural activity into good and poor episodic memory encoding states (Ezzyat et al., 2017). Specifically, we first train the classifier to discriminate between successful and unsuccessful word recall trials during a delayed free recall task using features based on spectral power of ECoG activity. We next evaluate the classifier on ECoG activity during the pre-stim epoch and the post-stim epoch of each stimulation trial and compute the change in probability of being in a good memory encoding state. We calculate the correlation between the modal controllability of the stimulation region and the average, within-session change in probability of being in a good memory state over stimulation trials (Figure 6D). We find a significant positive correlation (Pearson’s r(17) = 0.49, p = 0.03). These findings imply that the push towards better memory encoding states is associated with stimulation of strong modal controllers that are theoretically positioned to push the brain to more energetically unfavorable and distant brain states.

Collectively, our study reveals a link between brain structure and brain function that is grounded in network control theory. Using the network control framework, we uncover the important role of stimulation on reconfiguration of functional architecture that accounts for anatomical constraints on network dynamics via the topology of structural network connectivity. We find that brain networks may use a modal control strategy during transitions between difficult-to-reach dynamical states, which is associated with a reconfiguration in the localized coherence of individual network nodes to the broader functional brain network.

Here, we addressed the hypothesis that direct stimulation of cortical and subcortical structures alters the architecture of functional brain networks and shifts the dynamical brain state in accord with control strategies identified by applying tools from network control theory to electrophysiological and structural brain imaging data. In human epilepsy patients, we measured coherent patterns of ECoG activity thought to underlie coordinated functional interactions, and we mapped how these interactions vary with neurostimulation parameters. We observed that stimulation drives two modes of functional reconfiguration: The first mode involves distributed changes in the pattern of functional interactions across the network, and the second mode involves preferentially localized changes in the functional interactions associated with select brain regions. Notably, the mode of reconfiguration may be strategically selected based on the strength and location of stimulation. When we stimulated brain regions with weak structural connections to the rest of the network, we tended to invoke a modal control strategy marked by a modulation of the functional hubness of downstream brain regions and a large change in dynamical brain state.

The field of network neuroscience has long sought to understand how the rigid and interconnected anatomy of the structural brain network shapes interactions amongst functionally specialized brain areas, which may change from moment to moment and drive cognition and behavior (Hermundstad et al., 2013, 2014; Mišić et al., 2016; Park & Friston, 2013). Tools from network control theory (Liu, Slotine, & Barabasi, 2011; Pasqualetti et al., 2014; Tang & Bassett, 2017) have enabled researchers to identify controllability rules that prescribe how the dynamical state of a neural system can change based solely on the structural topology of the system (Betzel et al., 2016; Gu et al., 2017, 2015; Kim et al., 2018; Tang et al., 2017; Yan et al., 2017). Structural controllability rules account for network interactions that can occur, but they do not account for finer scale functional constraints that dictate whether these interactions will occur at a point in time (Liu et al., 2011; Tang & Bassett, 2017). Previous studies have incorporated these functional constraints into the study of network control by using neurophysiologically inspired dynamical mean-field models (Jirsa et al., 2017; Muldoon et al., 2016; Taylor et al., 2015), which require estimation of biologically plausible parameters.

In contrast to these studies, we used a data-driven, perturbative approach for inferring rules of functional network reconfiguration. By focally stimulating brain tissue at the millimeter scale, we mapped changes in the statistical interdependencies between brain regions. We note that there is a subtle distinction in the statistical methods used here to measure bidrectional, synchronized interactions and statistical methods used elsewhere to measure directed, effective functional interactions (Battaglia, Witt, Wolf, & Geisel, 2012; Friston, 2011; Lepage, Ching, & Kramer, 2013). We found that the observed change in statistical interdependency can be predicted by baseline levels measured before any stimulation is delivered. That is, a brain region that exhibits a strong, spontaneous functional interaction with the stimulated brain region is more likely to be modulated during stimulation than a brain region with a weak, spontaneous functional interaction with the stimulated brain region. In the long-standing debate regarding the validity of functional network models to explain causal dynamics (Jonas & Kording, 2017), these data provide compelling evidence in favor of a mechanism in which the functional network may causally convey the influence of stimulation on one brain region to other strongly interacting brain regions. Our findings implicate a strategy for the functional control of brain networks in which (a) stimulation of functionally isolated brain regions leads to spatially focal and strong downstream functional reconfiguration, and (b) stimulation of functionally hub-like brain regions leads to spatially diffuse and weak downstream functional reconfiguration.

To relate rules for functional control to theoretical predictions from structural controllability, we bridged electrophysiological data with structural brain imaging data. Confirming results from prior studies on controllability in healthy subjects (Gu et al., 2015; Tang et al., 2017), in epilepsy subjects we found that brain regions with structurally weak connections tend to be strong modal controllers, which facilitate more difficult transitions in brain state. Stimulations of structurally weak modal control points lead to strong, local change of functional architecture in a mean-field model of neuron population dynamics (Muldoon et al., 2016). In contrast, here we demonstrated that stimulation of modal control points leads to widespread, weak change in functional architecture because of the presence of strong functional hubs in the stimulated brain regions. We identified two potential explanations for the distinction between the model simulation and our empirical observations. First, the in silico model assumes identical biophysical parameters across different neuronal ensembles distributed across the brain network, which may limit the reproducibility of detailed spatial and temporal dynamics that would otherwise be expressed in vivo and measured by intracranial ECoG sensors. Second, the in silico model accounts for structural connectivity between neuronal ensembles using diffusion imaging, which measures white matter fiber pathways spanning long distances but does not capture gray matter pathways spanning short distances (Thomas et al., 2014) and synaptic microarchitecture responsible for plasticity over varying time scales (Song, Miller, & Abbott, 2000). When combined with additional data demonstrating a moderate relationship between white matter connectivity and correlated ECoG dynamics (Betzel, Medaglia, Kahn, et al., 2017), our findings suggest that non–white matter structural connectivity and other physiological factors may contribute to the reconfiguration of functional network architecture.

Neuronal synchronization is purported to play a crucial role in facilitating interareal communication between ensembles of neurons (Bonnefond et al., 2017; Buzsáki et al., 2012; Canolty & Knight, 2010; Fries, 2015; Schalk, 2015). Fries (2015) proposed that rhythmic oscillations in the local field potential give rise to states of excitability depending on the temporal position during an oscillatory cycle—two different ensembles of neurons are able to reliably transfer information between one another when they are mutually excitable, or exhibit oscillations that are in-phase. Equally important to communication is the frequency of the oscillation—higher frequency bands (γ) are thought to facilitate communication of bottom-up input over short distances, and lower frequency bands (θ, α, β) are thought to facilitate communication of top-down processes over long distances (Fries, 2015; Kopell et al., 2000).

Notably, we found that stimulation parameters may be tuned to selectively modulate different regions and spatial extents of the functional network. Stimulation intensity tends to provide greater control over reconfiguration of lower frequency networks, and stimulation location tends to provide greater control over reconfiguration of higher frequency networks. Specifically, we observed that the strength of the stimulation input has a greater effect on functional reconfiguration in lower frequency bands than in higher frequency bands, suggesting that stimulation parameters such as amplitude, pulse frequency, and duration may play an important role in the modulation of long-range, top-down functional interactions. We speculate that a stronger stimulation input may be more likely to penetrate wider spatial extent of cortex and heterogenously modulate network excitability at low frequencies (Fries, 2015). We also observed that the location of the stimulation input has a greater effect on functional reconfiguration in higher frequency bands than in lower frequency bands, suggesting that the hubness of the stimulated brain region may play an important role in the modulation of short-range, bottom-up functional interactions. Intuitively, we expect that brain regions involved in bottom-up communication associated with broadly conveying sensory input to higher order cortices might also be more sensitive to modulations via stimulation than brain regions involved in top-down communication.

We chose to perform a network analysis of intracranial data during neurostimulation, rather than a univariate analysis of individual activations. Our decision enabled us to examine the influence of neurostimulation on the distributed and interconnected physiology of the human brain, which is in line with previous in silico modeling work (Muldoon et al., 2016). Several notable studies have teased apart the effective connectivity between brain regions using neurostimulation to generate cortico-cortico evoked potentials (CCEPs; Keller, Honey, Entz, et al., 2014; Keller, Honey, Me, et al., 2014; Matsumoto et al., 2004), and have more recently predicted where CCEPs might occur using a combination of repetitive low-frequency stimulation and baseline connectivity (Keller et al., 2018). In a departure from these studies on effective connectivity, we test the novel hypothesis that neurostimulation can predictively perturb interareal statistical dependencies underlying distributed brain function. Our study is motivated by recent hypotheses regarding the role of coherent synchronization in neuronal communication. With the maturation of network analysis tools that enable simultaneous tracking of dynamic network architecture and dynamic activity (Murphy et al., 2016), we envision future studies where we investigate the potential for neurostimulation to selectively modulate brain activity or perturb functional network architecture.

We also chose to examine static changes in the functional network architecture aggregated over many repeated trials of neurostimulation. This approach allowed us to examine the ability of stimulation to drive network reconfiguration within each subject, and to evaluate the statistical robustness of effects across subjects. Generally, we find medium to large effect sizes on the order of 0.4–0.9. The advantage of our chosen approach is that the effects of stimulation are assessed separately for each patient, which enabled us to account for an individual’s stereotyped network topology while examining generalized rules that predict the network effects of stimulation. We can thus address questions such as, “Are some brain regions more likely to change their functional interactions than other brain regions?” and “How are these brain regions associated with the stimulated brain region?” In future studies, we aim to understand how stimulation influences the functional brain network from one moment in time to the next, as a function of brain state. For instance, does stimulation during the beginning of a coherent oscillatory cycle influence network architecture differently than stimulation during the middle of a coherent oscillatory cycle? This temporal mapping could inform control strategies to steer functional brain network reconfiguration in real time.

One important question regarding the neural adaptation effects of stimulation is beyond the scope of this analysis and remains unanswered: “How long do stimulation-induced changes to the network last?” An answer to this valuable question could support therapeutic applications of direct stimulation by reducing the need to stimulate the brain for long periods of time at high intensity—potentially improving safety, efficacy, and longevity of chronically implantable devices (Fisher, 2011; Stacey & Litt, 2008). The adaptation question will be addressed in a future study by introducing longer periods of rest between consecutive stimulation trials and quantifying the duration of connectivity changes. A related consideration of this study is that a sham condition was not included during the stimulation session and changes in functional network architecture were instead compared with passive changes in network architecture at baseline. Indeed, a dedicated sham condition during stimulation would be able to dissociate the effects of stimulation from natural fluctation in the brain’s internal dynamics. Previous studies have emphasized that sham stimulation may be more effective at assessing behavioral, “placebo” effects of stimulation than quantifying neural effects, given the biophysical differences between external input from a current source in comparison to spontaneous fluctuation of endogenous neural activity (Duecker & Sack, 2015). In this study, we enforced a fairer comparative condition to understand neural effects of stimulation by parameterically titrating the dose of stimulation via current amplitude, frequency, and duration—all likely inducing a more similar biophysical change to neural tissue than a sham condition with zero current delivery. Furthermore, stimulation trials were made blind to the patient, trial durations were kept short (below 1 s), and the start of a stimulation trial was randomly jittered to mitigate neural expectation of stimulation timing.

Approaches for recording and interrogating intracranial electrophysiology are inherently limited by spatial coverage of the ECoG sensors, which is determined during the management of a patient’s epilepsy. This sampling bias leads to a varied representation of the functional brain network between individual patients. While it is not yet possible to record from the entirety of the human brain using ECoG, we mitigated this shortcoming by taking key steps in our analysis. First, we used a statistically robust approach for characterizing the network impact of stimulation in individual patients. We computed separate measures of topology for each patient’s functional brain network, which enabled us to account for individual variability in sensor placement and physiological state. Our data demonstrated a set of functional rules for network reconfiguration that fundamentally depend on topological characteristics of the stimulated brain area that can vary within and between patients. Second, we used a large dataset consisting of 94 epilepsy patients, allowing us to account for a range of individual variability in functional brain network architecture that is often not possible in studies of human electrophysiology.

Here we mapped the impact of targeted neurostimulation on distributed functional architecture in the human brain. We demonstrated that network physiology can be predictably altered based on control theoretic rules that account for structural and functional organization of the brain network. Our results provide a causal, quantified description of the influence of structure and function on dynamical brain state. Our findings have significant translational implications in strategizing stimulation-based therapy based on a combination of behavioral biomarkers and neurophysiology.

Ankit N. Khambhati: Conceptualization; Formal analysis; Investigation; Methodology; Software; Validation; Visualization; Writing – Original Draft; Writing – Review & Editing. Ari E. Kahn: Methodology. Julia Costantini: Formal analysis; Software. Youssef Ezzyat: Conceptualization; Data curation. Ethan A. Solomon: Conceptualization; Data curation. Robert E. Gross: Data curation. Barbara C. Jobst: Data curation. Sameer A. Sheth: Data curation. Kareem A. Zaghloul: Data curation. Gregory Worrell: Data curation. Sarah Seger: Data curation. Bradley C. Lega: Data curation. Shennan Weiss: Data curation. Michael R. Sperling: Data curation. Richard Gorniak: Data curation. Sandhitsu R. Das: Data curation; Methodology; Software. Joel M. Stein: Data curation; Methodology; Software. Daniel S. Rizzuto: Data curation; Project administration. Michael J. Kahana: Data curation; Funding acquisition; Project administration; Resources. Timothy H. Lucas: Data curation. Kathryn A. Davis: Data curation. Joseph I. Tracy: Data curation. Danielle S. Bassett: Conceptualization; Funding acquisition; Methodology; Project administration; Resources; Supervision; Writing - Original Draft; Writing - Review & Editing.

Danielle S. Bassett, John D. and Catherine T. MacArthur Foundation (http://dx.doi.org/10.13039/100000870). Danielle S. Bassett, Alfred P. Sloan Foundation (http://dx.doi.org/10.13039/100000879). Danielle S. Bassett, Army Research Laboratory (http://dx.doi.org/10.13039/100006754), Award ID: W911NF-10-2-0022. Danielle S. Bassett, Army Research Laboratory (http://dx.doi.org/10.13039/100006754), Award ID: W911NF-14-1-0679. Danielle S. Bassett, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: 2-R01-DC-009209-11. Danielle S. Bassett, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: 1R01HD086888-01. Danielle S. Bassett, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: R01-MH107235. Danielle S. Bassett, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: R01-MH107703. Danielle S. Bassett, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: R21-M MH-106799. Danielle S. Bassett, U.S. Naval Research Laboratory (http://dx.doi.org/10.13039/100009917). Danielle S. Bassett, National Science Foundation (http://dx.doi.org/10.13039/100000001), Award ID: BCS-1441502. Danielle S. Bassett, National Science Foundation (http://dx.doi.org/10.13039/100000001), Award ID: CAREER PHY-1554488. Danielle S. Bassett, National Science Foundation (http://dx.doi.org/10.13039/100000001), Award ID: BCS-1631550. Michael J.Kahana, Defense Advanced Research Projects Agency (http://dx.doi.org/10.13039/100000185), Award ID: N66001-14-2-4032.

We thank Blackrock Microsystems for providing neural recording and stimulation equipment. The views, opinions, and/or findings contained in this material are those of the authors and should not be interpreted as representing the official views or policies of the Department of Defense, the U.S. government, or any of the funding agencies.

Direct neurostimulation:

Delivery of electrical pulses to electrodes located on the surface or in deep structures of the brain.

Control policy:

Rules for when, where, and how inputs should be delivered to a dynamical system to reach the desired target state.

Network control theory:

Framework to model changes in the state of a complex system under constraints prescribed by the structure of that system.

Graph:

Mathematical representation of a network consisting of nodes as brain regions and edges as structural or functional connections between regions.

Controllability:

Ability to steer a dynamical system from its current state towards a desired target state over a finite time horizon.

Coherence:

Functional communication capacity (connectivity) between two brain regions based on phase alignment of their oscillatory activity at a specific frequency.

Node strength:

Measure of global influence, or hubness, of a node based on its average connectivity to all other network nodes.

Memory encoding state:

Pattern of brain activity associated with the likelihood that a subject has successfully encoded a word into their episodic memory.

Edge reconfiguration:

Alteration in the spatial arrangement of connections in a dynamic network.

Modal controller:

Network node that facilitates energetically challenging state transitions and is defined based on its rigid structural connectivity with other network nodes.