Explosive Synchronization-Based Brain Modulation Reduces Hypersensitivity in the Brain Network: A Computational Model Study

Fibromyalgia (FM) is a chronic pain condition that is characterized by hypersensitivity to multimodal sensory stimuli, widespread pain, and fatigue. We have previously proposed explosive synchronization (ES), a phenomenon wherein a small perturbation to a network can lead to an abrupt state transition, as a potential mechanism of the hypersensitive FM brain. Therefore, we hypothesized that converting a brain network from ES to general synchronization (GS) may reduce the hypersensitivity of FM brain. To find an effective brain network modulation to convert ES into GS, we constructed a large-scale brain network model near criticality (i.e., an optimally balanced state between order and disorders), which reflects brain dynamics in conscious wakefulness, and adjusted two parameters: local structural connectivity and signal randomness of target brain regions. The network sensitivity to global stimuli was compared between the brain networks before and after the modulation. We found that only increasing the local connectivity of hubs (nodes with intense connections) changes ES to GS, reducing the sensitivity, whereas other types of modulation such as decreasing local connectivity, increasing and decreasing signal randomness are not effective. This study would help to develop a network mechanism-based brain modulation method to reduce the hypersensitivity in FM.


INTRODUCTION
Hypersensitive responses to external stimuli have been widely observed in various physical and biological systems such as cascading failures in power grids, abrupt state transitions in an electronic circuit, abrupt loss and recovery of consciousness in general anesthesia, and epileptic seizures in the brain (Chen et al., 2013;Boccaletti et al., 2016;Wang et al., 2017). Fibromyalgia (FM), a chronic pain disorder, is characterized by fatigue, poor memory, sleep problems, and mood disturbance (Hawkins, 2013;Menzies, 2016). Many of these individuals also present a hypersensitive response to external sensory stimuli, which is regarded to involve central sensitization associated with structural and functional changes in the brain (Harris et al., 2013;Harte et al., 2018). In our previous study, we found that explosive synchronization (ES; a first-order phase transition in a network) to be an underlying mechanism of the hypersensitivity in FM brain . A strong positive correlation was detected between the strength of ES conditions and chronic pain intensity in FM patients. This suggests that specific topological and functional network configurations of the brain could induce the first-order phase transition (abrupt state transition) against an external stimulus (Kim et al., 2016Lee et al., 2018).
Global synchronization in a network is initiated by the merging of synchronization clusters across nodes. In a brain network following a general synchronization (GS) path, hubs (nodes with intense connections) lead a global synchronization of the network, forming a large synchronization cluster (Gómez-Gardeñes et al., 2007). Once a large synchronization cluster has formed, it grows gradually, absorbing smaller local clusters. However, if the growth of the largest cluster is somehow suppressed by pharmacological, neurological, or pathological types of perturbation, it allows disconnected clusters to grow until the network reaches a critical threshold wherein a small stimulus triggers a singular explosive unification of all clusters termed ES. In these synchronization processes, hubs play a central role in creating a giant synchronization cluster in the brain network (Zhang et al., 2015;Kim et al., 2017). Therefore, we expected that modulating the hubs in the brain network can be the most effective and practical stimulation target to convert the type of state transition, from ES to GS, in brain networks.
A previous model study demonstrated the possibility to convert an ES network into a GS network by modulating the ratio of inhibitory to excitatory connectivity (Zhang et al., 2016). However, modulating the ratio of inhibitory connectivity in the global brain network may not be feasible with current stimulation techniques. Furthermore, it has not been tested whether the brain network state is near criticality, which presents complex brain dynamics of a conscious state, is convertible from ES into GS by a brain modulation. In principle, reversing the core mechanism of ES may convert an ES network into a GS network, but this would require reversing suppression of the giant synchronization cluster formation. In essence, this may change an abrupt state transition to a gradual state transition and also significantly reduce the network sensitivity. We hypothesized that enhancing the hub connectivity may be a way to facilitate the giant synchronization cluster formation, which in turn can convert an ES network to one of GS.
To justify the hypothesis, we constructed a large-scale human brain network model occurring ES, using a modified Stuart-Landau oscillator model on the anatomically informed human brain network structure . After constructing ES brain network models, two different network properties were modulated as control parameters for the brain network modulation: (1) the local structural connectivity of brain regions and (2) the randomness of node dynamics within a certain diameter centered on a target node. The correlation between node degrees and frequencies, which is one of the ES conditions, was measured at critical states of the brain networks before and after the modulation. To evaluate which is the most effective type of brain network modulation for reducing the brain network sensitivity, we directly applied an external stimulus to the brain network models and measured the change of brain network sensitivity with surrogate measures-responsivity and Lempel-Ziv complexity (Kim and Lee, 2020).

