Authors: Lia Papadopoulos, Christopher W. Lynn, Demian Battaglia, Danielle S. Bassett
Categories: Research Article, Computer and Information Sciences, Biology and Life Sciences, Medicine and Health Sciences, Social Sciences, Earth Sciences, Engineering and Technology
Source: PLOS Computational Biology
At the macroscale, the brain operates as a network of interconnected neuronal populations, which display coordinated rhythmic dynamics that support interareal communication. Understanding how stimulation of different brain areas impacts such activity is important for gaining basic insights into brain function and for further developing therapeutic neurmodulation. However, the complexity of brain structure and dynamics hinders predictions regarding the downstream effects of focal stimulation. More specifically, little is known about how the collective oscillatory regime of brain network activity—in concert with network structure—affects the outcomes of perturbations. Here, we combine human connectome data and biophysical modeling to begin filling these gaps. By tuning parameters that control collective system dynamics, we identify distinct states of simulated brain activity and investigate how the distributed effects of stimulation manifest at different dynamical working points. When baseline oscillations are weak, the stimulated area exhibits enhanced power and frequency, and due to network interactions, activity in this excited frequency band propagates to nearby regions. Notably, beyond these linear effects, we further find that focal stimulation causes more distributed modifications to interareal coherence in a band containing regions’ baseline oscillation frequencies. Importantly, depending on the dynamical state of the system, these broadband effects can be better predicted by functional rather than structural connectivity, emphasizing a complex interplay between anatomical organization, dynamics, and response to perturbation. In contrast, when the network operates in a regime of strong regional oscillations, stimulation causes only slight shifts in power and frequency, and structural connectivity becomes most predictive of stimulation-induced changes in network activity patterns. In sum, this work builds upon and extends previous computational studies investigating the impacts of stimulation, and underscores the fact that both the stimulation site, and, crucially, the regime of brain network dynamics, can influence the network-wide responses to local perturbations.
The brain is a multiscale system composed of many dynamical units that interact to produce a vast array of functions. At a large scale, macroscopic regions—each containing tens of thousands of neurons—are linked by a physical web of white matter tracts that facilitate the propagation of activity between distributed network elements. At the level of large neuronal ensembles or brain areas, collective activity is often rhythmic in nature [1], and these rhythms can become temporally coordinated between distant regions, giving rise to so-called functional interactions [2]. Importantly, oscillations have been implicated in a number of cognitive processes [3–9], and coherent activity is hypothesized to play an important role in interareal communication and information transfer among distributed brain areas [5, 6, 10]. Nonetheless, despite progress in mapping and characterizing the brain’s anatomical pathways and measuring neural oscillations, a number of questions remain as to how individual components in a brain network shape and modulate system-wide dynamics.
Among these questions, understanding how large-scale, oscillatory brain dynamics respond to localized perturbations is of critical importance [7, 11–14]. Because the brain is not a closed or static system, such activity changes could be induced by sensory inputs to primary sensory areas [15, 16], different tasks [17, 18], or other internal or regulatory processes [19–22]. In addition to naturally-induced changes, stimulation techniques such as transcranial magnetic stimulation [23], direct current stimulation [24], and alternating current stimulation [25] can also be employed to invoke modulations of dynamics in a specific brain area. By combining these techniques with imaging methods like EEG and MEG [26–31], it is possible to examine how the act of exciting a particular network component modifies rhythmic neural activity. Furthermore, in addition to its utility for basic science, neurostimulation has emerged as a promising approach for treating a number of neurological and psychiatric conditions [32–34].
Yet, while prior work has often focused on characterizing the proximal effects of local perturbations, a growing body of literature indicates that regional changes to neural activity can have more widespread consequences [11–14]. The realization that stimulation can have network-wide effects necessitates further investigations into the operating principles underlying such phenomena [35–42]. Furthermore, a crucial but seemingly understudied point is that the effects of perturbing a particular brain area can depend not only on the nature or location of the perturbation, but also on the intrinsic dynamical state of the system at baseline [43–45]. In particular, recent efforts have investigated the state-dependent effects of stimulation via precise experiments [46, 47]—focusing largely on alpha-band activity in single cortical areas—and via modeling [48–50]. These studies have uncovered robust relationships between the endogenous state of rhythmic activity and the capacity of external stimulation to modulate cortical oscillations in a given brain area. However, a pivotal next step is to extend the notion of state-dependence to the case of whole-brain networks, which acknowledge the fact that regions do not operate in isolation. Rather, in the case of large-scale brain networks, the macroscopic dynamical regime of the system arises from an interplay between units’ local activity and long-range anatomical coupling [51], leading to the emergence of collective oscillatory modes [40, 52]. Although it is reasonable to hypothesize that the global state of brain network activity should play a role in determining how a focal perturbation will manifest and influence distributed functional interactions, these ideas have yet to be systematically examined.
Thus, there is now a need to concurrently investigate and merge two outstanding (1) how regional stimulation spreads to induce distributed effects on brain network dynamics, and (2) how the global dynamical regime of the system impacts these effects. Here, we investigate these questions by constructing a biophysically-motivated model of large-scale, oscillatory brain activity, in which individual brain areas are modeled as Wilson-Cowan neural masses [53] coupled according to empirically-derived anatomical connectivity [51]. We first demonstrate that, in the absence of stimulation, the interareal coupling strength and the baseline excitation of the network transition the system between qualitatively distinct collective dynamical states. By providing additional excitation to a single brain area, we then systematically examine the consequences of such local stimulation on network activity. The primary contribution of this study is an exploration of how the effects of focal perturbations can depend not only on which area is stimulated, but also on the baseline dynamical regime of the non-linear model. Hence, this work builds upon previous whole-brain modeling efforts that have examined the effects of regional perturbations [35–37] with other work examining the state-dependent effects of stimulation in single cortical areas, but not large-scale networks [48, 49].
We find that in states of low baseline excitation, stimulation can significantly increase the frequency and power of regional activity, whereas in states of high background drive, local dynamics are less sensitive to perturbations. Importantly, these results show qualitative similarities and agreement with past work examining the focal effects of stimulation [48, 49]. We further find that, due to network interactions, regional perturbations can propagate and interact with brain areas’ ongoing rhythms. In particular, depending on the system working point, downstream areas that are strongly anatomically linked to the stimulated site also develop spectral components at the excited frequency of the stimulated region. Crucially, though, modifications to interareal phase-locking can additionally be induced in a broader frequency band comprising brain areas’ spontaneous, baseline oscillations, which may be well-separated from the excited frequency. Moreover, changing the dynamical regime of the system modulates the strength of associations between network-wide responses to perturbations in the baseline frequency band and structural or functional network connectivity. Hence, changing the collective oscillatory state of the system—which need not be entirely determined by the anatomical network—qualitatively changes the distributed effects of focal perturbations, and alters the relations between those effects and measures of either structural or dynamical organization. In sum, we use a simplfied, large-scale computational model to highlight that the effects of regional stimulation can depend both on the location of the perturbed site and on the global state of ongoing brain network dynamics. Though currently idealized, extending the reduced model to incorporate further biological realism and empirical constraints is an exciting direction for future work attempting to directly compare against experimental findings.
Human anatomical brain networks were reconstructed by applying deterministic tractography algorithms to diffusion-weighted MRI. In this study, we used a group-representative composite network assembled from 30 subject-level networks [54–56]. The mean age of participants was 26.2 years, the standard deviation was 5.7 years, and 14 of the subjects were female. To map anatomical networks, diffusion spectrum and T1-weighted anatomical images were acquired for each individual. For the DSI scans, 257 directions were sampled using a Q5 half-shell acquisition scheme with a maximum b-value of 5000 s/mm^2^ and an isotropic voxel size of 2.4 mm. We used an axial acquisition with repetition time TR = 5 seconds, echo time TE = 138 ms, 52 slices, and field of view of [231, 231, 125]mm. The T1 sequences used a voxel size of [0.9, 0.9, 1.0]mm, repetition time TR = 1.85 seconds, echo time TE = 4ms, and field of view of [240, 180, 160]mm. This data was initially collected for an earlier study [57], and was first published in [58]. The same data has also been used in several other prior investigations (e.g., [54, 56, 59]).
DSI Studio (www.dsi-studio.labsolver.org) was used to reconstruct DSI data using q-space diffeomorphic reconstruction (QSDR) [60], which reconstructs diffusion-weighted images in native space and computes the quantitative anisotropy (QA) of each voxel. Using the statistical parametric mapping nonlinear registration algorithm [61], the image is then warped to a template QA volume in Montreal Neurological Institute (MNI) space. Finally, spin-density functions were reconstructed with a mean diffusion distance of 1.25 mm with three fiber orientations per voxel. A modified FACT algorithm [62] was then used to perform deterministic fiber tracking with an angular cutoff of 55°, step size of 1.0 mm, minimum length of 10 mm, spin density function smoothing of 0.00, maximum length of 400 mm, and a QA threshold determined by DWI signal in the colony-stimulating factor [54–56, 58, 59, 63]. The algorithm terminated when 1,000,000 streamlines were reconstructed for each individual [54–56, 58, 59, 63] (Fig 1A).
Fig 1 Whole-brain imaging data, computational model of large-scale brain dynamics, and schematic of analysis.(A) An example of white matter streamlines reconstructed from diffusion imaging and tractography of a human brain. (B) Noninvasive magnetic resonance imaging scans of human brain anatomy are used to segment the cortex and subcortex into 82 regions. (C) Adjacency matrix for a group-averaged structural brain network. Individual brain areas are represented as network nodes, and normalized white matter streamline counts between region pairs are represented as weighted network edges. (D) Matrix of Euclidean distances between the centers of mass of all region pairs. (E) Left: Structural brain network representation; location of gray circles correspond to region centers of mass, and teal lines show the strongest 20% of interareal connections, with line thickness proportional to connection strength. The two encircled nodes correspond to an unperturbed region j and an excited region i in the large-scale brain network, with the perturbed region indicated by the yellow lightning bolt. Right: Schematic of the computational model of large-scale brain dynamics. The activity of a given brain region j is modeled as a Wilson-Cowan neural mass, composed of interacting populations of excitatory E and inhibitory I neurons. Neural masses are then coupled through their excitatory pools according to the structure of the anatomical brain network. A perturbation to region i (pictorially represented with the lightning bolt) is modeled as an increase in its excitatory input from P
E→ PE+ ΔPE. (F) The computational model generates oscillatory time-series of neural population activity for each brain region. These time-series can then be analyzed in Fourier space to determine relevant frequency bands for further analysis. After filtering time-series within the same frequency band of interest, functional interactions between brain region pairs are determined by extracting phase variables from each region’s filtered activity via the Hilbert transform, and then computing the phase-locking value to assess the consistency of phase relations over time and trials.
T1 anatomical scans were segmented using FreeSurfer [64] and parcellated using the Connectome Mapping Toolkit (http://www.connectomics.org) according to an N = 82 area atlas [65] of 68 cortical and 14 subcortical areas (Fig 1B; Table A in S1 Text). The parcellation was registered to the b0 volume of each subject’s DSI data, and region labels were mapped from native space to MNI coordinates using a b0-to-MNI voxel mapping [54–56, 58, 59, 63]. While we use a relatively coarse-grained atlas, it aligns with atlas sizes used in other computational modeling studies (e.g., [35, 66–69]),and was chosen to reduce the computational costs of numerical simulations. However, we do mention limitations involved with this choice in the Discussion.
All participants gave informed consent in writing and all protocols were authorized by the Institutional Review Board of the University of Pennsylvania.
To incorporate the structure of interareal connections into the model of large-scale brain activity, we represented the anatomical brain data as a network. This was achieved by first mapping each of the N = 82 regions to a unique node in a structural brain network A (see Table A in S1 Text for the mapping between node numbering and brain region labels). The structural edge weight Aij between two brain areas (nodes) i and j was then defined as the total number of streamlines between the two areas divided by the geometric mean of their volumes [54, 56]. Note that due to limitations in the non-invasive techniques available for constructing human connectomes [70], the resulting structural network is weighted but undirected. Because the brain is a spatially-embedded system [71], each region i also has a location ri = (xi, yi, zi) in real space. In the network representation, we defined the location of each node to be the center of mass of the corresponding region, allowing us to calculate matrix elements Dij representing the Euclidean distance between nodes i and j. Assuming a fixed conduction speed, these interareal distances are then used to approximate time delays for signal transmission in the computational model [35, 66, 67, 72, 73].
In this study, we report results using a group-representative structural brain network derived by combining individual brain networks across multiple subjects. We used a previously-established consensus method for constructing the group representative network that preserves both the average binary connection density of the individual brain networks, as well as the approximate edge-length distribution of intra- and inter-hemispheric connections [54]. More details on this pooling procedure can be found in [55]. A group-representative interareal distance matrix was constructed by averaging the pairwise Euclidean distance matrices across subjects. In what follows, we assume that A (or Aij) refers to the group-representative structural brain network, and that D (or Dij) refers to the group-averaged interareal distance matrix. We show the group-representative anatomical connectivity matrix in Fig 1C, and we show the group-averaged distance matrix in Fig 1D.
To model large-scale brain dynamics, we use a biophysically-motivated approach in which simulated activity is generated by a network of interacting neural masses [51]. In particular, the activity of each brain area is modeled as a Wilson-Cowan (WC) neural mass [53] and individual units are coupled according to the empirically-derived anatomical network. Importantly, these types of whole-brain computational models—which integrate non-linear, mean-field population dynamics with structural connectome architecture—have been utilized in a number of past efforts to gain insight into diverse neural phenomena [35, 66, 67, 69, 72–83].
Here, we employ such an approach to conduct a basic examination of how localized (regional) changes in neural activity affect dynamics across the brain. We offer a schematic of the model in Fig 1E. On the left, we show the structural brain network in real space. We focus on the two interconnected regions i and j encircled in black, of which the lower one (i) receives additional excitation (as denoted by the yellow lightning bolt). On the right, we show the setup of the coupled WC system for these two units. In the WC model, the activity of a particular brain region is defined by a coupled system of excitatory (E) and inhibitory (I) neuronal populations, and the dynamical variables are the mean firing rates of the E and I pools. The time-evolution of the average firing rates are in general governed by both intrinsic properties of the populations in a single region, as well as delayed, long-range input from other areas as dictated by the pattern of anatomical connectivity. In line with several previous studies [35, 66, 73–77, 81], we consider long-range connections to couple only the excitatory subpopulations of distinct brain areas.
The dynamics of the j^th^ brain area are governed by the following set of coupled differential
and
The variables Ej(t) and Ij(t) correspond to the firing rates of the excitatory and inhibitory populations of region j, and τE and τI are the excitatory and inhibitory time constants, respectively. The non-linear activation functions SE and SI of the excitatory and inhibitory pools are given by the sigmoidals
and
The quantities μE and μI give the mean firing thresholds for each subpopulation, and the gain parameters aE and aI set the spread of the firing thresholds for the two groups.
Dynamics of the excitatory ensemble are driven by (1) the local interaction strength within the excitatory population cEE, (2) the interaction strength from the inhibitory population to the excitatory population cIE, (3) constant, non-specific background drive PE,j, and also (4) interactions Wij corresponding to long-range excitatory inputs from different populations i that link to unit j via anatomical connectivity. Following [75–77, 84], we let Wij=Aij∑iAij, which is simply the connection weight from i to j, normalized by the total input to region j. Furthermore, C is a global coupling that tunes the overall interaction strength between different brain areas, and τij is a time delay between regions i and j that arises due to the spatial embedding of the brain network and the fact that signal transmission speeds are finite [35, 66, 67, 69, 72, 73]. We set τij=Dijv, where Dij is the Euclidean distance between regions i and j and v is a constant signal conduction speed. Activity in the inhibitory ensemble depends on (1) the interaction strength cEI from the excitatory population, (2) the local interaction strength within the inhibitory population cII, and (3) other possible non-specific inputs PI,j. Finally, to increase biological plausibility and incorporate the stochastic nature of neural dynamics, we add a term σE
ξ(t) to Eq 1 and a term σI
ξ(t) to Eq 2, which correspond to Gaussian white noise with zero mean and standard deviations σE and σI, respectively [35, 73]. In what follows, we will take the excitatory population activities Ej(t) of each brain area as the observables of interest [35, 66, 75, 76, 83].
Under appropriate parameter choices, the WC model can give rise to oscillatory dynamics [53]. Such rhythmic activity is ubiquitous in large-scale neural systems [1] and is the dynamical behavior of interest for this investigation. While oscillation frequencies observed in neural systems can span orders of magnitude [1], local neuronal populations often exhibit gamma band (30-90Hz) rhythms as a result of feedback between coupled excitatory and inhibitory neurons [4, 85, 86]. Furthermore, gamma oscillations and synchronization between distributed brain areas are associated with the flow of information between neuronal ensembles [10, 87, 88], are modulated by stimuli [15, 16], and are thought to underlie a number of cognitive processes [6]. Because gamma oscillations are robustly observed in excitatory-inhibitory circuits, we set parameters in the phenomenological WC model such that individual brain regions oscillate in the gamma band when coupled [73] (see Table 1). We also note that it may be interesting in future work to investigate other frequency bands or multiple frequency bands simultaneously [89].
As discussed further in Sec. SI of S1 Text, the non-specific background input PE is the typical control parameter used to tune the behavior of an isolated WC unit. At low values of PE, a single WC unit flows towards a low-activity steady-state (Fig. A, panel A in S1 Text), and at high values of PE, the system reaches a stable high-activity steady-state (Fig. A, panel C in S1 Text). At intermediate values of the excitatory drive, an isolated unit—with the parameters given in Table 1—will undergo a bifurcation and exhibit rhythmic activity in the gamma frequency band (Fig. A, panel B in S1 Text). Up to a point, increasing PE within this intermediate region leads to oscillations with increasing amplitude and frequency (Fig. A, panels D–F in S1 Text).
The situation becomes more complex when multiple WC units are coupled via the structural connectome. In this scenario, an individual region’s dynamics are determined by a combination of the constant drive PE and the strength of delayed inputs from other parts of the network, which are modulated by the coupling C and the structural connectivity A. To account for these two influences, we consider both PE and C as tuning parameters, and examine working points at which the combination of PE and C generate oscillatory activity in individual brain areas. Finally, we set the signal propagation speed to a fixed value of v = 10m/s, which is in the range of empirical observations and previous large-scale modeling efforts [35].
The baseline condition of the network corresponds to the situation in which all brain areas receive the same level of background drive, such that PE,j=PEbase for all j ∈ {1, …, N}. To investigate how regional perturbations affect brain-wide dynamics, we examine the effects of increased excitation to a single brain area. This is modeled as a selective increase in drive to the excitatory population of the perturbed neural mass i such that PE,i → PE,i + ΔPE, where ΔPE > 0 denotes the strength of the perturbation [35] (see Fig 1E for a schematic). The dynamics of the system in the baseline state can then be compared to the situation in which unit i receives additional input (i.e., where we have PE,i=PEbase+ΔPE and PE,j=PEbase for all j ≠ i).
We note that, phenomenologically, excitation of a given brain area could occur through a number of mechanisms, including sensory input to primary sensory regions, brain stimulation, or, alternatively, via internal processes that regulate inputs to or excitability levels of specific neuronal populations. The goal of this work is to study the effects of localized excitations generally, rather than to design a detailed model of a specific type of perturbation. For this reason, we choose to study the simplest case of constant excitation.
The equations governing the time evolution of the excitatory and inhibitory population activities form a system of coupled stochastic, delayed differential equations. We numerically integrate this system using the Euler-Mayurama method with a time step of dt = 5 × 10^−5^s. For the time delays, we round each τij to the nearest multiple of the integration time step dt, and for the initial conditions, we assume a constant history for each unit’s activity of length equal to the longest delay in the system. After running a simulation, we discard the first tburn = 1 second so that our analysis is not biased by transients or the specific choice of initial conditions. Each time-series is then downsampled to a resolution of dtds = 1 × 10^−3^s. The parameters for the numerical simulations are shown in Table 1.
Useful characteristics of the simulated activity are apparent in the frequency domain (see Fig 1F). Here, we use Welch’s method (as implemented in MATLAB R2019a) to estimate the power spectral density (psd) of the excitatory population activities. We use window sizes of 1 second with 50% overlap, and subtract the mean of each time-series before computing the psd.
To quantify the extent of temporal coordination between different brain areas, we use the phase-locking value (PLV) [90]. This measure is commonly used to assess the level of coherence between phases in a given frequency band. Importantly, because the state variables in the WC model are real-valued signals with possibly multiple spectral components, we compute PLVs for a given frequency band by (1) filtering all raw excitatory time-series within the same specified frequency range, and (2) extracting instantaneous phases for the given frequency band using the Hilbert transform (see Fig 1F). In the following two sections, we describe these steps in more detail.
Given a real-valued signal X(t), it is possible to define instantaneous phase and amplitude variables that describe the signal using the Hilbert transform. Importantly, although the Hilbert transform can theoretically be computed for an arbitrary signal X(t), the instantaneous amplitude A(t) and phase θ(t) are only physically meaningful for relatively narrowband signals [91]. It is therefore necessary to filter a signal before taking the Hilbert transform. Here, raw time-series were bandpass filtered in a frequency range fo ± Δf Hz using a 6th-order Butterworth filter in the forward and backward directions. In the results section, we describe how fo and Δf are determined during the presentation of various findings that depend on computing the Hilbert phase. Filtering was carried out in MATLAB using the ‘butter’ and ‘filtfilt’ functions. After filtering the simulated activity, the Hilbert transform was applied to extract instantaneous phases for the given frequency band. The Hilbert transform was implemented using the ‘hilbert’ function in MATLAB. More details on the Hilbert transform can be found in Sec. SXIII of S1 Text.
The outputs of the filtering and Hilbert transform processes described in the previous section are instantaneous phases θi(fo, t) derived from the excitatory activity Ei(t) of each brain region i at a given central frequency fo and time t (Fig 1F). From these phases, we can quantify the extent of phase-coherence between brain areas’ signals in a given frequency band using the phase-locking value (PLV); see Fig 1F. The PLV—here denoted symbolically as ρij—between two phase time-series θi(t) and θj(t) is given by
where Ts is the number of sample time points over which the phase-locking is computed. If the phase difference Δθij(t) = θi(t) − θj(t) is constant over a given time window, ρij will be equal to 1, whereas if the phase-differences are distributed uniformly, ρij will be approximately 0; in this way, ρij ∈ [0, 1].
We would also like to ensure that the PLV reflects the consistency of phase relations that arise from interactions (direct or indirect), and not locking arising from the fact that two regions happen to have the same frequency, but, possibly, a different phase relation in every trial. We therefore concatenate phase time-series from different trials before computing the PLV [92], where each trial is a simulation run with different random initial conditions and noise realizations. Accordingly, a high PLV indicates that across time and trials, the activity of the corresponding regions exhibits a consistent phase relationship within a particular frequency band.
As with structural connectivity, it is useful to think of a given N × N matrix of PLV values as a network where the element (edge) ρij is the phase-coherence between region (node) i and region (node) j. In contrast to the structural network, this PLV-based network represents the presence of functional associations between brain regions’ activity. Following common terminology, we will thus often refer to phase-locking as “functional connectivity” and phase-locking matrices as “functional networks”.
All data and statistical analysis was performed in MATLAB release R2019a. Statistical dependencies between two variables were assessed via the Spearman rank correlation, using the built-in MATLAB function ‘corr’. Throughout the text, we denote the Spearman correlation coefficient as rs. Rank correlations are considered statistically different from zero if the corresponding p-value is less than 0.05.
Throughout the text we compute a number of different measures to characterize the behavior of the system at baseline and under focal perturbation. To aid the readability of the manuscript, we list these quantities in Table 2 with a brief summary. The measures are listed according to the section in which they first appear.
Depending on the values of various parameters, the brain network model exhibits different qualitative behaviors. Importantly, different baseline states may in turn result in distinct modulations of brain-wide activity patterns in response to local perturbations. We thus begin by characterizing the behavior of the system at baseline (i.e. in the absence of regional stimulation). This initial study will provide context for our subsequent investigations examining how the effects of focal excitations depend upon the system’s baseline state.
We focus on two parameters of (1) the level of generic background input to the excitatory populations PEbase and (2) the global coupling strength C. Recall that for an isolated WC unit, PEbase is a bifurcation parameter that transitions population activity between a quiescent and an oscillatory state [53, 93]. However, when examining a network of coupled neural masses, the dynamics of each element are also dictated by inputs from other units in the system. The parameter C is a second control parameter that globally scales the interaction strength between brain areas by tuning how much input a given region receives from its neighbors in the network. The nature of both local and network-wide dynamical behaviors will thus change depending on the combination of PEbase and C, allowing the system to exist in markedly different states. Though the tuning parameters in the network model are phenomenological, from a biological standpoint, global changes in these parameters could represent, for example, the effects of neuromodulation [94–96]—which exerts widespread influences across the brain [97]—or changes in brain state more generically.
To quantify how model behavior varies as a function of PEbase and C, we perform a sweep over a broad range of these parameters, considering values of PEbase∈[0.5,0.85] in steps of ΔPEbase=0.05, and values of C ∈ [0, 5] in steps of ΔC = 0.1. These ranges were chosen to allow for the exploration of multiple oscillatory regimes of the system. For each parameter combination, we run five, 2-second-long simulations. The values of all other parameters are defined in Table 1, with the exception that, for these sweeps, we run noiseless simulations in order to more precisely demarcate the boundaries between different dynamical modes of the model.
We begin by computing two measures that quantify regional (1) the time-averaged firing rate E(t)¯, and (2) the frequency at maximum power (peak frequency) f^peak^ of a given region. To obtain summary measures characterizing the state of the system as a whole, we compute network-averages of these quantities, denoted by angled-brackets. In studying 〈E(t)¯〉 as a function of PEbase and C, we observe three principal regimes (Fig 2A). When both PEbase and C are low, the system settles to a state of low average firing rate (white region); this state corresponds to a non-oscillatory, low-activity equilibrium. In contrast, when PEbase and C are both high, the average firing rate saturates at a high level (dark green region); this state corresponds to a non-oscillatory, high-activity equilibrium. Finally, at intermediate values of these parameters, the mean firing rate varies between the low and high extremes, and the regional activity is oscillatory; because we wish to consider the rhythmic nature of brain activity, this is the relevant portion of parameter space.
Fig 2 Long-range coupling strength C and background drive PEbase modulate firing rates and oscillation frequencies at baseline.(A) The time- and network-averaged population firing rate 〈E(t)¯〉 as a function of C and PEbase (units are arbitrary). (B) The network-averaged peak frequency of regional activity 〈f^peak^〉 as a function of C and PEbase. (C) A segment of the activity of one brain area and (D) the corresponding power spectra of the same area at the working point denoted by the red dot in panels A and B (PEbase=0.553, C = 2.5). (E) A segment of the activity of one brain area and (F) the corresponding power spectra of the same area at the working point denoted by the orange dot in panels A and B (PEbase=0.57, C = 2.5). (G) A segment of the activity of one brain area and (H) the corresponding power spectra of the same area at the working point denoted by the yellow dot in panels A and B (PEbase=0.7, C = 2.5).
Next we seek to understand how 〈f^peak^〉 varies in the PEbase—C plane (Fig 2B). A clear wedge-shaped area marks parameter combinations that give rise to network-averaged peak frequencies in the gamma range. As with the firing rate, the peak frequency tends to increase (decrease) with either increasing (decreasing) background excitation or coupling strength. By comparing Fig 2B to 2A, we see that the white areas surrounding the purple wedge correspond to the regions of parameter space where the firing saturates at a fixed low or high value. In Sec. SII of S1 Text, we describe a systematic method for determining boundaries in the 2D space spanned by C and PEbase that indicate the onset or disappearance of oscillatory activity (see Fig. B in S1 Text). In what follows, we use PE*(C) to denote the level of background drive at which oscillations begin to emerge for a fixed coupling strength C. We refer the reader to Sec. SII of S1 Text for a detailed description of how this value is determined from the simulations. Furthermore, we often plot quantities as functions of the relative drive PEbase-PE*(C), such that PEbase-PE*(C)=0 indicates the transition point from a low-activity state to an oscillatory state at a coupling C.
To provide further intuition for how dynamics vary within this parameter space, we study example time-series and power spectra for three different baseline states (colored dots in Fig 2A and 2B). Note that these working points correspond to an intermediate coupling value of C = 2.5, but varying levels of the constant baseline input PEbase. We begin with the working point PEbase=0.553, which sits just beyond the boundary indicating the transition to sustained rhythmic activity. From the time-series, we observe that the activity is oscillating (Fig 2C), and the spectra indicates a peak frequency of ≈40Hz on a broadband background (Fig 2D). We next consider the working point PEbase=0.57. In this state, each unit receives slightly more drive, leading to higher-amplitude oscillations (Fig 2E and 2F). However, although peak spectral power increases, amplitude modulations can still be seen in the corresponding time-series (Fig 2E). Finally, we consider the working point PEbase=0.7. Here, the activity is characterized by regular, high-amplitude oscillations (Fig 2G). Furthermore, inspection of the power spectra indicates a single, narrow peak at a slightly higher frequency than the previous working point (Fig 2H).
Both the firing rate and the power spectra are measures that quantify the nature of individual regions’ activity. For networks, it is also imperative to define measures that capture information about the extent of dynamical order in the system as a whole. Indeed, for networks of coupled units, the system’s “state” is defined not only by the behavior of individual units, but also by how their dynamics are interrelated. Here, we are interested in the degree to which regional dynamics are coherent, which we quantify via the PLV between regions’ activities. To compute PLVs for baseline conditions, we begin by filtering the activity of each unit in the same, common frequency band. This band is determined by first finding the peak frequency of each unit at the given working point. Hence, we obtain a set of N values {fipeak} corresponding to the peak frequencies of all units i∈ {1, …, N} at baseline. We then filter the activity of every region in a frequency band spanning 10Hz above the maximum peak frequency in the network (max{fipeak}) and 10Hz below the minimum peak frequency in the network (min{fipeak}). After identically filtering each unit’s activity in this common band, we extract Hilbert phases from the filtered signals. Finally, PLVs between all pairs of brain areas are computed according to Eq 5, using 50 different simulations (trials) of 5 seconds each (with noise included).
To summarize how the overall level of coherence in the network varies as a function of the background drive and coupling strength, we defined a macroscopic order parameter as the average of the PLVs over all pairs of units in the ρ^global^ = 〈ρij〉. This quantity ranges between 0 and 1, where larger values indicate a more dynamically ordered state of the network. In general, we find that the background input and the coupling strength interdependently tune the level of coherence in the system (Fig 3A). At low coupling, brain areas cannot coordinate their dynamics and ρ^global^ remains at a relatively low value for a range of drives. In contrast, as the coupling is increased, we begin to see a qualitative change in behavior. For higher values of C, we observe that ρ^global^ varies non-monotonically (first increasing and then decreasing) as a function of the (relative) background drive. For a given coupling C, there appears to be a “critical” value at relatively small but non-zero PEbase-PE*(C) where the system develops a well-defined peak in global coherence. As the drive is increased further, ρ^global^ begins to decrease and then eventually plateaus, albeit with some fluctuations. More specifically, at levels of background drive well beyond the state of peak coherence, ρ^global^ relaxes to an intermediate value between its peak and its value at the lowest background input. In this regime, the system resides in a state of partial order. Increasing the coupling has the effect of amplifying the maximum value of ρ^global^ (although ρ^global^ remains well below 1 for all couplings considered), but does not appear to significantly affect the order parameter to the right side of the peak.
Fig 3 Long-range coupling strength C and background drive PEbase modulate network phase-coherence and relationships between structural and functional connectivity at baseline.(A) The global order parameter ρ^global^ vs. PEbase-PE*(C), for different fixed values of C. Error bars are estimated from 100 bootstrap samples of the simulations at each coupling and background drive, and correspond to ± one standard deviation of the bootstrap disribution of ρ^global^. (B) The difference between the global and local order parameters, ρ^global^ − ρ^local^, vs. PEbase-PE*(C), for different fixed values of C. Error bars are estimated from 100 bootstrap samples of the simulations at each coupling and background drive, and correspond to ± one standard deviation of the bootstrap distribution of ρ^global^ − ρ^local^. (C) Region-by-region PLV matrices for various values of PEbase-PE*(C) at fixed C = 2.5. The boxed matrices correspond to the red, orange, and yellow working points in Fig 2 and in panels A and B of this figure. (D) The Spearman correlation r
sbetween structural node strength s^struc^ and functional node strength s^func^ vs. PEbase-PE* at fixed C = 2.5. Empty circles indicate that the correlation was not statistically significant at the p = 0.05 level. The arrows mark three different working points—WP1, WP2, and WP3 (which correspond to the red, orange, and yellow dots/boxes in this figure)—that will be studied in detail.
To provide further intuition for this behavior, we focus on an intermediate coupling of C = 2.5 and examine the pairwise coherence patterns ρij for several values of the background drive (Fig 3C). At the lowest (relative) baseline input (PEbase-PE*(C)=0), some organization can be seen in the PLV matrix, but the system is weakly coherent overall. In this state, units exhibit relatively low amplitude oscillations, and are therefore more influenced by noise. It is thus reasonable to expect little phase-locking at low background drive. However, with only a small increase in the non-specific input (e.g., PEbase-PE*(C)=0.003), we observe distributed increases in coherence and a large spread of high, medium, and low coherence pairs dispersed throughout the network. Increasing the background drive slightly more (e.g., PEbase-PE*(C)=0.02) leads to the emergence of large, highly-coherent blocks that span the system. This working point sits near the peak of ρ^global^ and represents a highly ordered state of the system. As the background drive is increased further, though, phase-locking begins to decrease widely throughout the network and the coherence pattern markedly changes into a more segregated architecture. In particular, for high PEbase-PE*(C), we observe the emergence of smaller phase-locked clusters (Fig 3C, Row 3). To understand this shift in behavior, it is important to note that increasing PEbase increases the extent to which regional activity is independently generated in each area vs. driven by long-range network interactions. The strengthening of regional oscillations and enhanced influence of local dynamics with increasing PEbase seems to eventually hinder the ability of units to adjust their rhythms and achieve widespread coherence. Note that phase-locking is also made especially difficult by the large variance in the distribution of interareal delays imposed by the connectome’s spatial embedding, and indeed, for high background drive conditions, more strongly connected and spatially nearby units are those able to maintain stronger coherence.
In general, our observations point to complex behavior in which the macroscopic order parameter varies non-monotonically as a function of the baseline input and network coupling strength (Fig 3A and 3C). Therefore, a variety of qualitatively different regimes exist, beyond just a simple binary separation into a disordered and ordered state. To more quantitatively distinguish network states before and after the point of peak coherence, we also considered a local order parameter ρlocal=∑ijAijρij/∑ijAij, which is a weighted average of ρij with weights equal to the strength of structural network connections. In this way, ρ^local^ will be larger when more strongly connected brain areas are more phase-locked. In Fig 3B, we show ρ^local^ − ρ^global^
vs. PEbase-PE*(C) for different values of the coupling C. Beyond a certain point, the curves for all couplings exhibit a clear upward trend where the extent of local coherence increases relative to the extent of global coherence. This behavior indicates that the macroscopic state of the system becomes increasingly constrained by structure as the background drive increases. Hence, even though the level of global coherence can be similar to the left and right of peak ρ^global^, the system is in qualitatively different dynamical modes in the two regimes. Also note that for the higher couplings, ρ^local^ − ρ^global^ first decreases before consistently rising. This variation occurs because, for large enough coupling strengths, the level of global coherence is able to compete with the level of local coherence at background drives near peak ρ^global^.
As a final demonstration of the complexity of the structure-function landscape across operating points, we consider the relationship between brain areas’ structural and functional connectivity strengths as a function of PEbase-PE*(C) for a fixed coupling C = 2.5 (Fig 3D). The structural strength of node j, sjstruc=∑i=1NAij, is a common measure of a brain area’s importance in an anatomical network [98]. Similarly, the (baseline) functional strength of node j, sjfunc=∑i=1Nρij, quantifies how dynamically integrated that region is to the network as a whole. From Fig 3D, we observe that shifting the system’s working point can drastically alter how—and the extent to which—structural strength and functional strength are related. Specifically, while there tends to be a weak positive correlation between s^struc^ and s^func^ at high background drives (e.g. at WP3), the correlation disappears (e.g. at WP2) and then reverses in sign (e.g. at WP1) as the background drive is lowered. Critically, these transitions occur in the absence of any change to the anatomical connectome, and are instead driven by a global change in the behavior of brain areas’ dynamics (induced by changing the level of background input). Also note that when the correlations are significant, they are intermediately-valued. Together, these results indicate that while a given structural network may only be able to support specific patterns of coordinated activity, the relationships between the two are not trivial and are modulated by dynamic properties [99, 100]. In general, functional connectivity thus reflects a complex interplay between both anatomical connectivity and the system’s dynamical state.
It is crucial to remark that the behaviors seen here are more diverse than what tends to occur in simpler phase-oscillator models, where coupling strength is the main control parameter and typically induces a monotonic increase in synchrony. A critical difference between phase-based models and the more realistic WC model considered here is that, for the latter case, unit dynamics are described and coupled by real-valued signals that represent regional activity. Hence, widespread changes in the amplitude or stability of areas’ dynamics (in addition to changes in coupling strength) can affect the macroscopic state of the network. Indeed, the preceding analyses show that global modulations in the level of diffuse, constant input to the neural populations can push the system into very different oscillatory modes, beyond just a steady progression from an incoherent to a coherent state. In what follows, we will exploit this behavior to examine how the effects of focal perturbations depend not only on which region is targeted, but also on the baseline working point of brain network dynamics as a whole.
To investigate how local perturbations modulate brain network dynamics, and specifically how the effects may depend on the system’s collective state, we begin with an in-depth examination of three distinct working points. In particular, we focus on a fixed intermediate coupling strength C = 2.5 for which the system exhibits a clear peak in ρ^global^ (Fig 3A). We then examine two values of the background drive PEbase that place the system either in a state preceding (WP1) or following (WP3) the global coherence peak. In Sec. SIII of S1 Text, we also present results for a state in which the system is approximately at peak global coherence (WP2). We then proceed to more generally characterize the global impacts of stimulation as the background drive is varied across a wide range. Throughout the text, stimulation of a single brain area i is introduced by increasing its excitatory input by an amount ΔPE,i = 0.1, while keeping all other regions at their working-point-specific baseline drive. Finally, in S1 Text, we verify that results hold for different values of PEbase in the vicinity of those studied in the main text (Sec. SVIII), we examine the effects of varying the perturbation strength (Sec. SX), and we consider an alternative value of the global coupling (Sec. SXI). Note that our goal is not to exhaustively analyze all possible parameter combinations, but rather to demonstrate that the network response to stimulation qualititatively varies for different dynamical regimes.
We begin with the working point WP1 located at C = 2.5 and PEbase-PE*(C)=0.003, below peak coherence (Fig 3C, Row 1, Column 2). Here, the system is perched just past the boundary marking the transition between the quiescent state and the commencement of rhythmic dynamics. Hence, regional activity is oscillatory but of relatively low amplitude (see Fig 2C), and the power spectra is broad (see Fig 2D).
We first examine the effects of a regional perturbation on areas’ time series and power spectra (Fig 4A–4C). In agreement with past experimental and modeling studies [16, 89, 101, 102], increased drive to the excitatory pool of region i increases the amplitude and frequency of its oscillations (Fig 4C, Left). In particular, stimulation causes an increase in the power, narrowing of the spectra (associated with an increase in periodicity of regional activity), and a shift of the peak frequency from ≈ 40Hz at baseline to ≈ 50Hz when excited. We also note the appearance of modulation sidebands in the excited spectra to the left and right of the peak frequency, which arise due to the modulation of the excited region’s time-series by the lower-frequency input it receives from other areas in the network [92]. This modulation also results in a spectral peak at ≈ 16Hz—which is the difference between the new, excited frequency and the sideband peaks, and is a marker of quasiperiodic amplitude modulation in the time-series. To more carefully quantify the effects of an excitation to region i, we consider the shift in the peak frequency of unit i, Δfi,δipeak=fi,δipeak-fipeak, between its excited and baseline states (Fig 4D). Calculating these differences for all choices of the stimulated brain area, we find that they range from about 6Hz to 16Hz, with an average value of 〈Δfi,δipeak〉≈10.5 Hz. These perturbation-induced shifts thus yield excited peak frequencies that are well-separated from the range of peak frequencies in the baseline state (Fig 4E).
Fig 4 Regional excitation causes local and downstream changes to brain areas’ power spectra in different frequency bands at WP1.(A) Schematic of a brain network depicting the stimulated site i in brightest red. The black arrows point to two other regions j and k that lie at progressively further topological distances from the perturbed area in the structural network. In this figure, regions i, j, and k correspond to brain areas 1 (R–Lateral Orbitofrontal), 4 (R–Medial Orbitofrontal), and 10 (R–Precentral), respectively. (B) Left: A segment of region i’s activity time-course in the baseline condition. Right: A segment of region i’s activity time-course when it is stimulated. (C) Power spectra of area i and two other downstream regions j and k. In all three panels, the lighter curves correspond to the baseline condition, and the darker curves correspond to the state in which i is driven with additional input. The gray vertical lines indicate the peak frequency fi,δipeak of region i in the excited condition. (D) Histogram of the shift in peak frequency Δfi,δipeak induced by stimulating unit i, plotted over all choices of the perturbed area. (E) Distribution of peak frequencies of all units in the baseline condition {fipeak} (light gray) and distribution of the peak frequency units acquire when directly excited {fi,δipeak} (dark gray). (F) Average power spectra 〈psd〉
j≠iover all units j ≠ i at baseline (light gray) and when unit i is perturbed with additional input (dark gray). (G) Average difference 〈Δpsdj,δi〉j≠i of the spectra of unit j ≠ i when unit i is excited and in the baseline condition, where the average is over all units j ≠ i. For reference, the light gray vertical lines denote the minimum and maximum peak frequency across units in the baseline state, and the dark gray line indicates the peak frequency acquired by the stimulated region i. Shaded boxes denote two frequency bands of (1) the baseline band (purple) consisting of the main oscillation frequencies of brain areas under baseline conditions, and (2) the excited band (green) centered around the peak frequency that the stimulated region inherits. In subsequent analyses, we assess perturbation-induced changes in the PLV between brain areas in the baseline band, Δρδibase (purple), and in the excited band Δρδiexc (green).
We next consider the power spectra of two other units j and k located at increasing topological distances from the excited region, where a shorter topological distance indicates that two areas are linked by a path of stronger structural connections [98]. (Fig 4C, Middle, Right). We observe that unit j maintains its initial frequency content, but also develops new peaks centered at the frequency of the excited region and at the difference of the excited frequency and the baseline peak. In contrast, the spectra of unit k—which is more weakly structurally connected to the stimulated site—is relatively unchanged. Hence, depending on the network structure, stimulation of region i can also cause alterations to other regions’ spectra. In general, the power modulation of a downstream area’s spectra at the peak frequency of the stimulated site decays with increasing topological distance between the dowstream area and the perturbed region (see Sec. SV in S1 Text). To summarize how the spectra of other brain areas are altered by driving region i with additional input, we compare the average power spectral density 〈psd〉j≠i over all units j ≠ i at baseline and when unit i is stimulated (Fig 4F). At baseline, the network-averaged spectra is relatively broad and contains multiple peaks—a main one at 38Hz and a smaller peak around 34Hz. In addition, a local excitation produces complex and broadband alterations in power, as expected in a scenario of quasiperiodic entrainment between nonlinear oscillators [103]. For this example, we observe the appearance of an entirely new peak at 50Hz, but also an enhancement of the lowest baseline peak and a depression of the highest baseline peak. These changes are perhaps more apparent in Fig 4G, which shows the average difference 〈Δpsdj,δi〉j≠i in the spectra of unit j ≠ i between when unit i is excited and the baseline condition, where the average is over all units j ≠ i. In sum, we see that a regional enhancement of neural activity causes non-local modulations in power both at the frequency of the directly stimulated brain area, as well as at the system’s baseline oscillation frequencies. These analyses suggest that there are two relevant frequency bands to consider for subsequent (1) a relatively broad band containing the main frequencies of brain areas in the baseline state, and (2) a band centered around the peak frequency of the excited unit. In what follows, we will denote these two bands as “baseline” and “excited”, and consider changes in phase-locking, Δρδibase and Δρδiexc, in each band induced by local perturbations.
We are now prepared to study how focal perturbations alter the coordination of network-wide dynamics. Specifically, we examine changes in interareal phase-locking. We separate our analysis into two frequency bands—baseline and excited—by filtering regional activity in each band, extracting Hilbert phases from the filtered signals, and then calculating the PLV for each pair of regions within each band (see Fig 1E; Materials and methods). Since spectra are relatively broad at baseline, a single baseline frequency band for the network is determined by first finding the set of peak frequencies for each unit i in the baseline state, {fipeak}. Next, the lower frequency for the common baseline band is set to min{fipeak}-10 Hz, and the upper frequency is set to max{fipeak}+10 Hz. A region-by-region PLV matrix corresponding to the single baseline band is then computed after identically filtering each unit’s activity in this frequency range. To examine phase-locking between units within the much narrower excited band corresponding to a given stimulated region i, we first extracted the peak frequency of region i when it is stimulated, fi,δipeak. A PLV matrix corresponding to unit i’s excited frequency is then computed after filtering each region’s activity in the same frequency band ranging from fi,δipeak-1.5 Hz to fi,δipeak+1.5 Hz. This range was chosen to contain the majority of the excited band peak, while including as little of the original baseline band as possible. If the peak frequency of the stimulated area was not more than 3.5Hz above the largest baseline peak frequency, then we only examined PLV changes in the baseline frequency band. Our choices are motivated by the following the notion of an excited frequency band is only meaningful when a perturbation introduces a new spectral peak into the system that is separated from the frequencies present in the baseline condition. Also note that, unlike in the baseline band, areas exhibit little power at the excited frequency prior to stimulation; hence, we use changes in excited-band PLV as a measure of how effectively induced activity at the excited frequency spreads in the network.
To provide intuition about how phase-locking is altered upon a local perturbation, we consider the effect of stimulating two different brain areas (left and right panels in Fig 5A and 5C). These examples show that regional stimulation induces phase-locking at the excited frequency (Fig 5C), but can also cause changes in coherence in the frequency band containing the original oscillatory activity of the system (Fig 5A). Note that the excited band effects are mostly positive (due to the fact that power in the excited band is boosted by stimulation), whereas the baseline band effects can be both positive and negative. Furthermore, the patterns induced in the excited band (Fig 5C) are distinct from the modulations that occur in the baseline band (Fig 5A), and the phase-locking changes are markedly different between perturbation of region i (left panels) and perturbation of region j ≠ i (right panels). Thus, depending on the frequency band considered and the excited area’s location within the large-scale brain network, stimulation induces different responses across the system as a whole.
Fig 5 Phase-locking changes at WP1 are driven by local excitations of neural activity, differ between excited and baseline frequency bands, and are differentially related to structural and functional network properties.(A) Pairwise changes in the PLV inside the baseline band Δρ^base^ when region i (Left) or region j ≠ i (Right) is perturbed. In this figure, regions i and j correspond to regions 4 (R–Medial Orbitofrontal) and 23 (R–Lateral Occipital), respectively. (B) Network-averaged absolute PLV changes in the baseline band 〈|Δρδibase|〉 caused by stimulation of different brain areas. (C) Pairwise changes in the PLV inside the excited band Δρ^exc^ when region i (Left) or region j ≠ i (Right) is perturbed. (D) Network-averaged absolute PLV changes in the excited band 〈|Δρδiexc|〉 induced by stimulation of different brain areas. (E) The quantity 〈|Δρδibase|〉 vs. structural node strength sistruc (Left), and the quantity 〈|Δρδiexc|〉 vs. structural node strength sistruc (Right). (F) The quantity 〈|Δρδibase|〉 vs. functional node strength sifunc (Left), and the quantity 〈|Δρδiexc|〉 vs. functional node strength sifunc (Right). In panels (E) and (F), insets indicate Spearman correlation coefficients between the plotted quantities and their associated p-values).
To summarize the global effect of regional excitation, we calculate the average absolute change in PLV induced by driving each brain area with additional input. We use the notation 〈|Δρδibase|〉 and 〈|Δρδiexc|〉 to denote the network-average of the absolute PLV changes in the baseline and excited frequency bands, respectively, induced by stimulating region i. Note that since the PLV is always between 0 and 1, the maximum possible value of both quantities is 1. Furthermore, we use a phase-randomized null model (described in Sec. SXII of S1 Text) to assess whether pairwise PLV changes are significant. Prior to calculating the network-wide averages 〈|Δρδibase|〉 and 〈|Δρδiexc|〉, non-signicant changes are set to zero. We observe that the global responses exhibit a large degree of variability across different choices of the stimulated area (Fig 5B and 5D). That is, perturbation of some areas induces larger system-wide modulations of phase-locking than others. Furthermore, regions that induce the largest overall changes in PLV inside the excited frequency band are not necessarily those that cause the largest alterations of PLV in the baseline band. This observation suggests that distinct aspects of the network may be indicative of the overall effects generated at the baseline and excited frequencies.
What properties of the system drive or predict the diverse, distributed responses in the baseline and excited frequency bands brought about by focal stimulation? Because the network of anatomical connections couples different brain areas and allows them to directly interact, it is reasonable to hypothesize that the organization of this network should play a role in guiding the influence of a perturbation. To test this hypothesis, we study 〈|Δρδibase|〉 and 〈|Δρδiexc|〉 as functions of structural node strength sistruc (Fig 5E). Interestingly, the global PLV modulation induced in the network’s naturally-emergent frequency band is not well-predicted by the anatomical strength of the stimulated area (Fig 5E, Left). In contrast, though, we do observe a strong association (Spearman correlation rs = 0.96, p < 0.001) between structural strength and the PLV change elicited in the excited band (Fig 5E, Right). This result indicates that more structurally connected units generate larger overall effects at the enhanced frequency of the directly stimulated area. Because the excited band response is strongly constrained by structure, we also examined whether the effects differed between two broad, anatomically-defined classes of nodes. In particular, we compared the average excited band response for stimulation of cortical vs. subcortical areas (both of which are included in the anatomical parcellation). Given this breakdown, we find that the overall effect is significantly higher upon perturbation of subcortical regions (see Fig. E in S1 Text). This result is consistent with the findings of [35, 79], and reflects the notion that subcortical nodes make strong, distributed structural connections that may support large-scale network communication [104].
Although the brain’s structural connectivity plays a crucial role, macroscale activity patterns generally reflect an interplay between connectome architecture and the network’s dynamic regime. Indeed, in Fig 3D we observed that the correlation between structural and functional strength varies in intensity and sign with working point. Importantly, the presence of a functional connection between two brain areas implies an interdependence of their dynamics—enforced by the system’s oscillatory state—that can occur even in the absence of a direct structural connection. Intuitively, we may thus expect the organization of the system’s initial functional interactions (which could be non-trivially related to structure), to be indicative of how the coherence pattern is modulated under perturbation. Given this reasoning, we study 〈|Δρδibase|〉
vs. sifunc (Fig 5F). Relative to structural strength, we observe a strong positive relationship (Spearman correlation rs = 0.71, p < 0.001) between the average absolute change in baseline band coherence and functional strength. Thus, areas that are initially more coherent with other regions in the network tend to yield larger global modulations to baseline band interactions when perturbed. This should be contrasted to the results from structural node strength, for which there was not a strong relationship with absolute coherence changes at the baseline oscillation frequencies. Finally, we consider 〈|Δρδiexc|〉
vs. sifunc (Fig 5F, Right). Though they are correlated, (rs = −0.34; p < 0.05), the stimulation-induced responses in the excited band are much more strongly predicted by structural rather than functional strength.
The above results emphasize separate consideration of both anatomical network topology and the organization of emergent functional interactions, the latter of which is also driven by the dynamical regime of the system. In particular, for the working point considered here, structural and functional network properties relate to distinct types of perturbation-induced effects. First, phase-locking changes that arise in the excited band reflect the transmission and replication of oscillatory input from the directly excited area to and in downstream regions. If the structural connection between the stimulated site and a downstream area is strong enough, then the drive from the stimulated site will induce a new spectral component in the receiving area (see, e.g., Fig 4C); consequently, the two regions will exhibit phase-locking at the excited frequency. In addition, even two areas that are not directly linked can display a high PLV in the excited band due to strong common input from the stimulated region, or due to the propagation of the stimulated site’s signal along alternative paths in the network. In sum, because spreading of the perturbed area’s activity is highly constrained by the presence of structural connections, regions with stronger anatomical connectivity to other areas more forcefully drive downstream regions and lead to larger excited band effects. Perhaps more interesting are the modulations in coherence that occur in the baseline frequency band. These changes arise not due to a direct transmission of input, but rather via adjustments to the ongoing, mutual entrainment between units’ spontaneous rhythms. For WP1, the resulting alterations to the strength of coherent interactions are more related to the stimulated region’s initial functional connectivity rather than its anatomical strength. Intuitively, this may in part be due to the fact that perturbing a particular area tends to decouple it from other areas at the original oscillation frequencies, such that stimulating regions that are strongly coherent to begin with effectively reconfigures existing functional interactions in the baseline frequency band. Notably, the observed correlation between functional strength and baseline band coherence modulations does not uncover the deeper, precise mechanisms behind the effects. However, because the brain’s collective dynamics can reflect a complex interplay between its oscillatory state and its structural connectivity, the result highlights the importance of considering both aspects when trying to understand network-wide responses to perturbations.
Importantly, the collective state of the model can change even when anatomical connectivity is fixed (see Fig 3). To explore how this affects the impacts of regional perturbations on brain network dynamics, we next examine another working point—WP3—located at C = 2.5 and Pbase = 0.7. (Note that a working point WP2 between WP1 and WP3 is analyzed carefully in Sec. SIII of S1 Text). In this high background drive state, regional activity is characterized by more regular and higher-amplitude oscillations (see Fig 2G and 2H) relative to both WP1 and WP2. Furthermore, at WP3, the system resides well beyond the point of maximal coherence (Fig 3A), and the baseline PLV matrix is more constrained by anatomical connectivity (Fig 3C, Row 2, Column 4; Fig 3D).
Inspection of a region’s activity time-series at baseline and when stimulated with additional input indicates noticeable differences in how perturbations alter local activity at WP3 versus at either WP1 or WP2. Specifically, when the system operates in the high-drive state, a perturbation of the same strength has a much less drastic effect on the stimulated region’s activity, inducing only relatively small changes to its amplitude and frequency (Fig 6B).
Fig 6 Effects of local excitations on power spectra are more restricted at the high background drive working point (WP3).(A) Schematic of a brain network depicting the stimulated site i in brightest red. The black arrows point to two other regions j and k that lie at progressively further topological distances from the perturbed area in the structural network. In this figure, regions i, j, and k correspond to brain areas 1 (R–Lateral Orbitofrontal), 4 (R–Medial Orbitofrontal), and 10 (R–Precentral), respectively. (B) Left: A segment of region i’s activity time-course in the baseline condition. Right: A segment of region i’s activity time-course when it is stimulated. (C) Power spectra of area i and two other downstream regions j and k. In all three panels, the lighter curves correspond to the baseline condition, and the darker curves correspond to the state in which i is driven with additional input. The gray vertical lines indicate the peak frequency fi,δipeak of region i in the excited condition. (D) Histogram of the shift in peak frequency Δfi,δipeak induced by exciting unit i, plotted over all choices of the perturbed area. (E) The average shift in the peak frequency of the stimulated region 〈Δfi,δipeak〉 for WP1, WP2, and WP3 (error bars indicate the standard deviation over all choices of the excited unit). (F) Distribution of peak frequencies of all units in the baseline condition {fipeak} (light gray) and distribution of the peak frequency units acquire when directly excited {fi,δipeak} (dark gray). (G) Average power spectra 〈psd〉
j≠iover all units j ≠ i at baseline (light gray) and when unit i is perturbed with additional input (dark gray). Because stimulation does not induce a well-separated excited frequency band, we only assess perturbation-induced changes in the PLV between brain areas for a single baseline frequency band (purple area).
Consequences of regions’ enhanced baseline activity for focal stimulation are perhaps more evident from examples of areas’ spectra at baseline and under stimulation (Fig 6C). As for the lower-drive working points, the peak frequency and power of the stimulated region i shift to higher values (Fig 6C, Left). However, at WP3, the increase is modest relative to the shifts that occur at either WP1 or WP2, and no modulation sidebands arise in i’s spectra under excited conditions. As a result of the more unyielding nature of spontaneous dynamics, stimulation of unit i also has relatively little impact on the spectra of downstream regions (Fig 6C, Middle, Right), even if they are positioned topologically close to the perturbed site. To more generally quantify the effects of regional stimulation on areas’ power spectra, we examine the distribution of the shifts in peak frequency Δfi,δipeak that occur due to perturbation of each unit. The largest of these shifts is only about 3Hz (Fig 6D). Hence, relative to WP1 and WP2, the average shift in peak frequency 〈Δfi,δipeak〉 is greatly reduced at the high-drive working point (see Fig 6E). Furthermore, unlike the situation in the low-drive state, the distributions of peak frequencies at baseline and under focal stimulation begin to overlap at WP3 (Fig 6F), precluding the notion of separate baseline and excited frequency bands. For this reason, in our subsequent analyses we only consider phase-locking changes inside a single frequency band (Fig 6G). While we refer to this as the “baseline band”, we note that it still contains the peak frequency of the directly excited unit, since its frequency shift is so small.
The results presented in this section indicate that regional dynamics are more robust to perturbations at the high-drive working point. This can in part be understood by considering the effects of the background drive PEbase, which is the parameter tuned to move from WP1 → WP2 → WP3. In particular, the increased baseline input level at WP3 means that each unit, if disconnected from the network, would operate closer to the bifurcation separating the quiescent and oscillatory state than would be the case at WP1 or WP2. As a result, stronger oscillations emerge at WP3 when the network coupling is introduced, reflecting the increased influence of recurrent dynamics. The high-amplitude rhythms that arise in the high-drive state are more difficult to disrupt, leading to minimal changes in the power spectra under local perturbations. For the same reasons, it is also more difficult for a local change in activity to propagate and influence the dynamics of remote areas. In contrast, when the system operates at either WP1 or WP2, the baseline oscillations at each brain area are weaker. This lower-amplitude activity is easier to override, yielding the system more plastic and susceptible to local perturbations. This flexibility at WP1 and WP2 is reflected by clear modifications to regional spectra upon local stimulation and the signatures of the stimulation effect in downstream regions (Fig 4 and Fig. C in S1 Text).
While the rigidity of baseline rhythms at WP3 prevents the emergence of a well-defined excited frequency band, we can still assess the effects of selective perturbations on interareal phase-locking in the single, baseline frequency band. We carry out such an analysis in this section, focusing on contrasting the results obtained at WP3 to those obtained previously at WP1 (and WP2). To begin, we show examples of the pairwise changes in PLV induced by stimulation of two different brain areas i and j (Fig 7A), and the average absolute changes 〈|Δρδibase|〉 in phase-locking driven by perturbation of each brain area in the network (Fig 7B). As with WP1, the statistical significance of PLV changes is first assessed with a phase-randomized null model (see Sec. SXII of S1 Text) before computing network averages. Despite differences in the induced average responses across regions, visual comparison of the distribution at WP3 (Fig 7B) and at WP1 (Fig 5B) suggests that there may be less variation across regions when the system operates in the high-drive state. Indeed, although the mean of the average absolute changes across all brain areas is approximately the same at WP1 and WP3, when we compare the coefficient of variation (CoV) of the two distributions of {〈|Δρδibase|〉}, we find that the CoV for WP3 is 0.25, while for WP1 it is 0.45. This indicates that global coherence modulations elicited by regional stimulation are more homogenous across different choices of the stimulated brain area when the system operates in the high background drive regime. Because we are particularly interested in how perturbations may differentially affect network dynamics depending on baseline state, we also investigate if the average absolute PLV modulations are correlated between WP1 and WP3. This relationship is not significant (Fig. G in S1 Text), indicating that the rank-ordering of 〈|Δρδibase|〉 is quite different between the low and high drive regimes.
Fig 7 Phase-locking modulations induced by regional stimulation at WP3 and their associations with network properties.(A) Pairwise changes in the PLV inside the baseline band Δρ^base^ when region i (Left) or region j ≠ i (Right) is perturbed. In this figure, regions i and j correspond to regions 10 (R–Precentral) and 15 (R–Isthmus), respectively. (B) Network-averaged absolute PLV changes in the baseline band 〈|Δρδibase|〉 induced by stimulation of different brain areas. (C) The quantity 〈|Δρδibase|〉 vs. structural node strength sistruc (Left) and vs. functional node strength sistruc (Right). Insets indicate Spearman correlation coefficients between the plotted quantities and their associated p-values).
We next considered the associations between structural or functional properties of network nodes and the global phase-locking changes induced by regional stimulation (Fig 7C). At WP3, we find a strong positive correlation between 〈|Δρδibase|〉 and structural strength sistruc (rs = 0.82, p < 0.001; Fig 7C, Left). Thus, larger coherence modulations at the high-drive working point are dependent on the presence of strong anatomical connections emanating from the stimulated region. Note that functional strength sifunc also exhibits a positive correlation with 〈|Δρδibase|〉 (rs = 0.34, p < 0.05; Fig 7C, Right), which we might expect given the weak but positive correlation between s^struc^ and s^func^ at this working point (Fig 3D). However, structural connectivity is undoubtedly a more robust predictor of these effects.
A supplementary analysis (see Fig. H in S1 Text) indicates that the association between sistruc and 〈|Δρδibase|〉 is (largely) explained by a relationship between structural strength and the bulk decreases in coherence induced by perturbations. Hence, units with stronger anatomical connectivity to the network as a whole generate larger global breakdowns in coherence when stimulated with constant, additional input. To gain intuition for this result, recall that in the high-drive regime, brain areas exhibit relatively strong and inflexible oscillations at baseline. Furthermore, though not enough to bring about an entirely new excited band, stimulation still slightly accelerates the frequency of the perturbed unit. As a consequence, we may expect stimulation of regions with stronger structural connections—which have more direct influence on other areas—to drive more disruptions in the network’s ongoing dynamics. That is, a given node will tend to be driven out of coherence with other units when stimulated (since it acquires a slightly increased frequency), and high-strength nodes in particular will also be able to decouple other areas from their baseline functional assemblies, in turn causing larger dissociations of functional connectivity.
To conclude this section, we highlight the observation that depending on the oscillatory regime, different aspects of the system best predict global responses to focal stimulation. In the low-drive state, emergent coordination between units activities at baseline is more strongly related to coherence modulations at regions’ spontaneous frequencies, whereas structure is more strongly associated with the effects in the high-drive regime.
We have examined distinct operating points characterized by varying regional oscillation strengths, and in turn, different large-scale coordination patterns and relationships to anatomical connectivity. We found that focal perturbations had markedly different effects on brain network dynamics and associations to network connectivity depending on the system’s oscillatory mode. In this section, we more generally examine the effects of perturbations as the nature of regional activity is smoothly varied—by tuning the level of background drive PEbase—while holding network structure constant. Throughout this exposition, we set the coupling C = 2.5, as before.
To summarize the response in the baseline band, we first stimulate each region i to generate a set of N average absolute coherence modulations {〈|Δρδibase|〉}. We then consider the mean of this set across regions, 〈|Δρ^base^|〉 (which we refer to as the “grand average”), and study its behavior as the background drive is varied (Fig 8A,right axis). We find that the general shape of 〈|Δρ^base^|〉 as a function of the relative background drive PEbase-PE* tends to mimic that of the global baseline coherence ρ^global^ (Fig 8A,left axis), but with a slight shift towards lower PEbase-PE*. Importantly, the peak in the PLV modulation curve at intermediate background drive signifies a distinct state at which regional perturbations generate the largest overall changes to coherence in the system’s baseline frequency band. Of note is that this working point occurs just prior to the global PLV peak. Intuitively, the system may exhibit the largest response at this working point because there is heightened potential to both enhance and depress functional interactions when the system is perched at the transition into the maximally ordered state. At lower values of PEbase-PE*, stimulation of single brain areas has a smaller average effect, which likely occurs in part because there is less initial coherence for local excitations to disrupt. As ρ^global^ begins to decline with increasing PEbase-PE*, so too does the overall response to regional perturbation. The grand average of the baseline band PLV modulation approaches a local minimum at PEbase-PE*=0.07, after which it settles to intermediate values at high-drive working points well beyond peak ρ^global^ (e.g., at WP3).
Fig 8 Dependence of global phase-locking changes on the baseline state of the brain network model.(A) The left axis (gray) shows the network-averaged baseline PLV ρ^global^ as a function of background drive PEbase-PE* for a coupling C = 2.5. The right axis (purple) shows the grand average 〈|Δρ^base^|〉 of the perturbation-induced absolute changes in baseline band PLVs as a function of PEbase-PE*. (B) The left axis (gray) shows the network-averaged baseline PLV ρ^global^ as a function of background drive PEbase-PE* for a coupling C = 2.5. The right axis (purple) shows the coefficient of variation of the perturbation-induced average absolute changes in baseline band PLVs, CoV[〈|Δρδibase|〉], as a function of PEbase-PE*. (C) The left axis (gray) shows the network-averaged baseline PLV ρ^global^ as a function of background drive PEbase-PE* for a coupling C = 2.5. The right axis (green) shows the grand average 〈|Δρ^exc^|〉 of the perturbation-induced absolute changes in excited band PLVs as a function of PEbase-PE*. (D) The left axis (gray) shows the network-averaged baseline PLV ρ^global^ as a function of background drive PEbase-PE* for a coupling C = 2.5. The right axis (green) shows the standard deviation of the perturbation-induced average absolute changes in excited band PLVs, std[〈|Δρδiexc|〉], as a function of PEbase-PE*. (Note that here we consider the standard deviation rather than the coefficient of variation since the mean response in the excited band eventually drops to zero.).
In addition to the grand average, we also examined the coefficient of variation (CoV) of the distribution of average absolute coherence changes, which we denote as CoV[〈|Δρδibase|〉]. In general, CoV[〈|Δρδibase|〉] also exhibits a state dependence (Fig 8B). Specifically, CoV[〈|Δρδibase|〉] begins at an intermediate value for the lowest background drive, increases to a global maximum, and then declines to below its initial value for high-drive states. This behavior indicates that global responses exhibit less variability across different choices of the stimulated site when the system operates in states of strong baseline rhythms. In contrast, the network-wide impacts of regional stimulation display more dispersion in the low-drive regime, with variability peaking near maximum ρ^global^ (WP2). For these states, the overall response in the baseline band is thus more dependent on the precise location in the network that is perturbed.
We next assess how the dynamical state influences the global response induced in the excited frequency band. In this case, the grand average 〈|Δρ^exc^|〉 of the absolute excited-band PLV changes steadily decays and eventually vanishes with increasing PEbase-PE* (Fig 8C), as does the variability of the response distribution (Fig 8D). This behavior is due to the amplification of baseline oscillations with increasing background as regions’ spontaneous dynamics become more difficult for perturbations to override, the emergence of widespread phase-locking at regions’ excited frequencies declines. Eventually, stimulation no longer yields large enough frequency shifts to induce an excited band at all. This result thus indicates that activity from the stimulated area is most effectively propagated to downstream regions when the network operates in a state of weak baseline activity.
For our final analysis, we examine more generally how the relationships between stimulation-induced PLV modulations and structural or functional network connectivity depend on the system’s operating point. In particular, we consider the correlations between the average absolute baseline band coherence changes 〈|Δρδibase|〉 and the structural sistruc or functional sifunc strength of the perturbed unit. To highlight the differing levels of association between the global responses and either structural or functional connectivity, Fig 9A shows how the difference Δrs in the strength of these correlations varies with background drive PEbase-PE* (individual correlations are shown in Fig. J of S1 Text). We observe that the relationships between phase-coherence modulations and structural or functional strength can depend strongly on the system’s baseline state. For the bulk of the “low”-drive regime (working points below peak global coherence), there is a stronger relationship between 〈|Δρδibase|〉 and sifunc (rather than sistruc). This enhanced association between a region’s initial functional connectivity and the overall response it induces upon stimulation is largest at dynamical states near WP1, where there is an intermediate level of coherence in the network at baseline. In the “medium”-drive regime, stretching from just prior to peak global coherence up to the beginning of the plateau region in Fig 3A, the correlation between sistruc and 〈|Δρδibase|〉 starts to increase, but functional strength continues to remain more strongly associated with the global modulations elicited in the baseline band. Finally, we observe a clear “switch” at PEbase-PE*≈0.07, after which structural connectivity becomes relatively better at predicting the overall response to perturbation than functional connectivity (note that this transition occurs in tandem with baseline functional and structural strength becoming positively correlated; see Fig 3D). This behavior marks the onset of the “high”-drive regime, and persists across the remainder of the background drives examined.
Fig 9 Relationships between phase-locking modulations and the structural or functional connectivity of the stimulated site vary with working point.(A) Difference Δr
sin the strength of the correlation between the average absolute baseline band PLV changes 〈|Δρδibase|〉 and structural (sistruc) or functional (sifunc) node strength, plotted as a function of the baseline drive PEbase-PE* for a coupling C = 2.5. The difference is defined such that when the curve is positive, overall coherence modulations exhibit a stronger correlation with functional rather than structural strength. The arrows mark the locations of the different working points studied in detail in the main or SI Text. (B) A schematic summarizing how structural s^struc^ or functional s^func^ strength are related to either baseline or excited band PLV changes for different dynamical regimes. As the background drive varies from low (WP1) to medium (WP2) to high (WP3), the oscillatory state of the system changes, and so does the association of different phase-locking modulations to structural or functional network properties.
To conclude this section, we schematically summarize in Fig 9B the most robust associations between either structural or functional connectivity and the PLV modulations induced by focal perturbations in either the baseline or excited frequency bands. Although we have not comprehensively detailed the effects of perturbations for all possible working points of the model, the results discussed throughout the text highlight the critical influence of a brain network’s collective dynamical state in dictating the outcomes of localized stimulation.
In this study, we set out to explore relations between large-scale brain connectivity, dynamics, and the local and widespread impacts of regional perturbations to neural activity. Following the efforts of past work [35], we built a reduced computational model wherein brain areas were represented as Wilson-Cowan neural masses with long-range coupling between regions constrained by empirical diffusion tractography measurements. We then investigated how stimulation of a particular brain area affected network dynamics, and asked whether and how the collective dynamical state of the system plays a role in modulating such effects. Here, we chose to examine state-dependence by changing the combination of generic background drive and interareal coupling strength, which vary the extent to which local vs. network interactions drive regional dynamics, respectively. By tuning these parameters, we identified qualitatively distinct dynamical regimes of the model, and at these different working points, we assessed how local excitations of fixed strength induced or modulated interareal phase-locking. We found that, depending on the baseline regime of the system, the network exhibited different responses to regional perturbations. Furthermore, altering the working point of the model also qualitatively altered the relationships between stimulation-induced effects and properties of structural and functional network connectivity. To the best of our knowledge, these points have yet to be investigated at the whole-brain scale via computational modeling.
Before moving on to a further discussion of this study’s results and limitations, it is critical to state that it is not the first to examine the impacts of perturbations using computational models of brain dynamics. In fact, our work was motivated and inspired by a number of previous large-scale modeling studies that have uncovered key insights into either the network-wide [35–38] or state-dependent [48–50] influences of stimulation. Because our study builds on its methodology, we specifically emphasize the investigation in [35], where (as here) the authors used interconnected Wilson-Cowan units to simulate brain activity, and then utilized recent advances in network control theory [105] to make predictions about the overall functional effects of regional stimulation. We also bring particular attention to the study conducted in [48], where the authors built a simplified model to examine the state-dependent effects of alternating current stimulation on local cortical oscillations. As part of their study, they found that the strength of endogenous oscillations altered the susceptibility of local dynamics to external stimulation. In the present investigation, we vary the background drive (the main tuning parameter in the WC model), which also modulates the amplitude of regional activity. Therefore, although our study is focused on the network-wide effects of focal perturbations, it is necessary to heed and connect our results to those previously reported in [48].
While our work builds upon these and other past efforts, it is important to highlight some key distinctions and extensions of our analysis. First, in contrast to [36], we opted to use the more complex—but in some ways more biophysically-motivated—Wilson-Cowan system as the fundamental dynamical unit in the model, rather than a pure phase oscillator. Studying models that incorporate both phase and amplitude dynamics will likely be critical for a more complete understanding of oscillatory neural activity, since variations in amplitude can affect functional couplings in the brain, are modulated by behavioral conditions, and are also elicited by various forms of stimulation [68, 102, 106–108]. A second key difference between our study and those conducted in [35] and [37] is that the aforementioned studies analyzed a situation in which brain activity was assumed to be in a quiescent state prior to stimulation. Here, however, we were interested particularly in states corresponding to ongoing oscillatory activity, and the interaction between those baseline rhythms and increased excitation to a particular brain area. Finally, in an extension of prior work that has begun to examine the state-dependent impacts of stimulation on single cortical areas [48], we wished to consider a large-scale network model that allowed for an analysis of how stimulation can disseminate to induce or modulate dynamical interactions between widespread brain areas. Hence, the key contribution of this study is its simultaneous investigation of (1) not only the focal, but also the distributed impacts of regional stimulation, and (2) how the collective regime of system activity (i.e. that arising from a combination of local dynamical properties and network coupling) influences such network-wide effects. To the best of our knowledge, it remains an open question how the collective state of brain activity alters the way focal stimulation impacts functional relations between distant brain areas.
In the rate model implemented here, the general focal effect of stimulation was to increase the amplitude and frequency of activity in the perturbed area. Such changes in neural activity are broadly consistent with the effects of natural stimulation of certain cortical areas, such as visual stimuli impinging on visual cortex [16, 102] and the effects of increased excitatory drive to neuronal populations generally [101]. However, we also observed that an excitation of fixed strength had significantly stronger effects on local activity (i.e., the induced shift in power and frequency was larger) for states of lower background drive (e.g., at WP1 and WP2). Thus, for these working points, perturbations led to the emergence of a new, “excited” frequency in the system that was well-separated from the main oscillation frequencies at baseline. In contrast, the enhanced regional activity present in the high-drive regime (WP3) was significantly less responsive to perturbations (in that stimulation did not generate a well-separated excited frequency). Though the analogy is not perfect, this behavior is undoubtedly akin to and consistent with the findings reported in [48], wherein alternating current stimulation induced the strongest effects at the stimulation frequency in the absence of strong endogenous oscillations.
Beyond local responses, it is also critical to acknowledge that a given neuronal population exists in the context of a larger network of areas, such that local changes in activity can induce distributed effects [11–14]. Indeed, we found that in states of low background drive, areas topologically close to the stimulation site also developed spectral peaks at the frequency of the stimulated region. These downstream signatures—caused by strong oscillatory drive from the directly excited area—are thus qualitatively similar to those induced by rhythmic stimulation, where activity in the targeted site shows enhanced power at the stimulation frequency [48, 109]. In general, the network’s global response in the excited frequency band was tightly constrained by anatomical connectivity, demonstrating that these effects arise due to a direct propagation of the stimulated area’s strong rhythmic activity along structural pathways. Interestingly, we found that focal perturbations could also modify downstream regions’ spectra at frequencies other than the peak frequency of the directly excited area. Such generic broad-band alterations of spectra can occur in coupled neuronal populations with long-range excitatory interactions [110]. Indeed, in one of only a few possible scenarios, interacting nonlinear oscillators engage first into quasiperiodic dynamics, and then develop chaos, which is associated with spectral changes over continuous ranges [103]. We additionally saw that local excitations could alter levels of temporal coordination between brain areas’ activity in a frequency band containing the system’s spontaneous, “baseline” oscillations. Furthermore, in the low-drive regime, we observed that global absolute modulations in baseline-band coherence were actually more strongly related to the stimulated region’s functional strength than structural strength. This is in contrast to the excited-band responses, which were clearly mediated by anatomical connectivity. Importantly, changes in the baseline band arise from modifications of the dynamic interactions between brain regions’ ongoing rhythms, which can depend on a more complex and intricate interplay between both network structure and the collective oscillatory state of the system. Finally, in the high-drive regime, system-wide changes in phase-locking became more homogenous across different choices of the excited brain area, and were better predicted by the stimulated region’s structural strength. This reflects the idea that in a state of robust regional oscillations, focal stimulation will only yield a large global effect if the perturbed site has strong anatomical connectivity.
Although the field has not reached a clear consensus on the distributed effects of stimulation, the general finding in our study that increased drive to a single brain area can spread and reorganize functional interactions is consistent with previous modeling work [35, 37]. We also make some new, specific predictions about different types of changes that can occur and investigate potential mediators of these effects. Perhaps most interesting is the suggestion that both network structure and coordinated dynamical organization may play a role in guiding the non-local effects of perturbations. While network neuroscience has traditionally focused on how structure can predict function [74, 111–114], it is critical to acknowledge that the behavior of neural systems need not be completely constrained by structure alone [40, 88, 99, 115]. In particular, by modulating intrinsic properties of neural units, a single structural connectome can generate different patterns of coordinated activity [100, 116, 117]. Indeed, we observed that although structural connectivity may partially constrain the organization of functional interactions, the way and the extent to which it does so is dependent on the nature of regional dynamics. Our findings indicate that, as a consequence, the collective oscillatory state of the system may also be important in determining how a perturbation can lead to distributed changes in functional interactions. Specifically, we found that in certain cases, functional connectivity was actually a better predictor of the overall response to perturbations than structural connectivity. Although the presence of an association with functional strength should not be misinterpreted as a mechanistic explanation for how stimulation alters dynamics, such an observation can provide insight, especially when contrasted with the role of structure alone. In total, our results indicate that depending on the collective state of network activity, the widespread impacts of focal perturbations can differ, and may be driven by distinct processes.
Though we examined a simplified model, our results may be relevant to other work and a growing literature on brain state-dependent stimulation [43–50], which recognizes that the effects of exogenous perturbations can be conditional on the endogenous rhythmic or spontaneous activity of the system at the time of stimulation. For example, empirical studies have shown that outcomes of stimulation can differ as a function of cognitive state (i.e., task vs. rest) [44–47, 118]. Importantly, a few computational studies have investigated and provided explanations for such state-dependent responses to alternating current stimulation [48–50], though have only considered models of a single or a few coupled cortical and/or subcortical regions and have thus focused on stimulation’s local effects. Our goal here was to begin filling the need for larger-scale models that systematically consider the influence of whole-brain anatomical connectivity in conjunction with the collective dynamics of brain activity in shaping the widespread effects of stimulation. Importantly, stimulation also holds promise as a technique for moving the brain between specific, desired states [41, 42, 105]. The ability to control brain network dynamics has important applications in the treatment of neurological diseases such as Parkinson’s and epilepsy [32–34], in which neural activity goes awry. As they are refined, computational studies that put focus on how stimulation interacts with internal, ongoing dynamics in large-scale brain networks have the potential to inform future control and stimulation-based experiments and therapies [119, 120].
There are also a number of methodological considerations to comment on regarding this work. First, we used a relatively coarse-grained parcellation (N = 82 regions) to construct human structural brain networks. This resolution is roughly consistent with several other whole-brain modeling studies [35, 66–69], and is tractable for computationally intensive simulations. However, the employed parcellation is a simplification of the underlying anatomy, since the regions represent large pieces of neural tissue and hence remain agnostic to potentially important structural heterogeneities at finer scales. Moreover, we used Euclidean distances between region centers to approximate structural connection lengths [66, 69], but one could derive perhaps more accurate estimates based on fiber trajectories [82]. It is also important to note that the limitations of human brain imaging and tractography preclude a perfect reconstruction of interareal connections [121]. One primary drawback is that these methods cannot resolve the directedness of interareal connections, which could impact subsequent results [70]. Finally, we used a group-representative connectome in this study. On the one hand, this allowed us to focus on general trends and behaviors, but on the other, it leaves no room for examining how individual differences may affect certain findings. In future work, it will be interesting to explore how various results generalize to both higher-resolution and higher-quality brain data, and to understand how variability in brain structure across different human subjects [78, 122, 123] or even across different species [124] may relate to differences in how perturbations of neural activity are expressed in system dynamics.
Another limitation of this study concerns the model used to simulate brain activity. We opted to use the canonical Wilson-Cowan model [53]—which embodies a tradeoff between biological realism and tractability—but future work could examine other choices for the node dynamics. Futhermore, while the WC model constitutes an arguably realistic improvement over simpler linear models, phase-oscillators, or generic Hopf bifurcation models, it is still phenomenological in nature and our implementation of the model is highly idealized. For example, the simulated oscillations in this work tend to be much more regular than physiological rhythms, which limits the model’s biological plausibility. Similarly, the model yields unrealistically high levels of sustained, network-wide coherence. If observed in brain network dynamics, such states would be much more likely to occur in a transient manner as the brain transitions between segregated and integrated states dynamically [96, 125–127]. Though building whole-brain mean-field models that yield more realistic, transient and time-varying switches between synchronized or desychronized states is difficult, this is an exciting direction for future study.
In general, there are several ways the model could be improved upon to increase its biological validity. For instance, for simplicity and in line with past work [35, 66, 74–77, 81], we only considered long-range couplings between excitatory populations. However, interareal connections likely target both excitatory and inhibitory neurons, and it would therefore be relevant to examine the effects of long-range excitatory-to-inhibitory coupling. We also used a fixed value for the signal propagation velocity. Although we attempted to choose an empirically constrained value, it is known that delays can have significant consequences on brain dynamics [73, 128], and this could be explored further. Moreover, in the WC model, oscillations are generated via the interaction of excitatory and inhibitory neuronal populations, which is considered a biophysically-plausible mechanism underlying rhythmicity in the gamma and beta bands [4, 129]. However, lower frequency oscillations (e.g. alpha band activity), may be generated—at least in part—by thalamocortical loops [130]. Importantly, past models of small thalamocortical circuits (i.e. with only a few coupled regions) have shown that these interactions are important for explaining empirical results regarding the state-dependent responses of cortical alpha oscillations to AC stimulation [48]. Along these lines, an interesting direction for future work would thus be to incorporate more realistic dynamics of subcortical regions and biologically-motivated thalamocortical couplings [131]. Another assumption of the model is that each unit has identical parameters. While this is a reasonable and useful setup to analyze first, there is an important body of literature detailing specific heterogeneities across brain areas. For example, different regions may operate at different intrinsic time-scales and have different intrapopulation architectures (e.g., excitatory or inhibitory coupling strengths) that lead to functional specialization [132]. Recent modeling efforts have begun to incorporate some of these additional complexities, finding that doing so can lead to more realistic baseline dynamics and can explain certain empirically-observed behaviors not accounted for by simpler models [84]. The model implemented here also assumes oscillations in a single gamma frequency range. However, brain areas can exhibit rhythmic activity in multiple frequency bands, which potentially arise from different cortical layers [133]. In future work, it will be important to ask how the impacts of perturbations depend on collective state when additional details about the heterogeneity of the underlying anatomy or dynamics of brain areas are included in the model [36, 89]. Indeed, these steps will allow one to make more concrete and biologically-meaningful statements about how different regions influence large-scale activity patterns. That said, we also stress that the model considered in the present work does not prevent the ensuing dynamics from being complex and rich, and actually allows us to appreciate how dynamical—in addition to structural—complexity can be a key driver of stimulation-induced effects.
With respect to our analysis, it is also important to state that we examined only a few representative working points of the model, out of the many that exist. Though our goal in this work was to illustrate that the collective state of the system can influence the effects of local perturbations, we do not claim to have provided an exhaustive description of all possible behaviors. Furthermore, our analyses of the interareal coherence modulations induced by focal perturbations examined network-averaged responses and absolute changes. In other words, we focused on characterizing aggregate reorganizations in functional couplings. This focus leaves at least two clear and important directions for subsequent work. First, one could try to further disambiguate what leads to strengthening vs. weakening of interareal coherence upon focal perturbation. A deeper understanding of when each outcome occurs would likely require significantly more complex analyses of not only the structural and dynamical properties of the stimulated node, but also of other areas. Second, one could attempt to discern how focal perturbations alter dynamic interactions on finer scales, for example, between individual pairs of regions. Again, this would almost certainly require investigating which aspects of non-stimulated areas at baseline allow them to become engaged or disengaged when another region is perturbed. The complexity of both network structure and dynamics make these tasks challenging, but obtaining a more refined understanding of the underlying mechanisms is essential for making useful comparisons against empirical studies.
Perhaps the biggest limitation of the model presented here is that its baseline dynamics were not tuned according to empirical measurements of large-scale brain activity. To address this in future work, model parameters could be fit separately for different empirically-observed brain states [134], such that various aspects of simulated whole-brain activity match the experimental activity patterns in each case. Then, the effects of stimulation could be re-investigated in the context of the empirically-validated large-scale model to make more biologically-informed predictions. Because the current implementation lacks these constraints, it should be interpreted as a canonical, reduced example.
Another interesting direction for future work would be to systematically compare the impacts of perturbations when implemented on empirical brain network architecture versus other canonical network models. This may allow one to test which effects generalize across different network topologies and which are specific to brain network organization. For example, the recent study in [135] examines how human cortical networks may be structured to support ignition dynamics relative to other canonical network models, finding that the well-connected structural core present in actual brain networks plays a critical role. In extending this type of analysis, one could perhaps also consider whether the brain is especially well-organized for oscillatory control via localized perturbations, and how the brain’s controllability depends on its current state. In other words, does the brain exhibit special topological or dynamical features that enhance the ability of focal stimulation to push it into particular collective, oscillatory modes? Importantly, the controllability of large-scale brain networks is often studied using tools from linear network control theory [105]; however, understanding the control of oscillatory activity in the brain will require extending this framework to incorporate non-linear behaviors [136].
In addition to investigating alternative and/or more detailed models, forthcoming studies could also consider different ways of operationalizing perturbations to neural activity. As a first step, here we implemented “stimulation” as a constant increase in excitatory drive [35]. Given the already complex nature of network dynamics in the absence of stimulation, starting with simple, constant perturbations is critical for building intutions. However, it will be important for future whole-brain modeling studies to implement more realistic and experimentally viable forms of perturbations, which are often time-varying and transient in nature. One clear choice would be to model alternating current stimulation, for example. In this case, the stimulated area receives oscillatory input of a particular amplitude and frequency [137], both of which can be independently tuned. While this kind of stimulation has been examined in smaller thalamocortical circuits [48–50], it would be exciting to scale up to systems-level models. One could also investigate alternative ways of modulating brain state. In this study, we varied the baseline dynamical regime in perhaps the most straightforward way possible, by tuning the level of excitatory input globally for all network elements. This variation changed the local dynamics of each brain area, and in turn, the macroscopic state of the system as a whole. However, brain state could also be modulated by tuning a different physiologically-interpretable parameter, such as the gain in the sigmoidal activation functions. Indeed, recent modeling studies have shown that altering neural gain can lead to dynamical regimes in which functional integration and segregation are balanced [138]. It would also be interesting to understand how widespread changes of neural gain affect the way focal perturbations materialize.
In sum, we conducted an idealized investigation into the effects of focal stimulation on brain network dynamics, focusing on how the system’s collective state influences the distributed impacts of such perturbations. To the best of our knowledge, this latter point has only recently begun to be examined in the context of large-scale brain networks, and therefore warrants investigation via simplifed models. However, the results of this study must be interpreted cautiously as they are yet to be validated by data. Comparing conjectures based on reduced models against empirical findings is necessary to substantiate whether or not the model provides biologically meaningful insight. One testable prediction from the model presented here is that focal stimulation propagates to cause downstream modulations of power and phase-locking at both the dominant frequency acquired by the activated area, but also at regions’ baseline frequencies. Investigating this would require an experiment in which stimulation forces the perturbed area to oscillate at a well-defined frequency—which could perhaps be achieved by alternating current stimulation [7]—and a simultaneous measurement of other brain areas’ dynamics. Moreover, testing variations in the global effects of perturbing different regions would require multiple brain areas to be stimulated and responses from distributed regions simultaneously recorded. Recent advances combining non-invasive brain stimulation with measurement modalities like EEG and MEG [26–30] are making this increasingly possible. These types of experiments [26–30] will be crucial in validating results from computational models that posit how stimulation alters functional couplings across different parts of the brain. Furthermore, one could test how widespread changes of power and phase-locking depend on baseline brain activity, for example, during resting vs. task conditions. Notably, some experimental studies have indeed begun to test how brain state affects the local outcomes of stimulation [43–47], which has also been examined in biophysical models [48–50]. In order to use whole-brain simulation studies to better understand how the dynamic regime of the brain as a whole mediates the network-wide effects of stimulation, it will be necessary to inform these models with measurements of large-scale brain activity. Although here we implemented a more abstract model, working towards increasingly realistic and experimentally testable models is an exciting direction for forthcoming studies.
The authors would like to thank Richard Betzel and John Medaglia for providing assistance with the structural brain data. LP also thanks Panos Fotiadis and Keith Wiley for feedback on prior versions of this manuscript, as well as Jason Kim, Ann Sizemore Blevins, and Tanner Kaptanoglu for many helpful discussions.
The data used in this study is available from https://github.com/lia-papadopoulos/WilsonCowan_DynamicalState_Stim.
This work was primarily supported by National Science Foundation BCS-1631550, Army Research Office W911NF-18-1-0244, and National Institutes of Health R01-MH -116920. DSB, LP, and CWL would like to acknowledge additional support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the ISI Foundation, the Paul Allen Foundation, the Army Research Laboratory (W911NF-10-2-0022), the Army Research Office (Bassett-W911NF-14-1-0679, Grafton-W911NF-16-1-0474, DCIST-W911NF-17-2-0181), the Office of Naval Research, the National Institute of Mental Health (2-R01-DC-009209-11, R01-MH112847, R01-MH107235, R21-M MH-106799), the National Institute of Child Health and Human Development (1R01HD086888-01), National Institute of Neurological Disorders and Stroke (R01 NS099348), and the National Science Foundation (BCS-1441502, BCS-1430087, and NSF PHY-1554488). LP was also supported by a Graduate Research Fellowship from the National Science Foundation for part of this work. DB acknowledges support from the EU i-Innovative Training Network “i-CONN” (H2020 MSCA ITN 859937) and the French Agence Nationale pour la Recherche (“ERMUNDY”, ANR-18-CE37-0014-02). The content is solely the responsibility of the authors and does not necessarily represent the official views of any of the funding agencies. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
The data used in this study is available from https://github.com/lia-papadopoulos/WilsonCowan_DynamicalState_Stim.