Stress-induced changes in modular organizations of human brain functional networks.

Humans inevitably go through various stressful events, which initiates a chain of neuroendocrine reactions that may affect brain functions and lead to psychopathological symptoms. Previous studies have shown stress-induced changes in activation of individual brain regions or pairwise inter-regional connectivity. However, it remains unclear how large-scale brain network is reconfigured in response to stress. Using a within-subjects design, we combined the Trier Social Stress Test and graph theoretical method to characterize stress-induced topological alterations of brain functional network. Modularity analysis revealed that the brain network can be divided into frontoparietal, default mode, occipital, subcortical, and central-opercular modules under control and stress conditions, corresponding to several well-known functional systems underpinning cognitive control, self-referential mental processing, visual, salience processing, sensory and motor functions. While the frontoparietal module functioned as a connector module under stress, its within-module connectivity was weakened. The default mode module lost its connector function and its within-module connectivity was enhanced under stress. Moreover, stress altered the capacity to control over information flow in a few regions important for salience processing and self-referential metal processing. Furthermore, there was a trend of negative correlation between modularity and stress response magnitude. These findings demonstrate that acute stress prompts large-scale brain-wide reconfiguration involving multiple functional modules.

Introduction

From taking an examination, giving a public presentation, to almost being hit by a car, human beings are inevitably going through various stressful events in the daily life. These stressors - physical or psychological events that threaten homeostasis in organisms - initiates a stress response, which activates the autonomic nervous system and the hypothalamus-pituitary-adrenal axis, evoking a chain of neuroendocrine reactions that shift the brain into a state of sharpened vigilance and alertness and weakened higher-order cognitive capacity (Arnsten, 2009; van Marle et al., 2009; Qin et al., 2009; Roozendaal et al., 2009). Extensive evidence show that the amygdala plays a key role in evoking and regulating stress response through its efferent projections to a widely distributed set of brain regions (Arnsten, 2009; Roozendaal et al., 2009). For example, the amygdala's intrinsic connectivity (i.e., resting-state functional connectivity) with anterior cingulate cortex (ACC) and anterior insula was enhanced immediately following stress induction (van Marle et al., 2010), suggesting a state of hypervigilance that promotes salience processing. Intrinsic connectivity of amygdala with cortical midline structures including posterior cingulate cortex (PCC), precuneus, and medial prefrontal cortex (mPFC) were found to be enhanced 1 h after stress induction (Veer et al., 2011), suggesting that mPFC is engaged in top-down regulation of stress-induced amygdala activity and facilitation of self-referential mental processing. Enhanced intrinsic connectivity of amygdala with the hippocampus was found immediately after stress induction and this was sustained for as long as 2 h after stress induction (Vaisvaser et al., 2013), suggesting a long-lasting memory consolidation of stressful experiences. These findings together suggest that stress-induced changes involve alterations in multiple brain functional systems.

The human brain is intrinsically organized into distinct interconnected functional networks that orchestrate complex cognitive and affective processes (Sporns, 2011). There are three core neurocognitive systems crucial for cognitive and affective processing (Menon, 2011). The salience network (SN), which includes prominent nodes in the anterior insula, dorsal anterior cingulate cortex (dACC), and the amygdala, plays a key role in detecting salient internal or external events and assigning saliency to them (Menon, 2011; Menon and Uddin, 2010; Seeley et al., 2007); the default mode network (DMN), which includes important nodes in the mPFC and PCC, is critical for self-referential mental processing (Andrews-Hanna et al., 2010; Greicius et al., 2003); the central executive network (CEN), which includes pivotal nodes in the dorsolateral prefrontal cortex (DLPFC) and the posterior parietal cortex (PPC), is important for higher-order cognitive functions such as working memory, attention, and regulation (Dosenbach et al., 2007; Seeley et al., 2007). Deficits in the three networks are associated with a wide range of stress-related psychiatric disorders such as anxiety, depression, and posttraumatic stress disorder (Menon, 2011). It is noteworthy that brain regions affected by acute stress are largely overlapped with pivotal SN, DMN, and CEN nodes, suggesting that effects of acute stress on the brain may extend far beyond individual regions and pairwise inter-regional connectivity and into large-scale brain network reconfiguration.

