Alpha oscillation, criticality, and responsiveness in complex brain networks.

Brains in unconsciousness are characterized by significantly limited responsiveness to stimuli. Even during conscious wakefulness, responsiveness is highly dependent on ongoing brain activity, specifically, of alpha oscillations (∼10 Hz). We hypothesized that the variety of brain responses to external stimuli result from the interaction between state-specific and transient alpha oscillations and stimuli. To justify this hypothesis, we simulated various alpha oscillations in the human brain network, modulating criticality (a balanced state between order and disorder), and investigated specific alpha oscillation properties (instantaneous amplitude, phase, and global synchronization) that induce a large or small response. As a result, we found that the alpha oscillations near a critical state show a more complex and long-lasting response, which is more prominent when stimuli are given to globally desynchronized and low-amplitude oscillations. We also found specific phases of alpha oscillation that barely respond to stimuli, which implies the presence of temporal windows in the alpha cycle for a large or small response. The results explain why brain responses are so variable across conscious and unconscious states and across time windows even during conscious wakefulness, and emphasize the importance of considering ongoing alpha oscillations for effective brain stimulation.

INTRODUCTION

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.

 

Schematic diagram of the study design. In this study, we simulated alpha oscillations, applying a coupled Stuart-Landau model to a human brain network structure consisting of 82 brain regions. Brain states were classified as below, near, or above a critical state based on the pair correlation function (PCF, the variance of the instantaneous global synchronization level). The critical state C was defined with the largest PCF. Global pulsatile stimuli u were applied to the 82 alpha oscillation nodes in three brain states, and the following perturbation responses of phase synchronization (PR(t)) for node jwere binarized by comparing them to the prestimulus period. Using PR, we first defined a responsivity R, an average of PR during certain epochs after stimulation, to measure the magnitude of the response. Second, the Lempel-Ziv complexity (LZc) of PR was calculated to measure the perturbational complexity of the response. Then we investigated the relationship between the ongoing alpha oscillation properties (global synchronization, amplitude, and phase) at stimulus onset and brain network responsiveness to find specific alpha oscillation conditions that induce large (or small) responsiveness. Finally, we stimulated the occipital region and compared effective, random, and less effective stimulation conditions.

RESULTS

Alpha Oscillations Near and Far from a Critical State

