Network-level macroscale structural connectivity predicts propagation of transcranial magnetic stimulation.

Information processing in the brain is mediated by structural white matter pathways and is highly dependent on topological brain properties. Here we combined transcranial magnetic stimulation (TMS) with high-density electroencephalography (EEG) and Diffusion Weighted Imaging (DWI), specifically looking at macroscale connectivity to understand whether regional, network-level or whole-brain structural properties are more responsible for stimulus propagation. Neuronavigated TMS pulses were delivered over two individually defined nodes of the default mode (DMN) and dorsal attention (DAN) networks in a group of healthy subjects, with test-retest reliability assessed 1-month apart. TMS-evoked activity was predicted by the modularity and structural integrity of the stimulated network rather than the targeted region(s) or the whole-brain connectivity, suggesting network-level structural connectivity as more relevant than local and global brain properties in shaping TMS signal propagation. The importance of network structural connectome was unveiled only by evoked activity, but not resting-state data. Future clinicals interventions might enhance target engagement by adopting DWI-guided, network-focused TMS.

Introduction

In the last two decades, transcranial magnetic stimulation (TMS) has been widely employed to study functional connectivity and plasticity in the human cortex ( Ziemann, 2004 ), as well as transiently manipulate brain spontaneous dynamics ( Santarnecchi et al., 2018 ) and human cognition ( Momi et al., 2019 ). Indeed, TMS allows a relatively focal stimulation (~2–2.5 cm diameter) of the human brain using very brief duration magnetic waves administered by an electromagnetic coil positioned on the scalp ( Hallett, 2007 ). However, the effect of the TMS pulse is not limited to the underlying stimulated tissue, but it can possibly reverberate within interconnected cortical networks, inducing synchronization of distant cortical areas and changes in functional connectivity dynamics ( Di Lazzaro, 2004; Massimini et al., 2005; Siebner et al., 2001 ). A compelling way to quantify these cortico-cortical connectivity patterns is to combine TMS with electroencephalography (EEG) ( Voineskos et al., 2010 ), thus looking at millisecond-level brain activity propagation throughout the brain and potentially capture individual variability in response to external perturbation, as well as disease-specific alterations of brain activity ( Benussi et al., 2020; Massimini et al., 2012 ). In a recent TMS-EEG study ( Ozdemir et al., 2020 ), we have stimulated two resting-state networks (RSNs) ( Fox et al., 2005 ) (i.e. Default Mode Network, DMN, and Dorsal Attention Network, DAN) showing how the external-induced signal follows network dynamics usually measured via functional magnetic resonance imaging (fMRI). Changes in RSNs connectivity emerge when pathological states arise ( Anderson et al., 2011 ), suggesting resting-state fMRI (rs-fMRI) as a potential biomarker of disease symptomatology and progression ( Ferguson et al., 2019 ). The possibility of indexing network activity via TMS-EEG would open new possibilities for disease tracking and early detection of network-level dysfunction, but more information about network engagement and propagation are needed.

While the understanding of TMS propagation pathways and principles is still incomplete, the vast majority of studies suggest that TMS pulses propagate along white matter tracts ( O’Shea et al., 2008 ), with the conductivity of white matter bundles potentially shaping the propagation of action potentials ( McCann et al., 2019 ). In this context, individual differences in the complexity and organization of the structural connectome might play a critical role in determining the response to TMS, and such information could be used to guide TMS targeting and achieve maximal target engagement. Here we used image-guided TMS-EEG to selectively perturb one of two neighboring nodes of the DMN and DAN ( Fox et al., 2005 ) to then extrapolate individual network engagement values by computing the area under the curve (AUC) for significant post-TMS time points, thus providing higher temporal resolution. Afterwards, we correlated individual response to TMS with macroscale structural properties of the brain. In particular, we combined TMS, EEG and Diffusion Weighted Imaging (DWI) to test an ad-hoc five-layer hierarchical framework ranging from local to network-level up to whole-brain structural connectivity, therefore moving along a gradient including local, regional patterns on one side and very unspecific whole-brain connections on the other end. Specifically, the response to TMS targeting the DMN and DAN was correlated with (a) the structural connectivity of the stimulated region with the rest of the brain, (b) the connectivity of the stimulated region with the other nodes of the same network, (c) the intrinsic connectivity of the targeted network, (d) the modularity of the brain as an index of between and within networks connectivity ( Sporns and Betzel, 2016 ), and (e) the overall average structural connectivity of the entire brain. Given the relevance of local activation induced by TMS ( Krieg et al., 2013 ), we hypothesized that layers of the hierarchical model related to local connectivity (e.g. a and b) would constitute important predictors of the response to TMS, followed by network-level ones. To further investigate the specificity of the results, the structural connectivity of non-stimulated RSNs was also correlated as a control condition. Additionally, to verify whether measuring the response to TMS-based perturbation provides an advantage in capturing the relevance of individual structural connectivity properties, the same analysis was repeated by looking at network activity during resting-state EEG recording. Finally, considering the quest for data reproducibility, the same analyses were repeated on data collected on the same sample of healthy individuals across two separate study visits one month apart. For details on the study design please see Fig. 1 and the Methods section of the manuscript.

Fig. 1.

