Alpha oscillation, criticality, and responsiveness in complex brain networks

The brain’s responsiveness to external stimuli is significantly dependent on its state. Previous studies with transcranial magnetic stimulation (TMS) have demonstrated that stimulation effects during sleep (NREM), anesthesia, and vegetative states are simple and constrained to local areas, whereas they are more complex and diffuse in the conscious state (Casali et al., 2013; Ferrarelli et al., 2010; Sarasso et al., 2015, 2014). In the conscious state, perceptual responses to sensory stimuli are associated with ongoing brain activity at the moment when the stimulus is applied. In visual perception tasks, for example, target detectability is dependent on the power and phase of ongoing alpha oscillations (∼10 Hz) (Benedetto, Lozano-Soldevilla, & VanRullen, 2018; Busch, Dubois, & VanRullen, 2009; Dugue, Marque, & VanRullen, 2011; Ergenoglu et al., 2004; Mathewson, Gratton, Fabiani, Beck, & Ro, 2009; Mathewson, Lleras, Beck, Fabiani, Ro, & Gratton, 2011; Milton & Pleydell-Pearce, 2016; Nunn & Osselton, 1974; Romei et al., 2008; van Dijk, Schoffelen, Oostenveld, & Jensen, 2008). Recent repetitive TMS (rTMS) studies have suggested that the instantaneous phase of the sensorimotor μ-rhythm in the alpha band (8–12 Hz) is an important factor in determining the stimulation effect (Schaworonkow, Triesch, Ziemann, & Zrenner, 2018; Zrenner, Desideri, Belardinelli, & Ziemann,2017).

However, despite many relevant experimental studies, researchers have not elucidated how external stimulation interacts with specific amplitudes and phases of alpha oscillation in ways that result in a large or small perceptual response. In this study, we hypothesized that perceptual binding associates with the brain’s network response, and the network response largely depends on the ongoing coupled oscillation properties at stimulus onset. Recent studies have revealed that alpha oscillations transfer information through traveling waves in the cortex, and the global functional connectivity of alpha oscillations reflect conscious and unconscious states (Deco et al., 2018; Deco, Woolrich, Stevner, van Hartevelt, & Kringelbach, 2017a; Jobst et al., 2017; H. Kim, Moon, Mashour, & Lee, 2018; M. Kim, Kim, Mashour, & Lee, 2017; Moon et al., 2017; Moon, Lee, Blain-Moraes, & Mashour, 2015; Zhang, Watrous, Patel, & Jacobs, 2018). Therefore, we assumed that the propagation and complexity of perturbed responses in the brain network may be determined by the alpha oscillations at stimulus onset. In this modeling study, we examined the specific amplitude, phase, and synchronization level of alpha oscillations that evoke a large or small network response.

To identify such a condition, we constructed a large-scale human brain network model with the Stuart-Landau oscillators of the alpha band frequencies (∼10 Hz) that are coupled with each other in an anatomically informed human brain connection structure. We tested three dynamic brain states, modulating criticality (below, near, and above critical state). The three brain states produced various transient alpha oscillations in terms of global synchronization, amplitude, and phase. The critical state was defined with the largest pair correlation function (PCF), measuring the variance of global synchronization fluctuation across time (Yoon, Sorbaro Sindaci, Goltsev, & Mendes, 2015), which is equivalent to susceptibility in statistical physics models. Theoretically, near the critical point, a large dynamic fluctuation would be expected. (Chialvo, 2010; Cocchi, Gollo, Zalesky, & Breakspear, 2017). Empirically, the brain during conscious wakefulness presents a large dynamic fluctuation, especially in synchronization (Haimovici, Tagliazucchi, Balenzuela, & Chialvo, 2013; Kitzbichler, Smith, Christensen, & Bullmore, 2009); however, various brain perturbations (e.g., sleep, anesthesia, traumatic injuries) mitigate this fluctuation, which implies that the brain deviates from criticality (Hutt, Lefebvre, Hight, & Sleigh, 2018).