Various alpha oscillations in the human brain network were simulated in three brain states, C, C, and C (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, C, and C. Twenty different frequency configurations of the brain network were tested. The average levels of synchronization at C (orange triangle) are widely distributed between the small and large average levels of synchronization at C (blue circle) and C (yellow square). The alpha oscillations at C (orange triangle) show a larger PCF (mean ±SD = 1.39 ± 0.68) than those of C (blue circle) and C (yellow square), which suggests a large variance of brain network dynamics at C. Figure 2C and2D demonstrate the instantaneous synchronization r and instantaneous amplitude of the alpha oscillations at stimulus onset (s = 1, 2,…, 600, and j = 1, 2, …, 82). At C (orange), both the instantaneous synchronization levels rand amplitudes at are variable between the two small and large rand of C (blue) and C (yellow), respectively (in Figure 2C, median ± SD, 0.11 ± 0.04 for C, 0.32 ± 0.14 for C, and 0.78 ± 0.05 for C; in Figure 2D, 0.21 ± 0.04 for C, 0.48 ± 0.57 for C, and 1.75 ± 1.65 for C). 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.

 

Characteristics of alpha oscillations near and far from a critical state, C, C, and C. (A) Examples of global synchronization fluctuations at C (blue), C (red), and C (yellow) for one initial frequency configuration in the brain network. (B) The average instantaneous global synchronization r(t) and the pair correlation functions (PCF) for C (blue), C (red), and C (yellow) for 20 different initial frequency configurations in the brain network. Without a stimulus, the critical state C is characterized by the largest PCF (mean ± SD = 1.39 ± 0.68), compared with the PCF of C and C. (C) The estimated density distributions of instantaneous global synchronization (r) at the 600 stimulus onsets for C (blue), C (red), and C (yellow). The small white circles indicate the median of r, and the thick black lines indicate standard deviation. The estimated Gaussian density distribution for the 600 stimuli is plotted with a gray line. The alpha oscillations at C have a broad and balanced distribution (median ± SD = 0.32 ± 0.14). (D) The estimated density distributions of the alpha oscillation amplitudes at the 600 stimulus onsets for C, C, and C. The values at C and C are biased with small and large amplitudes, respectively, while the at C is relatively balanced between the others. The distinct alpha oscillation properties at stimulus onset may result in different responses of the brain network to the stimulus.

Distinct Responses of the Alpha Oscillations Near and Far from a Critical State

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, C, and C. 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(t)) for node j (j = 1, 2, …, 82) at time t, setting PR(t) as 1, if the synchronization value of node j was significantly increased from the prestimulus period (1 s); otherwise PR(t) was set as 0 (see Materials and Methods for details). Figure 3B presents examples of PR. The black dot represents PR(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(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 persists longer than the PRs at C and C (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.

 

Quantifying the perturbation responses (PRs). (A) Examples of the global synchronization level r(t) for the three brain states, C (left), C (middle), and C (right). A triangle indicates the timing of a stimulus. The response to a single pulsatile stimulus was stimulated repeatedly 600 times with random stimulus timings. (B) The PR of a pulsatile stimulus is presented for the three brain states. The red and blue shaded areas indicate the stimulus period (50 ms) and response period (500 ms), respectively. The black dot at each node indicates the time point at which phase synchronization significantly increased after the stimulus compared with the prestimulus period; PR = 1. (C) The average PR(t) was calculated over 600 stimulation trials, using a new time axis in which the timing of the stimulus coincides with zero. The thick lines indicate the median of the average PR, and the shaded areas cover the 25% to 75% quantiles. The PR persists longer at C than at C and C (sign test, p < 0.05 after 190 ms).

Magnitude and Complexity of Perturbation Responses

In this study, we defined two responsiveness indexes to quantify the PRs. First, we defined the responsivity R of a node j, counting the total amount of 1s in PR (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 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).

Correlation Between Responsiveness and Synchronization of Alpha Oscillations

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

 

Correlations between instantaneous synchronization (r) at stimulus onset and the responsiveness of the brain network. (A) Spearman correlation coefficient ρbetween r and the average responsivity ¡R¿ and (B) Spearman correlation coefficient ρ between r and the average complexity of perturbation response ¡LZc¿ for C (blue), C (red), and C (yellow). The Spearman correlations were calculated across 600 stimulation trials. The average Spearman correlation coefficients at C for both responsiveness indexes are −0.49 and −0.51, respectively, with the significance level ***p < 0.001. At C, applying a stimulus at a lower level of global synchronization induces a larger magnitude and complexity of perturbation response. A Kruskal-Wallis test with a Tuckey-Kramer multiple comparison test was performed to compare the correlation coefficients across the three states (***p < 0.001).

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 and < LZc > for the three brain states. Only at C is < LZc > significantly correlated with r (Spearman correlation ρ = −0.52,p < 0.001 for C). At C 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, C, and C).

We found that near a critical state, stimulating a globally more desynchronized alpha oscillation (small r) 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.

Specific Phases of Alpha Oscillation Unresponsive to Stimulus

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 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 (instantaneous global synchronization level, s = 1, 2, …, 600) and (instantaneous node amplitude, j = 1, 2, …, 82)—at stimulus onset. These nodes were classified into high/low r (HS/LS) and high/low (HA/LA) by comparing them to the average r and , respectively. We further separated the phases of node 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 and C, 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 and C; 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 and C; however, in this study we focused on the response behavior at critical state C to interpret the empirical results during conscious wakefulness.

 

Phase dependence of responsiveness at C. (A) For each stimulation trial, the responsivity and complexity (the number of trials, s = 1, …, 600, the number of nodes, j = 1, …, 82) were calculated for each trial and each node; each node has the instantaneous alpha oscillation properties at stimulus onset: and . These oscillation properties were separated into four groups. In the main text, we show the results of two groups: low synchronization and low amplitude (LS and LA) and high synchronization and high amplitude (HS and HA). The phases of alpha oscillation were segmented into 30 phase bins of 12 intervals for each group. Results for other groups are presented in Supporting Information Figure S2. (B) The responsivity R and (C) the complexity LZc of perturbation responses for the HS and HA group significantly depend on the phase of alpha oscillation. The R and LZc of specific phases (around 180°) show a lower responsiveness to stimuli. A thick red (blue) line indicates the median value; the shaded areas cover the 25% to 75% quantiles. (D and E) For statistical evaluation, the phases of alpha oscillation were separated into two regimes (60°− 240° and 240°− 60°). In the LS and LA group, stimulation induced greater responsiveness than that in the HS and HA group for both (D) R and (E) LZc (p < 0.001, multiple comparison test using the Tukey-Kramer method). Phase dependence is more prominent in the HS and HA condition than the LS and LA condition, showing larger responsiveness (R and LZc) in the 60°− 240° regime (***p < 0.001, Wilcoxon rank sum test). In the case of high-amplitude, highly synchronized alpha oscillations, the results imply the existence of temporal windows (with an alpha frequency cycle) that barely respond to external stimuli. The error bar indicates the standard deviation.

Testing Global Stimulus Conditions with a Local Stimulus

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.

 

Testing the responsiveness of a local stimulus under (effective, random, and less effective) global stimulation conditions. Three global stimulation conditions were tested to see if they also hold for local stimulation in the brain network. (A) R and (B) Lempel-Ziv complexity (LZc) of the occipital region stimulation were compared under three global stimulation conditions. Under the effective stimulation condition, specific phases of alpha oscillation (240° to 60°) were stimulated with low synchronization and low amplitude (red). Under the random stimulation condition, stimulation of alpha oscillations was random (gray). In the less effective stimulation condition, specific phases of alpha oscillation (60°− 240°) were stimulated with high synchronization and high amplitude (blue). For occipital region stimulation, the effective (less effective) stimulation conditions produce larger (smaller) responsiveness. The responsiveness of the three conditions is significantly different in both R and LZc (Kruskal-Wallis test, p < 0.001). Multiple comparison tests were performed with the Tukey-Kramer method (**p < 0.005 and ***p < 0.001). The results showed that the effective and less effective stimulation conditions of global stimulation also hold for local stimulation.

DISCUSSION

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, C, and C) 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), low amplitude was associated with an incoherent state; near the critical point (C), high and low amplitudes were balanced with a large degree of synchronization fluctuation; above the critical point (C), 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.