In fact, Hermans et al. (2011) have reported enhanced responsiveness and interconnectivity within a salience network including a set of cortical (e.g., frontoinsula, dACC) and subcortical (e.g., amygdala, thalamus) regions during exposure to aversive movies. They further proposed that exposure to acute stress prompts large-scale brain network reconfiguration during which resources are reallocated to a salience network at the cost of a central executive network (Hermans et al., 2014). A recent review indirectly tested this hypothesis by evaluating findings from a fair amount of stress studies (van Oort et al., 2017). They found that studies focusing on brain activity or connectivity changes during exposure to acute stress show varying and inconsistent results, with some studies reported increased activity in SN regions (Gianaros et al., 2008; Koric et al., 2012; Seo et al., 2011; Sinha et al., 2004) whereas others reported decreased activity in SN regions (Khalili-Mahani et al., 2009; Pruessner et al., 2008), and some studies reported no changes in activity in CEN regions (Boehringer et al., 2015; Khalili-Mahani et al., 2009; Lord et al., 2012; Pruessner et al., 2008; Seo et al., 2011; Sinha et al., 2004) whereas others reported increased activity in CEN regions (Fechir et al., 2010; Gianaros et al., 2008; Koric et al., 2012). Most of these studies also identified stress-induced changes in DMN, including increased (Boehringer et al., 2015; Fechir et al., 2010; Koric et al., 2012; Lord et al., 2012; Pruessner et al., 2008; Seo et al., 2011; Sinha et al., 2004) or decreased (Albert et al., 2015) activity in DMN regions, increased mPFC (DMN region) connectivity with anterior insula (SN region), and decreased mPFC connectivity with DLPFC (CEN region), PCC and IPL (DMN regions) (Vaisvaser et al., 2016). The discrepancy may be due to differences in experimental design and contrast (van Oort et al., 2017). Moreover, van Oort et al. (2017) found that few studies have investigated the effects of stress on intrinsic connectivity after stress induction. To date, the most consistent findings are increased connectivity between SN regions and DMN regions under stress (van Marle et al., 2010; Maron-Katz et al., 2016; Quaedflieg et al., 2015; Vaisvaser et al., 2013; Veer et al., 2011). Taken together, these findings suggest that stress prompts large-scale changes involving all three core neurocognitive networks and highlight the need to directly examine large-scale network reconfiguration using a method beyond regional activity and pairwise inter-regional connectivity.

In the current study, we aimed to directly investigate how large-scale brain networks are reconfigured in response to stress with the aid of graph theoretical analysis. Specifically, we focused on the period immediately after stress induction. As a natural framework for mathematical representation of complex networks, graph theory provides a unique and powerful way to quantitatively describe the segregation and integration of the brain network from the perspective of topological organization (Bullmore and Sporns, 2009; He and Evans, 2010). In particular, we are interested in modularity which refers to the formation of local modules that nodes in the same module are densely connected while those in different modules are sparsely connected (Newman and Girvan, 2004; Radicchi et al., 2003). The modular organization provided a template to further characterize stress-induced changes in the topology of brain regions and modules. We hypothesized that acute stress prompts large-scale network reconfiguration through altering interactions among modules and key SN, DMN, and CEN regions.

Materials & methods

Participants

Seven graduate and twenty-three undergraduate students participated in the resting state fMRI scanning. Two were excluded because of excessive head movements, leaving 28 participants for data analysis (age: 20.6 ± 1.9 years (18–25 years); 14 females). Participants were instructed to refrain from heavy exercise and caffeine for at least 12 h before the experiment. None of the participants had previously enrolled in stress-related research or had psychology majors. None of them reported any history of smoking, use of illicit substances or psychoactive medications, or psychiatric, neurological, or cognitive disorders. Informed written consent was obtained from each participant before the experiments. The study was carried out in accordance with the Declaration of Helsinki and was approved by the South China Normal University Institutional Review Board.

Design and procedures

After acclimation period of 20 min following participants’ arrival (T1), saliva samples and affect ratings were collected. Then participants were informed that they had 5 min to prepare for the following stress or control task (see stress induction and Fig. 1). After the preparation (T2), affect ratings were collected again. Then participants went through the formal cognitive task, with saliva samples and affect ratings collected upon task completion (T3). Next, participants were required to lie still in the fMRI scanner with eyes closed and resting-state fMRI data were collected. After completion of resting state fMRI scan, saliva samples and affect ratings were collected every 12 min for three times (T4, T5, and T6; Fig. 1). Note, participants received both stress and control manipulations which were separated by at least 30 days, with sequence of manipulation counter-balanced across participants. All experiments were conducted between 1300 h and 1800 h to control diurnal variation of cortisol secretion.

Fig. 1