Study design and conceptual framework. (A) TMS targets were individualized based on resting-state fMRI data. A high degree of variability in functional connectivity of the TMS targets was present. (B) fMRI-guided TMS was applied to two neighboring parietal nodes corresponding to the DMN and DAN. Anatomical MRI were used for the neuronavigation of the TMS spots while hd-EEG with 64 channels was simultaneously recorded. (C) TMS-induced electric field was modelled with SimNIBS (Thielscher et al., 2015). The final map was overlapped with the 7 Network parcellation (Schaefer et al., 2018) in order to assure network engagement specificity. (D) The EEG signal was projected at source level using dynamic statistical parametric mapping (dSPM) and constraining source dipoles to the cortical surface. The RSNs time series were extracted for both DMN and DAN. The raw time series were first rectified (Cheng et al., 2013 ) and then a baseline bootstrapping procedure (Lv et al., 2007 ) was applied. Then, 1000 permutation t -test were performed in which the surrogated post-TMS vs pre-TMS difference was computed after each iteration and statistically compared with the real difference ( Pernet et al., 2015 ). Finally, the cluster threshold was determined as the 95th percentile of the cluster’s surrogate distribution and the area under the curve (AUC) of the significant clusters was extracted. (E) The individual whole brain structural connectome was computed. A five-layers hierarchical model was created where several connectivity metrics were extracted, ranging from local to whole-brain measures. Note: V/m: Volt per meter; DAN: Dorsal attention network; DMN: Default mode network. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Material and methods

Participants

The study was approved by the Institutional Review Board of the Beth Israel Deaconess Medical Center. Each participant provided written informed consent conformed to the Declaration of Helsinki and was compensated for the entire study. Twenty-one right-handed ( Oldfield, 1971 ) healthy volunteers (mean age = 32 ± 10 years, ranging from 19 to 49 years) with normal neurological and psychiatric evaluation and no history of drugs acting on the central nervous system were recruited through flyers and on-line advertisement. Participants carried out a pre-TMS assessment comprehensive of structural (e.g. T1-weighted (T1w) and DWI) and functional (e.g. fMRI) MRI. After that, two TMS visits, separated by one month, were performed where 120 single pulses were delivered in two neighboring parietal nodes corresponding to the DMN and DAN. At the beginning of each TMS visit, resting motor threshold (RMT) was identified for each participant, targeting the left motor cortex (M1) in order to define TMS intensity for each participant according to international TMS guidelines ( Rossi et al., 2020 ). Methods for data acquisition are presented in the following paragraph and have been further described in ( Ozdemir et al., 2020 ).

Conceptual framework

To explore the relationship between local/global brain features and response to network-targeted TMS, we delineated a granular, ad-hoc five-levels hierarchical model ranging from high modularity to high generality where several connectivity measures were extracted (see Fig. 1 E):

Stimulated region to brain connectivity (Stim2Brain): a metric representing the number of structural connections between the stimulated region and any other brain region. This index was calculated by averaging the number of streamlines which originate from the stimulated area and terminate in every other brain region. This metrics captures global connectivity of the stimulated region regardless of the stimulated region/network;

Stimulated region to network connectivity (Stim2Network): an index expressing the connectivity between the stimulated region and any other node of its network. This index was calculated by averaging the number of streamlines which originate from the stimulated region and terminate in any other node of the same network. This metric captures local connectivity of the stimulated region with respect to the targeted network;

Stimulated network connectivity (Network): this metric incorporates the intrinsic connectivity of the stimulated network extracted by considering the streamlines that connect each pair of network nodes. This index was calculated by averaging the number of streamlines which connect every brain region belonging to the same network. This represents the stimulated network connectivity without considering global property of the brain;

Within/between network connectivity (Modularity): a general property of the individual connectome expressing how much each structural connectivity matrix was arranged in sub-modules ( Rubinov and Sporns, 2010 ). This index was computed using the Louvain algorithm ( Blondel et al., 2008 ), implemented in the Brain Connectivity Toolbox ( Rubinov and Sporns, 2010 ) as follows:

where m is the sum of all the nodes in the network; Aij is the adjacency matrix representing the edge weight between node i and node j; k i and k j are the sum of the weights of the edges attached to nodes i and j, respectively; δ(ci,cj ) are the communities of the nodes and is 1 if nodes i and j belong to the same subset of the maximized partition of brain nodes, and 0 otherwise.

Whole-brain connectivity (Brain): the most generic level where an index of global structural connectivity was computed regardless of network differentiation. This index was calculated by averaging the number of streamlines which traverse every brain region. This level represents global properties of each individual, that could affect TMS propagation regardless of the stimulated region/network;

For each subject, the structural connectivity metrics were accordingly extracted following this model and then used to predict TMS-EEG response. Specifically, the number of streamlines within each level was correlated with the average source-level induced activity in the target network (see DWI analysis section for more details on streamlines calculation). A multiple regression analysis was implemented to identify which level of the hierarchical model better explain the TMS-EEG response. For a summary of the study design, analytical framework and structural hierarchy please see Fig. 1. For details on the specific methods see dedicated sections below.

MRI data acquisition

A T1w anatomical MRI scan was obtained in all participants and used for neuronavigation. MRI data was acquired on a 3T scanner (GE Health-care, Ltd., United Kingdom) using a 3D spoiled gradient echo sequence: 166 axial-oriented slices for whole-brain coverage; 240-mm isotropic field-of-view; 0.937-mm × 0.937-mm × 1-mm native resolution; flip angle = 15°; TE/TR ≥ 2.9/6.9 ms; duration ≥ 432 s. DWI sequence were also acquired using a single-shot echo planar imaging (slices = 71; matrix size = 256 × 256 × 71; voxel size = 0.8 mm × 0.8 × 2.2; repetition time = 8500 ms, time echo = 79 ms; 30 non-colinear directions, b-value = 1000s/mm 2 ).