MATERIALS AND METHODS
The whole computational modeling procedures are summarized as follows (see Figure 1 for a schematic diagram of the study design).
1. We constructed a large-scale human brain network model containing ES as in the FM brain. We assumed that the large network sensitivity of ES is the underlying mechanism of hypersensitivity of FM brain, and the degree of ES condition correlates with the pain intensity of FM patients according to our previous empirical study . 2. We then modulated the degree of ES condition in the human brain network model with four different types of network modulation: increasing or decreasing the local structural connectivity and increasing or decreasing the randomness of node dynamics within a certain diameter centered on a target node. 3. The critical states for the brain network models were identified by searching through parameter space, at which the simulated brain activities reflect the characteristic brain dynamics during a conscious resting state. The brain networks with different ES and GS conditions had their own distinctive critical states. 4. We evaluated the ES condition (correlation between node degrees and frequencies) and the change of brain network sensitivity before and after a brain network modulation. We directly induced a global pulsatile stimulus to the brain network models before and after the modulations and quantified the changes of brain network sensitivity with surrogate measures-responsivity and Lempel-Ziv Complexity. 5. Finally, we identified the type of network modulation and effective stimulation sites that can reduce the ES condition and the brain network sensitivity.

Construction of Explosive Synchronization Human Brain Network Model
We used a coupled Stuart-Landau model to simulate the oscillatory dynamics of the human brain network. The coupled Stuart-Landau model with an anatomically informed brain network structure has been widely applied to simulate the signals from various types of imaging modalities including electroencephalogram (EEG), magnetoencephalogram (MEG), and functional magnetic resonance imaging (fMRI) (Cabral et al., 2014(Cabral et al., , 2017Moon et al., 2015;Deco et al., 2017Deco et al., , 2018 FIGURE 1 | Schematic diagram of the study design. We simulated the fibromyalgia (FM) brain network using a model with the explosive synchronization (ES) mechanism. The ES brain network model was constructed using a modified coupled Stuart-Landau oscillator on an anatomically informed human brain network structure. We found a critical state of ES brain network model by calculating autocorrelation function (ACF) to simulate the FM brain during conscious wakefulness. Four different types of network modulation (local structural connectivity increase and decrease, signal randomness increase and decrease) with thirty different target brain regions were applied to the model to investigate which modulation types and which target brain regions can convert an ES network to one of general synchronization (GS), thereby reducing the sensitivity. Then, we induced external stimuli to the brain networks at a critical state before and after the modulation, evaluated sensitivity (responsivity and complexity) of the brain network responses, and compared the sensitivity between the networks before and after the modulation. Kim et al., 2018;Kim and Lee, 2019). Since the goal of this study is to convert ES brain network into GS brain network, we used a modified Stuart-Landau model including an adaptive feedback term. Our previous study introduced the model to simulate a hysteresis phenomenon during general anesthesia . With this model, we can adjust the ES strengths of the brain network by modulating the adaptive feedback term R Z j , that is, a recursive interaction process between a node j and the other nodes. The modified coupled Stuart-Landau model, composed of the N number of oscillators, is defined as follows: Here, r j (t) is the amplitude of oscillator (node) j at time t. j is a parameter modulating the randomness of the amplitude dynamics, and S is a coupling strength among anatomically connected oscillatory dynamics. Changing these parameters induces various amplitude and phase dynamics through the competition of independent node dynamics and its topological connections with other oscillators upon an anatomical brain network structure (Cavanna et al., 2018). Each node shows supercritical Hopf bifurcation, and the dynamics of the oscillator settle on a limit cycle if λ j > 0, and on a stable fixed point if λ j < 0. A jk is the anatomical connection weight between oscillator j and k. The connection matrix A consists of 82 brain regions including cortical and subcortical regions which were constructed from diffusion tensor imaging (DTI) (van den Heuvel and Sporns, 2011). A certain value from 0 to 1 was assigned to A jk based on the white matter connection weight between brain region j and k. The τ jk = D jk /s is a time delay between region j and k, where D jk is a distance between brain regions with s = 7 ms, an average speed of axons across brain regions (Caminiti et al., 2013). The node j interacts with a connected node k after the time delay τ jk . j (t) is a phase of the oscillator j at time t, and j is a natural frequency of the oscillator j. R j is defined as 1/2(e iθ j + 1/N N k = 1 e iθ k ), which is the extent of synchrony of node j with the other nodes. The R j enhances a heterogeneity in the amplitude and phase dynamics of the brain network, incorporating the memory of a given synchronized state into the dynamics. The Z is a scale term for the adaptive feedback process R j , and it plays a role of a control parameter for adjusting ES strengths of the brain network model, which indicates it inhibits a giant synchronization cluster formation centered around a node j with large synchrony, which is a core mechanism of ES . The effect of Z was first introduced by Filatrella et al. (2007) with Kuramoto model that has phase dynamics only. We then applied the Z term with some modifications to Stuart-Landau model that has both phase and amplitude dynamics to study the hysteresis phenomenon during loss and recovery of consciousness in general anesthesia . Furthermore, we showed that the adaptive feedback term is equivalent to the Hill coefficient (the dose-response slope of anesthetics) in the dose-response equation. The relationship between the adaptive feedback term and the hysteresis size was analytically studied (Khanra et al., 2020), which also shows the relationship between the adaptive feedback term and ES.
We used Gaussian distribution for the natural frequency ω with a mean frequency of 10 Hz and a standard deviation of 0.4 Hz to simulate the dominant frequency bandwidth of human EEG alpha activity (Moon et al., 2015(Moon et al., , 2017Kim et al., 2017;Kim et al., 2018;Lee et al., 2018). We first fixed λ j λ as 0. The coupling strength among the oscillators S was modulated from 0 to 30 with δS = 0.2, which yields the change of brain network from a fully incoherent state to a fully synchronized state. We set Z = 3, which is large enough to simulate ES in the brain network. We numerically solved differential equations of the Stuart-Landau model using the Runge-Kutta 4th method with 1,000 discretization steps and resampled the data with 250 Hz. For every coupling strength, the last 40 s of spontaneous oscillatory dynamics was used for the analysis after 20 s of saturation periods. Thirty different initial frequency configurations were simulated in each parameter set, and the results were averaged over all configurations.

Four Types of Brain Network Modulation to Convert Explosive Synchronization Into General Synchronization Network
According to our hypothesis, enhancing hub connectivity between brain network nodes may facilitate giant synchronization cluster formation, which in turn converts ES network into GS network, reducing the brain network sensitivity significantly. In addition, it is known that repetitive transcranial direct current stimulation (tDCS) targeting a local brain region can induce a long-lasting change in the brain network activity through neuroplasticity and cortical excitability shifts (Nitsche et al., 2007). Therefore, we tested four types of brain network modulation: increasing or decreasing the regional structural connectivity centered on a target node and increasing or decreasing the randomness of the regional activities centered on a target node in the ES human brain network model. These types of modulation may reflect neuroplasticity and cortical excitability shifts that are empirically observed in tDCS experiments. For changing the regional connectivity, we increased or decreased the local structural connectivity by three times within a radius of 2 cm centered on the targeted brain region, roughly simulating the effective range for connectivity change covered by a tDCS patch. For changing the randomness of regional brain activities, we increased or decreased the bifurcation parameter λ of nodes, respectively λ = 2 and λ = − 2, within a radius of 2 cm centered on the target brain region. The parameters were chosen for simulating the most realistic results after the parameter test. The 30 highest-degree nodes in the human brain network were selected as the target brain regions and tested the modulation effect for each target brain region (left and right putamen, amygdala, pallidum, thalamus, hippocampus, caudate, temporal pole, insula, precuneus, entorhinal cortex, precentral cortex, superior parietal cortex, accumbens area, isthmus cingulate cortex, left superior frontal cortex, and left postcentral cortex). The initial frequency configurations of the brain network models are all the same for the four types of modulation.

Identification of Critical States for the Brain Network Model
Recent computational modeling and empirical studies suggest that brain dynamics in conscious states are near criticality (i.e., an optimally balanced state between order and disorder), and losing criticality is related to altered states of consciousness (Kitzbichler et al., 2009;Tagliazucchi et al., 2012;Haimovici et al., 2013;Muñoz, 2018;Kim and Lee, 2019). The network dynamics at a critical state reflect the characteristics of the brain activities in the conscious resting state: an optimal balance between stability and instability, optimal information processing, large flexibility to adapt to a changing environment, and wide repertoires of brain states (Beggs, 2008;Kitzbichler et al., 2009;Tagliazucchi et al., 2012;Haimovici et al., 2013;Cocchi et al., 2017;Kim and Lee, 2019). Therefore, we assumed that the critical states of a brain network before and after modulation of ES condition might reflect the distinctive brain activities in conscious states before and after the brain modulation. In this modeling study, we searched the parameter space of various coupling strengths and bifurcation parameters and determined the parameter set having a maximal autocorrelation function (ACF) as a critical state. The maximal ACF is one of the characteristics of a system approaching a critical transition, which is called "critical slowing down" and refers to the tendency of a system to take longer to return to equilibrium after a perturbation (Scheffer et al., 2009). To calculate the ACF of a simulated brain signal, we selected the time lag of 50 ms, which catches the dynamics of the alpha oscillations in the simulated brain signals.

Evaluation of Brain Network Sensitivity for the Four Types of Brain Network Modulation
To test which type of brain modulation effectively reduces the brain network sensitivity and changes the type of state transition from ES (first-order phase transition) to GS (second-order phase transition), we investigated three different properties. (1) We examined the typical phase transitions of ES and GS. Converting ES into GS should present the typical change in transition pattern from a discrete synchronization transition (first-order phase transition) to a continuous synchronization transition (secondorder phase transition). (2) We examined whether a brain network modulation induces the typical network configuration of ES, a positive correlation between node degrees and node frequencies in the brain network. (3) Then, we directly measured the change of network sensitivity after the modulation directly implementing an external stimulus to the brain network.
To investigate the shape of synchronization transition, the instantaneous network synchronization r(t) at time t was measured by the order parameter of the oscillators, where θ j (t) is a phase of j th oscillator and ψ(t) is the average global phase at time t. Here, r(t) equals 0 if phases of oscillators are uniformly distributed and 1 if all oscillators have the same phase. The global network synchronization R is calculated by taking the average of the instantaneous network synchronization over time. As the coupling strengths increase in the brain networks before and after the four types of modulation, we investigated the shape of R to determine whether the transition is discrete or continuous. A positive correlation between node degrees and frequencies was introduced as one of the network configurations that can induce ES in a heterogeneous network (Gómez-Gardeñes et al., 2011). Therefore, we tested which type of brain network modulation and which target node effectively mitigate the ES condition. The degree of j th node was calculated by G j = N k = 1 A jk and the frequency of the j th node was calculated using the instantaneous phases of the j th node in the brain network at a critical state. A Spearman's correlation between the degrees and frequencies, ρ deg−freq , was calculated for all the brain networks before and after the modulation of different target nodes.
Finally, we directly measured the network sensitivity against an external stimulus and investigated which type of brain network modulation significantly reduces the brain network sensitivity. For the brain networks before and after the four types of modulation, we induced a global pulsatile stimulus to the brain networks and quantified the change of the network sensitivity. The global pulsatile stimulus was induced with u(t) as follows: Here, p is the intensity of the stimuli during a period T = t 2 − t 1 . We chose a stimulus strength p = 20 and a duration T = 100 ms based on our previous study (Kim and Lee, 2020).
In the previous study using the same model, we demonstrated a significant dependence of brain responsiveness on spontaneous brain activities at the stimulus onset and found that the stimulus strength p = 20 and the duration T = 100 ms were strong enough to rule out the potential of such dependence (Kim and Lee, 2020). The global stimulus was given at a randomly selected timing t 1 for each trial. Twenty different trials were performed for each 30 different initial frequency configurations. Therefore, the same external stimulus was given 600 times to all the brain networks before and after modulation. The sensitivity was obtained by measuring the responsivity and Lempel-Ziv complexity of the brain network's response after the stimulation (Kim and Lee, 2020). To measure the responsivity and complexity, we first calculated an instantaneous amplitude of the nodes for each stimulation trial. The instantaneous amplitude value for the j th node of one trial was normalized by the mean and standard deviation of the baseline amplitude values of the j th node. The baseline amplitude values were obtained from a total of 200 sec, consisting of 20 different stimulation trials of 10sec prestimulus segment for each trial. We considered (1 − α) * 100 th percentile with α = 0.05 as a significantly changed amplitude after a stimulus onset. A perturbation response (PR) of the amplitude for the j th node at time t was defined in a binary fashion: PR j (t) = 1, if the amplitude after the stimulus onset is significantly changed, and PR j (t) = 0, otherwise.
The responsivity P and Lempel-Ziv complexity LZc were calculated with the binary PR during 300 ms after the stimulus onset. The responsivity P was calculated by taking the average of PR(t) for all nodes over 300 ms. The LZc(t) of the PR(t) was defined as below: Here, c(t) is the nonnormalized Lempel-Ziv complexity calculated by LZc76-algorithm (Lempel and Ziv, 1976), which is the number of unique patterns in the PR(t) at time t, and N is the number of nodes. The LZc(t) is the normalized c(t).
The LZc was calculated by taking the average of LZc(t) over 300 ms after the stimulus onset. The P and LZc were calculated for each stimulation trial. The sensitivity results for each target brain regional modulation were the average of the 30 frequency configurations with the 20 stimulation trials.

Increasing Regional Brain Connectivity Converts Abrupt Transition Into Gradual Transition
Four types of brain network modulation [local structural connectivity increase: CI; local structural connectivity decrease: CD; local randomness (bifurcation parameter) increase: RI; local randomness (bifurcation parameter) decrease: RD] were tested to convert an ES network into a GS network. For each type of modulation, we applied the modulation to the brain networks centered on thirty different highest-degree nodes to investigate the effect of regional modulation. Figure 2 presents representative examples for the shape of synchronization transition before (blue line) and after (red line) the modulation. The examples are the shapes of synchronization transition of brain network modulation targeting the right precuneus for each type of modulation. The simulated brain network before the modulation showed an abrupt synchronization transition near the critical state (blue lines in Figures 2A-D), which reflects ES (the first-order phase transition). Only CI modulation converted the type of synchronization transition from ES to GS (red line in Figure 2A), whereas the other types of modulation (CD, RI, and RD) did not (red lines in Figures 2B-D). The shapes of synchronization transition targeting other nodes of all types of modulation were presented in Supplementary Figures 1-4. The CI modulation targeting other nodes mostly showed the shape change of synchronization transition from abrupt to gradual (Supplementary Figure 1). Other types of modulation (CD, RI, and RD) did not show the shape change in most of the nodes (Supplementary Figures 2-4).

Increasing Regional Brain Connectivity Mitigates Explosive Synchronization Condition
We examined a correlation between node degree and node frequency (ρ deg−freq ) in the brain network to investigate whether the four types of modulation disrupt a typical ES condition. In generic network models, a larger positive ρ deg−freq induces ES with a higher probability. Figures 3A-D, respectively, show the ρ deg−freq values before and after the four types of modulation for the thirty target nodes. The gray shaded areas in Figure 3 indicate the mean standard error of ρ deg−freq before the modulation and the square dots with error bars indicate the mean standard error of ρ deg−freq after the modulation. The name of the target brain regions is presented on the x-axis.
We found that only CI modulation gives rise to significant decreases of ρ deg−freq for the most target node modulations (left pallidum, thalamus, hippocampus, amygdala, temporal pole, insula, entorhinal cortex, precentral cortex, postcentral cortex, isthmus cingulate, right thalamus, temporal pole, cortex, right precuneus, entorhinal cortex, and accumbens area; t-test, * p < 0.05, Figure 3A), whereas the CD modulation showed an increase of ρ deg−freq for some target nodes (left and right putamen, right amygdala, left and right pallidum, right thalamus, left amygdala, left and right caudate, right insula, and right entorhinal cortex, Figure 3B). The RI and RD modulations did not change ρ deg−freq or only a few target node modulations induced alternation (Figures 3C,D). The comparison among ρ deg−freq values for four types of modulation is presented in Supplementary Figure 5. The CI modulation was significantly lower than other types of modulation. In addition, we investigated the change of signal properties such as a global synchronization and a signal amplitude, which is associated with the responsiveness of a network against stimulation (Kim and Lee, 2020). Only CI modulation significantly increased the global synchronization and the amplitude of the signal (Supplementary Figure 6), which implies the decreased responsiveness and the mitigated ES condition.

Increasing Regional Brain Connectivity Decreases Network Sensitivity
Next, we implemented an external stimulus and compared the changes in brain network responses before and after the modulation to directly evaluate the brain network sensitivity change. Responsivity P and complexity LZc were calculated to evaluate the brain network sensitivity with perturbed signals during 300 ms after the stimulation. Only CI modulation with specific target brain regions showed decreases in both responsivity and complexity. Figure 4 shows the responsivity (left) and complexity (right) for the CI modulation. The names of the target brain regions are presented in the x-axis. The gray area covers the 25-75% values of P and LZc before the CI modulation. The blue and green shaded areas in Figures 4A,B, respectively, denote the 25-75% quantiles of P and LZc after the CI modulation for 30 different target brain regions. The red star markers indicate the effective target brain regions where the network sensitivity is significantly decreased (ttest, * p < 0.05). Increasing the regional brain connectivity centered on subcortical regions and some cortical hub regions (left putamen, pallidum, amygdala, insula, entorhinal cortex, accumbens, isthmus cingulate cortex, right putamen, amygdala, pallidum, hippocampus, temporal pole, caudate nucleus, insula, precuneus, and entorhinal cortex) significantly decreases the brain network sensitivity. The responsivity of the brain networks in the other types of modulation (RI and RD) was also decreased significantly; however, the changes in responsivity were not (Supplementary Figure 7). Some extreme cases such as λ = 5 or λ = − 5 are also presented in Supplementary Figure 8.

Hub Nodes Are Effective Target Sites to Reduce Brain Network Sensitivity
We hypothesized that the hub nodes in the brain network can be effective target sites to reduce brain sensitivity. Therefore, we evaluated the relationship between the node degree of target brain regions and the sensitivity of the node modulation. Figures 5A,B presents the relationship between the average node degree of target brain regions and the responsivity (complexity) of the target node CI modulation. Both responsivity and complexity showed significant negative Spearman's correlations with the node degree of target brain regions (-0.52 and -0.43, The gray area indicates the mean standard errors of ρ deg−freq over 30 different initial conditions of brain networks before the modulation. The colored square with error bars indicates the mean standard errors of ρ deg−freq over 30 different initial conditions of brain networks after the modulation. Each marker indicates the centered target node for the modulation. A target node presenting a significant change after the modulation is marked with "*" (t-test, *p < 0.05). The CI modulation induces a significant decrease of ρ deg−freq for most of the target nodes. In contrast, other types of modulation induce significant changes for the much smaller number of nodes (CD and RI). The RD modulation shows no change for all nodes.
FIGURE 4 | Changes in brain network responsivity and complexity after the I modulation. The (A) responsivity and (B) complexity of the brain networks before and after the CI modulation are presented. The gray area covers 25-75% values of network responsivity (sensitivity) before the CI modulation over different initial conditions. The blue (green) colored area covers 25-75% values of network responsivity (complexity) after modulation over different initial conditions. A target node presenting significantly different responsivity (complexity) is marked as " " (t-test, *p < 0.01). The CI modulation to the right and left insula, right precuneus, and left isthmus cingulate cortex results in decreased responsivity and complexity in the brain network.
respectively), which suggests that the CI modulation targeting hub nodes gives rise to lower responsivity and complexity than targeting peripheral nodes. The large correlation between network sensitivity and node degree implies the existence of effective target sites in the brain network, which also reveals the possibility of developing a systematic brain modulation method with better outcomes considering the network topology.

DISCUSSION
Hypersensitivity of ES may provide an opportunity for a system to have flexible adaptation to external stimuli with high susceptibility and fast response, but also it implies high risks of unwanted excessive responses that might cause a neurologic problem in the brain which may explain previous neuroimaging FIGURE 5 | The relationship between the node degree of target brain regions and the sensitivity after the CI modulation. The sensitivity was evaluated by responsivity and complexity. The relationships between the node degree of target brain regions and the (A) responsivity and (B) complexity are presented. The Spearman's correlation coefficient between the network responsivity (complexity) and the average degree of target nodes is -0.52 (-0.43). The CI modulation to larger degree nodes (i.e., hubs) induces more decrease in the brain network sensitivity.
findings in FM. Empirical evidence of ES in FM brain network and the possibility for converting ES into GS network in previous model studies motivated us to investigate brain modulation methods that could reduce hypersensitivity of the brain network at the fundamental level. We found that increasing the regional brain connectivity, especially centered on hubs, changed the ES brain network to GS brain network, disrupted the ES condition, and reduced network sensitivity against an external stimulus. Considering a central role of hubs in FM (Kaplan et al., 2019) and giant synchronization cluster formation, changing hub connectivity in the brain network could potentially reverse the key mechanism of ES, facilitating giant synchronization cluster formation in the network, so that in principle, the network is converted from ES (a first-order phase transition) into GS (a second-order phase transition) with reduced brain network sensitivity.

Converting Explosive Synchronization Into General Synchronization in the Brain Network
Previous computational modeling studies have identified network conditions to prevent the onset of ES (Zhang et al., 2016;Dai et al., 2020). The authors have found that ES could be changed to GS if the population of inhibitory connections was larger than a threshold (about 10%) in the networks consisting of inhibitory and excitatory connections. The ratio of the inhibitory connections required for converting ES into GS depended on the network size (the total number of nodes) and the average node degree (the average number of connections among nodes). Such existence of the threshold implies that converting ES into GS is controllable by modulating the population of inhibitory connections. However, it is difficult to directly apply this method to a large-scale human brain network due to its complex hierarchical modular structure and complex dynamics. In addition, it has not been tested in network dynamics near a critical state, presumably reflecting brain network dynamics in a conscious state. Therefore, we instead searched a practical way to convert ES into GS in the brain network, modulating local brain connectivity and local brain signal randomness that can facilitate a giant synchronization cluster formation, that is, reversing the core mechanism of ES. In addition, targeting local brain regions with intense connectivity such as hubs in the brain network is feasible with current noninvasive brain stimulation tools.

Hub as an Effective Target Site
We selected hubs in the brain network as the most promising target site for the brain network modulation because brain network hubs are the most powerful candidate that can control the whole brain network with local modulation. Hubs have been known that they show the largest network controllability in the model (Gu et al., 2015) and play a central role in chronic pain (Kaplan et al., 2019). In addition, hub structure in the brain network plays as a major determinant in global synchronization and desynchronization (Schmidt et al., 2015;Vlasov and Bifone, 2017). Topologically, the hierarchical hub structure in the brain network called "rich club, " is a highly connected and centralized collection of nodes that occupies only 10% of the brain network but facilitates 70% of the communication pathways across the brain regions (van den Heuvel and Sporns, 2011;Schmidt et al., 2015), which can be one of the reasons why hubs are the most effective target sites to promote or suppress global synchronization in the brain network. Recent brain imaging studies have presented that the hub or the membership of the rich club is varied with the intensity of clinical pain. The posterior insula, primary somatosensory, and motor cortices (S1/M1) belonged to the rich club only in FM patients with the highest clinical pain intensity, which was not observed in patients with low pain intensity and healthy controls. The altered hub topology suggests that pronociceptive regions such as insula appear to have acquired new hub status that would influence the quality and quantity of information processing in the brain network (Kaplan et al., 2019). Moreover, the altered hub topology with new hubs may enhance heterogeneity of the brain network, suppressing the onset of global synchronization and changing the brain network closer to the ES network. Therefore, we propose that generating a dominant hub node whose local connectivity is enhanced via CI modulation can successfully convert ES to GS. The enhanced hub connectivity can lead to a gradual synchronization propagation starting from the dominant hub node to other nodes, instead of a sudden singular unification near a critical point.

Increasing Local Structural Connectivity in Brain Network
In this study, we showed that only increasing local structural connectivity of target nodes can mitigate the ES condition with a more negative correlation between node degree and frequency (Figure 3). A positive correlation between node degree and frequency is a well-known network configuration for ES, which suppresses a giant synchronization cluster formation (Gómez-Gardeñes et al., 2011). If higher frequencies were distributed to hubs in a network (i.e., a positive correlation between node degree and frequency), hubs would suppress the overall network synchronization until the network reaches a critical point because higher frequency oscillators have small or irregular amplitudes and incoherent phases in a network (Honey et al., 2007;Gollo et al., 2015;Moon et al., 2015Moon et al., , 2017. On the contrary, if slower frequencies were distributed to hubs (i.e., a negative correlation between node degree and frequency), hubs would tend to increase the local synchronization around them and progressively expand the synchronization throughout the whole network (Moon et al., 2015).
We evaluated how the network sensitivity is significantly changed by a direct brain stimulus after the four types of brain network modulation (increasing or decreasing brain regional connectivity around target nodes and increasing or decreasing randomness of target node dynamics). Only enhancing local structural connectivity around target nodes significantly reduced the brain network sensitivity (Figure 4). In addition, the reduced network sensitivity was highly correlated with the average node degree of target brain regions, which verifies our hypothesis on the central role of the hubs in reducing the network sensitivity. In contrast, decreasing the local connectivity and increasing or decreasing the randomness of node dynamics gave rise to the opposite results because they made the target sites functionally more separated from other nodes, hampering a hub-dominated progressive synchronization growth. Notably, the results imply that despite the importance of hub nodes as effective target sites, the type of modulation is also crucial to converting ES into GS network. These computational modeling results provide a theoretical criterion for developing an effective brain modulation method to reduce the hypersensitivity in FM. Irrespective of brain modulation methods, for instance, pharmacologic, electric, and magnetic stimulation, the outcomes of brain modulation should avoid significant decreases in the hub connectivity and significant change in the hub dynamics (i.e., decreasing/increasing randomness of local brain activity around hubs) that can enhance the ES condition and also the network sensitivity in the brain.

Modulating the Brain Network at Criticality and Its Relationship With Fibromyalgia
We implemented criticality to the brain network model to simulate spontaneous brain activities in a conscious resting state and examined the network sensitivity change between the critical states of the networks before and after the modulation. Brain activities in the conscious resting state share various signal properties of a system near criticality: large variance, large autocorrelation, long-range correlation, high susceptibility, and power law (Beggs and Plenz, 2003;Shew et al., 2011;Tagliazucchi et al., 2016;Kim and Lee, 2019). As far as we know, modulating the brain network at criticality and examining the dependence of the network sensitivity on different types of network modulation, especially for the brain network occurring ES, have not been studied yet. In addition, this computational modeling study and the empirical evidence of ES in FM brain give us insights into the potential application of other important characteristics of ES and its relationship to chronic pain. In general, a system with ES conditions can be characterized by a hysteresis phenomenon, that is, differential pathways between forward and backward state transitions (i.e., transitions from incoherent to synchronized states and vice versa) (Gómez-Gardeñes et al., 2011;Boccaletti et al., 2016). Hysteresis is a universal phenomenon observed in nature and has been investigated in various research fields such as physics, biology, and engineering, and also state transitions in the brain such as sleep and general anesthesia (Hopkinson and Williams, 1912;Chikazumi et al., 1997;Steyn-Ross et al., 2004;Bertotti and Mayergoyz, 2006;Friedman et al., 2010;Rempe et al., 2010;Voss et al., 2012;Joiner et al., 2013). In particular, it has been found in general anesthesia of diverse species (drosophila, murine, mouse, rat, and human) that the anesthetic concentration required for inducing unconsciousness is higher than the concentration where consciousness is regained during general anesthesia (Sepúlveda et al., 2019;McKinstry-Wu et al., 2020;Luppi et al., 2021). In addition, our previous study has demonstrated that human subjects who have larger ES conditions in EEG networks show larger hysteresis in state transitions during the loss and recovery of consciousness , which implies individuals with larger ES conditions take longer time to recover from general anesthesia. Therefore, we expect that the brain with larger ES conditions such as FM brain experiences a large hysteresis phenomenon during state transition, for instance, faster loss of consciousness and more prolonged recovery in anesthesia.
In addition, another representative characteristic of ES at a critical point is the high instability of functional connectivity, which is originated from the large variance of network synchronization at a critical point (Kim and Lee, 2020). In a GS network, the structural connectivity significantly shapes functional connectivity patterns at a critical point with a large correlation Kim et al., 2018). We have mathematically proved the underlying mechanism of this large correlation, empirically confirmed it with brain networks of different species (mouse, monkey, and human), and suggested the large correlation between functional and structural connectivity as an indicator of brain criticality . However, ES network has an aberrant hub network organization because the hub connectivity of the ES network is directly or indirectly disrupted. Considering the essential role of hubs in the organization of higher-order brain functions such as cognition (Vatansever et al., 2020;Zdanovskis et al., 2020;Wang et al., 2021), such aberrant functional brain organization may be associated with the cognitive deficit in FM. However, the detailed association between the cognitive deficit and the aberrant hub structure in the brain remains to be tested.

Implication for Effective Brain Modulation in Fibromyalgia
The identified effective target sites in the human brain network model with the ES condition provide many implications for developing an effective brain modulation method for FM. The simulation results showed that increasing the local structural connectivity centered on hub regions including the insula, left isthmus cingulate cortex, and right precuneus is the most effective way to reduce the brain network sensitivity changing ES to GS. Interestingly, the brain regions such as the insula and the cingulate cortex have been already known as hubs in FM brain and currently take into account as the target sites for brain modulation to reduce the pain intensity (Kaplan et al., 2020). Recently, the motor cortex (M1) has been also considered as a target site for tDCS for pain treatment (Foerster et al., 2015;Cummiford et al., 2016). Here, we suggest the right precuneus, which is a hub region in the cortex, as a feasible target site for noninvasive brain stimulation such as tDCS to treat FM because of the limitation of direct stimulation to subcortical regions.

LIMITATION
This study contains several limitations. First, the Stuart-Landau model has limitations in interpreting the simulation results because the model is not based on a biophysical mechanism. Instead of simulating the realistic brain activities, we focused on identifying a general network principle that could be associated with the emergence of hypersensitivity in the brain networks and effectively reduce the hypersensitivity at a global brain network level. This general network principle is applicable to diverse brain networks with different genders, ages, and diseases. In addition, it would also need to test the robustness of these simulation results with other types of network structures such as multiplex or multilayer networks (De Domenico, 2017;Buldú and Porter, 2018). Second, the target brain sites that we suggested to reduce the brain network sensitivity were selected from an anatomical brain network structure averaged over individuals. Thus, acquiring individual brain network structures would enable us to identify individualized target sites that might result in better performance to convert ES into GS. Third, the modulation we induced can reflect neuroplasticity and cortical excitability changes that are empirically observed in tDCS trials but cannot explicitly show that it changes the excitatory or inhibitory connectivity. Since the excitatory and inhibitory connectivity is important in the cortical network (Soriano et al., 2008;Vogels and Abbott, 2009), further studies will be needed in various types of models including excitatory and inhibitory connections (Chowdhury et al., 2020(Chowdhury et al., , 2021. Fourth, the effects of the four types of brain network modulation are not independent of each other. However, despite the theoretic limitation of separating the intertwined effects, the different types of brain modulation (increasing/decreasing local connectivity or randomness of local node dynamics) provide us measurable aims that can be achieved with the experimental design. We expect that future studies with a fine model setting can address some of these limitations.

CONCLUSION
This computational model study suggests a network mechanismbased brain modulation method that can significantly reduce brain network sensitivity. Increasing the local structural connectivity centered on hubs such as insular, isthmus cingulate cortex, and precuneus fundamentally changes the type of state transition from a first-order synchronization transition (ES) to a second-order synchronization transition (GS), mitigates the ES condition, and reduces network sensitivity of the brain. The ES mechanism-based brain network modulation will provide a theoretic framework to design an effective brain stimulation method to systematically control the hypersensitivity in chronic pain patients.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
UL and MK designed the research and wrote the original draft. MK simulated the model. UL, MK, RH, and AD reviewed and edited the manuscript. All authors have contributed to the interpretation of the results and approved the final manuscript. FUNDING UL, RH, and AD were supported by National Institute of Health (NIH) grant R01AT010060. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.