Relations between large-scale brain connectivity and effects of regional stimulation depend on collective dynamical state

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.


SI Dynamics of an isolated Wilson-Cowan unit
In this section we briefly describe the behavior of a single Wilson-Cowan (WC) unit, which forms the fundamental dynamical component of the whole-brain computational model. The excitatory and inhibitory pools of an isolated WC unit evolve according to Eqs. ?? and ??, with the coupling term set to C = 0. Here we take all other model parameters to be those displayed in Table ??, with the exception that we consider noiseless simulations for the purpose of demonstration.
A typical bifurcation parameter for the WC model is the drive P E to the excitatory population; when other parameters are appropriately tuned, varying P E can induce oscillatory activity. The top row of Fig. A shows phase plane representations of a single WC unit for three different levels of the input P E , and the bottom row of each panel shows the time-evolution of the excitatory activity E(t) for the given parameter value. In the phase planes, the blue lines correspond to the excitatory variable nullcline (dE/dt = 0), the red lines correspond to the inhibitory variable nullcline (dI/dt = 0), and the purple lines show an example trajectory that begins at the point denoted by the star. For a low input level of P E = 0.6 (panel A), the system has a single stable fixed-point corresponding to a low activity steady-state. For an intermediate drive of P E = 1.25 (panel B), the system exhibits a stable limit cycle and the firing-rate activity oscillates in time. For a high input of P E = 3 (panel C ), the system again exhibits a single stable fixed-point, but corresponding to a high activity steady-state.
To summarize the effect of the drive P E on the behavior of an isolated WC unit, we first plot the time-average of the excitatory firing-rate E(t) as a function of the input P E (Fig. A, panel D). Note that E(t) increases with increasing drive. Second, we consider how the peak frequency f peak of the excitatory activity varies with the input level P E (Fig. A, panel E), where the peak frequency is that for which the Welch's power spectral density of the excitatory time-series is maximum (see ?? for details). For low inputs, f peak is approximately zero; the system resides in the low activity steady-state and there are no intrinsic oscillations. As the input is increased, though, oscillations emerge with frequencies in the gamma range (30 -70 Hz). Increasing the input P E in the oscillatory regime first raises the peak frequency from its initial value up to ∼65Hz. Such an increase in oscillatory frequency with increasing drive has also been found experimentally [1][2][3][4]. However, beyond a certain point, further increasing the excitatory drive causes the peak frequency to decline back to zero as the system approaches the high activity fixed point where oscillations again cease completely. In Fig. A, panel F, we show excitatory activity time-series for three different values of the input P E that place the system in the oscillatory regime. It is clear by eye that for these parameters, increasing the excitatory drive increases the amplitude and frequency of the oscillations.

SII Determining the onset of oscillatory activity in the whole-brain model
In order to systematically determine the boundary marking the onset of oscillatory activity as a function of the background drive P base E and the coupling C, we examine the network-average of the standard deviation of the excitatory activity across time, std[E i (t)] . This quantity measures the strength of fluctuations of the excitatory population activities around their mean values. Thus, by noting when std[E i (t)] jumps from a value near zero to a higher, positive value, we can qualitatively determine the transition from the state of low, non-oscillatory firing-rates to the onset of rhythmic dynamics in regional activity. Here, we are interested in finding the level of background excitation P * E (C) that is needed to induce oscillations at each brain area for a given interareal coupling C. To determine these "boundary" points P * E (C), we thus hold C fixed, and consider the difference in std[E i (t)] between consecutive values of P base E . We plot this difference ∆ std[E i (t)] as a function of P base E and C in Fig. B, where we indeed observe a clear boundary separating the low-activity and oscillatory regimes. In particular, we define the border point for each coupling C as the value P * E (C) for which the difference ∆ std[E i (t)] is maximized. We mark the boundary corresponding to the onset of oscillatory activity with red squares in Fig. B.

