Abstract
The brain is an inherently dynamic system, and much work has focused on the ability to modify neural activity through both local perturbations and changes in the function of global network ensembles. Network controllability is a recent concept in network neuroscience that purports to predict the influence of individual cortical sites on global network states and state changes, thereby creating a unifying account of local influences on global brain dynamics. While this notion is accepted in engineering science, it is subject to ongoing debates in neuroscience as empirical evidence linking network controllability to brain activity and human behavior remains scarce. Here, we present an integrated set of multimodal brain–behavior relationships derived from fMRI, diffusion tensor imaging, and online repetitive transcranial magnetic stimulation (rTMS) applied during an individually calibrated working memory task performed by individuals of both sexes. The modes describing the structural network system dynamics showed direct relationships to brain activity associated with task difficulty, with difficult-to-reach modes contributing to functional brain states in the hard task condition. Modal controllability (a measure quantifying the contribution of difficult-to-reach modes) at the stimulated site predicted both fMRI activations associated with increasing task difficulty and rTMS benefits on task performance. Furthermore, fMRI explained 64% of the variance between modal controllability and the working memory benefit associated with 5 Hz online rTMS. These results therefore provide evidence toward the functional validity of network control theory, and outline a clear technique for integrating structural network topology and functional activity to predict the influence of stimulation on subsequent behavior.
SIGNIFICANCE STATEMENT The network controllability concept proposes that specific cortical nodes are able to steer the brain into certain physiological states. By applying external perturbation to these control nodes, it is theorized that brain stimulation is able to selectively target difficult-to-reach states, potentially aiding processing and improving performance on cognitive tasks. The current study used rTMS and fMRI during a working memory task to test this hypothesis. We demonstrate that network controllability correlates with fMRI modulation because of working memory load and with the behavioral improvements that result from a multivisit intervention using 5 Hz rTMS. This study demonstrates the validity of network controllability and offers a new targeting approach to improve efficacy.
Introduction
A dominant narrative in contemporary neuroscience is that the brain is a vastly dynamic system of interacting networks called into action to achieve different mental states and accomplish cognition. An emerging notion in such network neuroscience views is that global brain states might be controllable from a single node within a cortical network (Betzel et al., 2016), thereby imparting certain privileged roles for different brain areas in directing neural activity, such a conception is supported by growing evidence for network hubs, such as the dorso-lateral prefrontal cortex (DLPFC). Within this context, brain states can be altered through exogenous perturbations via noninvasive neuromodulation, such as repetitive transcranial magnetic stimulation (rTMS). This approach is attractive because it raises the possibility of causally influencing the brain in a controlledmanner to test specific hypotheses about network dynamics. However, questions remain about the manner in which control is exerted in the brain (Tu et al., 2018; Pasqualetti et al., 2019) because of limited empirical evidence supporting the role of controllability in human behavior. Of particular interest for cognitive neuroscience is the concept of modal controllability (MC), a measure of the ability of a single region to drive the network system (brain) into difficult-to-reach states potentially associated with high cognitive performance (Pasqualetti et al., 2014). Modal control was originally defined within engineering frameworks to develop systematic means of adjusting input matrices representing large-scale linear time-invariant multivariable systems (Simon and Mitter, 1968). In the present context, MC may be defined as control that changes the modes (i.e., the eigenvalues of the system matrix) to achieve the desired control objectives. Defining “brain state” as the magnitude of neurophysiological activity across brain regions, we can operationalize MC as reflecting a region's contribution to difficult-to-reach global activity patterns. The functional relevance of MC may therefore be evidenced by relationships to cognitive demand, metabolic cost (e.g., changes in BOLD activation), or rTMS-induced behavioral changes.
In light of these questions, we sought to address the validity of the controllability account by linking MC to fMRI-based brain states and assessing the value of MC in predicting behavioral performance after exogenous perturbation was used to manipulate these states. Specifically, we examined the relationship between MC and the behavioral benefit associated with online rTMS during a working memory (WM) task. We chose one of the most reliable brain patterns in fMRI-based cognitive neuroscience research to manipulate: the increased activity in frontoparietal areas associated with greater WM difficulty (Blumenfeld and Ranganath, 2007). Our choice was influenced by the relatively high estimation of MC in frontoparietal cortex (Gu et al., 2015) and a number of computational simulations demonstrating the capacity of MC to predict functional patterns (Betzel et al., 2016; Muldoon et al., 2016), although there are currently little empirical data supporting this interpretation. Further, the effects of neuromodulation are typically observed not only in the stimulated site, but also in distally connected regions (Bestmann et al., 2004; Beynel et al., 2020). Therefore, the controllability framework offers an opportunity to observe how the network properties of a single node may predict rTMS effects dependent on a more global WM network.
We tested our empirical data against the predictions set out by network control theory. To bridge the gap between theoretical and observed brain states, we first investigated the relationship between global brain activation pattern and the eigenmodes of the structural brain networks, from which MC is derived. Specifically, we predicted that the activation pattern in the hardest task condition should align with eigenmodes representing theoretically difficult-to-reach states. We related MC at the stimulation site to the parametric increase of functional activation as a function of WM difficulty. A higher parametric increase in activity would suggest that additional neural resources are being successfully recruited in the more difficult task conditions, and that state transitions can be relatively easy to achieve. Within this context, if MC indicates a transition into a difficult-to-reach state, we predicted a negative relationship between MC and difficulty-driven BOLD activation (i.e., as a function of WM load). Last, if difficult state transitions can be augmented by exogenous stimulation with excitatory rTMS to a region with high MC, then we predicted that subjects with higher MC and lower parametric activation at the stimulated site would show the greatest benefit from rTMS. Such a set of findings will help validate the value of network controllability in predicting brain states and help predict rTMS effects foundational to many FDA-approved treatments.
Materials and Methods
Subjects
Twenty-nine young adults (mean age 23.38 ± 5.13 years; 14 females) completed this 6-day study (Fig. 1A), approved by the Duke Medical School Institutional Review Board (#Pro00065334). All procedures and analyses were performed in accordance with Institutional Review Board guidelines and regulations for experimental testing. Subjects had no history of psychiatric or neurologic disorders and were not using psychoactive drugs or medications known to lower seizure threshold. Five subjects were excluded because of incomplete imaging data (because of excessive movement).
Study protocol. A, Timeline of the experimental protocol. B, The DRAT, with examples of the four levels of difficulty (levels were individually titrated); the easiest difficulty level was dropped for Visits 3-6.
Experimental design
The following description outlines relevant design details and parameters for our 6-day protocol; additional information can be found in Davis et al. (2018) and Beynel et al. (2019). This study is a reanalysis of data previously reported by Beynel et al. (2019) that is now extended to report on controllability, including diffusion weighted imaging (DWI) data, which were not reported on in these previous articles. During the first visit, subjects consented to participate and were screened to make sure they did not have any contraindications to TMS or MRI. Resting motor threshold (RMT) was then assessed and subjects performed the WM task. We used a modified Sternberg task in which arrays of letters were manipulated within WM (D'Esposito et al., 1999), noted here as the delayed-response alphabetization task (DRAT; Fig. 1B). For each trial, letters were presented, followed by a delay period during which subjects mentally reorganized letters into alphabetical order. A letter with a number above it was then presented, and subjects were asked whether (1) the probe letter was or was not in the original set; (2) the letter was in the original set and the number matched the serial position of the letter once the sequence was alphabetized; or (3) the letter was in the original set but the number did not match the serial position of the letter once alphabetized. These conditions are referred to as New, Valid, and Invalid, and consisted of 20%, 40%, and 40% of the trials, respectively.
During the second visit, MRI was performed in a 3 T gradient-echo scanner (General Electric 3.0 Tesla Sigma Excite HD short bore scanner), equipped with an 8-channel head coil at the Duke Brain Imaging Analysis Center. A structural MRI and a DWI scans were acquired, as well as four blocks of functional acquisitions while subjects performed the DRAT with four difficulty levels (Very Easy, Easy, Medium, Hard) that were individualized based on performance from the first visit. The anatomic MRI was acquired using a 3D T1-weighted echo-planar sequence (matrix = 2562, TR = 12 ms, TE = 5 ms, FOV = 24 cm, slices = 68, slice thickness = 1.9 mm, sections = 248). DWI data were collected using a single-shot EPI sequence (TR = 1700 ms, slices = 50, thickness = 2.0 mm, FOV = 256 × 256 mm2, matrix size 128 × 128, voxel size = 2 mm3, b value = 1000 s/mm2, diffusion-sensitizing directions = 36, total images = 960, total scan time = 5 min). Finally, in the fMRI runs, coplanar functional images were acquired using an inverse spiral sequence (64 × 64 matrix, TR = 2000 ms, TE = 31 ms, FOV = 240 mm, 37 slices, 3.8 mm slice thickness, 254 images). The total scan time was ∼1 h 40 min. fMRI analysis was then performed (see fMRI analysis) to define a target for the subsequent rTMS visits.
During Visits 3-6, active and sham rTMS were performed with an active/placebo figure-8 coil (A/P Cool-B65) and a MagPro X100 stimulator with MagOption (MagVenture). Sham stimulation was applied using the same coil in placebo mode, which produces clicking sounds and somatosensory sensation via electrical stimulation with scalp electrodes mimicking those occurring in the active mode, but with a substantially attenuated E-field penetrating the head. The coil position over the individualized fMRI target within the left medial frontal gyrus was continually monitored through a stereotaxic neuronavigation system (Brainsight, Rogue Research) and maintained at a high level of precision throughout the session with real-time robotic guidance (Smart Move Robot, Advanced Neuro Technology).
rTMS was applied “online” (i.e., while subjects performed the WM task at three difficulty levels: Easy, Medium, Hard). Online stimulation was chosen over offline rTMS to observe the immediate and specific effects of rTMS rather than cumulative effects over time; and to promote Hebbian-like plasticity through the concurrent activation of the frontoparietal network both endogenously through the WM task and exogenously with rTMS. Although in many studies online rTMS has been used to disrupt neurocognitive processes, online rTMS has also been shown to enhance cognitive processes (Hamidi et al., 2008). One explanation is that rTMS can entrain endogenous task-related oscillatory dynamics (Thut et al., 2011; Luber and Lisanby, 2014). In our study, online rTMS was applied at 5 Hz to entrain theta oscillations because (1) the theta frequency has been shown to play a crucial role in WM (Roux and Uhlhaas, 2014), since it increases with the number of items being maintained (Boonstra et al., 2013), and increases at the onset of a WM task, continuing its activation until the trial's end (Raghavachari et al., 2001); and (2) because there is evidence that theta-band rTMS can enhance WM (Luber et al., 2007; Yamanaka et al., 2010; Albouy et al., 2018; Riddle et al., 2020). During each visit, 10 blocks of the DRAT were performed (25 trials per block, 20% of those trials were New conditions and were therefore excluded from the analysis, resulting in 20 trials of interest per block): four blocks with active stimulation, and four blocks with sham stimulation, both preceded by a block without stimulation. The order of the stimulation type (active or sham) was presented according to an ABBA schedule on Visits 1 and 2, and a BAAB schedule on Visits 3 and 4. Twenty-five rTMS pulses were applied either before encoding or during the delay period of the DRAT, and this stimulation timing was alternated from block to block. Subjects' responses were collected on each trial, and feedback was displayed at the end of each block to inform the participant about their overall accuracy on that block of trials. The protocol included four visits, each with 8 blocks with stimulation, and 20 trials of interest per block, resulting in 640 trials that were divided across 24 conditions: Valid/Invalid, Active/Sham, rTMS before the encoding/during the delay, and Easy/Medium/Hard difficulty levels. This led to a total of ∼27 trials per condition.
Statistical analysis
fMRI analysis
Functional images were preprocessed using image processing tools, including FLIRT (FMRIB's. Linear Image Registration Tool) and FEAT (FMRIB Expert Analysis Tool) from FMRIB's Software Library (FSL, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/). Images were corrected for slice acquisition timing, motion, and linear trend; motion correction was performed using FSL's MCFLIRT, and six motion parameters estimated from the step were then regressed out of each functional voxel using standard linear regression. Images were then temporally smoothed with a high-pass filter using a 190 s cutoff and normalized to the MNI stereotaxic space. White matter and CSF signals were also removed from the data and regressed from the functional data using the same method as the motion parameters. Spatial filtering with a Gaussian kernel of FWHM of 6 mm was applied. BOLD activations during the array presentation, the delay period, and the response period were entered in a standard GLM, using HRF-convolved trial regressors and the temporal derivatives. Since the delay period has been associated with WM processing, the delay period was modeled by four independent regressors (Very Easy/Easy/Medium/Hard) to capture the whole-brain activation pattern as a function of task difficulty. To quantify the difficulty-driven state transition, the functional modulation (FM) of task difficulty was modeled in an additional GLM analysis by a single delay-period parametric regressor weighted linearly by the four levels of difficulty in the DRAT (orthogonalized to a nonweighted regressor for the same event), as well as its temporal derivative. Output from this model thus reflects the parametric increase in brain activation associated with the increase in task difficulty (for more conventional univariate models of these same data, see Davis et al., 2018 and Beynel et al., 2019).
Electric-field modeling
For rTMS targeting, the peak activation within left medial frontal gyrus was selected in each participant, and electric field (E-field) modeling was used to refine the optimal coil position and orientation. The E-field was simulated using T1, T2, and DWI images and the finite element method in the SimNIBS software package (Thielscher et al., 2015). The models featured five compartments with distinct electrical conductivities: skin, skull, CSF, gray matter (GM), and white matter. The DWI information was used to generate anisotropic conductivities for the white matter using the volume-normalized approach. The E-field was simulated at 54 candidate coil placements: 9 positions generated by placing a 3 × 3 grid with 1 cm spacing above the peak fMRI activation and 6 orientations per position corresponding to 30° rotation increments in a 180° semicircle. The coil position and orientation showing the highest correlation between the E-field magnitude and the fMRI activation pattern were selected as the rTMS target. This final target was then entered into the neuronavigation system (BrainSight, Rogue Research) and maintained at a high level of precision throughout the session with real-time robotic guidance (SmartMove Robot, Advanced Neuro Technology).
Structural network construction
Information on the structural connections based on diffusion tractography between each pair of regions in our data was assessed with a standard DWI processing pipeline used previously in our group. DWI data were analyzed using FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) and MRtrix (http://mrtrix.org) software packages. Data were denoised, corrected with eddy current correction, and bias-field corrected using MRtrix (Tournier et al., 2007) and FSL. Constrained spherical deconvolution was used in calculating the fiber orientation distribution. This fiber orientation distribution was used along with the brain mask to generate whole-brain tractography. Relevant parameters regarding track generation were as follows: seed = at random within mask; step size = 0.2 mm; 10 million tracts. After tracts were generated, they were filtered using spherical-deconvolution informed filtering of tractograms (Smith et al., 2012) to improve the quantitative nature of the whole-brain streamline reconstructions used here. This process utilizes an algorithm that determines whether a streamline should be removed or not based on information obtained from the fiber orientation distribution, which improves the selectivity of structural connectomes by using a cost function to eliminate false-positive tracts. Tracts were SIFTed until 1 million tracts remained. Connectomes were then generated by using FLIRT to apply a linear registration to the atlas used in the present study (described below) to register them to native diffusion space. Subsequent connectomes describe the number of streamlines connecting any pair of regions within this atlas, the summary network density measure describes the total number of streamlines within a subject's structural connectome.
Cortical parcellation and ROI extraction
A consistent parcellation scheme across all subjects and all modalities (DWI, fMRI) was used. Subjects' T1-weighted images were segmented using SPM12 (www.fil.ion.ucl.ac.uk/spm/software/spm12/), yielding a GM and white matter mask in the T1 native space for each subject. The entire GM was then parcellated into 471 ROIs, each representing a network node by using a sub-parcellated version of the Harvard-Oxford Atlas (HOA), defined originally in MNI space. This parcellation scheme was built by transforming the standard HOA of 110 anatomically defined regions. The advantage of the standard HOA over other atlases is that the results are easier to link to the cognitive neuroscience literature on various brain regions than with other atlases, whereas the disadvantage is that the ROIs could have different sizes. This issue was addressed by the sub-parcellation into 471 ROIs (Fornito et al., 2010), which has four important strengths for the current study: (1) the cortical surface areas (∼2 cm3) are consistent with the estimated area of maximal effect for standard rTMS targeting (Deng et al., 2013) (in contrast, most widely used parcellations that are comprised of much larger frontal ROIs); (2) it provides better coverage of many subcortical structures, such as the medial temporal lobes known to facilitate memory functioning; (3) it produces nearly isometric volume across all ROIs, which reduces bias in both DWI tractography and in the estimation of BOLD-related FM; and (4) it results in a level of dimensionality comparable to other investigations of cortical control (Gu et al., 2015; Medaglia et al., 2018). The T1-weighted image was then nonlinearly normalized to the ICBM152 template in MNI space using fMRIB's Non-linear Image Registration Tool (FNIRT, FSL, www.fmrib.ox.ac.uk/fsl/). The inverse transformations were applied to the HOA atlas in the MNI space, resulting in native-T1-space GM parcellations for each subject. T1-weighted images were then coregistered to native diffusion space using the subjects' unweighted diffusion image as a target; this transformation matrix was then applied to the GM parcellations above, using FSL's FLIRT linear registration tool, resulting in a native-diffusion-space parcellation for each subject. For each subject, the peak fMRI activation coordinates were transformed back into MNI coordinates. The atlas ROI corresponding to the peak coordinates was then extracted as the individualized stimulated site ROI; the same ROI was used to define the MC values at the targeted site.
Brain states and controllability measures
We adopted a mathematical definition of controllability from Pasqualetti et al. (2014) where brain system dynamics are modeled with a linear discrete-time and time-invariant model as follows:
Here, the N-by-N weighted adjacency matrix A describes the structural connectivity between each pair of nodes (N = 471). The N-by-1 column vector 
describes the brain state (i.e., functional activities in all brain regions). The N-by-1 column vector 
indicates external input, where the only non-zero entry 
denotes rTMS applied at brain region i. While modeling neural dynamics that occurs in a highly complex system, such as the brain would require nonlinear modeling, a necessary simplification of this approach is to use canonical linear dynamics (Kenett et al., 2018), since the controllability of nonlinear systems is structurally similar to the ones of linear systems. This approach is now commonly used (Liu et al., 2011), and the predictions of this structurally derived framework show a strong correlation to in vivo output of observed brain state dynamics (e.g., Kenett et al., 2018; Stiso et al., 2019). To establish a link between actual brain state and this system dynamic model, we may represent the brain state using the eigenmodes of the brain network system as follows:
Here, 
represents the unit-length eigenvector (eigenmode) with corresponding eigenvalue 
, both of which are intrinsic to the structural brain network A. 
denotes the loadings of each eigenmode, and larger 
indicates greater contribution of 
to brain state x. Combining (1) and (2), without external input, the brain state at t can be expressed as 
, which theoretically predicts that an eigenmode 
with larger |
would be more robust (i.e., slow-decaying) and easier to reach, whereas another eigenmode 
with smaller |
would be less robust and harder to reach (Stiso et al., 2019). In addition, a parallel interpretation of the eigenmodes comes from previous work by Medaglia et al. (2018), which suggests that easy-to-reach eigenmodes resemble possible activation patterns that are aligned to the anatomic brain network, whereas difficult-to-reach eigenmodes resemble activation patterns that are liberal (i.e., misaligned to the network). We offer an intuitive illustration to such interpretation by showing the spatial patterns of several eigenmodes in one participant (see Fig. 3B).
To test whether this theoretical model predicts the influence of task difficulty on brain states, we investigated how the brain states consist of the eigenmodes differently because of the change in task difficulty for each participant, according to Equation 2. Here, A, {
}, and {
} are participant-wise measurements, and {
} are participant- and condition-wise measurements (see Fig. 3). To quantify the contribution of each eigenmode state to the brain state, we scaled the loadings {
} with the reciprocal of 2-norm of the brain state, such that 
. Greater 
indicates greater contribution of the eigenmode 
to the brain state.
The eigenmodes of the brain network are directly linked to the MC measure (Kalman, 1964), interpreted here and elsewhere (Gu et al., 2015) as the ability of one region to steer the brain toward difficult-to-reach states. The MC of brain region i is defined as follows:
where 
denotes the i-th entry in 
's corresponding eigenmode 
, that is, the value of this eigenmode at the i-th brain region. By extension of the Popov-Belovich-Hautus test (Lathi, 2005), the controllability of thej-th eigenmode in region i is proportional to 
; in other words, this eigenmode would be uncontrollable from region i if it does not contribute to the value in this region. Meanwhile, the coefficient 
favors the controllability of the respective eigenmode state with low robustness (i.e., harder to reach). Therefore, as a composite measurement, the MC theoretically predicts whether stimulating a brain region can movethe brain into difficult-to-reach states. Last, while baseline differences between MC values across subjects were modest (mean MC = 0.994; SD = 0.005); the rank of MC within each individual was used rather than the raw MC values. This allowed us to: (1) remain consistent with previous investigations of MC (Gu et al., 2015), (2) control for cross-subject noise, and (3) indicate to what extent rTMS has optimally targeted at the difficult-to-reach state for each participant. Furthermore, a principle reason for using rank instead of absolute value is that relative value (e.g., the coordinate of a peak voxel, an ROI that ranks the highest in a certain measure) is more helpful in determining optimal TMS target in individual subjects, and such practice is consistent with a few recent TMS studies using peak values of neuroimaging measures for individualized targeting (Beynel et al., 2019).
In this study, since our rTMS effects were found only in the most difficult task conditions, we focused our analyses on MC. Nonetheless, other controllability measures are often discussed in the context of neuroscience-based applications of network control theory. For example, average controllability purports to describe the capacity of a network node to move the brain into easy-to-reach states with a little energy, and is often discussed in the context of resting-state-based analysis since they have been mainly reported in the default mode network (Gu et al., 2015). Boundary controllability is purported to control the dynamics between different networks in hierarchical modular networks (Pasqualetti et al., 2014). This measure therefore describes the capacity of a network node to integrate or segregate the information from different networks, and would be associated with integration of sensory-motor areas that underlie creativity processes (Kenett et al., 2018) as well as into cognitive functions (Medaglia et al., 2018). Such approaches therefore lie orthogonal to our principle goal relating MC to difficulty-to-reach brain states characterized by increasing WM load.
Mediation analysis
Mediation analysis was performed with structural equation modeling, using the lavaan package (version 0.5) (Rosseel, 2012) in R version 3.3.3 (R Core Team, 2016). The standard errors in the mediation model were estimated using robust structural equation modeling option. In addition, age, gender, and starting set size were modeled as three independent covariates to eliminate potential confounding effect to each of the three variables of interest (i.e., rTMS effect, MC, and FM). To ensure the reliability of the mediation analysis, we repeated the analysis on all possible alternative mediation models, where the directions of mediation and the mediator were reversed.
Results
Modal controllability, Functional modulation, and eigenmode analysis
DWI-based streamline counts were used to construct structural networks for each participant, and MC values were computed across 471 cortical sites. Figure 2A illustrates the rank of group mean MC values for all regions across the brain (lighter warm colors representing greater MC). While the frontal regions stimulated in the current study (Fig. 2A, inset; Fig. 4) were not selected on the basis of their MC value, the global distribution of MC values is nonetheless consistent with earlier investigations based on other cortical parcellations (Gu et al., 2015; Khambhati et al., 2018; Medaglia et al., 2018). Also consistent with Gu et al. (2015), regions with higher MC were characterized by a lower structural connectivity strength with the rest of the brain (r(470) = –0.88; Fig. 2C).
MC and FM and their correlation across the brain. (A) Mean MC and (B) FM across all ROIs. Inset, Stimulation site overlap (see Fig. 3). As expected, MC showed a consistent linear relationship with node strength (C). While frontoparietal regions appear to show some overlap in both MC and FM associated with increasing WM load, we observed no global relationship between these measures (D).
The current analysis focused on whole-brain evidence for the importance of functional transitions between active brain states, based on the modulation across levels of WM difficulty. Targeting focused on the cortical ROI, which demonstrated the greatest FM in WM load within each individual (for full imaging results, see Davis et al., 2018); thus, while the group map presented below represents a general pattern across subjects, each subject was stimulated at a unique cortical site determined by their individual parametric map (see Fig. 4A). As such, a group summary of the rank value based on our univariate fMRI analysis revealed a typical pattern of BOLD activity associated with increasing WM load (lighter cool colors represent greater FM; Fig. 2B). The spatial pattern of FM is independent of the spatial pattern of MC, as evidenced by a weak correlation between the ranks of mean MC and mean FM across regions (r(470) = 0.11; Fig. 2D). Such independence between MC and functional activity further suggests the possibility that the interplay between MC, FM, and the behavioral rTMS effect should be highly specific to the rTMS-targeted region.
In order to address the biological plausibility of the state dynamic model from which MC is derived, we sought to relate the eigenmodes of the structural network to the brain states (i.e., BOLD-related activity) assessed at four levels of task difficulty. The actual brain state can be decomposed into eigenmodes, whose loadings reflect their contribution to the brain state. Figure 3A displays a clear trend indicating that theoretically easy-to-reach eigenmodes contributed primarily to the brain states. With the example eigenmodes provided in Figure 3B, this suggests that whole-brain activation tends to follow connectivity constraints (Medaglia et al., 2018). We then grouped the eigenmodes into three subspaces according to their absolute eigenvalues, such that Subspace 1 contained the most easy-to-reach eigenmodes, Subspace 3 contained the most difficult-to-reach eigenmodes, Subspace 2 contained intermediate eigenmodes, and each subspace contributed almost equally to the actual brain state (∼33%) (collapsed across subjects and task difficulty). A Task Difficulty-by-Subspace repeated-measures ANOVA exhibited a trend interaction (Task Difficulty × Subspace: F(6,132) = 1.99, p = 0.07), suggesting that the activation patterns at different task difficulty level may load onto the three subspaces to different degrees. Intriguingly, the total contribution of Subspace 1 eigenmodes to the brain states was greater in the easiest task condition (paired t test for Very Easy vs Hard: t(22) = 1.87; p = 0.07), although this was also a trend. In contrast, the total contribution of Subspace 3 eigenmodes to the brain states was reliably greater for the most difficult task condition (main effect of Task Difficulty: F(3,66) = 2.94, p = 0.04; paired t test for Very Easy vs Hard: t(22) = −3.06 p = 0.01; Fig. 3C). The effect of difficulty in the most difficult-to-reach subspace was robust, as it also existed when the eigenmodes were grouped into 4 or 5 equal subspaces with regard to their contribution to the brain states (F(3,66) = 2.92, p = 0.04; F(3,66) = 2.75, p = 0.05, respectively). These results align with the conceptualization of MC as quantifying controllability of difficult-to-reach states, since the eigenmodes, as the constituents of MC, demonstrated differential association with brain states in varying levels of task difficulty.
A, The contribution of eigenmodes to brain state. The eigenmodes are arranged according to the order of their absolute eigenvalue (left to right, high → low); examples of the contribution of individual regions to these eigenmodes are displayed below the x axis (although these are distinct for each subject's structural matrix). B, Example of the spatial distribution of loadings for various eigenmodes v (here, the 1st, 4th, 176th, and 471st eigenmodes) for a representative subject. The spatial continuity diminishes as the eigenmode moves from early (easy-to-reach) to late (difficult-to-reach). C, Average factor loadings for eigenmodes grouped into three ordered subspaces demonstrate a clear relationship with task difficulty, such that eigenmodes in early subspaces contributed more to the brain state in easiest conditions (i.e., low WM load), whereas the subspaces characterizing the later eigenmodes (reflecting more difficulty-to-reach states) contributed more to the brain state in more difficult task conditions (i.e., high WM load). **p < 0.005. *p < 0.05. +p < 0.1.
Effects of rTMS on behavior
Having characterized the spatial distribution of MC measures, the distribution of functional states associated with WM function, and the eigenmode state synthesis of these two systems, we next sought to address the relevance of these metrics in predicting rTMS-related behavioral outcomes. To clarify the behavioral benefit associated with 5 Hz online rTMS, a repeated-measures ANOVA conducted on WM accuracy to evaluate effects of difficulty (Easy, Medium, Hard) and stimulation type (Active, Sham) revealed a significant interaction between these two factors (F(2,44) = 3.04, p = 0.05). Decomposition of this interaction demonstrated that subjects were significantly more accurate with active compared with sham rTMS, but only in the most difficult task condition (Fig. 4B). Accuracy in this hardest condition improved from 49.8 ± 15.5% with sham stimulation to 53.3 ± 17.5% with active rTMS. Subsequent analyses used a within-subject metric to describe the percent increase in WM accuracy for active relative to sham rTMS, hereafter referred to as the “rTMS effect” (i.e., [Active rTMS Accuracy – Sham Accuracy]/Sham Accuracy × 100); subjects in the current study demonstrated an average percentage increase of 7.32% in the most difficult condition. While the effect size was small (Cohen's d = 0.21), the number of subjects showing performance improvement with active rTMS is higher in the hard condition (n = 17) compared with the easy and medium conditions (n = 10 for both easy and medium difficulties). One participant has been excluded because their improvement was beyond 2 SDs from the group mean. More details on initial findings from this dataset can be found in Beynel et al. (2019).
Stimulation sites and rTMS effect across all subjects. A, Overlap of stimulation targets for all subjects, based on region with the highest BOLD activity during a WM localizer task. B, Bean plots represent the effects of online rTMS on WM accuracy (percentage increase in accuracy, e.g., Active – Sham/Sham × 20), with significant rTMS effect only in the hardest difficulty level (*t22 = 2.19, p = 0.013) (for full analysis, see Beynel et al., 2019). Within each level, tick marks represent the range across subjects in the benefit in WM function because of TMS.
Relationships between modal controllability, functional modulation, and behavior
To test the hypothesis that exogenous input to a node with high MC will facilitate brain-state transitions to difficult-to-reach states, and hence a greater rTMS-induced effect, we investigated therelationship between MC, FM, and the benefit of rTMS in the most difficult WM condition. Correlations examining functional and structural properties at the targeted stimulation location showed that MC and fMRI-based FM at the stimulated site were negatively correlated (Fig. 5A; r(22) = –0.53, p = 0.01), providing evidence that additional recruitment of neural resources was impeded by the local presence of a difficult-to-reach eigenmode (i.e., high MC). Moreover, the average fMRI activation across these levels of difficulty was unrelated to MC (r(22) = 0.04, p = 0.85), confirming that this result was indeed specific to task difficulty. Second, MC was significantly and positively related (r(22) = 0.46, p = 0.03) to the benefit from rTMS in the difficult WM condition (Fig. 5B). This relationship was not present at Easy (r(22) = 0.14, p = 0.51) or Medium (r(22) = −0.05, p = 0.81) difficulty levels, confirming a selective rTMS effect only at a difficult-to-reach cognitive state. Finally, the FM associated with difficulty at the stimulated ROI was also significantly related to the rTMS effect (Fig. 5C; r(22) = –0.49, p = 0.02). This result supports the theory that rTMS effects depend on brain state at the time of stimulation (Silvanto and Pascual-Leone, 2008), and suggests an interpretation that the benefit of rTMS is stronger when the stimulated region faces difficulty in effectively modulating its own state.
Correlations between FM, rTMS effect, and MC. A, Relationship between MC and FM (i.e., the parametric increase in BOLD-related activity associated with increasing WM load). B, The rTMS effect (i.e., the percentage increase in WM accuracy after 5 Hz online rTMS). C, Correlation between FM and the rTMS effect. All correlations represent values derived from the stimulated ROI (see Fig. 4).
Finally, a mediation analysis was performed to test the hypothesis that the MC-rTMS relationship was mediated by the FM associated with more difficult trials (Fig. 6). This analysis yielded a significant total effect (TMS ∼ MC: zc = 2.43, p = 0.03), as well as significant relationships between the predictor and mediator (TMS ∼ FM: za = –3.36, p = 0.0008) and the mediator and outcome variable (FM ∼ MC: zb = –2.69, p = 0.007). Crucially, a significant mediation by FM (zab = 2.12, p = 0.04) explained a major proportion of the total effect (64%), resulting in a nonsignificant direct contribution of the rTMS effect to MC (zc' = 0.68, p = 0.53). To ensure the reliability and specificity of the mediation by FM, we repeated the analysis on all possible alternative mediation models (i.e., reversing the direction and order of the three principle variables in Fig. 6), where the directions of mediation and the mediator were changed. Standard practice in the assessment of mediation effects is to reverse the directions of causal inference to determine the selectivity of effect; in our data, no strong mediation effects were observed among all alternative mediation models, although a model in which MC was the predictor and TMS effect was the outcome was the closest to reaching significance (zab = 1.71, p = 0.058). As such, these results help to confirm that our observations are best explained by the model in Figure 6, as functional activity mediating the relationship between rTMS effect and structural controllability.
Mediation analysis investigating the capacity for fMRI activation to explain the relationship between MC and rTMS effect. Age, gender, and baseline task performance were controlled as confounding covariates (Cov).
Discussion
The current study sought to assess the validity of network control theory, and specifically the element of MC, in characterizing whether regional control properties may predict the modulatory effects of rTMS. In this study, rTMS has been found to improve performance only in the most difficult condition of our task, and the effect size was small. While this task-difficulty specificity has already been observed in former studies (Luber et al., 2007; Violante et al., 2017), a potential way to improve rTMS effect size is to find which factors are associated with stronger effect. We used fMRI and DWI to investigate the functional and structural information associated with brain-state transitions, and rTMS to test whether providing exogenous perturbation to a specific node could help the brain transition to difficult-to-reach states. The notion that brain activity could be controlled by exogenous input to a single cortical node is still highly debated (Tu et al., 2018; Pasqualetti et al., 2019), and only a handful of in vivo studies have investigated its validity in human models (Medaglia et al., 2018; Stiso et al., 2019). The present analysis sought both to test these assumptions by examining the relationships between empirical brain and behavioral data, as well as to use those empirical relationships to clarify the expectations of what the MC measure can predict. We discuss these results and the implications of these findings below.
At its heart, network controllability theory assumes that the state of a system at a given time depends on the functional state of the network, the structural links connecting the nodes, and any additional perturbation (Ruths and Ruths, 2014). As such, given the metabolic costs, difficult-to-reach states would only be expressed during demanding cognitive states, making it valuable to promote these states through exogenous means. We examined the interactions between MC, functional brain states, and behavior consequences resulting from neuromodulation, and found significant relationships between these measures. This can be interpreted as the cost associated with transitions needed to reach a cognitive state, which was further resolved in our single mediation model. These results help to show that brain network information, based on control theory, is supported by both functional and neuromodulatory data, and helps support the use of a new family of targeting approaches based on the dynamic control of global brain states.
The principle finding from this study was the observation that MC was reliably associated with the FM of BOLD activity associated with parametric increases in WM load (Fig. 5A), and that the global distribution of these functional patterns aligned well with the higher, more disordered eigenmodes (i.e., moredifficult-to-reach states; Fig. 3) as defined by the spectral eigenanalysis of the structural connectome. This latter finding is particularly important because it directly relates a decomposition of the MC measure (a linear combination of the eigenmodes for a given region) to task-related functional imaging data at the level of an individual. As such, the clustering of task activations representing more difficult WM task conditions with eigenmodes of a higher value (and a generally more disordered distribution across the cortex) bears direct relevance for the inference that MC relates some capacity of the network to “transition to difficult-to-reach states” (Gu et al., 2015; Pasqualetti et al., 2019). Furthermore, this finding suggests a benefit in applying more advanced network science measures to structural and functional imaging data to individualize TMS targets to enhance the efficacy of stimulation. Nonetheless, more basic graph measures, such as node strength and overall network density, must also be considered, especially when exploring TMS targeting options in different patient populations or age groups in which these more basic graph measures vary systematically.
A second important finding from this study was that MC at the stimulated site was reliably predicted by the behavioral benefits caused by rTMS in the WM task. This rTMS benefit was seen as increases in response accuracy and was only found in the most difficult task condition (Fig. 4B). Correspondingly, the association between MC and this performance benefit was selective only to this difficult condition (Fig. 5B). This result, therefore, supports existing computational models of network dynamics based on structural connectivity information, which purport to describe the energy associated with transitioning between brain states.
Predicting the effect of brain stimulation by using brain network topology is of theoretical and applied significance, and the current results provide a mechanistic application of these principles as set forth by Pasqualetti et al. (2014). According to these authors, the brain constitutes a controllable network that can be moved into different states with input perturbation, such as brain stimulation. Recent theoretical (Muldoon et al., 2016) and empirical work (Betzel et al., 2016) addressing neuronal control theory have focused on how the brain transitions from one state of activation to another, and what anatomic substrates facilitate these transitions (Stiso et al., 2019). Within such a framework, controllability is the mathematical formulation of how a system can move through state-space along a desired trajectory (Gu et al., 2015). A recent TMS study provided behavioral evidence supporting the theory that the controllability is a factor accounting for the influence of brain stimulation (Medaglia et al., 2018); the authors observed that a median split of subjects based on MC values identified subtle performance differences in a language production task, providing some support for the conclusion that MC may influence task-based response to rTMS. Nonetheless, the observed effects did not supply evidence that higher MC allows the brain to more easily traverse functional brain states, as the influence of MC was found only in the easiest language task condition. Furthermore, this condition did not show a significant benefit of active rTMS over sham, and so it is difficult to draw conclusions concerning the functional relevance of MC in that study.
Finally, we found that the core relationship between rTMS effect and MC was significantly mediated by the FM associated with increasing difficulty. The mediation analysis reported above adds to this finding to suggest MC reflects not only the structural value, but also the functional capacity for a given node to modify the derived eigenmatrix. This result helps to advance the notion that all three metrics examined here (the behavioral benefit from rTMS, FM, and MC) reflect distinct, but related, constructs of “difficulty.” In the current approach, we posited that the networks identified by parametric changes in functional activation can be influenced by modulation of a single site. Following Gu et al. (2015), we treated this influence as an eigenvalue problem, such that regions with high MC have the greatest influence on less persistent eigenmodes. MC would thus provide a measure of the capacity to drive the dynamics of a network toward difficult-to-reach configurations. This mathematical formulation helps to explain the more basic mediation finding in that, while MC of the network reflects the capacity for a node to shift to a new brain state, FM reflects the degree to which regions effectively do traverse such a brain state. These cooperative contributions in turn lead to greater performance improvements after rTMS, and presumably a successful shift in brain state. Even so, we agree with Tu et al. (2018) that this is an interpretation at the level of regions and may not capture multivariate interactions inherent in a broader, and perhaps less tenable, question: is the brain controllable? The current analysis nonetheless seeks to ground this interpretation by applying the control theory metrics to real human brain data. Notably, the direction of the best-fitting model suggested MC as the outcome and TMS effects as a predictor, whereas the mediation effect (zab) failed to reach significance in the reverse model (zab = 1.71, p = 0.058). Nonetheless, given the lack of longitudinal data or postscan responses in this study, it is difficult to build inferences on this asymmetry. Future studies using both pre- and post-TMS neuroimaging may help to disambiguate this causality, by assessing measures of longitudinal change in both FM and structural connectivity, which may shift the underlying controllability topology. Of note, the critical triad relationships between MC, FM, and TMS effects were not present when using the more basic measurement of network strength (i.e., weighted degree). While strongly correlated (Fig. 2C), MC differs from strength in that strength is uncorrelated with parametric activation (r = −0.035, p = 0.88), and furthermore lacks a strong theoretical basis: strength values merely represent the sum of all weighted connections for a given node. We therefore see MC as a refinement of the streamline values (of which strength is a simple sum), and endeavored here to show the utility of this particular structural network metric.
The approach used here has clear implications for theindividualization of TMS targeting, and advances targeting approaches that have improved from scalp-based, to anatomic, to functionally guided techniques (Herwig et al., 2001). Both clinical and basic science targeting approaches for TMS should consider the consequences of regional brain stimulation in manipulating more global functional and structural networks connected to the stimulation site (J. X. Wang and Voss, 2015; Davis et al., 2017; W. C. Wang et al., 2018). When regional stimulation is applied to models that differ only in their structural connectivity, the resulting activity patterns differ between individuals (Jbabdi et al., 2015; Spiegler et al., 2016), suggesting that functional responses are sensitive to subject-specific differences in white matter structure. Computational models informed by both forms of brain information are therefore necessary to provide realistic, individualized models of immediate TMS outcomes. The utility of the MC value is clear: when using rTMS to improve performance in one specific cognitive task, it is straightforward to derive an fMRI target from this task; however, when it comes to clinical application, the choice of the optimal task to use becomes more complicated. Furthermore, as task-related fMRI activations suffer from intrasubject and intersubject variability and depend largely on task performance, using MC, based on structural acquisition, could represent a way to derive a more efficient TMS target for clinical applications. Our findings relating the structural and functional components of a controllable system suggest that MC may be used as an alternative when functional localizers for TMS are unfeasible.
In conclusion, by combining information from structural network topology, whole-brain functional activity associated with a WM task, and brain stimulation, we confirmed the validity of network controllability by showing a strong interdependence of these measures. This work helps to outline a clear technique for integrating structural network topology and functional activity to inform network-based approaches in selecting an individualized specific cortical target for neurostimulation that could help increase rTMS efficacy.
Footnotes
This work was supported by National Institute of Aging Grants U01AG050618 and K01AG053539, and National Institute of Mental Health Intramural Research Program. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. We thank Duy Nguyen, Alexandra Brito, and Joyce Wang for assistance with data collection.
The authors declare no competing financial interests.
- Correspondence should be addressed to Simon W. Davis at simon.davis{at}duke.edu