Experimental procedure and manipulation check. (A) After acclimation period of 20 min following arrival, participants were required to go through the Trier Social Stress Test (TSST) which consists of preparation (5 min) and formal tasks (10 min). The formal tasks were performed either with (stress condition) or without (control condition) social evaluative processes. After the formal tasks, resting-state fMRI were collected. Saliva samples were collected at T1, T3, T4, T5, and T6. Affective ratings were collected at T1, T2, T3, T4, T5, and T6. (B) Cortisol and positive/negative emotional responses under control and stress condition. Compared with control condition, stress induced higher cortisol responses at T3, T4, T5, and T6, as well as lower positive and higher negative emotional responses at T2 and T3. *p < 0.05, **p < 0.01, ***p < 0.001.

Stress induction

Participants were required to go through the Trier Social Stress Test (TSST) (Kirschbaum et al., 1993), which included a preparation period (5 min) and a formal task period (5-minute public speaking and 5-minute mental arithmetic). In the stress condition, participants were instructed to prepare for a job interview during which they had 5 min to introduce themselves to the selection committee in the form of public speaking. They were required to convince the committee that they were the most suitable candidates for the job. To increase task engagement, participants were asked to write down their dream job before the preparation period. After the preparation, they delivered a 5-minute speech (without notes) in front of a video camera and the selected committee (one woman and one man) who were trained to remain neutral and unresponsive. Unbeknown to participants, upon completion of the speech, they would be required to complete a 5-minute mental arithmetic task, that is, serially subtracting the number 13 from 1022 in English as fast and accurately as possible in front of the same committee. On every failure, the committee asked the participant to restart from 1022. In the control condition, to ensure a comparable cognitive load, participants went through the same tasks with less social evaluative processes (i.e., no video camera and no committee). The instruction to start over for math and number of observers were the active differences between groups. Admittedly, although the control condition is not stress-free, we expected the experimental condition to induce higher level of stress than the control condition.

Endocrine and psychological measures

Salivary cortisol was assessed at multiple time point throughout the procedure (Fig. 1). Saliva samples were collected with the Salivettes (Sarstedt, Germany) and were stored at −15 °C until assayed. Cortisol concentrations in saliva (in ng/mL) were measured by ELISA (DRG, Germany). Three participants’ cortisol cannot be assayed due to insufficient saliva, leaving 25 participants for data analysis. Positive (i.e., calm, relaxed, peaceful, confident, and energetic) and negative (i.e., nervous, anxious, scared, tired, and upset) emotion were rated on a four-point scale (1 = not at all; 4 = extremely) at multiple time points throughout the procedure (Fig. 1).

Neuroimaging data acquisition

Images were acquired using Siemens Trio 3T scanner at the Brain Imaging Center of South China Normal University. T2*-weighted echo-planar images (EPI) were obtained with blood oxygenation level-dependent (BOLD) contrast. Thirty-two slices of 4 mm thickness (with no gap in-between) that covered the whole brain were acquired in an inferior to superior direction in interleaved order (repetition time = 2000 ms, echo time = 25 ms, field of view = 220 × 220 mm2, flip angle = 85°, matrix size = 64 × 64). During the data acquisition, participants were required to lie still in the scanner with their eyes closed. Two hundred and forty images were acquired for the resting state scan. T1-weighted images were acquired with the following parameters: repetition time = 1900 ms, echo time = 2.52 ms, field of view = 256 × 256 mm2, flip angle = 9°, matrix size = 256 × 256.

Neuroimaging data preprocessing