DWI data preprocessing and analysis

A customize pipeline running in Ubuntu 18.04 LTS was used for the preprocessing of DWI images using tools in FMRIB Software Library (FSL 5.0.3; www.fmrib.ox.ac.uk/fsl ) ( Jenkinson et al., 2012 ), MRtrix3 ( mrtrix.readthedocs.io/en/latest/ ) ( Tournier et al., 2012 ), FreeSurfer ( Fischl et al., 2004 ) and ANTs ( stnava.github.io/ANTs/ ) ( Avants et al., 2011 ). All images were denoised ( Veraart et al., 2016 ), preprocessed via FSL’s EDDY ( Andersson and Sotiropoulos, 2016 ), and bias field corrected ( Zhang et al., 2001 ). The response function for a single fiber population was estimated using spherical deconvolution Tournier algorithm ( Tournier et al., 2007 ).

Simultaneously, the T1w images were coregistered to the b0 volume and then segmented using FAST algorithm ( Zhang et al., 2001 ). Following the anatomically constrained tractography was employed to generate the initial tractogram with 50 million streamlines using second-order integration over fiber orientation distributions ( Tournier et al., 2010 ). Then, spherical-deconvolution Informed Filtering of Tractograms (SIFT2) methodology ( Smith et al., 2015 ) was applied in order to provide more biologically accurate measures of fiber connectivity. The Schaefer’s atlas ( Schaefer et al., 2018 ) which divided the brain into 100 regions and 7 Networks was then mapped to the individual’s FreeSurfer parcellation and then used to construct the final structural connectome calculating the number of estimated tracts between any two brain regions. The final connectivity matrices were then normalized based on the size of each ROIs ( Bonilha et al., 2015 ). For further details on grand mean average structural connectome please see Supplementary Results and Figure S5.

TMS

TMS was delivered using a figure-of-eight shaped coil with dynamic fluid cooling (Magspro 75 mm cool B-65, Magpro A/S., Denmark) attached to a MagPro X-100 stimulator (MagVenture A/S, Denmark). Individual high-resolution T1w images were imported into the Brainsight ™ TMS Frameless Navigation system (Rogue Research Inc., Montreal, Canada), and co-registered to digitized anatomical landmarks for online monitoring of coil positioning. Motor evoked potentials (MEPs) were recorded from the right first dorsal interosseous (FDI) and the abductor pollicis brevis (APB) muscles. Ag-AgCl surface electrode-pairs placed on the belly and tendon of the muscles and a ground on the right ulnar styloid process. EMG data were amplified and digitized using a Power-lab 4/25T data acquisition system (ADInstruments) at a sampling rate of 4000 Hz (bandpass filtered at 10 Hz to 2000 Hz). EMG signals were continuously streamed by using LabChart software (LabChart 8.0) to monitor MEPs and epochs were recorded with a 150 ms window length covering from 50 ms before to 100 ms after TMS pulse.

EEG

Whole scalp 64-channel EEG data was collected with a TMS-compatible amplifier system (actiCHamp system, Brain Products GmbH, Munich, Germany) and labeled in accordance with the extended 10–20 international system. EEG data were online referenced to Fp1 electrode. Electrode impedances were maintained below 5k Ω at a sampling rate of 1000 Hz. EEG signals were digitized using a BrainCHamp DC amplifier and linked to BrainVision Recorder software (version 1.21) for online monitoring. Digitized EEG electrode locations on the scalp are also co-registered to individual MRI scans using Brainsight ™ TMS Frameless Navigation system.

TMS targets

In order to identify individualized TMS targets, group-level resting-state functional networks maps were used, based on a 7 networks parcellation covering cortical and subcortical structures ( Yeo et al., 2011 ). The 7 networks correspond to visual (VIS), somatosensory (SM), limbic (LIM), dorsal attention (DAN), anterior salience (AS), default mode (DMN), and fronto-parietal (FPN) RSNs. Confidence maps for each RSN were used, representing the confidence of each vertex belonging to its assigned network across a sample of 1000 healthy subjects (expressed as valued between −1 and 1), with larger values indicating higher confidence. By using group-level functional parcellations and confidence maps, we were able to target the most consistent and reliable regions within each network, therefore increasing the generalizability of TMS-EEG findings. We first projected the 7-network functional cortical atlas and the confidence maps onto subject’s cortical surface using the spherical registration implemented in Freesurfer software ( Fig. 1 A). The resulting maps were then resampled to native structural T1w MRIs. Voxels within each network were weighted by the confidence map and the voxels with the highest confidence in angular gyrus and superior parietal in the right hemisphere were chosen for DMN and DAN stimulations respectively.

TMS-EEG data collection

TMS stimulation intensity was determined based on individual resting motor threshold (RMT), defying as the lowest stimulation intensity necessary to elicit an MEP of at least 50 μV in 5 out of 10 trials ( Rossini et al., 2015 ). The hotspot of stimulation was therefore defined as based on the cortical hand region where MEPs were the largest and more consistent, as recorded in the first dorsal interosseous (FDI) ( Rothwell et al., 1999 ). Throughout the stimulation visit, participants wore earplugs to protect their hearing ( Rossi et al., 2009 ), and auditory white noise masking was used to minimize the impact of the TMS click ( ter Braack et al., 2015 ). A thin layer of foam was placed under the TMS coil to minimize somatosensory contamination of the TMS-evoked EEG potentials. A total of 120 single TMS pulses were delivered to each stimulation target (DMN target site within the angular gyrus, DAN target site within the superior parietal lobule)at an intensity of 120% RMT with randomly jittered (3000–5000 ms) inter stimulus intervals over two repeated blocks each consisting of 60 trials ( Fig. 1 B). Each participant completed two identical experimental sessions 4 weeks apart.