Dependence on Network Criticality

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 and C was larger than that at C (Figure 3D). However, the persistence of responsivity is longer at C, which implies restricted propagation at the states far from criticality, C and C (Figure 3D). The brain network at C 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.

Dependence on Instantaneous Global Synchronization

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.

Dependence on the Amplitude and Phase of Ongoing Oscillations

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.

Multi-Scale Mechanisms of Large and Small Responsiveness

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.

Conclusion

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.

Limitations

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.

MATERIALS AND METHODS

Simulation of Networked Alpha Oscillations in a Human Brain Network Structure

We constructed a large-scale functional brain network using coupled Stuart-Landau models on the human anatomical brain structure to generate spontaneous oscillations. The Stuart-Landau model with brain topology has been used to replicate the oscillatory dynamics of different types of electromagnetic brain signals such as electroencephalography (EEG), magnetoencephalography (MEG), and functional magnetic resonance imaging (fMRI) (Deco et al., 2018, 2017a; Deco, Kringelbach, Jirsa, & Ritter, 2017b; H. Kim et al., 2018; Moon et al., 2015). The coupled Stuart-Landau model is defined as the following:

Here, the state of the j oscillator, j = 1, 2, …, N, is determined by a complex variable z(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 = 1 if there is a connection between the j and k oscillators, and A = 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 λ > 0 and on a stable focus if λ < 0. Here, we fixed all λ as 1. The ω = 2πf and is an initial angular natural frequency of each j 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 = K between the j and k oscillators from 0 to 0.4 with δK = 0.002, which determines the global connection strength among brain regions. The time delay between oscillators, τ = D/s, was introduced with the average speed of axons in brain areas, s = 7ms (Caminiti et al., 2013), and the distance D between brain regions. The node j received input from connected node k after the time delay τ. The time delays are various, but as long as the time delay is smaller than a quarter of the period of oscillation, here τ <∼ 25ms, results are not qualitatively different (Moon et al., 2015). A Gaussian white noise ξ(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.

Three Brain States: Below, Near, and Above Critical State

We selected three representative brain states at different coupling strengths based on the level of global synchronization and the pair correlation function (PCF) to understand the responsiveness of different systems. We first calculated an instantaneous global synchronization level r(t) at time t using phase difference Δθ(t) = θ(t) − θ(t) at each coupling strength K. Here r(t) = 1 if all phases are equal, but r(t) is nearly zero if all phases are randomly distributed. Then we calculated the variance of r(t) for each coupling strength as a PCF to define the critical state C (Shanahan, 2010). The PCF is equivalent to susceptibility in statistical physics models. Coupled phase oscillators in complex networks showed the largest susceptibility and PCF at critical state (Yoon et al., 2015). We considered generated signals at certain coupling strengths with the largest PCF to be equivalent to spontaneous alpha oscillations in conscious wakefulness. There has been emerging evidence in computational model studies that brain dynamics at critical point show the largest spatiotemporal diversity (Haimovici et al., 2013; Hudetz, Humphries, & Binder, 2014). Notably, a large repertoire of metastable states exists in brain dynamics during the conscious state (Deco et al., 2017b). With the critical state, we selected two other representative states, one below (C) and one above (C) critical state. We defined C and C as the states at the 10th and 90th percentiles of the averaged order parameters over all coupling strengths, respectively. These states are the representative states far from the critical states that show small PCF. Generated signals in these three states were used for the analysis. We used 20 different initial frequency distributions and selected three representative states for each frequency distribution.

Brain Network Stimulation Procedure

Global pulsatile stimuli were induced as a stimulation with u(t) as the following: Here p is the strength of the stimulus during a period T=t2 − t1. We fixed the p = 10 for the analysis after testing the effects of various stimulus strengths (Supporting Information Figure S11). We also fixed the duration of the stimulus, T = 50ms, which is about half of the period of the alpha frequency range cycle. We first gave the same stimulus to the whole-brain network to understand the dependence of the brain’s response derived only from the network dynamics of the system (as opposed to its network structure). The stimuli were induced at randomly selected timing t1 for one trial. In total, 30 different timings, each with 10 iterations, were selected for 20 different frequency distributions. Therefore, C, C, and C each underwent 600 stimulation trials of pulsatile stimuli. Each stimulus was applied at a unique instantaneous brain state (i.e., a different one was used for each trial). We decomposed the instantaneous brain states into the level of instantaneous global synchronization and the amplitude and phase of alpha waves to identify the relationship of each of these separate factors with responsiveness.

Calculation of Significant Response After Stimulation

First, we calculated the instantaneous phase synchronization of the j node, , for each trial. The was obtained by taking the average of the pairwise synchronization between node j and node k at time t, Here N is the number of connections of the j node. The instantaneous synchronization value for the j node of one trial was normalized by the mean and standard deviation of the baseline synchronization values of the j node. Baseline values were obtained by using a total of 10 s, consisting of 10 iterations of a 1-s prestimulus segment for each trial. We considered the one tail (1 − α) * 100 quantile with α = 0.05 as a significantly increased synchronization after stimulus onset. A perturbation response (PR) of the j node at time t was defined in a binary fashion: PR(t) = 1, for the significantly increased synchronization of node j, and PR(t) = 0, otherwise.

Responsiveness Measures: Responsivity and Perturbational Complexity

We defined two different measures to quantify the responsiveness of stimulation at different instantaneous brain states. First, a responsivity R was defined by taking the average PR 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.

Second, we also calculated the Lempel-Ziv complexity (LZc) of the PR to understand how complex the response was. The LZc of the PR was defined as the following: Here (c) is the nonnormalized Lempel-Ziv complexity of PR calculated by a LZ76-algorithm (Lempel & Ziv, 1976), that is, the number of unique patterns in the PR during time t. We defined a perturbational complexity, LZc, as the normalized c.

We first investigated the relationships between the level of instantaneous global synchronization at the stimulus onset r, s = 1, 2, …, 600, and the average R and LZc over all nodes at C, C, and C (Figure 4). Here we calculated spatial LZc, which is , N = 82, to see the global effect of the stimulation. Spearman correlations with p value were calculated between r and the average R and LZc at C, C, and C. 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. 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).

Local Stimulation Effect

At C, 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, 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

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

ROLE INFORMATION

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.

FUNDING INFORMATION

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

Click here for additional data file.