SIII Working point 2: Global coherence peak
In the main text, we examined in detail the effects of focal perturbations for two distinct working points -WP1 and WP3 -corresponding to low and high background drive states situated below and above peak global coherence, respectively. In this section, we also consider the effects of perturbations at WP2, for which P base E = 0.57 (the coupling is kept at C = 2.5). Here, the background input is intermediately valued, and the system resides at approximately the state of peak ρ global (see Fig. ??A). We thus observe blocks of relatively strong phase-locking in the baseline PLV matrix ( Fig. ??C, Row 2, Column 1) at this working point. However, note that the network-averaged PLV is only ≈ 0.57, which is still significantly less than the maximum possible value of ρ global = 1. ranges between approximately 5 and 15 Hz (Fig. C, panel D). Consequently, there is a clear separation between the distribution of units' peak frequencies in the baseline condition and when excited with additional input (Fig. C, panel E), indicating that there are again two frequency ranges of interest for further analysis. Due to interactions with the network, the excited area's spectra also develops sidebands to the left and right of its main peak, and a second bump at a frequency equal to the difference in the sideband and peak frequencies. This latter feature reflects the enhanced amplitude modulations that emerge in the time-series of unit i under stimulation (Fig. C, panel B, Right) and is indicative of quasiperiodic dynamics in the network. We also examine the power spectra of two downstream regions j and k located at increasing topological distances (hence receiving progressively weaker structural input) from the perturbed area i (Fig. C, panel C, Middle, Right). As for WP1, downstream region j, which receives strong input from i, develops a new spectral component at the main frequency of the excited region, and also at a lower frequency approximately equal to the difference of its baseline peak and the excited peak. On the other hand, unit k, which is more weakly connected to i, does not exhibit these same modulations. In Sec. SV of ??, we show that the relationship between the strength of the power modulation at the peak frequency of the stimulated site i and the topological distance from i to the dowstream unit holds across different choices of the excited area.
To see more generally how focal stimulation can modulate downstream spectra, we compute the average spectra psd j =i over all units j = i in the baseline state and in the state when unit i is selectively excited (Fig. C, panel F) and the average difference ∆psd j,δi j =i of the spectra of unit j = i between when unit i is stimulated and in the baseline condition (Fig. C, panel G). We again observe power modulations in both the baseline frequency band and in an excited frequency band centered around the peak frequency of the directly stimulated region. However, for the perturbed area studied in the example shown at WP1 (Fig. ??F-G) and again here at WP2 (Fig. C, panels F-G), the maximum relative power modulations are weaker in the latter case. In general, it is important to note that the strength of the observable changes to the power spectra induced by focal perturbations are dependent on the choice of the excited area and on the system's working point. In the next section, we proceed to examine changes in the PLV in the baseline and excited bands -∆ρ base δi and ∆ρ exc δi -induced by regional perturbations.