The various alpha oscillations of the three brain states were decomposed into basic oscillation properties: instantaneous global synchronization, amplitude, and phase. We applied a pulsatile stimulus to the brain network in order to investigate the relationship between these alpha oscillation properties at stimulus onset and brain network responsiveness. We found specific phases of alpha oscillation that barely respond to the pulsatile stimulus, which suggests the presence of periodic temporal windows in the alpha frequency cycle for a large/small responsiveness, and which is consistent with the empirical studies (Callaway & Alexander, 1960; Ergenoglu et al., 2004; Nunn & Osselton, 1974; VanRullen, 2016, 2018; VanRullen & Koch, 2003). A schematic diagram of the study design is presented in Figure 1.

Various alpha oscillations in the human brain network were simulated in three brain states, C_b, C_p, and C_a (below, near, and above a critical state, see Materials and Methods for details). To simplify our model, we applied a simple form of stimulus (a pulsatile stimulus) to all nodes for 50 ms. In total, 600 stimulation trials were performed for each state. First, we examined the characteristics of alpha oscillations near and far from a critical state. Figure 2A provides examples of distinct global synchronization fluctuations for the three states. Synchronization fluctuation was measured with the PCF. Figure 2B presents the average synchronization < r(t > and PCF of the alpha oscillations at C_b, C_p, and C_a. Twenty different frequency configurations of the brain network were tested. The average levels of synchronization at C_p (orange triangle) are widely distributed between the small and large average levels of synchronization at C_b (blue circle) and C_a (yellow square). The alpha oscillations at C_p (orange triangle) show a larger PCF (mean ±SD = 1.39 ± 0.68) than those of C_b (blue circle) and C_a (yellow square), which suggests a large variance of brain network dynamics at C_p. Figure 2C and2D demonstrate the instantaneous synchronization r_s and instantaneous amplitude $Zsj$ of the alpha oscillations at stimulus onset (s = 1, 2,…, 600, and j = 1, 2, …, 82). At C_p (orange), both the instantaneous synchronization levels r_sand amplitudes $Zsj$ at are variable between the two small and large r_sand $Zsj$ of C_b (blue) and C_a (yellow), respectively (in Figure 2C, median ± SD, 0.11 ± 0.04 for C_b, 0.32 ± 0.14 for C_p, and 0.78 ± 0.05 for C_a; in Figure 2D, 0.21 ± 0.04 for C_b, 0.48 ± 0.57 for C_p, and 1.75 ± 1.65 for C_a). In the next sections, we will show how these state-specific alpha oscillations respond to the same stimulus and identify specific conditions of alpha oscillations that produce a large or small response in the brain network.

Alpha oscillations at diverse states of criticality show distinct responses to the stimulus. Figure 3 shows examples of the instantaneous global synchronization level r(t) at C_b, C_p, and C_a. Each triangle indicates a discrete application of the stimulus. One pulsatile stimulus was applied at random per simulation; a total of 600 such simulations were performed for each state. To calculate the responsiveness of each stimulation trial, we first defined a perturbation response (PR) of phase synchronization (PR_j(t)) for node j (j = 1, 2, …, 82) at time t, setting PR_j(t) as 1, if the synchronization value of node j was significantly increased from the prestimulus period (1 s); otherwise PR_j(t) was set as 0 (see Materials and Methods for details). Figure 3B presents examples of PR. The black dot represents PR_j(t) equals 1 and the red and blue shaded areas indicate the pulsatile stimulus of 50 ms and the response period, respectively. The average PR_j(t) was calculated over 600 stimulation trials. Figure 3C shows the temporal changes of the median (thick line) and of the 25% and 75% quantiles (thin lines) of the average PR for the three states. The results show that the PR at C_p persists longer than the PRs at C_b and C_a (sign test, p < 0.05 after 190 ms, Figure 2C). We used the PRs during the first 500 ms to quantify responsiveness for further analysis. In addition to the PR of phase synchronization, we also calculated the PR of the amplitudes of the alpha oscillations. These results are presented in the Supporting Information.

In this study, we defined two responsiveness indexes to quantify the PRs. First, we defined the responsivity R_j of a node j, counting the total amount of 1s in PR_j (significantly changed phase synchronization) during the first 500 ms after the stimulus (see Materials and Methods for details), which measures the magnitude of PRs. We also defined the Lempel-Ziv complexity (LZc) (Lempel & Ziv, 1976) of the PR_j for a node j, which measures the perturbational complexity (see Materials and Methods for details). The LZc measures the amount of information contained in the spatiotemporal pattern of PR, which is similar to the perturbational complexity index (PCI) that has been successfully applied to quantify the level of consciousness in sleep, anesthesia, and pathologic unconsciousness (Casali et al., 2013; Sarasso et al., 2015).

First, we investigated the relationship between instantaneous synchronization at stimulus onset (r_s, s = 1, 2, …, 600) and the responsiveness indexes. The average responsivity < R > was calculated by taking the average of R_j over all nodes. Figure 4 presents the Spearman correlation coefficients of the three states (C_b, C_p, and C_a) between <R> and the instantaneous synchronization r_s for 600 stimulation trials. Only at C_p does < R > show a significant negative correlation with r_s (Spearman correlation ρ_R = −0.48, p < 0.001 for C_p), suggesting a more desynchronized alpha oscillation induces a larger magnitude of PR at C_p.

Next, we calculated the LZc at each time point and averaged it over 500 ms, defining < LZc > (see Materials and Methods for details). Figure 4B presents the Spearman correlation coefficients between r_s and < LZc > for the three brain states. Only at C_p is < LZc > significantly correlated with r_s (Spearman correlation ρ_LZc = −0.52,p < 0.001 for C_p). At C_p a more desynchronized alpha oscillation induces a larger complexity of PR, and the complexity itself is the largest among the three states (Supporting Information Figure S1, mean ± SD = 0.57 ± 0.05, 0.67 ± 0.17, and 0.39 ± 0.18 for C_b, C_p, and C_a).

We found that near a critical state, stimulating a globally more desynchronized alpha oscillation (small r_s) gave rise to greater brain responsiveness in terms of both magnitude and complexity of PR, whereas far from a critical state, responsiveness was not correlated with the global synchronization of alpha oscillations.

In the previous section, we found a significant correlation between the synchronization of alpha oscillations and responsiveness. In this section, we focus on the critical state C_p and explore whether the amplitude and phase of alpha oscillation also play a role in responsiveness. Interestingly, we found specific phases of alpha oscillation that barely respond to stimuli. We first classified all nodes into two groups according to their oscillation properties—r_s (instantaneous global synchronization level, s = 1, 2, …, 600) and $|Zsj|$ (instantaneous node amplitude, j = 1, 2, …, 82)—at stimulus onset. These nodes were classified into high/low r_s (HS/LS) and high/low$|Zsj|$ (HA/LA) by comparing them to the average r_s and $|Zsj|$⁠, respectively. We further separated the phases of node $θsj$ in each group into 30 phase bins of 12° intervals. The graphical description of classification is presented in Figure 5A. Therefore, responses were classified by phase groups with intersections of ongoing oscillation properties (global synchronization level, node amplitude, and phase). Figure 5B and5C present responsiveness (R and LZc) in each phase bin for the low synchronization and low amplitude (LS and LA) group and high synchronization and high amplitude (HS and HA) group. The thick lines indicate the median value of responsiveness for each phase bin, and the shaded area covers the 25% to 75% quantiles of the values of each phase bin for each group. We found that the average responsiveness of LS and LA is larger than that of HS and HA (p < 0.001; Kruskal-Wallis test). Notably, the HS and HA group showed significant phase dependence. For the purpose of statistical evaluation, we separated these phases into two phase regimes, 60° to 240° and 240° to 60°. The HS and HA group showed significant phase dependence in both responsiveness measures (R and LZc) (Wilcoxon rank sum test, Figure 5D and5E); in particular, stimuli applied to specific phases from 60° to 240° failed to evoke significant brain response (***p < 0.001). These results suggest the existence of a periodic temporal window that does not significantly respond to external stimuli. We also tested various other conditions: LS and HA, HS and LA, different stimulus strengths, the two other states, C_b and C_a, and different model parameter λ = −1 (see Supporting Information Figure S2 for LS and HA and HS and LA; Supporting Information Figures S3–S7 for the responses to different stimulus strengths; Supporting Information Figures S8 and S9 for the two other states, C_b and C_a; and Supporting Information Figures S10 for the different model parameter). An excessively small stimulus failed to evoke a response, whereas an excessively strong stimulus produced large responses that eliminated the difference between LS and LA and HS and HA. Phase dependence was also observed in C_b and C_a; however, in this study we focused on the response behavior at critical state C_p to interpret the empirical results during conscious wakefulness.

In the previous sections, we examined the various degrees of responsiveness of alpha oscillations to a global stimulus. For this section, we tested if stimulation conditions for large/small responsiveness to the global stimulus also held for a local stimulus. Here we applied a pulsatile stimulus to nodes within a radius of 50 mm centered on the left cuneus in the occipital region to roughly simulate stimulation of the visual area (13 regions; (left) cuneus, inferior parietal, isthmus cingulate, lateral occipital, lingual, pericalcarine, precuneus, superior parietal, (right) cuneus, isthmus cingulate, lingual, pericalcarine, precuneus). We considered three stimulation conditions: effective, noneffective, and random stimulation. The effective stimulation condition was defined as the phases from 240° to 60° of LS and LA (red bar in LS and LA of Figure 5D and 5E) that induce a large degree of responsiveness. The noneffective stimulation conditions were defined as the phases from 60° to 240° of HS and HA (green bar in HS and HA of Figure 5D and 5E) that induce a small degree of responsiveness. For the random stimulation condition, we randomly selected the stimulus onset, which may correspond to conventional stimulation methods. In comparisons of responsiveness (R and LZc), we found that the effective (less effective) stimulation condition induces larger (smaller) responsiveness (Figure 6A and 6B), while the responsiveness of the random stimulation condition fell in the middle. The R and LZc of the PRs were significantly different across the three conditions (Kruskal-Wallis test, p < 0.001; in addition, a multicomparison test was performed using the Tukey-Kramer method for R and LZc, respectively; **p < 0.005 and ***p < 0.001). The results demonstrated that the effective and less effective stimulation conditions we identified for global stimulation can be also applied to local stimulation.

Empirical studies have suggested that brain responsiveness is associated with ongoing brain activities at specific frequencies (e.g., alpha oscillation, ∼10 Hz) when a stimulus is applied. In this study, we hypothesized that brain responsiveness is determined by the interaction between a stimulus and coupled alpha oscillations in the brain network. We simulated various alpha oscillations for three brain states (C_b, C_p, and C_a) and investigated the conditions of alpha oscillation that facilitate large and small responses to a stimulus. Each brain state presented characteristic alpha oscillation features: below the critical state (C_b), low amplitude was associated with an incoherent state; near the critical point (C_p), high and low amplitudes were balanced with a large degree of synchronization fluctuation; above the critical point (C_a), high amplitude was associated with a highly synchronized state. The alpha oscillation response was the most complex and persisted longer at critical state; in this state, responses were larger and more complex when global pulsatile stimuli were applied to a globally desynchronized, low-amplitude state. Importantly, we found there are specific phases of high-amplitude, highly synchronized alpha oscillations that barely respond to stimuli.

Empirical and computational model studies have proposed that the brain in conscious wakefulness resides near a critical state, whereas brain states under significant alterations such as anesthesia and traumatic injury are distant from a critical state (Fekete et al., 2018; Haimovici et al., 2013; Hutt et al., 2018; Kitzbichler et al., 2009). In this model study, we observed that the responsivity (the magnitude of the significant response) at C_b and C_a was larger than that at C_p (Figure 3D). However, the persistence of responsivity is longer at C_p, which implies restricted propagation at the states far from criticality, C_b and C_a (Figure 3D). The brain network at C_p also exhibited the most complex response when we measured perturbational complexity (Supporting Information Figure S1). These characteristics are consistent with TMS outcomes in empirical studies. TMS-evoked activities during unconsciousness show large but short-lasting simple response patterns (Ferrarelli et al., 2010; Sarasso et al., 2014), whereas EEG responses induced by TMS during conscious wakefulness show more complex perturbational patterns compared with anesthesia and NREM sleep (Casali et al., 2013; Ferrarelli et al., 2010; Sarasso et al., 2015, 2014). The similarity between the results of our model study and empirical results from previous studies suggests that state-dependent responses of the human brain to TMS may be due to the distinct response behaviors of the characteristic states (synchronous/incoherent or high/low amplitude) of networked oscillators below, near, and above a critical state.

In our modeling study, we found that brain network responsiveness correlates with the level of instantaneous global synchronization at the critical state (Figure 4). Stimulation at lower levels of instantaneous global synchronization produced larger responses for both amplitude and phase synchronization, whereas higher levels of instantaneous global synchronization induced smaller responses (Supporting Information Figures S12 and S13). Therefore, considering the large synchronization fluctuation at a critical state, it might be important to stimulate the brain network at appropriate times in order to induce large network perturbation. Not considering the timing of stimulation may produce greater variability in outcomes. In addition, synchronization dependence provides novel insight into the role of synchronization fluctuation in brain information processing. In an epoch of lower levels of synchronization, the brain might not be able to integrate distributed information within the brain network, but may be highly susceptible to external stimuli. In contrast, in an epoch of higher levels of synchronization, the brain might easily integrate distributed information within the brain network, but may not be able to respond to external stimuli. Regarding the typical large synchronization fluctuation at a critical state, the response behaviors of a network at high and low levels of instantaneous synchronization seem to create temporal windows that can integrate internal and external information, respectively, while systematically separating internal and external information processing in the brain network.

In our modeling study, we found that responsiveness also associates with the amplitude and phase of the ongoing alpha oscillations. We decomposed the alpha oscillations into high and low instantaneous synchronization levels, amplitudes, and phases from 0° to 360°. The average responsiveness of lower instantaneous synchronization and low amplitude (LS and LA) was larger than that of higher instantaneous synchronization and high amplitude (HS and HA) (Figure 5). Importantly, considering the alpha oscillations of high synchronization and high amplitude, we found that there are specific phases (60° to 240°) that exhibit significantly smaller responsiveness in both responsiveness indexes. The results suggest periodic temporal windows (phases from 60° to 240° in each cycle of an alpha wave) in which the responsiveness of the brain network was largely inhibited.

The amplitude and phase of alpha waves (8–13 Hz) are associated with sensory stimuli processing, especially in regard to visual perception (Benedetto et al., 2018; Busch et al., 2009; Dugue et al., 2011; Ergenoglu et al., 2004; Mathewson et al., 2009, 2011; Milton & Pleydell-Pearce, 2016; Nunn & Osselton, 1974; Romei et al., 2008; van Dijk et al., 2008). Stimulation at lower amplitudes and at a phase near the trough of an alpha wave show increased target detectability (Brüers & VanRullen, 2018; Busch et al., 2009; Ergenoglu et al., 2004; Mathewson et al., 2009, 2011; Milton & Pleydell-Pearce, 2016; van Dijk et al., 2008). Phosphenes artificially induced by TMS are more detectable in lower prestimulus alpha power and specific alpha phase ranges (Dugue et al., 2011; Romei et al., 2008; Samaha, Gosseries, & Postle, 2017). Recent TMS studies have also shown that motor-evoked potential (MEP) amplitude is dependent on the prestimulus cortical μ-rhythm phase (Schaworonkow et al., 2018; Zrenner et al., 2017). The established dependence of the brain’s responsiveness on the amplitude and specific phase of alpha waves is consistent regardless of stimulus type. Our human brain network model showed that stimuli applied at specific phases of alpha oscillation resulted in larger responsivity and perturbational complexity at the critical state. Interestingly, phase dependence was more pronounced when the stimulus was given at higher amplitudes of alpha oscillation (Figure 5C and 5D), which is consistent with the visual target detectability experiments (Mathewson et al., 2009). This result emphasizes the importance of considering effective stimulation conditions, especially when the stimulus is applied to alpha oscillations with high amplitudes and high levels of synchronization. We also note that phase-dependent responses are diminished when stimulus strength is too small or too large (Supporting Information Figures S3–S7).

We also tested if the conditions identified under global stimulation would still hold for local stimulation. We applied the same stimulus to the occipital region across three different conditions: effective (low instantaneous synchronization, low amplitude, and phases of 240°− °60), less effective (high instantaneous synchronization, high amplitude, and phases of 60°− °240), and random (conventional stimulation method without considering ongoing brain activity). Relative to the random condition, the effective and less effective conditions produced significantly larger and smaller brain responsiveness. These results may explain why the outcomes of conventional brain stimulation studies show such large intrasubject variability (Huang et al., 2017), which degrades the reliability of brain stimulation studies. In addition, finding the specific phases of alpha oscillations for large or small responsiveness reilluminates the historical concept of the “neuronic shutter” described in the 1960s (Callaway & Alexander, 1960), which implies the existence of temporal windows that periodically prohibit sensory information processing during conscious states (Callaway & Alexander, 1960; Cecere, Rees, & Romei, 2015; Harris & Thiele, 2011; Nunn & Osselton, 1974; VanRullen, 2016, 2018; VanRullen & Koch, 2003). Further study is required to identify the relevance of temporal windows of sensory information processing to the response behavior of networked alpha oscillations to external stimuli.

A desynchronized EEG, as characterized by low amplitude, noise, and fast-frequency activity, is associated with arousal and increased awareness (Zerlaut & Destexhe, 2017). This desynchronized state emerges as the product of highly recurrent and stochastic interactions within, and between, the excitatory-inhibitory neuronal populations at both the cellular and network level (Zerlaut & Destexhe, 2017). At the cellular level, the stochastic interplay of excitatory and inhibitory conductance sets neurons into a high-conductance state with enhanced responsiveness. At the population level, neuronal ensembles also become highly responsive to afferent inputs (Zerlaut & Destexhe, 2017). At the large-scale brain network level, responsiveness is determined by the characteristics of ongoing networked oscillators. However, it remains to be answered how the specific oscillation properties induce such large or small responsiveness. One of the potential mechanisms is the phase response curve (PRC) that has been studied in physics. The PRC characterizes the way that a system with a collective periodic behavior responds to external stimuli. The response of a periodically driven oscillating system is measured by the phase shift from the original phase, and the phase shift (advancing or delaying the original phase) is an inherent characteristic of any oscillatory system. This method has been applied to many rhythmic biological systems such as circadian rhythms, cardiac rhythms, and spiking neurons to study how external stimuli perturb the original rhythms (Granada, Hennig, Ronacher, Kramer, & Herzel, 2009; Hannay, Booth, & Forger, 2015; Ko & Ermentrout, 2008). The previous analytic studies discovered influential factors causing phase shift and phase synchronization perturbation: a low (high) phase coherence induces a large (small) phase shift (Hannay et al., 2015). Dynamically, stimulation to the phases around a stable fixed point of the PRC increases phase coherence, whereas stimulation to the phases around an unstable fixed point decreases phase coherence. These properties also hold for networks with different coupling functions, network structure, connectivity, and even for the amplitude dynamics of Stuart-Landau oscillators (Levnajić & Pikovsky, 2010). If the alpha oscillations are globally coupled in the human brain network (Zhang et al., 2018), the response behavior may be governed by the same mechanism of PRC. Furthermore, future study needs to answer how responsiveness across multiple scales, including cellular, neural population, and large-scale brain networks, are related to each other.

In summary, this modeling study demonstrated for the first time a relationship between brain network responsiveness and the ongoing alpha oscillations at different states (below, near, and above a critical state). We found properties of alpha oscillation that induce a large or small responsiveness with the same stimulus. The results explain why brain responsiveness is so variable across brain states in consciousness and unconsciousness and across time windows even during conscious wakefulness. They will also provide a theoretical foundation for developing an effective brain stimulation method that considers state-specific and transient brain activity.

This modeling study has several limitations. First, despite the potential contribution of other frequency bands to periodic brain responsiveness (Fiebelkorn, Pinsk, & Kastner, 2018; Helfrich et al., 2018), our model simulated only alpha oscillations. Further studies are required to confirm periodic responsiveness for other frequency bands. Our modeling results suggested that the periodic prohibition of sensory information processing in the human brain is due to the general response behavior of networked oscillations. However, to justify this argument, the association of the global synchronization of alpha oscillations (HS and HA and LS and LA) with responsiveness, which our modeling study discovered, needs to be tested empirically. Second, we only studied a pulsatile stimulus, which is similar to the stimulus type applied in previous experiments (Brüers & VanRullen, 2018; Busch et al., 2009; Ergenoglu et al., 2004; Mathewson et al., 2009, 2011; Milton & Pleydell-Pearce, 2016; van Dijk et al., 2008). Even though we tested various stimulus strengths, we cannot assure that other types of stimulus, for instance, a continuous sinusoidal stimulus, would give the same result. Third, in this study, we mainly focused on the effect of a global stimulus and in addition, tested a local stimulus targeting the occipital region. However, considering the significant influence of network topology on local brain dynamics, further study is required to identify the effects of other regional stimuli.

Here, the state of the j_th oscillator, j = 1, 2, …, N, is determined by a complex variable z_j(t) at time t. N is the total number of brain regions acquired from group-averaged diffusion tensor imaging (DTI) with 82 nodes (van den Heuvel & Sporns, 2011). The A_jk = 1 if there is a connection between the j_th and k_th oscillators, and A_jk = 0 if there is not, based on the structural brain network. Each node undergoes supercritical Hopf bifurcation, and the dynamics of the oscillator settle on a limit cycle if λ_j > 0 and on a stable focus if λ_j < 0. Here, we fixed all λ_j as 1. The ω_j = 2πf_j and is an initial angular natural frequency of each j_th oscillator. We used a Gaussian distribution for natural frequency with a mean frequency of 10 Hz and standard deviation of 0.5 Hz to simulate the alpha bandwidth of human EEG activity (H. Kim et al., 2018; M. Kim et al., 2017; Lee et al., 2018; Moon et al., 2017, 2015). We also used the homogeneous coupling term K_jk = K between the j_th and k_th oscillators from 0 to 0.4 with δK = 0.002, which determines the global connection strength among brain regions. The time delay between oscillators, τ_jk = D_jk/s, was introduced with the average speed of axons in brain areas, s = 7ms (Caminiti et al., 2013), and the distance D_jk between brain regions. The node j received input from connected node k after the time delay τ_jk. The time delays are various, but as long as the time delay is smaller than a quarter of the period of oscillation, here τ_jk <∼ 25ms, results are not qualitatively different (Moon et al., 2015). A Gaussian white noise ξ_j(t) for each node was added with the standard deviation β = 0.05. We numerically solved the differential equations of the Stuart-Landau model by using the Stratonovich-Heun method with 1,000 discretization steps. We also tested the Runge-Kutta fourth-order method and the results were qualitatively same. The first 10 s of time series were discarded and last 25 s were used for the analysis of each simulation. Therefore, each brain region generates its own spontaneous oscillatory dynamics within the alpha bandwidth at each coupling strength K for one simulation.

We defined two different measures to quantify the responsiveness of stimulation at different instantaneous brain states. First, a responsivity R_j was defined by taking the average PR_j during a specific epoch after stimulation. It quantifies the total amount of significant responses after stimulation. Here we used t = 500 ms, which covers the maximum response changes after stimulation.

We first investigated the relationships between the level of instantaneous global synchronization at the stimulus onset r_s, s = 1, 2, …, 600, and the average R and LZc over all nodes at C_b, C_p, and C_a (Figure 4). Here we calculated spatial LZc, which is $LZct=ctN/log2N$⁠, N = 82, to see the global effect of the stimulation. Spearman correlations with p value were calculated between r_s and the average R and LZc at C_b, C_p, and C_a. A Kruskal-Wallis test with a multiple comparison test using the Tukey-Kramer method was used to see the statistical differences across brain states (***p < 0.001).

We also examined the instantaneous alpha amplitude/phase dependences of responsiveness (R and LZc) to the stimuli (Figure 5). Here we found the certain oscillatory conditions of alpha oscillations that can induce large or small responsiveness at C_p. We considered the oscillation properties with certain levels of instantaneous global synchronization, amplitudes, and phases that induce large (small) responsiveness as effective (less effective) stimulation conditions and applied these conditions to stimulate the local areas described below.

The same procedures were performed for the amplitude response; these results are shown in the supplementary materials (Supporting Information Figures S11–S13).

At C_p, we induced stimuli to 13 regions in the occipital area within a 50-mm radius of the left cuneus to confirm that the specific stimulation conditions we found from global stimulation could be applied to local stimulation (Figure 6). On the left and right hemispheres of the brain, respectively, the stimulated regions included: (left) cuneus, inferior parietal, isthmus cingulate, lateral occipital, lingual, pericalcarine, precuneus, superior parietal; (right) cuneus, isthmus cingulate, lingual, pericalcarine, precuneus. We selected the occipital region because abundant previous sensory stimulation studies mostly targeted the occipital region for visual tasks within the alpha frequency band. We compared the responsiveness measures with effective, less effective, and random stimulation conditions. For the effective (less effective) stimulation, we applied the stimulus to the alpha phases from 240° to 60° (60° to 240°) with low (high) levels of instantaneous global synchronization and low (high) amplitude. The stimuli were applied to the random instantaneous brain states of the occipital region for the random stimulation condition. The number of nodes for random stimulation was determined by a mean value of the number of nodes for effective and noneffective conditions. We performed 600 different trials at C_p, in which each trial consisted of stimuli with p = 30 for 50 ms, and classified the target nodes with specific conditions. We used a Kruskal-Wallis test to compare the significant differences in responsiveness for three conditions. A multiple comparison test was also performed using the Tukey-Kramer method (**p < 0.005 and ***p < 0.001).

Supporting Information for this article is available at https://doi.org/10.1162/netn_a_00113.

Minkyung Kim: Conceptualization; Formal Analysis; Investigation; Methodology; Visualization; Writing - Original Draft; Writing - Review & Editing. Uncheol Lee: Conceptualization; Funding Acquisition; Project Administration; Supervision; Writing - Original Draft; Writing - Review & Editing.

Uncheol Lee, Foundation for the National Institutes of Health (http://dx.doi.org/10.13039/100000009), Award ID: R01 GM098578.

Responsiveness: 

Here it is defined as responsivity and complexity of perturbed response (magnitude and complexity of synchronization and amplitude changes) after stimulation.

Ongoing brain activity: 

Spontaneous brain activity without any explicit tasks. In electroencephalography (EEG), it is referred to as the signal components that are not induced by external stimuli or other events.

Alpha oscillation: 

Neural oscillations in the frequency range of 8–13 Hz in brain.

Criticality: 

In physics, a system at criticality exhibits a state between ordered and disordered with the largest fluctuations it its dynamics. Emerging evidence shows that human brain in conscious wakefulness operates near criticality.

Global synchronization: 

A global quantity measuring consistent phase relationship among the population of oscillators.

Stuart-Landau model: 

A mathematical model to describe a nonlinear oscillating system near the Hopf bifurcation. It shows noisy to oscillating system with changes of parameters of the model.