The fMRI data were preprocessed and analyzed using Statistical Parametric Mapping software SPM12 (http://www.fil.lon.ucl.ac.uk/software/SPM12/; Wellcome Trust Department of Cognitive Neurology, London, UK), resting-state fMRI data analysis toolkit (REST) (Song et al., 2011), and Data Processing Assistant for Resting-State fMRI (DPARSF) (Yan and Zang, 2010). After discarding the first 10 vol, the remaining fMRI volumes were first slice-time corrected and then motion corrected using a least squares approach and a 24-parameter autoregressive model (Friston 24-parameter model) (Friston et al., 1996). Nuisance signal corrections were also performed on the 24 head-motion parameters, cerebrospinal fluid (CSF), white matter and global signals (Fox et al., 2009; Power et al., 2014). Next, corrected functional data were normalized to MNI space, re-sampled to 3 mm × 3 mm × 3 mm isotropic voxel, spatially smoothed with a 4 mm FWHM Gaussian filter, and temporally band-pass filtered into 0.01–0.1 Hz to reduce the effect of very low frequency drift and the high frequency physiological noise. In addition to the abovementioned motion correction, we further addressed the residual effects of head motion in group analysis through frame-wise displacement (FD) derived from Jenkinson's relative root mean square algorithm as nuisance covariate (Jenkinson et al., 2002). Participants with mean Jenkinson FD greater than group mean motion + 2 × S.D. were excluded for following graph theoretical analysis (Yan et al., 2013a). The mean FD in stress and control conditions were 0.113 ± 0.027 and 0.110 ± 0.027, respectively. Paired t-test showed no difference in FD between the two conditions (p = 0.34).

Network construction

To construct the functional connectivity matrix, we employed Automated Anatomical Labeling (AAL) template (Tzourio-Mazoyer et al., 2002) (Table S1) to parcellate the brain into 90 regions of interest (ROIs). Representative time series of each ROI was then obtained by averaging the time series of each voxel within that region. Pearson correlation coefficient, which stands for the functional connectivity strength, was calculated for all pairs of ROIs. A Fisher's-Z transformation was then applied to the correlation matrices to improve the normality and to further reduce relationships between motion and connectivity strength (Yan et al., 2013b). Given the ongoing debate on physiological meaning of negative correlations (Anderson et al., 2011; Chang and Glover, 2009), we set all negative correlations to zero. These matrices (one for each participant) served as the input for modularity analysis and the corresponding results were reported as main findings. However, we performed additional validation analysis with negative correlations preserved and found the main findings intact (see supplementary materials).

Modularity analysis

In graph theoretical analysis, modularity reflects how well the network can be delineated into communities and modularity analysis aims to find a specific partition that yields the maximum modularity. Modularity is one of the most fundamental and intriguing properties and has been repeatedly revealed in both functional and structural brain networks. The maximum modularity value Q quantifies how well a partition differentiates subsets of tightly connected nodes. Q of a network G with N nodes for a given partition is defined as the average difference between present and chance-expected within-module connection weights (Newman, 2004, 2006; Rubinov and Sporns, 2011):where is the connection weight between node i and j, is the sum of all weights in the network, computed as , is the strength of node i, computed as , is the module containing node i, and is 1 if node i and j are in the same module and 0 otherwise to ensure only intra-modular edges are added to the sum. Q value ranges from 0 to 1, with ~1 corresponding to a strongly modular network whereas ~0 corresponding to a non-modular organization. A high modularity value indicates more balanced topology between brain functional segregation and integration while conserving wiring cost and high resilience to network node or edge damages (Bullmore and Sporns, 2009). In previous studies, Q ≥ 0.3 is generally accepted as a criterion indicating non-random community structure (Newman and Girvan, 2004). Here, modularity analysis was conducted with Brain Connectivity Toolbox (Rubinov and Sporns, 2010), in which a modified greedy optimization algorithm (Clauset et al., 2004) was adopted for community detection.

Prior to modularity analysis, a sparsity thresholding approach was applied to the obtained individual functional brain networks to ensure that the wiring cost of each participant was equivalent. For a given network, sparsity is defined as the ratio of the number of existing edges to the maximum possible edge number in the network. Previous studies have indicated that a high sparsity threshold would lead to low-modularity graphs that is equivalent to a random graph whereas a low sparsity threshold would lead to disconnected graphs that contain isolated brain regions (Meunier et al., 2009). Therefore, we estimated modular structure over a sparsity range of 10%–40% with 1% stepwise for each participant. To avoid multiple comparisons at individual sparsity threshold and to reduce the dependency of any significant topological differences on the arbitrary choice of a single sparsity, the AUC of Q values within the modular structure regime of each participant were fed into statistical analysis to quantify modularity difference between the control and the stress conditions.

Modular organization and regional role assignment

To characterize modular organization in each condition, fine-grained analysis was further performed on functional brain networks averaged across individuals for each condition and thresholded at a sparsity of 13%. This sparsity allows us to capture the network backbone underlying the modular organization and maintain the network fully-connected (the ability for each node to reach other nodes in the network) in most participants (Meunier et al., 2009; Wu et al., 2012). Once obtaining a modular organization, we can further assess topological roles of nodes depending on their intra- and inter-modular connectivity measured by the within-module betweenness centrality (sBCi) and participation coefficient (PCi) (Guimerà and Amaral, 2005; Wu et al., 2012).

Betweenness centrality (bci) of a node i is defined as the number of shortest paths between pairs of other nodes that pass through node i (Freeman, 1977). The sBCi is the betweenness centrality obtained within the module where node i belongs. Therefore, sBCi quantifies the influence of a node over information flow between other nodes within the same module.

Betweenness centrality (bci) of a node i is defined as the number of shortest paths between pairs of other nodes that pass through node i (Freeman, 1977). The sBCi is the betweenness centrality obtained within the module where node i belongs. Therefore, sBCi quantifies the influence of a node over information flow between other nodes within the same module. PCi measures the inter-modular connectivity of a node i:where is the total number of edges linking to node i, is the number of edges linking node i to other nodes in the module s, and is the total number of modules. If node i has a homogeneous connection distribution with all modules, PCi will be close to 1; otherwise, if node i links exclusively to other nodes within the same module, PCi will be close to 0.

In terms of sBCi and PCi, a node could be characterized as connector hub (R1; sBC > mean + S.D. and PC > 0.25), connector non-hub (R2; sBC mean +S.D. and PC > 0.25), provincial hub (R3; sBC > mean + S.D. and PC 0.25), or provincial non-hub (R4; sBC mean +S.D. and PC 0.25) (He et al., 2009; van den Heuvel and Sporns, 2013; Wu et al., 2012). Based on topological roles of nodes, we further characterized topological roles of modules and investigated their differences between conditions. Specifically, we calculated connector coefficient for each module, which refers to the ratio between the number of connector nodes (R1 + R2) and the total number of nodes (R1 + R2 + R3 + R4) within the module. If a module has a high connector coefficient (>0.6) and a high ratio of inter-modular connections (>1/number of modules), it could be defined as connector module (Meunier et al., 2009; Wu et al., 2012). Finally, we calculated intra-modular connectivity strength (i.e., the sum of all connectional weights within a module) and inter-modular connectivity strength (i.e., the sum of connection weights between two modules) for each individual and each condition to assess differences in module significance between conditions (Shin et al., 2014).

Statistical analysis

Repeated-measures ANOVA and two-tailed paired-samples t-test were used to assess differences between conditions. Correlation analysis was conducted to examine whether stress-induced changes in brain measures (i.e., modularity, intra- and inter-module connectivity) are associated with individual difference in stress response and recovery as well as emotional response. Specifically, for each condition, T4 (the point when cortisol reached its maximum under stress; see Fig. 1B) and T6 (the point when cortisol decreased to baseline level under stress) cortisol was baseline corrected by subtracting corresponding T1 cortisol. Then cortisol difference in T4 between stress and control conditions was used as a measure of stress response whereas cortisol difference in T6 – T4 (i.e., habituation from cortisol peak) between stress and control conditions was used as a measure of stress recovery. Similarly, for each condition, T3 (the point when emotional response reached its peak) emotional rating was baseline corrected by subtracting corresponding T1 rating. The positive or negative rating difference in T3 between stress and control conditions was used as a measure of positive or negative emotional response. Family Discovery Rate (FDR) correction was used to correct for multiple comparisons. Results with FDR corrected p < 0.05 were considered as significant. In addition, for significant findings, we further examined if they were affected by experimental randomization order using multiple regression analysis.

Results

Endocrine and psychological measures

To examine time-dependent cortisol variation during stress and control conditions, we conducted repeated-measures ANOVA on cortisol with Manipulation (Control vs. Stress) and Time (T1, T3, T4, T5, and T6) as within-subjects factors. Results revealed a significant main effect of manipulation, F (1, 24) = 13.14, p = 0.001, and a significant main effect of time, F (4, 96) = 3.63, p = 0.008 (Fig. 1B). Manipulation by Time interaction was also significant, F (4, 96) = 3.40, p = 0.012. Further pairwise comparison revealed higher cortisol in stress condition than that in control condition at T3 (F (1, 24) = 4.74, p = 0.04), T4 (F (1, 24) = 20.01, p < 0.001), T5 (F (1, 24) = 17.17, p < 0.001), and T6 (F (1, 24) = 12.62, p = 0.002). The intra- and inter-assay coefficients of variation (CVs) for cortisol analysis were 3% and 5%, respectively.

To examine time-dependent emotion variation during stress and control conditions, we conducted repeated-measures ANOVAs on positive and negative emotion ratings separately, with Manipulation (Control vs. Stress) and Time (T1, T2, T3, T4, T5, and T6) as within-subjects factors. For positive emotion ratings, results revealed a marginally significant main effect of manipulation, F (1, 27) = 3.73, p = 0.064, and a significant main effect of time, F (5, 135) = 6.76, p < 0.001 (Fig. 1B). Manipulation by Time interaction was also significant, F (5, 135) = 8.71, p < 0.001. Further pairwise comparison revealed lower positive emotion ratings in stress condition than that in control condition at T2 (F (1, 27) = 8.80, p = 0.006) and T3 (F (1, 27) = 24.55, p < 0.001). For negative emotion ratings, results revealed a significant main effect of manipulation, F (1, 27) = 7.20, p = 0.012, and a significant main effect of time, F (5, 135) = 3.50, p = 0.005. Manipulation by Time interaction was also significant, F (5, 135) = 7.83, p < 0.001. Pairwise comparison revealed higher negative emotion ratings in stress condition than that in control condition at T2 (F (1, 27) = 17.92, p < 0.001) and T3 (F (1, 27) = 21.52, p < 0.001).

Modularity and modular reorganization

For both conditions, maximum modularity decreased monotonically with increase in sparsity, and the functional brain networks consistently exhibited a community structure (Q 0.3) over the predefined sparsity range. No condition difference was found in modularity. Similar results were found when negative correlations were preserved.

Moreover, the fine-grained analysis, which was conducted based on the functional brain network averaged across individuals for each condition and thresholded at a sparsity of 13%, revealed 5 connected modules for both conditions. Specifically, in the control condition, the 5 connected modules were designated as the frontoparietal module (Module I, 20 regions), the default mode module (Module II, 22 regions), the occipital module (Module III, 14 regions), the subcortical module (Module IV, 14 regions), and the central-opercular module (Module V, 20 regions) (Fig. 2A). In the stress condition, the 5 modules were similar to those revealed in the control condition (Fig. 2B), with a few differences in relative size and topological role of brain regions (see Fig. 2C and D, Table 1, and Table S2; modular organization derived from analysis with negative correlations preserved can be found in supplementary results). Specifically, six regions were identified as connector-hub under control condition, including SN regions right amygdala, right insula, and left putamen as well as DMN regions right SFGmed and bilateral PCG. In contrast, only right frontoinsula area (SN regions) were identified as connector-hub under stress. Moreover, bilateral inferior orbitofrontal gyrus switched allegiance from the frontoparietal module to default mode module and bilateral precuneus switched allegiance from the default mode module to occipital module.

Fig. 2

Modular structure of condition-average brain functional network for both control and stress condition at sparsity level of 13%. Five connected modules were identified in both control (A) and stress (B) conditions, including the frontoparietal module (Module I), the default mode module (Module II), the occipital module (Module III), the subcortical module (Module IV), and the central-opercular module (Module V). Topological roles of brain regions (i.e., connector hub, connector non-hub, provincial hub, or provincial non-hub) in each module were presented for control (C) and stress (D) conditions. Shape stands for connector (square) or provincial (circle) nodes whereas size stands for hub (large) or non-hub (small) nodes.

Table 1

The distribution of connector nodes and the topological roles of modules in the control and stress conditions.

Module	Control			Stress
	Regions	Connectors	Inter-modular	Regions	Connectors	Inter-modular
I (frontoparietal)	20	12 (0.60)	60 (0.22)	18	11 (0.61)	56 (0.21)
II (default mode)	22	16 (0.73)	75 (0.28)	23	13 (0.57)	66 (0.25)
III (occipital)	14	4 (0.29)	21 (0.08)	16	6 (0.38)	29 (0.11)
IV (subcortical)	14	6 (0.43)	42 (0.16)	14	7 (0.50)	39 (0.15)
V (central-opercular)	20	10 (0.50)	70 (0.26)	19	12 (0.63)	72 (0.27)
Total	90	48	134	90	49	131

Note: Module I – V stand for the frontoparietal, default mode, occipital, subcortical, and central-opercular module, respectively. Region column indicates the number of regions in each module. Connector column indicates the number of connector nodes in each module and the ratio between the number of connector nodes and the total number of nodes in each module. Inter-modular column indicates the number of connections between a module and the other four modules as well as the ratio between a module's inter-modular connections and the total number of inter-modular connections across all modules. Bold number indicates connector module defined as a module with high connector coefficient (>0.6) and high inter-modular connections (>1/number of modules).

Furthermore, the frontoparietal module was identified as connector module for both conditions as it had a high connector coefficient (>0.6) and a high ratio of inter-modular connections (>1/number of modules; Table 1). The default mode module was identified as connector module only in the control condition and central-opercular module was identified as connector module only in the stress condition (Table 1). In addition, compared with control condition, intra-frontoparietal module connectivity was significantly weakened under stress (t (27) = −3.96, p = 0.001, FDR corrected) while intra-default-mode module connectivity was significantly enhanced under stress (t (27) = 4.58, p < 0.001, FDR corrected). Similar results were found after controlling for randomization order and no significant effect of order or order-related interaction were found.

Brain-behavioral associations

Although we did not find condition difference in modularity, we found a trend of negative correlation between stress-induced changes in modularity and individual differences in stress response (t (24) = −1.80, p = 0.084), suggesting that the stronger the stress response, the more severe the balance between brain functional segregation and integration was disrupted. We did not find any other significant associations between stress-induced changes in brain measures and individual differences in stress response. Moreover, we did not find any significant associations between stress-induced changes in brain measures and individual differences in stress recovery. In addition, we found that stress-induced changes in intra-occipital module connectivity negatively correlated with individual differences in positive emotional response (t (26) = −2.87, p = 0.040, FDR corrected) and positively correlated with individual differences in negative emotional response (t (26) = 3.71, p = 0.005, FDR corrected). Similar results were found after controlling for randomization order or when negative correlations were preserved, and no significant effect of order or order-related interaction was found.

Discussion

Combining the TSST and modularity analysis, we demonstrated significant endocrine and psychological stress responses (Campbell and Ehlert, 2012) immediately after stress induction, along with brain-wide network reconfiguration involving multiple functional systems.

The optimal modular architecture was preserved under stress in healthy individuals but there was a trend of negative correlation between modularity and stress response magnitude, suggesting that the stronger the stress response the more disturbance in the balance between brain functional segregation and integration. Moreover, we found five distinct modules encompassing frontoparietal, medial prefrontal and posterior cingulate cortices, occipital, subcortical, and sensorimotor cortices under both control and stress conditions. These modules corresponded to several well-known functional systems (He et al., 2009; Meunier et al., 2009), with frontoparietal module largely overlapping with the CEN, default mode module largely overlapping with the DMN, and subcortical module largely overlapping with the SN. Among these modules, the frontoparietal and default mode modules were further identified as connector module under control condition because of the high proportion of connector nodes and inter-modular connections, consistent with a recent meta-analysis study which aggregated data from 9208 experiments and revealed that connector nodes mainly located in the frontal and parietal regions in which brain activity were associated with many different cognitive functions (Bertolero et al., 2015). In fact, brain regions within the frontoparietal cortices (or CEN) can flexibly and rapidly adapt to ongoing tasks through shifting their patterns of interaction with other functional systems (Cole et al., 2013). Similarly, brain regions within the medial prefrontal and posterior cingulate cortices (or DMN) were proved to constitute of connector hubs that link all major structural models (Raichle, 2015). Thus, our findings provide further evidence to support the notion that the CEN and the DMN play a crucial role in coordinating information communication and integration across distributed regions and systems.

Although the frontoparietal module remained as connector module under stress, interconnectivity within this module was weakened in response to stress, which may indicate stress-induced disturbance in the CEN and suggest dampened higher-level cognitive functions under stress (Hermans et al., 2014; Menon, 2011). In fact, Hermans et al. (2014) has proposed that stress prompts a reallocation of resources to the salience network at the cost of the executive control network. However, most stress studies have shown no changes in the CEN (van Oort et al., 2017). Our findings provide new evidence to support a large-scale network reconfiguration involving the CEN.

Moreover, the default mode module lost its role of connector module and interconnectivity within this module was enhanced under stress, suggesting that the DMN became more segregated and its communication with other functional systems was hindered by stress. In fact, most stress studies have identified stress-induced changes in DMN regional activity (Albert et al., 2015; Boehringer et al., 2015; Fechir et al., 2010; Koric et al., 2012; Lord et al., 2012; Pruessner et al., 2008; Seo et al., 2011; Sinha et al., 2004) or pairwise connectivity between DMN region and other brain regions (van Marle et al., 2010; Maron-Katz et al., 2016; Quaedflieg et al., 2015; Vaisvaser et al., 2013; Veer et al., 2011). We speculate that enhanced connectivity within DMN may reflect involvement of neural processes supporting self-referential mental processes related to preceding stressful experience (Fransson, 2005; Vaisvaser et al., 2013). While our findings are consistent with previous studies in that stress elicited changes in the DMN, we provide novel evidence demonstrating that stress altered DMN topological profile, which is far beyond regional activity and pairwise connectivity.

Furthermore, connectivity between midline posterior DMN nodes and other DMN nodes was interrupted by stress. Specifically, precuneus switched allegiance from the DMN to occipital module, and PCC, which was believed to directly interact with all other DMN nodes (de Pasquale et al., 2012; Fransson and Marrelec, 2008), lost its hub position within the module under stress. Weakened connectivity between PCC/precuneus and other DMN nodes has been previously reported in patients with posttraumatic stress disorder (Bluhm et al., 2009), and dysfunction in the DMN has been consistently identified as a major feature of mental disorders related to altered episodic/autobiographical memory (Menon, 2011). Given that the PCC/Precuneus is anatomically interconnected with the medial temporal memory system (Buckner et al., 2008), we speculate that disconnect between PCC/Precuneus and other DMN regions may affect memory integration/update of stressful experience. Additionally, bilateral inferior orbitofrontal gyrus (ORBinf) switched allegiance from the frontoparietal module to default mode module, suggesting increased connectivity between ORBinf and other DMN nodes under stress. This finding may indicate enhanced involvement of ORBinf in evaluating and regulating internal representation of stressful experience, given that the inferior orbitofrontal area has been found playing an important role in appraisal of internally generated events (Dixon et al., 2017) and regulation of stress and emotion (Dedovic et al., 2009).

In addition, previous stress studies have consistently identified stress-induced changes in the SN, including increased (Gianaros et al., 2008; Koric et al., 2012; Seo et al., 2011; Sinha et al., 2004) or decreased (Khalili-Mahani et al., 2009; Pruessner et al., 2008) activity in SN regions, and increased connectivity within SN (Hermans et al., 2011) or between SN and DMN regions (van Marle et al., 2010; Maron-Katz et al., 2016; Quaedflieg et al., 2015; Vaisvaser et al., 2013; Veer et al., 2011). While we did not identify stress-induced changes in subcortical module or its connectivity with other modules, we found topological changes in pivotal SN regions right amygdala and bilateral ACC under stress. Specifically, amygdala lost its central place within the subcortical module whereas bilateral ACC lost its capacity of communication with other modules, suggesting stress-induced alterations in SN regional control over information flow in the brain network.

We also found stress-induced changes in the central-opercular and occipital modules associated with sensorimotor and visual systems, respectively. The central-opercular became a connector module under stress, indicating enhanced capacity of sensorimotor system in communication with other functional systems. Previous studies have reported enhanced activation in sensorimotor and visual areas in stressed participants (Soares et al., 2013). We speculate that the alteration in sensorimotor topology may suggest a hypersensitized perception-action system, which was crucial to making fight-or-flight reaction in response to stress. Moreover, stress-induced changes in intra-occipital connectivity were correlated with stress-induced changes in emotional responses. This is not surprising given that occipital lobe is important for imagery of emotional events via retrieving and modifying sensory information from memory system, which can activate the amygdala and autonomic nervous system and result in physiological changes (Kosslyn et al., 2001).

Finally, two issues should be taken into consideration while interpreting the current results. First, the automated anatomical labeling parcellation was used for the network construction. The limitation of anatomical parcellation is that the anatomical boundaries of nodes may not match well functional boundaries and that a large anatomical region might contain heterogeneous resting state fMRI signals, which further influence network construction (Smith et al., 2011). Functional parcellation approaches may be a better option. However, the widely used functional parcellation templates (Power et al., 2011; Yeo et al., 2011) do not include subcortical regions which play a crucial role in stress-related responses. Given this restriction, we finally chose the AAL parcellation which is widely used and proved valid in the literature. Second, the aftermath of stress is not a homogeneous period. Dynamic interactions between autonomic nervous system and hypothalamus-pituitary-adrenal axis lead to dynamic (time-dependent) changes at the brain level (Hermans et al., 2014). However, previous stress studies and the current study characterized brain functional organization in a static way that assume constant functional connectivity between distinct brain regions throughout resting-stats period or a specific task condition. Actually, a few recent studies have demonstrated distinct dynamic interactions between large-scale brain networks (Chen et al., 2016; Taghia et al., 2017, 2018), or significant changes in dynamic properties of functional connectivity in clinical population (Demirtaş et al., 2016; Kaiser et al., 2016; Wise et al., 2017). Future stress studies are needed to elucidate the dynamic neural mechanism underlying stress response and recovery.

In conclusion, our study revealed a five-module architecture of the brain functional network under stress, which corresponded to several well-known functional systems underpinning cognitive control, self-referential mental processing, visual, salience processing, sensory and motor functions. Acute stress weakened the intra-frontoparietal module connectivity and strengthened the intra-default-mode module connectivity, and the stress response tended to be negatively correlated with modularity. These findings demonstrate that acute stress prompts large-scale brain-wide reconfiguration involving the frontoparietal, default mode, occipital, subcortical, and sensorimotor networks.