SIII.2 Structural and functional network connectivity continue to predict overall changes in excited and baseline band phase-locking at the peak-coherence working point
Akin to WP1, local stimulation at WP2 induces phase-locking changes that differ depending on which part of the network is perturbed and which frequency band is examined (Fig. D, panels A-D). To appreciate this fact, it is helpful to study examples that display the changes in the PLV in both the baseline band and the excited band when either region i or region j = i is excited (Fig. D, panels A,C). To characterize the network-wide impact of stimulating each region, we set non-significant PLV changes to zero (determined by comparing the observed changes against phase-randomized surrogate data; see Sec. SXII) and then compute network-averages of the absolute PLV modulations. Upon examination of the global responses |∆ρ base δi | and |∆ρ exc δi | (Figs. D, panels B,D), we again find significant variability across the choice of the perturbed area. Actually, in comparing the overall baseline band changes at WP1 (Fig. ??B) to the changes at WP2 (Fig. D, panel B), we find that the dispersion (as quantified by the coefficient of variation of the set of changes { |∆ρ base δi | } and the mean of the global responses are both larger at the second working point. We additionally note, however, that the mean response induced in the excited frequency band, |∆ρ exc | , decreases from WP1 to WP2. Hence, when the system operates around the state of peak ρ global , the coherence modulations in the baseline frequency band are larger -on average -and also more heterogeneously distributed relative to the responses at the working point below the state of peak ρ global , while the average of the global responses in the excited band decreases. We consider these points further in the final section of the ??, where we more generally examine the state-dependence of perturbation-induced changes in phase-locking. To conclude this section, we consider the relationships between changes in phase-locking induced by stimulating a given region and the structural or functional node strength of the stimulated area. A number of the associations that exist at WP2 (Fig. D, panels E,F) were also observed at WP1 (see Fig. ? (F) Average power spectra psd j =i over all units j = i at baseline (light gray) and when unit i is perturbed with additional input (dark gray). (G) Average difference ∆psd j,δ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 interest: (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, ∆ρ base δi (purple), and in the excited band ∆ρ exc δi (green).
the baseline frequency band continue to be best predicted by the strength of regions' initial coherence with the system as a whole. The main difference between WP1 and WP2 is that for the second working point, a positive correlation

SIV Average excited-band responses differ between stimulation of cortical and subcortical areas
In this section we show that the overall responses induced in the excited frequency band (at WP1 and WP2) areon average -larger for stimulation of subcortical vs. cortical brain areas. To do so, we first compute the average absolute change in excited-band phase-locking |∆ρ exc δi | for stimulation of each region i ∈ {1, ..., N }. Note that for these computations, only statistically significant PLV changes are retained as non-zero (see Sec. SXII). We then collect the sets of responses { |∆ρ exc δi | } c and { |∆ρ exc δi | } s corresponding to perturbation of cortical and subcortical areas, respectively. Fig. E shows the mean and spread of these two groups at WP1 (panel A) and at WP2 (panel B), from which we observe that the mean excited-band response is larger for stimulation of subcortical regions at both working points. To determine whether this effect is statistically significant, we perform a non-parametric permutation test of the null hypothesis that there is no difference in the mean excited-band response between cortical and subcortical areas. Using 1000 randomizations, we find that the mean excited-band responses are indeed significantly larger for stimulation of subcortical regions (p = 0.001) at both WP1 and WP2. Error bars indicate ± one standard deviation. The mean subcortical response is significantly greater than the mean cortical response (p = 0.001).

SV Topological distances in structural network predict power modulations at frequency of stimulated unit
In the main and supplementary text, we show examples of how focal stimulation of one region in the network affects the spectra of downstream areas (see Figs. ??C, ??C, and C, panel C). At WP1 ( Fig. ??C) and WP2 ( C, panel C), we observed a propagation effect in which a downstream area located at a short topological distance from the perturbed site developed an increase in power at the excited frequency of the directly stimulated area. In constrast, a unit located topologically further from the stimulated area exhibited much weaker power modulations. In this section, we show quantitatively that the topological distance from the perturbed site to downstream regions is a relatively good predictor of the downstream power modulation at the peak frequency of the stimulated unit. Hence, the effects observed in the examples shown in Figs. ??C and C, panel C generalize to other choices of the stimulated unit.
To begin, we more formally define the measure of topological distance that we employ. In general, a topological distance between any two nodes i and j is the length of the shortest path between those nodes, where a path is a non-intersecting sequence of edges that share a common node. Hence, to compute a topological distance, we first need to assign lengths to each edge in the network. Here, we define the length of an edge to be the inverse of the corresponding structural connectivity edge weight [5]. (For this analysis, we use the edge weights from the normalized adjacency matrix W ij ; see ??). With this definition, the topological distance from node i to node j will thus be shorter when node i links to node j via a path of stronger structural connections. Now, letting i denote the stimulated region and j = i denote a downstream region, we define D t i,j as the topological distance from node i to node j, and∆P j (f peak i,δi ) as the relative change in power of unit j's spectra at the excited peak frequency of the stimulated unit i. In particular, the relative change is computed between baseline and the condition in which i is driven with additional input.
With these quantities defined, we now study their relationship at WP1. Fig. F, panel A shows a plot of log |∆P j (f peak 1,δ1 )| vs. D t 1,j (i.e. we consider the case that the stimulated unit i = 1, which was the example in Fig. ??C.) This scatter plot exhibits a relatively clear negative trend, indicating that units located topologically nearer to the stimulated site (in the structural network) tend to exhibit stronger spectral modulations at the excited peak frequency of the stimulated region, whereas units that are further away (hence more weakly structurally connected) show little change. Furthermore, by examining the Spearman correlation between∆P j (f peak i,δi ) and D t i,j for all choices of the stimulated unit i, we see that this relationship holds more generally, regardless of which site is given the perturbation (Fig. F, panel B). Specifically, the correlation between the relative power modulation∆P j (f peak i,δi ) and the topological distance D t i,j is consistently negative and statistically significant across all choices of the stimulated area. Fig. F

SVI Average responses to perturbations in the baseline frequency band at WP1 vs. WP3
To investigate how the dynamical state of the brain network model influences the effects of focal stimulation, we consider the relationship between the global responses to stimulation at two different working points. In particular, we examine the average absolute change in phase-locking in the baseline frequency band, |∆ρ base δi | , at WP3 vs. at WP1 (Fig. G). Recall that WP1 corresponds to a low background drive working point preceding peak baseline coherence, whereas WP3 corresponds to a high background drive working point following peak baseline coherence. It is clear upon visual inspection of Fig. G that there is no consistent relationship between these two quantities. Hence, regions that induce a large response at the system's spontaneous frequencies at WP1 are not necessarily those that induce a large response at WP3.

Fig G. Network-averaged absolute change in baseline band phase-coherence |∆ρ base
δi | plotted at WP3 vs. WP1. Each point corresponds to a different choice of the stimulated region. Also note that PLV changes that were not statistically significant (see Sec. SXII) were set to zero before computing these network-wide averages.

SVII Negative PLV changes drive correlation with structural connectivity at WP3
In the main text we observed a strong positive correlation between a node's structural strength s struc i and the absolute coherence modulation |∆ρ base δi | it induces upon perturbation at WP3 (see Fig. ??). Here, Fig. H demonstrates that this association is largely driven by a strong relationship between s struc i and the network-averaged absolute decreases | ↓ ∆ρ base δi | in baseline band coherence induced by the focal stimulation. In particular, s struc i and | ↓ ∆ρ base δi | have a Spearman correlation of r s = 0.72, which is almost as strong as the correlation between structural strength s struc i and the absolute change | ↓ ∆ρ base δi | (r s = 0.82).

SVIII Verification of relationships between phase-locking modulations and structural or functional connectivity at alternate working points in the low, medium, and high background drive regimes
At WP1, WP2, and WP3, we considered the associations between the average changes in phase-locking induced by regional perturbations (within both the baseline and excited frequency bands) and structural or functional network properties of the stimulated region. We found that, depending on the baseline state of the system, different In this section, we verify that qualitatively similar relationships hold for other working points in the immediate vicinity of those studied in the main text. Note that for each alternative working point, we consider the same excitation strength used originally (i.e., ∆P E,i = 0.1). Moreover, throughout this section, we set non-significant PLV changes to zero (determined by comparing the observed changes against phase-randomized surrogate data; see Sec. SXII) prior to computing network-averages of the absolute PLV modulations.
We begin by analyzing an alternative working point near WP1, which we term WP1 alt . WP1 was located at P base E = 0.553 and C = 2.5; for WP1 alt , we consider parameters P base E = 0.555 and C = 2.5. Note that because peak global baseline coherence ρ global is reached rapidly as a function of P base E (see Fig. ??A), in order to consider a second working point located prior to ρ global but still near WP1, we can only shift P base E slightly from its value at WP1. We find the same set of relationships between phase-locking modulations and structural or functional strength at WP1 alt as we did at WP1 (see We next analyze an alternative working point near WP2, WP2 alt . Recall that WP2 was located at P base E = 0.57 and C = 2.5; for WP2 alt , we consider P base E = 0.572 and C = 2.5. In order to examine a second working point in close vicinity of the peak in global baseline coherence -which was the condition used to determine parameters for WP2 -we again must consider only a small change in P base E away from its value at WP2. As before, this is because the dynamical state of the system changes quickly as a function of P base E in this regime (see Fig. ??A). Using the specified parameter choices, we find consistent relationships at WP2 and WP2 alt in terms of how perturbation-induced phaselocking modulations are related to structural and functional node strength. First, the average absolute PLV changes that arise in the baseline frequency band |∆ρ base δi | are significantly correlated with the structural strength of the perturbed region s struc i (Fig. I, panel B,Top), but remain most strongly related to functional strength s func i (Fig. I,  panel E,Top). Second, the average absolute phase-locking modulations in the excited frequency band |∆ρ exc δi | are strongly associated with the structural strength s struc i of the perturbed region (Fig. I, panel B,Bottom), and are not significantly correlated with functional strength (Fig. I, panel E,Bottom).
Lastly, we analyze an alternative working point near WP3, WP3 alt . WP3 was located at P base E = 0.7 and C = 2.5; for WP3 alt , we consider P base E = 0.68 and C = 2.5. We once more find that the relationships between phase-locking modulations induced by regional stimulation and structural or functional node strength are consistent across WP3 and WP3 alt (see Fig. I, panel C,F). In particular, there is a strong positive correlation between the average absolute phase-locking modulations induced in the baseline frequency band |∆ρ base δi | and the structural strength of the perturbed region s struc i (Fig. I, panel C  weaker (Fig. I, panel F,Top). Finally, note that there is no excited frequency band at WP3 or WP3 alt .

SIX General dependence of associations between phase-locking modulations and structural or functional connectivity as a function of background drive
In Fig. 9A of the main text, we showed the difference ∆r s in the strength of the correlation between the average absolute baseline band PLV changes |∆ρ base δi | and structural (s struc  (Fig. J, panel A). Although the precise values of the correlations can vary in a somewhat complex manner as a function of the baseline input, we reiterate the key point that functional strength is more strongly related to baseline band coherence modulations in the low-and medium-drive regimes, whereas structural strength dominates in the high-drive regime.
We also examine the relationships between structural or functional node strength and the network-averaged absolute change |∆ρ exc δi | in excited-band PLV induced by focal perturbations. Because stimulation fails to induce an excited band when the background drive is too high, we compute Spearman correlations between these quantities only for values of the drive where at least half of the regions yield an excited band upon perturbation. Furthermore, at each working point, correlations are computed only between regions that induce an excited band. We find that areas' structural node strength robustly predicts the global response in the excited band (see Fig. J, panel B). In particular, the Spearman correlation coefficient r s between s struc i and |∆ρ exc δi | is greater than 0.89 across all baseline inputs for which at least half of the units in the network induce excited frequency bands. Furthermore, the strength of the correlation with structural connectivity is consistently much higher than the strength of the correlation using functional connectivity (compare yellow and blue curves in Fig. J, panel B). These results are consistent with the conclusions drawn in the main text positing that network structure mediates the excited-band effects.

SX Effects of perturbation strength
In the main text we studied a single perturbation strength of ∆P E,i = 0.1. In this section, we assess the dependence of various results on the level of additional excitatory input ∆P E,i received by the stimulated unit (Fig. K). In particular, for both WP1 and WP3, we vary ∆P E,i between 0.01 and 0.15 in steps of 0.02.
We first analyze how the perturbation strength affects the shift in the peak frequency of the stimulated area. As a summary measure, we consider the change in peak frequency averaged across all choices of the stimulated region, ∆f peak i,δi . As expected, this quantity increases with increasing perturbation strength for both WP1 (Fig. K, panel A) and WP3 (Fig. K, panel B). We next study the stimulation-induced changes in phase-locking in the baseline frequency band as a function of the stimulation strength. In particular, we examine the global response (grand average) |∆ρ base | , where the mean change is computed first over all pairs of brain areas for a given stimulation site, and then across all choices of the perturbed region. For both working points, this measure also increases monotonically as a function of ∆P E,i (Fig. K, panels C,D). Hence, as the strength of the stimulation increases, so does the overall amount of functional reconfiguration at the system's baseline frequencies. For WP1, we find that the grand average PLV change |∆ρ exc | in the excited band also grows as a function of the perturbation strength (Fig. K, panel E). Note, however, that at WP3 (Fig. K, panel F), no units generate an excited frequency band for ∆P E,i < 0.15, and for ∆P E,i = 0.15, only about 10% of units do so.
In the main text, we found that at WP1, the absolute change in baseline band phase-coherence induced by stimulating region i, |∆ρ base δi | , was strongly correlated with the functional strength of region i, s func i ( Fig. ??F,  Left). In contrast, there was not a strong association between |∆ρ base δi | and the structural strength s struc i at WP1 (Fig. ??E, Left). Here, we observe that the nature of these two relationships remains qualitatively the same across the considered range of stimulation strengths ∆P E,i (Fig. K, panel G). A second result from the main text was a strong positive correlation between the PLV change induced in the excited frequency band |∆ρ exc δi | , at WP1, and the structural strength s struc i of the stimulated unit ( Fig. ??E, Right). The present analysis reveals that the relationship between |∆ρ exc δi | and s struc i holds across a range of perturbation strengths ∆P E,i > 0.05 (Fig. K, panel I), for which at least half of the units in the network yield an excited frequency band. Note that for ∆P E,i < 0.05, local excitations do not induce an excited frequency band at all (so the correlations are undefined) and for ∆P E,i = 0.05, fewer than half of the units generate excited frequency bands (so we do not consider correlations with structural or functional strength). Finally, for the case of ∆P E,i = 0.1 studied in the main text, there was a strong positive correlation between |∆ρ base δi | and s struc i at WP3 (Fig. ??C, Left)  ( Fig. ??C, Left). These trends also hold across the range of stimulation strengths examined in this section (Fig. K,  panel H). Note that we do not consider correlations between structural and functional node strength and excited band PLV changes, because at WP3, less than half of the network generates an excited frequency band (even for the strongest stimulation strength).  In conclusion, we note that while an in-depth examination of the effects of the stimulation strength is beyond the scope of the present study, it is important that the main relationships between interareal phase-locking modulations and network properties hold over a range of values for this parameter.

SXI Results for an alternative value of the global coupling
In the main text, we examined the effects of focal excitatory stimulation for a global coupling of C = 2.5. Here, we analyze an alternative (but relatively nearby) coupling value of C = 2.0, and show that qualitatively similar results are found. As for C = 2.5, we consider three different working points by varying the level of background drive P base E , while holding the coupling fixed. Specifically, we consider P base E = 0.60 (WP1), P base E = 0.615 (WP2), and P base E = 0.745 (WP3), which place the system below, at, or above the state of peak global coherence (see Fig. ??A of the main text), respectively. As before, these working points represent three distinct dynamical states of the system.
We begin by considering the effects of regional stimulation on the power spectra of the perturbed area. To summarize this, we examine the average shift in the peak frequency of the stimulated region, ∆f peak i,δi , for each of the three working points (Fig. L). For all three states, additional excitation has the effect of increasing the peak frequency of the stimulated region. However, for WP1 and WP2, the peak frequency shifts by a noticeably larger amount ( ∆f peak i,δi = 9.5Hz for WP1 and ∆f peak i,δi = 10.8Hz for WP2) relative to the more modest effect at WP3 ( ∆f peak i,δi = 2.2Hz). These general trends are consistent with the results in the main text, and again demonstrate that individual areas are most responsive to additional excitation in states of lower background drive (WP1 and WP2). In contrast, given the same excitation strength, regional dynamics are relatively imperturbable when the system operates in the high background drive state (WP3).
We next examine how regional stimulation affects interareal phase-locking at each of the three working points. For WP1 and WP2 we analyze separate "baseline" and "excited" frequency bands, since the peak frequency of the stimulated area becomes separated from the peak frequencies of the system at baseline. For WP3, we consider a single "baseline" band, as the peak frequency of the excited area shifts only slightly and can overlap with the main frequencies at baseline. For the present analysis, we use the same protocol described in the main text to define baseline and excited frequency bands. In general, we refer the reader to the ?? section of the main document for further details and discussion regarding the results presented below.
We first show -for WP1 -examples of the phase-locking modulations within the baseline and excited frequency bands for two different choices of the stimulated area (Figs. M, panels A,B). As in the main text (see Figs. ??A,C), we see that the network response to a local perturbation differs between the two frequency bands, and for different choices of the stimulated region. We next study the associations between the network-wide average of the phaselocking modulations induced by regional stimulation and structural or functional strength (Figs. M, panels C,D). To characterize the network-wide impact of stimulating each region, we set non-significant PLV changes to zero (determined by comparing the observed changes against phase-randomized surrogate data; see Sec. SXII) and then compute network-averages of the absolute PLV modulations. In comparing the results presented here for a coupling of C = 2 to those in Fig. ??E We next conduct the same analyses regarding changes to interareal phase-locking, but for WP2. Here, Figs. N, panels A,B show examples of the phase-locking modulations within the baseline and excited frequency bands for two different choices of the stimulated area, and Figs. N, panels C,D depict relationships between the global PLV changes induced by regional stimulation and structural or functional node strength. We again find qualitatively similar behavior between the results shown here and those depicted in the main text. Note that for both values of the coupling (C = 2 here and C = 2.5), the main difference between WP1 and WP2 is that structural strength s struc i also exhibits a positive correlation with the average absolute change in phase-coherence for the baseline frequency band |∆ρ base δi | (Fig. N,   spontaneous frequencies continue to be most strongly associated with the stimulated region's baseline functional strength (Fig. N, panel D, Left).
For completeness, we lastly consider phase-locking changes induced by regional perturbations at WP3. Fig. O, panel A shows the change in PLV between each pair of regions (for the single, baseline frequency band) for two different choices of the stimulated area. As found in the main paper (e.g., Fig. ??A), the response, in general, differs across the choice of the excited region. The relationships between the phase-coherence modulations and structural or functional network strength found here for C = 2 (Fig. O, panel B) are also consistent with the analysis performed in the primary text for C = 2.5 ( Fig. ??C)

SXII Determining statistical significance of PLV changes using phaserandomized surrogates
Throughout the text, we compute changes in interareal phase-locking between baseline conditions and the case of focal stimulation. In particular, we subtract the PLV matrix computed from baseline simulations from the PLV matrix computed from simulations corresponding to regional stimulation. Of note is that this method does not account for the possibility that observed changes in phase-locking arise only from differences in the autocorrelation structures of units' time-series in the baseline vs. stimulation conditions, which could affect results in the case of finite sample sizes. To determine when observed changes in phase-locking are different than the changes expected from a potential change in the autocorrelations of units' dynamics alone, we can conduct statistical significance testing using autocorrelation-preserving surrogate data. We describe this methodology in what follows. To begin, we will use the terminology ∆ρ base jk,δi and ∆ρ exc jk,δi to denote the change in baseline band and excited band PLV, respectively, between regions j and k induced by perturbation of unit i. Note that i, j, k are all ∈ {1, ..., N }. To determine whether these changes are statistically significant, we need to build null distributions {∆ρ base jk,δi } and {∆ρ exc jk,δi } against which the observed changes are compared to. Furthermore, we want these null distributions to represent the expected changes in phase-locking for time-series that have the same autocorrelations as those from the baseline-and stimulation-condition simulations, but constructed such that the dependencies between surrogate time-series from different units are destroyed. One well-known autocorrelation-preserving surrogate method is that of phase-randomization [6]. This surrogate method maintains the power spectrum of an original signal s(t), but randomizes its Fourier phases. In short, this is accomplished by computing the discrete Fourier transform of the original signal, adding phases drawn independently and at random from the interval [0, 2π] to the phase for each frequency (while preserving the fact that the signal must be real), and then transforming back to the time-domain to obtain a "phase-randomized" surrogates(t) [7]. If a different phase-randomization is applied to the surrogates corresponding to different units in the network, then each unit's autocorrelation will be preserved, but interdependencies between different units will be destroyed. In this way, we can test the null hypothesis that the observed PLV changes ∆ρ base jk,δi or ∆ρ exc jk,δi are due only to finite sample size bias and the differing autocorrelation structures of units' time-series between the baseline and stimulation simulations. If we can reject this null hypothesis, then we conclude that the observed modulations in phase-locking reflect actual changes in interareal coherence, beyond what is expected due to differing autocorrelations.
To generate the null distributions {∆ρ base jk,δi } and {∆ρ exc jk,δi }, we follow the same steps used to generate the true changes ∆ρ base jk,δi or ∆ρ exc jk,δi , with one exception: instead of using units' actual excitatory time-series in the PLV calculations, we use corresponding phase-randomized surrogates. In particular, to generate one instance of {∆ρ base jk,δi } or {∆ρ exc jk,δi }, we generate a phase-randomized surrogate from each unit's activity in every trial, using different randomizations for all units and trials. Then, using the surrogate time-series, we compute one instance of∆ρ base jk,δi and∆ρ base jk,δi following the same steps used to compute the original PLV changes in the baseline and excited frequency bands (see the ?? section of the main text). This process is then repeated 50 times -each time using a different set of surrogate realizations -to generate null distributions {∆ρ base jk,δi } and {∆ρ exc jk,δi }. Once we have generated null distributions {∆ρ base jk,δi } and {∆ρ exc jk,δi } for a given choice of the excited unit i, the next step is to compare them to the observed differences ∆ρ base jk,δi and ∆ρ exc jk,δi . To determine if an observed ∆ρ jk is statistically different from a null distribution {∆ρ jk }, we first check if ∆ρ jk is positive or negative. If ∆ρ jk > 0, then we compute the fraction of surrogates p + for which∆ρ jk > ∆ρ jk . If p + < 0.05, then we conclude that ∆ρ jk is statistically greater than expected under the null distribution, and the increase in phase-locking is significant. If ∆ρ jk < 0, then we compute the fraction of surrogates p − for which∆ρ jk < ∆ρ jk . If p − < 0.05, then we conclude that ∆ρ jk is statistically less than expected under the null distribution, and the decrease in phase-locking is significant. If the change in phase-locking ∆ρ jk is not found to be statistically different from the null distribution, then we set its value to zero prior to computing network-averaged changes in the PLV. Using this procedure, we can check the statistical significance of all PLV changes in the baseline and excited frequency bands, and for each choice of the stimulated node. The statistical significance testing described above is performed for WP1, WP2, and WP3 in the analyses where we consider network-averaged absolute modulations in phase-locking and their relationships to the structural or functional strength of the perturbed node (e.g., Fig. ??B,D-F, Fig. D, panels B-F, and Fig. ??B,C). In particular, at these working points, PLV changes that are not significant (relative to the phase-randomized null model) are set to zero before computing the network-wide averages of PLV modulations in the baseline |∆ρ base δi | and excited |∆ρ exc δi | frequency bands. We do not repeat these comparisons against surrogate data for the parameter sweeps over the background drive and stimulation strength. However, findings from WP1, WP2, and WP3 (as well as WP1 alt , WP2 alt , and WP3 alt ) indicate that results are largely unaffected by whether or not one sets non-significant PLV changes to zero prior to computing the global response. Throughout the text, we indicate when the statistical significance testing was performed on the PLV changes.

SXIII Details on the Hilbert Transform
A common way to extract an instantaneous phase variable from a real-valued, oscillatory signal is with the Hilbert transform. To begin, one writes the analytic (complex-valued) signal representation X A (t) of the real-valued timeseries X(t) as where X H (t) is the Hilbert transform of X(t), A(t) is the instantaneous amplitude of X(t), and θ(t) is the instantaneous phase of X(t). Once one has computed X H (t) and thus X A (t), it is apparent from Eq. 1 that the phase θ(t) can be computed as The Hilbert transform of a signal X(t) is defined as where the integral is evaluated as a Cauchy principal value. From Eq. 3, one observes that the Hilbert transform is the convolution of X(t) and 1/πt: X H (t) = X(t) * 1/πt, so the Fourier transform (FT) of X H (t),X H (f ), is just the product of the FTs of X(t) and 1/πt. For frequencies f > 0, we thus have thatX H (f ) = −iX(f ), from which it becomes clear that the Hilbert transform just induces a phase shift of π/2 to each frequency component in the signal.
In this study, we computed Hilbert transforms of the simulated neural activity using the 'hilbert' function in MATLAB. As described in the ?? section, the Hilbert Transform was applied after first filtering the raw time-series in a specified frequency band, in order to ensure that the corresponding phase variable is well-defined [8]. Table A. Brain region ID numbers with their corresponding hemisphere (L = left hemisphere, R = right hemisphere) and anatomical labels. The ID numbers match the node numbering used in the adjacency matrix (see Fig. ??C), which also applies to all other region-by-region figures in the main text and in the supplementary text.

Region ID
Hemisphere-Label Region ID Hemisphere-Label

SXV Citation diversity statement
Recent work in neuroscience and other fields has identified a bias in citation practices such that papers from women and other minorities are under-cited relative to the number of such papers in the field [9][10][11][12][13][14]. Here we sought to proactively consider choosing references that reflect the diversity of the field in thought, form of contribution, gender, race, geographic location, and other factors. We used automatic classification of gender based on the first names of the first and last authors [9], with possible combinations including male/male, male/female, female/male, female/female. Excluding self-citations to the senior authors of our current paper, the main text references contain 61% male/male, 8% male/female, 24% female/male, 7% female/female. We look forward to future work that could help us to better understand how to support equitable practices in science.