EEG data processing

All EEG data pre-processing was performed offline using EEGLAB 14.1 ( Delorme and Makeig, 2004 ), and customized script running in Matlab R2017b (Math-Works Inc., USA). The two single blocks of 60 trials each were merged into a single block of 120 trials and then segmented into epochs of 1500 ms each (from −500 ms (pre-pulse) to 1000 ms (post-pulse)). An amplitude of the mean pre-pulse ( −500 ms to −100 ms) signal was used to perform baseline correction. Noisy channels were removed following visual inspection, with an average of 3 ± 2 channels removed out of 63. Early TMS pulse artefact was removed by performing zero-padding on a window of −2 ms to 14 ms. Noisy epochs were then rejected based on the voltage ( ≥ 100 μV), kurtosis ( ≥ 3), joint probability (single channel-based threshold ≥ 3.5sd) and visual inspection. The data were reduced into 60 components via principal component analyses (PCA) to minimize overfitting and noise components. Then, a first round of fast independent component analysis (fICA) ( Hyvärinen and Oja, 1997 ), further aimed at removing remaining early TMS-evoked and EMG artefacts (1 ± 1 component was removed; range 0–3 out of 60). The EEG data were then interpolated for previously zero-padded time window around TMS pulse using linear interpolation, band pass filtered using a forward-backward 4th order Butterworth filter from 1 to 100 Hz, notch filtered between 57 and 63 Hz, and referenced to global average. Data were further reduced into 57 dimensions using a second PCA followed by a second round of fICA to remove any remaining artefact ( Rogasch et al., 2017 ) including eye movement/blink, muscle noise (EMG), single electrode noise, TMS evoked muscle, cardiac beats (EKG) and auditory evoked potentials (22 ± 6 components were removed; range 18–28 out of 57). A semi-automated artefact detection algorithm incorporated into the open source TMS-EEG Signal Analyzer (TESA v0.1.0-beta; https://nigelrogasch.github.io/TESA/ ) was used during both fICA ( Rogasch et al., 2017 ). Finally, the data were low pass filtered with a 4th order Butterworth filter at 50 Hz and previously removed channels were spherically interpolated. Two participants dropped out before performing the re-test visit and were not included in the data analysis and two others were excluded because of bad signal to noise ratio in the EEG. For a schematic representation of the preprocessing steps please see Supplementary Figure S2.

EEG source reconstruction

All TMS evoked EEG source reconstruction was performed using Brainstorm ( Tadel et al., 2019 ). First, digitized EEG channel locations and anatomical landmarks of each subject were extracted from Brainsight ™ (nasion ‘NAS’, left pre-auricular ‘LPA’, and right pre-auricular ‘RPA’ points), and registered onto individual MRI scans in brainstorm. Next, the EEG epochs, from −500 ms to 1000 ms with respect to TMS pulse, for each TMS trial were uploaded and average epoch time series was generated for each subject. Forward model solution of neuro electric fields was performed using the open MEEG symmetric boundary element method ( Gramfort et al., 2010 ), all with default parameter settings ( Tadel et al., 2019 ). Noise covariance was estimated from individual trials using the pre TMS (from −500 ms to −100 ms) time window as baseline. The inverse model solution of the cortical sources was performed using dynamic statistical parametric mapping (dSPM) and by constraining source dipoles to the cortical surface. The resulting output of EEG source reconstruction was the current density time series for each cortical vertices.

Source-level metrics

In order to determinate the network engagement, the average current density timeseries were extracted (after flipping the sign of sources with opposite directions) from the all 7 RSNs maps projected on surface space for each individual both for DMN and DAN stimulation ( Fig. 1 D). The final timeseries were normalized (z-score) and rectified ( Cheng et al., 2013 ). Then, to assess the threshold for significance of timeseries, a bootstrap method ( Lv et al., 2007 ) which does not assume normal distribution of the observations, was applied by shuffling the time samples of pre TMS activity (from −500 ms to 0 ms). Specifically, 500 surrogated pre-stimulus timeseries were computed to obtain a maximum distribution (control for type I error) and significance level was set at p < 0.01. Finally, both a time-wise permutation testing and a cluster-based thresholding ( Pernet et al., 2015 ) as a correction for multiple comparisons were performed. Specifically, the permutation test transformed the difference between the TMS condition (from 0 ms to 500 ms) and the baseline (from −500 ms to 0 ms) into a z value with respect to a null distribution of surrogate conditions difference values, obtained by swapping condition labels at each of 1000 permutations. The resulting z-scores were thresholded at p < 0.05. With an additional 1000 iterations permutation test, a distribution of cluster sizes of contiguous significant time points under the null hypothesis of no condition difference was computed, and only clusters that exceeded the 95th percentile of this distribution were retained. Finally, the AUC was extracted from the significant clusters ( Fig. 1 D).

It is important to mention that for the EEG timeseries it was not possible to create a five-levels hierarchical model equivalent to the DWI metrics. Indeed, the current density extracted from the stimulated region does not reflect neither the Stim2Brain nor the Stim2Network level of the DWI framework. Specifically, both these metrics express the connectivity of the stimulated region with either the rest of the brain (Stim2Brain) or the network (Stim2Network) of interest. Such connectivity information is not addressed by EEG-source current density extracted from the stimulated region which simply reflect the engagement of the such region following a TMS-pulse. Therefore, in this paper we focused on the amount of network engagement following an external perturbation of a single node.

Control analyses

To control for the specificity of the stimulated network, intrinsic connectivity measures were extracted for the other RSNs by considering the streamlines that connect each network’s node. We hypothesized that TMS-evoked network engagement was only predicted by the intrinsic connectivity of the stimulated network instead of the connectivity of the other RSNs. In addition, EEG source metrics were also computed for the non-stimulated RSNs and correlated with the modularity of the structural connectome. Moreover, similarly to a recent work published by our group ( Ozdemir et al., 2020 ), AUC of the cortical activations in the resting-state period (from −500 ms to −100 ms) preceding TMS were also extracted and correlated with DWI measures. We hypothesized that structural connectome profiles only predicted TMS-evoked activity instead of conventional resting-state EEG. Finally, to assure the specificity of the network engagement, electric field induced by TMS pulses was mathematically modelled ( Fig. 1 C) and overlapped with the individual the RSNs. For further details on method and results please see Supplementary Information and Fig. 2 B&C.

Fig. 2.

Specificity of network engagement and reproducibility of TMS-EEG measures. (A) Individualized DMN and DAN targets mapped to MNI space are provided to show variability of TMS sites across individuals. (B) Average TMS-induced electric field as modelled with simNIBS ( Thielscher et al., 2015 ). The normalized electric field (Efield) distribution was thresholded considering only the 83% of the maximal E-field. For each site/session, the thresholded cluster was overlapped with the RSNs parcellation by Yeo et al. (2011) ( Schaefer et al., 2018 ) in order to quantify network engagement. (C) Quantitative spatial overlap analysis ( Dice, 1945 ) between thesholded E-field maps and stimulated networks. High overlap was found between the stimulated network and the E-field maps both for DMN (top, 91.3%) and DAN (bottom, 86.9%). (D) The percentage of network engagement as measured via EEG source analysis for each network is shown, demonstrating high propagation specificity after TMS of DMN and DAN. (E) Test-retest reliability of source-level network engagement showing high reproducibility across visits for both DMN (top, red lines) and DAN (bottom, green lines). (F) Subject’s topographical maps for visit 1 (top, magenta) and visit 2 (bottom, cyan) show the reproducibility of TMS-evoked measures. (G) Evoked activity map (left) and EEG time series for electrodes F6, C6 and P6 (right) for visit 1 (magenta line) and visit 2 (cyan line). Note: V/m: Volt per meter; AUC: Area under the curve; DAN: Dorsal attention network, SM: Sensorimotor network; VIS: Visual network; DMN: Default mode network; FPN: Fronto-parietal network; LIM: Limbic network; AS: Anterior salience network. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Results

Accuracy of TMS targeting, target engagement and reliability

A set of analyses were conducted to verify the goodness of our dataset before proceeding with the exploration of the hierarchical framework, including accuracy of targeting procedures, location of individual stimulation sites, amount of induced activity in each targeted network and reliability of such activity across sessions.

To verify the goodness of TMS targeting procedures across participants, a quantitative spatial overlap analysis (DICE coefficient) ( Dice, 1945 ) between E-field maps and stimulated networks was computed. Results suggest high targeting accuracy for both networks, with high overlap between the stimulated network and the E-field maps, both for DMN (top, 91.3%) and DAN (bottom, 86.9%) ( Fig. 2 B&C). For details on biophysical modeling please see Methods section and Supplementary Materials.

The percentage of network engagement as measured via EEG source analysis for each network was also calculated, indexing the amount of activity generated in the target network as compared to the rest of the brain/other RSNs. Importantly, this metric includes permutation testing and cluster correction ( Pernet et al., 2015 ) and partially differs from what was previously published by our group ( Ozdemir et al., 2020 ), where only the stimulated network sources were extracted without taking into account other RNSs. Stimulation of the DMN resulted in greater activity of the DMN network (30.60%) compared all other RSNs (DAN = 10.29%, VIS = 11.01%, SM = 12.83%, AS = 12.34%, LIM = 11.40%, FPN = 11.50%). Conversely, stimulation of DAN node engaged selectively the DAN network (32.08%) more than all other RSNs (DMN = 14.14%, VIS = 12.16%, SM = 12.52%, AS = 8.95%, LIM = 9.45%, FPN = 10.67%) (see Fig. 2 D).

To test the reproducibility of TMS evoked cortical activation dynamics, TMS-evoked potentials (TEPs) time series were computed for each EEG channel across both stimulation conditions and visits. AUC extracted from significant clusters were highly reproducible ( Fig. 2 E) within each participant and across visits (DMN stimulation: Visit 1 AUC M = 242, SD = 67.93 - Visit 2 AUC M = 233.66, SD = 56.96; DAN stimulation: Visit 1 AUC M = 213.45, SD = 61.51 - Visit 2 AUC M = 206.46, SD = 49.62). For further details on test re-test reliability of TEPs please refer to Supplementary Results.

Hierarchical framework and TMS-evoked network engagement

Considering our hierarchical model, a significant positive correlation was found between the structural connectivity of the stimulated network (see level c of the hierarchical model in Methods/ Fig. 1 ) and the TMS-EEG response ( Fig. 3 A) for both DMN (R 2 = 23%, p < 0.02) and DAN (R 2 = 32%, p < 0.006). Results were highly replicable at the retest visit (DMN: R 2 = 23%, p < 0.02, DAN: R 2 = 25%, p < 0.01). Moreover, brain structural modularity (Q, level d of the hierarchical model; M = 0.65, SD = 0.2) was also significantly related to the TMS-EEG response ( Fig. 3 B). Specifically, a positive correlation was found between modularity and the magnitude of the evoked response within the stimulated network for both DAN and DMN stimulation, and this correlation was reproducible across the two TMS visits (DMN Visit 1: R2 = 18%, p = 0.04; Visit 2: R2 = 27%, p = 0.01; DAN Visit 1: R2 = 21%, p = 0.03; Visit 2: R2 = 23%, p = 0.02). Interestingly, correlation between modularity and evoked TMS-EEG activity was present only for the stimulated network; there was no significant correlation between modularity and evoked activity within the DAN when the DMN was stimulated (Visit 1: R2 = 0.002%, p = 0.83; Visit 2: R2 = 0.005%, p = 0.73), and between modularity and evoked activity in the DMN when the DAN was stimulated (Visit 1: R2 = 0.02%, p = 0.45; Visit 2: R2 = 0.01%, p = 0.60). For further details on the correlation between brain modularity and other RSNs please see Supplementary Results and Figure S3.

Fig 3.

structural connectivity predictors of TMS-EEG propagation. (A) At visit 1 (left panel) a significant positive correlation was found between the structural connectivity of the stimulated network and the TMS-EEG response for both DMN (R 2 = 23%, p = 0.02) and DAN (R 2 = 32%, p = 0.006). The same pattern was observed at visit 2 (right panel) for both DMN (R 2 = 23%, p = 0.02) and DAN (R 2 = 25%, p = 0.01). (B) With DMN stimulation (top), a significant positive correlation was found between brain modularity and DMN response (red dots) for both visit 1 (R 2 = 18%, p = 0.04) and visit 2 (R 2 = 27%, p = 0.01), while no significant correlation was found for DAN (green dots) for both visit 1 (R 2 = 0.002%, p = 0.83) and visit 2 (R 2 = 0.005%, p = 0.73). With DAN stimulation (bottom), a significant positive correlation was found between the brain modularity and the DAN response (green dots) for both visit 1 (R 2 = 21%, p = 0.03) and visit 2 (R 2 = 23%, p < 0.02), while no significant correlation was found for DMN (red dots) for both visit 1 (R 2 = 0.02%, p = 0.45) and visit 2 (R 2 = 0.01%, p = 0.60). Note: AUC: Area under the curve; DAN: Dorsal attention network; DMN: Default mode network. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Conversely, no significant correlation was reported between the TMS-EEG response and Brain (Visit 1: DMN: R2 = 0.009%, p < 0.65; DAN: R2 = 0.08%, p < 0.18; Visit 2: DMN: R2 = 0.005%, p < 0.73; DAN: R2 = 0.11%, p < 0.14), Stim2Network (Visit 1: DMN: R2 = 0.001%, p < 0.85; DAN: R2 = 0.01%, p < 0.49; Visit 2: DMN: R2 = 0.08%, p < 0.21; DAN: R2 = 0.04%, p < 0.26) or Stim2Brain (Visit 1: DMN: R2 = 0.03%, p < 0.39; DAN: R2 = 0.008%, p < 0.69; Visit 2: DMN: R2 = 0.03%, p < 0.38; DAN: R2 = 0.0001%, p < 0.99) connectivity ( Fig. 4 A).

Fig 4.

Control analyses. (A) No significant correlation was found between the TMS-EEG response and Bra in (Visit 1: DMN: R2 = 0.009%, p < 0.65; DAN: R2 = 0.08%, p < 0.18; Visit 2: DMN: R2 = 0.005%, p < 0.73; DAN: R2 = 0.11%, p < 0.14), Stim2Network (Visit 1: DMN: R2 = 0.001%, p < 0.85; DAN: R2 = 0.01%, p < 0.49; Visit 2: DMN: R2 = 0.08%, p < 0.21; DAN: R2 = 0.04%, p < 0.26) or Stim2Brain (Visit 1: DMN: R2 = 0.03%, p < 0.39; DAN: R2 = 0.008%, p < 0.69; Visit 2: DMN: R2 = 0.03%, p < 0.38; DAN: R2 = 0.0001%, p < 0.99) connectivity. (B) FPN (top) and VIS (bottom) intrinsic connectivity do not correlate with visit 1 and visit 2 of TMS-EEG response both for DMN (left panel: visit 1 FPN: R2 = 0.006%, p = 0.72; VIS: R2 = 0.004%, p = 0.75; visit 2: FPN: R2 = 0.02%, p = 0.44; VIS: R 2 = 0.00001%, p = 0.98) and DAN (right panel: visit 1: FPN: R2 = 0.06%, p = 0.25; VIS: R2 = 0.04%, p = 0.34; visit 2: FPN: R2 = 0.03%, p = 0.37; VIS: R2 = 0.05%, p = 0.28). Note: AUC: Area under the curve; DAN: Dorsal attention network; DMN: Default mode network.

TMS-evoked network engagement in non-targeted RSNs

Connectivity metrics were also extracted for other RSNs as a control analysis. Considering TMS-EEG response following DMN stimulation ( Fig. 4 B, left panel; Figure S4), no significant correlations were found with the structural connectivity of DAN (Visit 1: R 2 = 0.001%, p = 0.85; Visit 2: R 2 = 0.07%, p = 0.21), VIS (Visit 1: R 2 = 0.004%, p = 0.75; Visit 2: R 2 = 0.00001%, p = 0.98), SM (Visit 1: R 2 = 0.08%, p = 0.19; Visit 2: R 2 = 0.11%, p = 0.09), AS (Visit 1: R 2 = 0.005%, p = 0.74; Visit 2: R 2 = 0.11%, p = 0.09), LIM (Visit 1: R 2 = 0.02%, p = 0.48; Visit 2: R 2 = 0.07%, p = 0.21) and FPN (Visit 1: R 2 = 0.006%, p = 0.72; Visit 2: R 2 = 0.02%, p = 0.44).

As for TMS-EEG response following DAN stimulation ( Fig. 4 B, right panel, Figure S4), no significant correlations were found with the structural connectivity of DMN (Visit 1: R 2 = 0.06%, p = 0.25; Visit 2: R 2 = 0.03%, p = 0.37), VIS (Visit 1: R 2 =0.04%, p = 0.34; Visit 2: R 2 = 0.05%, p = 0.28), SM (Visit 1: R 2 = 0.01%, p = 0.58; Visit 2: R 2 = 0.16%, p = 0.07), AS (Visit 1: R 2 = 0.04%, p = 0.32; Visit 2: R 2 = 0.08%, p = 0.19), LIM (Visit 1: R 2 = 0.02%, p = 0.45; Visit 2: R 2 = 0.02%, p = 0.51), FPN (Visit 1: R 2 = 0.006%, p = 0.72; Visit 2: R 2 = 0.03%, p = 0.37).

Relationship between structural connectivity and resting-state EEG

In order to demonstrate that the aforementioned results were not a function of resting-state brain oscillations, but due to specific activity elicited by TMS instead, source-level activity was computed (i.e. AUC) for both networks using baseline resting-state EEG data, and further correlated with the structural connectivity of the stimulated network and modularity. As shown in Fig. 5 A, no significant correlation was found between the baseline resting-state EEG data and the structural connectivity of the stimulated network for both DMN (left panel: visit 1: R2 = 0.002%, p = 0.83; visit 2: R2 = 0.01%, p = 0.54) and DAN (right panel: visit 1: R2 = 0.09%, p = 0.16; visit 2: R2 = 0.01%, p = 0.64).

Fig. 5.

Spontaneous vs Evoked EEG activity. (A) While TMS-evoked activity display correlations with network-level structural connectivity (red, green), resting-state EEG is not correlated with structural connectivity of the stimulated network, for both DMN (left panel: visit 1: R2 = 0.002%, p = 0.83; visit 2: R2 = 0.01%, p = 0.54) and DAN (right panel: visit 1: R2 = 0.09%, p = 0.16; visit 2: R2 = 0.01%, p = 0.64). (B) TMS-evoked activity was correlated with brain structural modularity (red, green) whereas no relationship was found considering resting-state EEG baseline for both DMN (left panel: visit 1: R2 = 0.005%, p = 0.74; visit 2: R2 = 0.0003%, p = 0.91) and DAN (right panel: visit 1: R2 = 0.002%, p = 0.82; visit 2: R2 = 0.03%, p = 0.39). Note: AUC: Area under the curve; DAN: Dorsal attention network; DMN: Default mode network. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

As shown in Fig. 5 B, no significant correlation was found between baseline resting-state EEG data and the structural modularity for both DMN (left panel: visit 1: R2 = 0.005%, p = 0.74; visit 2: R2 = 0.0003%, p = 0.91) and DAN (right panel: visit 1: R2 = 0.002%, p = 0.82; visit 2: R 2 = 0.03%, p = 0.39).

Discussion

A recent study by our group has used a network-perturbation approach to characterize individual brain dynamics within discrete brain networks with high temporal resolution ( Ozdemir et al., 2020 ), showing how source reconstruction of individual TEPs was highly reliable and propagated within the stimulated network. Here, we further expanded this concept showing how the structural modularity of the whole-brain and the connectivity of the stimulated network explain individual variability in propagation of TMS pulses. Results were replicable across two separated visits and selective for the stimulated network, supporting the relevance of macroscale structural connectivity in predicting network-level response to perturbation. Finally, brain structural wiring also seems more related to propagation of activity after perturbation rather than spontaneous resting-state activity, suggesting perturbation-based approaches as a valuable tool to investigate structure/function relationship in the human brain.

Contrary to our original hypothesis, the structural connectivity of the stimulated network revealed to be more relevant in predicting TMS-evoked activity than the connectivity of either the stimulated region or the whole brain. Instead, we found that signal propagation induced by single-region stimulation was highly related to the structural connectivity throughout the stimulated network, showing how the amount of fibers within the perturbed network might play a role in the propagation of action potentials. The idea that the underlying anatomical architecture of the cerebral cortex shapes brain functioning on multiple time scales has been already amply demonstrated by combining DWI and rs-fMRI techniques ( Honey et al., 2009, 2007 ). Recent studies have also evaluated the relationship between functional and structural connectivity using source estimated resting-state EEG recording and probabilistic tractography ( Chu et al., 2015 ). Moreover, a study reported a temporary decrease in the whole-brain correlation between source EEG activity and DWI connectivity measures after TMS in the α, β, and γ frequency bands ( Amico et al., 2017 ). In this context, our result represents the first set of empirical evidence that the efficacy of TMS, as a tool to induce network changes, is dependent upon the white matter connectivity of the same network, therefore offering crucial insight for the selection of optimal stimulation targets for maximal target engagement (e.g. stimulate networks with higher intrinsic connectivity over highly connected regions) as well as to predict response to TMS in clinical populations.

Our results also stress the importance of brain modularity for the prediction of individual propagation trajectories after TMS. Researches from various fields, such as physics or computational biology, have shown how modular network organization is associated with several leverages, including greater robustness ( Kirschner and Gerhart, 1998 ), the minimization of wiring costs ( Clune et al., 2013 ) as well as resilience to damage ( Wig, 2017 ). In the context of brain functioning, previous studies have shown how modular organization of brain networks accounts for better cognitive functions, such as working memory ( Stevens et al., 2012 ) and general intelligence ( Hilger et al., 2017 ). It has been suggested that information processing within segregated networks sub-serves specific cognitive functions, while the exchange of information between modules is assumed to be responsible for the coordination and integration of cognitive processes ( Gratton et al., 2012 ). The fact that we observe a relationship between network engagement induced by an external perturbation and the amount of brain modularity is in line with the aforementioned studies stating that a more segregated system might support more organized network-level information flow ( Baum et al., 2017 ). Interestingly, such relationship was not present for the non-stimulated network, suggesting how network engagement following an external perturbation may constitute a valuable biomarker of network dynamics.

Among the five hierarchical levels, non-significant correlations between TMS-evoked network activity and the connectivity of stimulated region with the rest of the brain (level a) is less expected than the negative result for connectivity of the entire brain (level e). One would assume that the specific region directly targeted by TMS always carries more information than either the whole network or the entire brain, instead our data suggest network structural connectivity as the best predictor. This finding highlights the need for investigation of whole-network dynamics, rather than local ones, when trying to unveil TMS propagation patterns. Related to this, a previous study showed that the size of TEPs response was related to local property of the stimulated tract (i.e., fractional anisotropy), but it did not control for the same measure at the network or at the whole brain level ( Kearney-Ramos et al., 2018 ). Contextually, recent evidences demonstrated that, in order to control a given brain network, a single node modulation is insufficient and a multi-node controllability model is needed instead ( Tu et al., 2018 ), which might partially explain why the stimulated region connectivity by itself is not able to predict network-level response.

Furthermore, our control analysis revealed that network structural connectivity and overall modularity were not related to resting-state EEG data, suggesting that the characteristics of network engagement following perturbation may significantly index the integrity of structural networks (outperforming spontaneous brain activity). So far, the majority of studies have focused on correlating individual measures with data acquired during unconstrained spontaneous activity showing how resting-state correlation patterns between various networks can be used to predict individual variability in several cognitive functions ( Finn et al., 2015 ), personality traits ( Adelstein et al., 2011 ) and behavior ( Fox et al., 2007 ). Nevertheless, recent evidences reported that the association between cognitive abilities and individual connectivity patterns could be better highlighted by task-based neuroimaging ( Finn et al., 2017; Greene et al., 2018 ). In this framework, a recent study published by our group ( Ozdemir et al., 2020 ) reported a significant positive correlations between TMS-evoked network engagement and high-order cognitive abilities (i.e. abstract reasoning, IQ) which was not observed when considering resting-state EEG. Present results demonstrate that brain structural wiring is more related to the propagation of activity after network perturbation than to spontaneous activity, implying perturbation-based approaches as a promising tool to investigate structure/function relationship in the human brain. Future studies should systematically compare intrinsic network structural connectivity with resting-state, task-based and perturbation-based activity data.

Finally, recent evidences have reported that patients with neurological and psychiatric diseases might have reduced structural network configuration compared to healthy subjects ( Alexander-Bloch et al., 2010 ). The modulation of rs-fMRI networks has become a relevant topic especially in clinical settings, but requires extreme precision because of individual differences in fMRI patterns ( Fox et al., 2012 ). In order to modulate pathological network interactions, the majority of TMS studies have used intrinsic fMRI connectivity to identify TMS targets ( Eldaief et al., 2011 ). We propose that structural connectivity might be a suitable way to better engage network dynamics, as the future clinical utility of TMS is dependent on considering structural architectures as a tool for patient-specific dosing.

In regard of the limitations of the study, a possible caveat is represented by the unfeasibility of applying the ad-hoc hierarchical framework used for structural connectivity to TMS-EEG data. In fact, only intrinsic network dynamics were extracted from source EEG data (i.e. average current density in a given network mask), and even if it would be possible to extrapolate the current density just from the stimulated region this local metric would be partially different from its “equivalent ” structural connectivity layers (i.e. layer a and b), representing connection from the stimulated region and not just activity or fibers related to the stimulated area. Moreover, source reconstructed TMS-EEG metrics might be contaminated by peripherally evoked artifacts, such as somatosensory and auditory potentials ( Conde et al., 2019 ), therefore challenging reliability. However, our control analysis reported no significant correlation between TMS-EEG measures and the intrinsic connectivity of the others RSNs, especially those that might be related with artifactual activity/propagation (e.g., somatosensory network related to muscular artifacts), suggesting our findings reflect transcranial evoked cortical activity alone.

In conclusion, intrinsic network structural connectivity provides valuable information to estimate network engagement following controlled, functional networks perturbation. Moreover, brain structural wiring is more related to propagation of activity after perturbation than spontaneous activity, implying perturbation-based approaches as a valuable tool to investigate structure/function relationship in the healthy and pathological brain.