Reduced spatiotemporal brain dynamics are associated with increased depressive symptoms after a relationship breakup.

Depressive symptoms following a stressful life event, such as a relationship breakup, are common, and constitute a potent risk factor for the onset of a major depressive episode. Resting-state neuroimaging studies have increasingly identified abnormal whole-brain communication in patients with depression, but it is currently unclear whether depressive symptoms in individuals without a clinical diagnosis have reliable neural underpinnings. We investigated to what extent the severity of depressive symptoms in a non-clinical sample was associated with imbalances in the complex dynamics of the brain during rest. To this end, a novel intrinsic ignition approach was applied to resting-state neuroimaging data from sixty-nine participants with varying degrees of depressive symptoms following a relationship breakup. Ignition-based measures of integration, hierarchy, and metastability were calculated for each participant, revealing a negative correlation between these measures and depressive ratings. We found that the severity of depressive symptoms was associated with deficits in the brain's capacity to globally integrate and process information over time. Furthermore, we found that increased depressive symptoms were associated with reduced spatial diversity (i.e., hierarchy) and reduced temporal variability (i.e., metastability) in the functional organization of the brain. These findings suggest the merit of investigating constrained dynamical complexity as it is sensitive to the level of depressive symptoms even in a non-clinical sample.

Introduction

The experience of a stressful life event, such as a relationship breakup, can lead to the development of depressive symptoms (Field et al., 2009, Kendler et al., 1999, Najib et al., 2004, Verhallen et al., 2019). Even in the non-clinical population, these can have serious effects on the quality of life (World Health Organization, 2017), and increase the risk for the onset of a major depressive disorder (MDD) (Cuijpers and Smit, 2004, Karsten et al., 2011). The main method to understand the neural basis of depression involves highlighting the differences between healthy and patient populations. Therefore, it remains unclear whether the depressive symptoms in individuals without a clinical diagnosis have reliable neural underpinnings. In this study, we investigated whether individual differences in the severity of depressive symptoms after a relationship breakup were associated with changes in the resting-state whole-brain dynamics.

Resting-state functional magnetic resonance imaging (fMRI) has been increasingly applied to investigate alterations in the intrinsic form of brain functional connectivity across discrete neural systems. Substantial evidence suggests that depressive symptoms may evolve due to imbalanced communication of large-scale brain networks (for a review see Mulders et al., 2015), and some of these alterations have also been reported in mild levels of depression. For example, previous findings reported reduced functional connectivity of the cognitive control network, suggesting deficits in cognitive processing and cognitive biases in subclinical depression (Hwang et al., 2015, Schultz et al., 2019, Wei et al., 2014). Similarly, abnormal communication between the frontoparietal and the default network may underlie increased attention to internal thoughts and ongoing rumination, which are commonly linked to increased depressive symptoms (Grimm et al., 2009, Kaiser et al., 2015). Likewise, weaker functional dominance in the dorsal attention network may reflect less goal-directed attention to the external world and may contribute to sadness and depressive symptoms (Petrican et al., 2015). Alterations involving the frontolimbic network, which plays a key role in emotion regulation and memory, have also been linked to the poor cognitive performance observed in depressive disorders (Qi et al., 2018). Together, these studies indicate that abnormal functional network connectivity during rest may contribute to altered cognitive and affective functioning, considered to underlie depressive symptoms. However, it is still unclear how exactly these network alterations relate to the functional deficits at the global level.

Advances in the field of connectomes have provided further insight into the mechanisms implicated in global neural communication and how these mechanisms shape brain function in health and disease. An essential property for brain function involves the integration of neural information across specialized brain regions within distributed networks (Lord et al., 2017). A widely used approach to explore integration relies on graph theoretical analysis. In graph theory, topological integration relates to the efficiency in which distributed information is incorporated to enable global communication (Rubinov and Sporns, 2010). In the context of depression, numerous studies have reported changes in different properties of the network, such as increased characteristic path length (i.e., it takes more steps to go from a brain region to any other) and reduced global efficiency (i.e., less efficient exchange of information at the global level). Another relevant property of the neural networks, which is also involved in depression (Gong and He, 2015), is the so-called rich club. This is a set of highly interconnected regions (van den Heuvel and Sporns, 2011) whose centralization of connectivity plays an important role in information integration (Senden et al., 2014, van den Heuvel et al., 2012).

While neuroimaging research on network integration using graph theory, has yielded valuable insights into the topological aspects of brain organization, it might also be beneficial to explore the temporal variability underlying information integration. Accumulating evidence indicates that the levels of integration dynamically change over time (Cohen and D’Esposito, 2016, Deco et al., 2015, Lord et al., 2017, Shine et al., 2016). Furthermore, there is increasing awareness of the relevance of these time-varying properties in clinical settings (Ponce-Alvarez et al., 2015) as they potentially play a central role in supporting cognitive function (Allen et al., 2014, Bassett et al., 2011, Calhoun et al., 2014, Hutchison et al., 2013). Indeed, patients with depression demonstrated aberrant dynamic functional connectivity, particularly among systems involved in attention regulation and self-referential thinking such as the frontoparietal network and the default-mode network (Demirtaş et al., 2016, Kaiser et al., 2016, Wise et al., 2017). Moreover, recent findings have indicated alterations in the relative expression of specific functional networks over time. For example, Zhi et al. (2018) showed that patients with depression engaged more frequently in a state characterized by weaker between- and within-network connectivity. Another study by Figueroa et al. (2019) demonstrated that remitted patients with recurrent MDD had more difficulties to engage in a state involving the cognitive control system.

As a result, several methods have been developed to characterize the spatiotemporal unfolding of whole-brain connectivity (Allen et al., 2014, Deco et al., 2017a, Hansen et al., 2015), including the novel intrinsic ignition method described by Deco and Kringelbach (2017). The advantage of this local-to-global approach is that it describes whole-brain communication in terms of how the broadcast of information spreads across the brain over time. Specifically, the concept of intrinsic ignition measures the contribution of each brain region to propagate neural activity to other regions inducing varying degrees of global integration.

In the present study, we set out to investigate the dynamical complexity of the brain at rest by applying the intrinsic ignition framework to a dataset of 69 participants with varying degrees of depressive symptoms following a relationship breakup. We hypothesized that greater levels of self-reported depressive symptoms are associated with reduced global integration and reduced spatiotemporal variability in the functional organization of the brain.

Materials and methods

Participants

We analyzed a dataset of healthy participants (without a disorder diagnosis) who had recently experienced a relationship breakup. Recent work by colleagues in our lab (Verhallen et al., 2019) on this dataset revealed varying degrees of self-reported depressive symptoms among these individuals. A full description of the participants, acquisition, and analysis of the Major (ICD-10) Depression Inventory (MDI) is provided in Verhallen et al. (2019).

Recruitment

Seventy-one participants who went through a recent relationship breakup participated in the study. Participants were recruited by email, online forums, word of mouth, and local advertisement with posters and flyers. Participants aged 18–26 years were eligible for enrolment if they had a relationship for at least six months and ended in the preceding six months, to ensure variability in depressive symptoms, in line with other studies in the field (Aron et al., 2005, Fisher et al., 2010, Lewandowski and Bizzoco, 2007, Stoessel et al., 2011). Additional eligibility criteria were heterosexuality, right-handedness, a Western background, and oral contraceptive use in women. Participants were excluded if they had any MRI contraindications or any neurological or psychiatric disorder. Participants received a full explanation of the study procedure, gave written informed consent, and received a financial recompense for their participation. The present study was approved by the local ethics committee of the University Medical Center Groningen and was in accordance with the Declaration of Helsinki.

Final participants included in the analysis

Two participants were excluded from the analysis after neuroimaging data pre-processing, as we will explain in Section 2.3.1. Pre-processing. The considered data consisted of 37 women and 32 men, between 18 and 25 years old (mean = 22).

Depressive symptom scores

The presence and severity of depressive symptoms were assessed with the MDI. The MDI is a brief self-report mood questionnaire based on the DSM-IV and the ICD-10 symptoms of depression (Bech et al., 2001, Olsen et al., 2003). This inventory can be used as a continuous scale, providing an estimate of symptom severity. Moreover, it is considered a reliable tool for assessing depressive symptoms in the general population (Cuijpers et al., 2007). The MDI consists of 12 questions asking how frequently, ranging from never (1) to very often (5), participants have experienced symptoms related to depression over the last 14 days. In our sample, MDI scores (mean = 14.3, IQR = 7–21) varied from a minimum of 1 to a maximum of 45 (out of a possible 50). This indicates that there is a wide range of depressive symptoms in our population sample. Supplementary Fig. S1A shows the distribution of MDI scores.

Imaging acquisition

Neuroimaging data were acquired on a 3 Tesla Philips Intera MRI scanner (Philips Medical Systems, Best, the Netherlands) with a 32-channel SENSE head coil. Structural images were collected using a T1-Turbo Field Echo sequence with the following parameters: 0 gap; 1 mm slice thickness; repetition time (TR) = 9000 ms; echo time (TE) = 3.59 ms; field of view (FOV) in mm = 256 × 170 × 232 (RL × FH × AP); voxel size in mm = 1 × 1 × 1; 8° flip angle. Functional images were acquired with a T2* weighted fast field echo-planar imaging (EPI) pulse sequence with 0.3 mm gap; 3.5 mm slice thickness; TR = 2000 ms; TE = 29 ms; FOV in mm = 220 × 220 × 121.8 (RL × FH × AP); voxel size in mm = 3.44 × 3.44 × 3.30; 70° flip angle. Scanning time = 308 s. A total of 150 volumes (37 slices per volume) were collected for each participant. Slices were oriented parallel to the AC-PC transverse plane and acquired in descending order. Participants used earplugs to minimize the noise of the scanner and foam pads to reduce head motion. They were instructed to hold their head still, close their eyes, let their mind flow, and not fall asleep.

Imaging analysis

Image data analysis was performed using Statistical Parametric Mapping, SPM12 (http://www.fil.ion.ucl.ac.uk/spm), implemented in Matlab 2015b (The MathWorks Inc., Natick, MA).

Pre-processing

Functional images were realigned to the first image using rigid body transformations and the mean EPI image, created during the former step, was co-registered to the anatomical T1 image. Then, images were spatially normalized to MNI T1-template and resampled to 2 × 2 × 2 mm voxel size. Bounding box parameters were changed to −90:90, −126:90, −72:108 to ascertain overlap of all regions of interest (ROIs) with the bounding box for extraction of time series data. Subsequently, a set of eight nuisance variables and their first-order temporal derivatives were regressed out: six head motion parameters (three translation and three rotation), white matter (WM) signal, and cerebrospinal fluid (CSF) signal. The extraction of these two signals involved the creation of a WM and CSF masks by segmenting the T1-weighted image and then extracting the first eigenvariate from the time series of the overlapping voxels. Smoothing was applied using 8 mm full width at half maximum Gaussian kernel. After data pre-processing, two participants were excluded; one due to excessive head motion (more than 3 mm of displacement or 3 degrees of rotation in any direction), and another due to data quality issues. This resulted in a total of 69 participants for the fMRI data analysis.

Time series extraction and filtering

To obtain a fine-grained parcellation of the whole cerebral cortex, brain regions were created by drawing a 10-mm-diameter sphere around two different sets of ROI coordinates; the 264-region system proposed by Power et al. (2011) and the Harvard-Oxford Subcortical Structural Atlas. The latter was used to define three missing subcortical structures: bilateral amygdala, hippocampus and caudate, which are relevant in depression-related research (Servaas et al., 2013). These resulted in a total of no-overlapping 270 ROIs. A group whole-brain mask was generated in MNI space based on the EPI images of all participants to locate the parts of the brain that were free from susceptibility artifacts in all participants. When a region overlapped less than 50% with the group mask at 90% mean signal intensity, that region was excluded from further analysis. This procedure resulted in the exclusion of 47 Power ROIs. For the remaining 223 regions (Supplementary Table S1), the BOLD time series were extracted as the average signal across voxels within each region. For each participant, the 223 ROI time series were filtered using a ninth-order Butterworth filter with passband 0.02 to 0.1 Hz, discarding low-frequency drifts (Smith et al., 1999) and high-frequency components associated with cardiac and respiratory signals (>0.1 Hz) and covering the most meaningful frequency range of resting-state fluctuations (Biswal et al., 1995). All analyses were performed on this pre-processed data.

Pre-defined functional networks

The brain was mapped into 13 predefined functional networks from the Power atlas, which was derived from resting-state fMRI data (Supplementary Table S1). These include the sensorimotor (25 ROIs), default mode (45 ROIs), cingulo-opercular task control (12 ROIs), frontoparietal task control (23 ROIs), subcortical (19 ROIs), salience (18 ROIs), auditory (13 ROIs), visual (29 ROIs), ventral attention (8 ROIs), dorsal attention (11 ROIs), memory retrieval (5 ROIs), cerebellum (4 ROIs), and uncertain (11 ROIs) networks. For the explorative analysis of specific functional networks, 12 networks were included; the uncertain network was excluded because it does not represent a network with explicable function.

Ignition-based measures of dynamical complexity

The intrinsic ignition analysis, proposed by Deco and Kringelbach (2017), allows characterizing the degree of integration in the brain that results from spontaneous events arising over time. These events reveal the ability of a given region to start the propagation of neural activity (i.e., ignition) to other regions eliciting varying degrees of integration in the brain. In turn, integration reflects the capacity of the brain to become interconnected and exchange information. This process is illustrated in Fig. 1.

Fig. 1

Calculation of ignition-based measures of dynamical complexity. (A) For each filtered time series, a threshold (i.e., mean plus standard deviation of the signal) is used to define events as those crossing the threshold from below. These intrinsic ignition events reflect the capability of a given brain region to start the propagation of neuronal activity to other regions in the brain. For each event elicited (green rectangle), the activity in the rest of the brain (red rectangle) is measured in the time window of 3 TRs (grey area). (B) In each time window, co-occurring events are represented in a binary matrix of size 223 × 223 regions. (C) In this matrix, we find the length of the largest connected component as a measure of integration. This process is repeated for each driving event. (D) For each brain region, local measures of ignition and ignition variability are calculated, respectively, as the average and standard deviation of the integration across events. Ignition is a dynamic spatiotemporal measure that enables the investigation of how global integration is elicited from a single brain region over time. In turn, integration reflects the overall capacity of the brain to become interconnected and exchange information. At the whole-brain level, global integration and spatial diversity are calculated as the mean and standard deviation of the ignition scores across regions, respectively. Spatial diversity indicates the level of hierarchy concerning the contribution of each region to global integration. Temporal variability (i.e., the degree of dynamic flexibility or metastability) is calculated as the ignition variability averaged across regions. Ignition variability quantifies how local activity in one region changes over time, which is closely related to its local metastability. Adapted from Deco and Kringelbach (2017). (For color figure the reader is referred to the web version of this article). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For a given brain region (n), driving events are captured following the established procedures of Tagliazucchi et al. (2012) and set as a binary signal as follows. The filtered BOLD time series is first transformed into z-scores, z. Subsequently, a threshold θ (i.e., the mean plus the standard deviation of the signal) is applied, where the binary sequence σ = 1 if z > θ and crosses the threshold from below, otherwise σ = 0. This method has been shown to be threshold-independent (Tagliazucchi et al., 2012) given that the binarization results from the so-called Poincaré section (i.e., an approach that aim to reduce the dimensionality of a dynamical system). When a brain region triggers an event (green rectangle in Fig. 1A), the level of integration in the rest of the brain is calculated for a time window of 3 TRs –from the time when the trigger occurs (grey area in Fig. 1A), and the highest integration level within that window is used. Note that the number of TRs in the window should be shorter than the interevent distances and large enough to integrate information (Deco et al., 2017b, Deco and Kringelbach, 2017). This window width was selected given the time that the integration took to return to base level. For each time window, co-occurring events (i.e., two regions triggering an event) can be represented in a binary matrix (Fig. 1B). As introduced by Deco et al. (2015), the amount of integration in the binarized matrix (considered as an adjacency matrix) can be calculated using the size (i.e., the number of brain regions) of the largest connected component –i.e., the largest subgraph in which paths are connecting any pair of vertices and which connects to no additional vertices in the supergraph (Fig. 1C). The integration value is then normalized by the maximal number of connected brain regions (i.e., N = 223 regions). This process is repeated for each driving event. Subsequently, the average and the standard deviation of the elicited integration across events are calculated to characterize the ignition and ignition variability profiles of a given brain region.

At the global level, ignition and ignition variability scores can be averaged across all brain regions to produce a global measure of integration and temporal variability, respectively. Temporal variability indicates the degree of dynamic flexibility, also referred to as metastability. We also computed a measure of spatial diversity by calculating the standard deviation of the ignition values across all brain regions. This measure represents a good indicator of the hierarchy of information processing (Deco and Kringelbach, 2017). The underlying idea is that some regions are more relevant than others for the broadcasting of information across the whole brain.

Phase-based metastability

Furthermore, using an approach from dynamical systems, we characterized the degree of metastability in the brain in terms of how much the synchronization between regions fluctuates over time (Wildie and Shanahan, 2012). The overall metastability was calculated as the standard deviation of the Kuramoto order parameter (OP) across time. The OP measures the degree of synchronization in the brain and is expressed by the following equation:where n = 1, …, 223 brain regions and φn(t) refers to the instantaneous phase of the BOLD signal at a given time point. Instantaneous phases were obtained by applying the Hilbert transform to the bandpass-filtered time series. Constant levels of the OP over time mean that the system is in stable equilibrium whereas an increased variance of the OP indicates that the system is systematically switching from one state to the other, and therefore displays a more dynamical configuration.

Statistical analysis

Statistical analysis was conducted in R (R Core Development Team, 2013), scatterplots were produced with the R package ggplot2 (Wickham, 2016) and ggpubr (Kassambara, 2019) and brain networks were displayed using the BrainNet Viewer toolbox (Xia et al., 2013), implemented in Matlab 2015b (The MathWorks Inc., Natick, MA).

For our main analysis, we investigated the relationship between depressive ratings, as measured by MDI, and three global measures of dynamical complexity: global integration, spatial diversity (i.e., hierarchy), and temporal variability (i.e., metastability). We observed that the MDI scores were not normally distributed (Shapiro-Wilk test for normality; W = 0.878, p < 0.001). Due to the skewness (s = 0.98) of the data, non-parametric statistics were employed for all analyses. Specifically, Spearman correlation was used as a non-parametric alternative to Pearson correlation, to test the relationship between MDI scores and measures of brain dynamics. The significance level of the correlation coefficient was determined by means of permutation testing. We used gender-restricted permutations to account for the Kolmogorov-Smirnov distance (D = 0.34, p = 0.039) between the distribution of MDI scores in women and men (Supplementary Fig. S1B). Specifically, a 10,000-permutation vector was generated by randomly reordering MDI scores within the female and male categories. Within the permutation test, the significance threshold was set to alpha = 0.05. This approach was chosen after having examined the potential effect of age or gender on the outcome measures. Spearman correlation tests revealed no significant effect of age, nor Mann–Whitney U test revealed a significant effect of gender in any of the measures, except for MDI. As expected, women reported higher MDI scores than men (U = 813.5, p = 0.008). Supplementary Table S1 shows the results of these tests.

In addition to our main analysis, we explored the relationship between MDI scores and the measures of integration and metastability within each of the 12 functional networks (see Section 2.3.3. Pre-defined functional networks). We calculated the network-driven integration and the network temporal variability (i.e., metastability) as the average ignition and ignition variability scores across the regions within a given network, respectively. The relationship between these two measures and MDI scores was likewise examined using Spearman correlation, whereby the significance of the correlation coefficients was determined using gender-restricted permutations (N = 10,000). In absence of any pre-established hypotheses about specific networks involved, Bonferroni corrected p-values were applied to account for multiple comparisons –i.e., 12 hypothesis tests, corresponding to the 12 functional networks being tested–, by controlling the family-wise error rate (FWER). Only networks for which significant results were found after Bonferroni corrections (alpha = 0.05/12) are reported in the main text. All computed correlations, corrected and uncorrected p-values are shown in Supplementary Table S3 and S4.

Data/code availability statement

The dataset is not publicly available due to restrictions from the informed consent, as participants have not given consent to have their data publicly stored, and European Data Privacy Regulations (GDPR). De-identified behavioral and imaging data (NIfTI format) are available upon request via j.b.c.marsman@umcg.nl.

The code for computing these measures of dynamical complexity is publicly available at github.com/sonsolesalonsomartinez/ignition

Results

Measures of integration, spatial diversity, and temporal variability were calculated to characterize the dynamic spatiotemporal organization of resting-state whole-brain connectivity. A Spearman correlation was computed to assess the relationship between these measures and the degree of self-reported depressive symptoms as measured by MDI. Bivariate outliers were identified for all correlation analyses and depicted in all scatterplots by a red circle. The results presented here included the entire dataset (N = 69) since there were no theoretical reasons to exclude the possible outliers, and given that removing the outliers did not change our results (Supplementary Table S6).

Global integration

Global integration (mean = 0.102, SD = 0.009) was calculated as the average ignition values across all brain regions. As described in the methods Section 2.4. Ignition-based measures of dynamical complexity, ignition denotes the contribution of a given brain region to starting the propagation of neural activity that ultimately leads to different levels of global integration. In turn, integration reflects the brain’s capacity to incorporate and process spatially distributed information. Individual differences in global integration negatively correlated with MDI scores, such that lower levels of integration were associated with individuals reporting more frequent depressive symptoms (r = -0.32; p = 0.007). The scatterplot in Fig. 2A summarizes the results.

Fig. 2

Association of decreased ignition-driven integration with increased depressive ratings. (A) displays the negative correlation between MDI scores and whole-brain integration (i.e., ignition scores averaged across all 223 regions). (B) displays the negative correlation between MDI scores and network-driven integration (i.e., ignition scores averaged across the 18 regions of the salience network). Brain regions were rendered on a brain surface using the BrainNet Viewer. Only functional networks showing a significant correlation after Bonferroni corrections are displayed. In all scatterplots, individual dots represent scores for 69 participants. The circled points are outliers identified using the bagplot rule (Rousseeuw et al., 1999; see bagplots in Supplementary Fig. S3). R values are from the Spearman correlation coefficient, the solid line represents the line of best fit and the shadowed area indicates the 95% confidence intervals. The significance level of a correlation coefficient resulted from a 10,000 (gender-restricted) permutation test at alpha = 0.05 (in Fig. 2A) and alpha = 0.05/12 (in Fig. 2B).

We further explored the importance of particular functional networks for the broadcasting of information across the brain. For each of the 12 functional networks, we calculated the network-driven integration as the average ignition scores across the regions within a given network. After Bonferroni corrections (alpha = 0.05/12), a negative relationship (r = -0.36; p = 0.024) was found between MDI scores and averaged ignition scores within regions of the salience network (Fig. 2B). Supplementary Table S3 displays the results for the 12 functional networks.

Spatial diversity

Spatial diversity (mean = 0.019, SD = 0.007), computed as the standard deviation of the ignition values across regions, defines the level of heterogeneity in the brain concerning the contribution of each region to the global integration. This heterogeneity provides further information about the level of hierarchy in the functional organization of the brain. There was a negative correlation (r = -0.25; p = 0.048) between MDI scores and spatial diversity (Fig. 3A), indicating that higher ratings of depressive symptoms were associated with a decreased hierarchy of information processing. To provide a visual representation of this hierarchy (Fig. 3B), participant-specific ignition values of all regions were sorted from the largest scores and color-coded according to their MDI scores. The depicted inverted S-curve, or more formally, the inverse sigmoid function, revealed that the organization of brain activity is (non-uniformly) hierarchical (Deco and Kringelbach, 2017). Moreover, this curve became less steep in individuals reporting more frequent depressive symptoms, which indicates a decrease of hierarchy whereby the contributions to global integration were more evenly distributed across regions.

Fig. 3

Association of decreased hierarchy with increased depressive ratings. (A) displays a negative correlation between MDI scores and spatial diversity (i.e., the variability of the ignition scores across all 223 regions). Individual dots represent scores for 69 participants. No outliers were identified using the bagplot rule (Rousseeuw et al., 1999; see the corresponding bagplot in Supplementary Fig. S3). R values are from the Spearman correlation coefficient, the solid line represents the line of best fit and the shadowed area indicates the 95% confidence intervals. The significance level of a correlation coefficient resulted from a 10,000 (gender-restricted) permutation test at alpha = 0.05. (B) participant-specific ignition values of all regions were sorted from the largest scores and color-coded according to their MDI scores. The color bar indicates the MDI scores from 0 to 45. (For color figure the reader is referred to the web version of this article).

Temporal variability

Ignition-based metastability

Temporal variability (mean = 0.046, SD = 0.013) was computed as the average ignition variability scores across all regions in the brain. Ignition variability quantifies how local activity in one region changes over time, which is closely related to its local metastability (Deco and Kringelbach, 2017). A higher ignition variability indicates that the propagation of neural activity fluctuates more over time, whereas lower ignition variability represents more stability –less dynamism. We found a significant negative association (r = −0.34; p = 0.005) between temporal variability and MDI scores (Fig. 4A). This indicates that the higher the depressive ratings, the less dynamic is the pattern of activity over time.

Fig. 4

Association of decreased metastability with increased depressive ratings. (A) displays a negative correlation between MDI scores and temporal variability (i.e., ignition variability scores averaged across all 223 regions). (B) displays negative correlations between MDI scores and network temporal variability (i.e., ignition variability scores averaged across the regions within the COC, DAN, MR and salience networks). Brain regions were rendered on a brain surface using the BrainNet Viewer. Only functional networks showing a significant correlation after Bonferroni corrections are displayed. In all scatterplots, individual dots represent scores for 69 participants. The circled points are outliers identified using the bagplot rule (Rousseeuw et al., 1999; see bagplots in Supplementary Fig. S3). R values are from the Spearman correlation coefficient, the solid line represents the line of best fit and the shadowed area indicates the 95% confidence intervals. The significance level of a correlation coefficient resulted from a 10,000 (gender-restricted) permutation test at alpha = 0.05 (in Fig. 4A) and alpha = 0.05/12 (in Fig. 4B). COC = cingulo-opercular task control; DAN = dorsal attention; MR = memory retrieval.

We also explored the relationship between MDI scores and network temporal variability for each of the 12 functional networks. A network temporal variability score was calculated as the average ignition variability across regions within a given network. After Bonferroni corrections (alpha = 0.05/12), a significant negative correlation between temporal variability and MDI scores appeared in the cingulo-opercular task control (r = −0.35; p = 0.035), dorsal attention (r = −0.36; p = 0.044), memory retrieval (r = −0.39; p = 0.005), and salience (r = −0.35; p = 0.025) networks. Scatterplots of these results are displayed in Fig. 4B. Results of the correlation analysis for all 12 functional networks are presented in Supplementary Table S4.

Phased-based metastability

Additionally, we investigated temporal variability using a phase-based measure of metastability, which measures how much the level of synchrony changes over time. Phased-based metastability was calculated as the standard deviation of the Kuramoto parameter. Similar to the results obtained using the ignition-based measure of metastability, we found that phase-based metastability (mean = 0.090, SD = 0.018) gradually decreased as MDI scores increased (r = −0.33; p = 0.003). A scatterplot in Supplementary Figure S2A summarizes the results.

To ensure that the observed measure can be ascribed to the dynamics, we used surrogate time series under the assumption of linearity and stationarity. First, 10,000 surrogate time series were generated for each participant, by randomly permuting the phases of the original time series across time. Subsequently, we calculated the Spearman correlation coefficients between MDI scores and phase-metastability values for all the surrogate set. The p-value to reject the null hypothesis was computed as the probability of the observed correlation coefficient given the null‐distribution. The null-hypothesis of linearity and stationarity using the surrogate time series was rejected at a significance level of alpha = 0.05 (p = 0.001). Supplementary Fig. S2B shows the null distribution of the correlation coefficients calculated using the surrogate time series.

Effect of motion on brain dynamics

We conducted a follow-up analysis to confirm that the observed relationship between depressive symptoms and the measures of dynamical complexity were unrelated to head movement. To assess motion, we calculated the mean Framewise displacement (FD) for each participant (mean = 0.15; SD = 0.04). FD estimates the participant́s head movement of any given frame relative to the previous frame. Supplementary Fig. S4 displays for each participant, the mean and standard deviation FD values across frames. Subsequently, we computed the Spearman correlation coefficient between mean FD and the global and network measures of dynamical complexity (i.e., integration, spatial diversity, and temporal variability). To avoid false-negative effects of motion we did not correct for multiple comparisons. No significant Spearman correlations were found between any of these measures and mean FD (Supplementary Table S5).

Discussion

In the present study, we applied a novel intrinsic ignition analysis to investigate whether the severity of depressive symptoms in non-clinical individuals was associated with changes in the dynamical complexity of the brain at rest. Our approach was motivated by the emerging discipline of whole-brain computational connectomes, which suggests the importance of a dynamical regime in cognitive, behavioral, and social abilities (Deco et al., 2017a, Deco and Kringelbach, 2016). Dynamical complexity was characterized using ignition-based measures of integration, hierarchy, and metastability. As expected, decreased global integration (Fig. 2), decreased spatial diversity (Fig. 3), and decreased temporal variability (Fig. 4) correlated with increased depressive ratings as measured by MDI.

Reduced global integration was associated with increased depressive ratings

Increased depressive symptoms were associated with diminished global integration, implying a reduced capacity of the brain to regulate and merge the incoming of information from distributed regions. This finding is consistent with other studies reporting reduced global integrity in patients with depression (Manelis et al., 2016, Wang et al., 2017b). However, also increases in global integration have been documented in the MDD literature (for a review see Suo et al., 2018). These conflicting results may be partially explained by the heterogeneity of the population sample, such as individual differences in the number of previous episodes (Meng et al., 2014) and the use of different definitions of altered network properties (i.e., degree centrality, global efficiency, and characteristic path length). Unlike these studies, we estimated the level of integration across the whole-brain network by considering the diversity of computation both, in space and time. In other words, global integration resulted from the distinct and heterogenous contributions across brain regions over time (i.e., ignition profiles). Global integration is crucial for efficient global communication, which enables adaptive cognitive and behavioral responses (Deco et al., 2015). Accumulating evidence suggests that depressive symptoms are linked to suboptimal cognitive and emotional processes (Joormann and Quinn, 2014, Nolen-hoeksema et al., 1993, Rood et al., 2009). It is possible that the breakdown of information processing underlying the observed reduced global integration accounts for these maladaptive processes. Impaired functional integration has also been identified in individuals scoring high in neuroticism (Servaas et al., 2015). This finding is of particular interest given that this personality trait is considered a strong risk factor for the onset of depression (Jeronimus et al., 2016). In general, integration deficits have been linked to other diseases, including schizophrenia (Damaraju et al., 2014, Lynall et al., 2010), bipolar disorder (Rashid et al., 2014) and Alzheimer’s disease (Seo et al., 2013), highlighting the cross-diagnostic relevance of integration in neuropsychiatric disorders (van den Heuvel and Sporns, 2019).

Reduced hierarchy was associated with increased depressive ratings

Via the intrinsic ignition analysis, it was possible to estimate the importance of each brain region for the propagation of neural information (i.e., ignition). Then, the assessment of the variability of these ignition levels across brain regions provides information about the degree of spatial diversity (i.e., hierarchy) in the whole brain. A reduced hierarchy of information processing was associated with increasing self-reported depressive symptoms. Higher hierarchy implies higher functional specialization, whereby distinct brain regions perform distinct roles. In particular, regions with low levels of ignition are considered to be relevant for sensory processing, whereas regions with high levels of ignition are more computationally relevant for the integration of information processing (Deco and Kringelbach, 2017). Notably, our results showed that depression-related interindividual variability was maximal in the regions with the highest ignition scores. Prior studies have noted the importance of a hierarchical organization for whole-brain integration (Deco et al., 2015) as it enables the complex brain functions that are required for normal behavior and cognition (Deco and Kringelbach, 2017, Lord et al., 2017, Park and Friston, 2013, Sporns, 2014, Zamora-López et al., 2010). Thus, a possible interpretation for the association between increased depressive ratings and reduced degree of hierarchy might be related to the altered affective and cognitive processing styles underlying depressive symptoms (Nolen-hoeksema et al., 1993).

Reduced metastability was associated with increased depressive ratings

Along with reduced hierarchy, diminished metastability (i.e., temporal variability) was also associated with the severity of depressive symptoms, suggesting that individuals with more symptoms of depression exhibited a narrower brain dynamic regime (Deco et al., 2017a). This finding agrees with the argument of Ghosh et al. (2008) that a reduced dynamic repertoire is implicated in disorders characterized by emotional and cognitive dysfunction. Recent computational models have revealed that the brain dynamics at rest operate at a maximal level of metastability (Deco et al., 2017a, Deco and Kringelbach, 2016, Jobst et al., 2017, Lee and Frangou, 2017). In other words, they propose that a healthy brain tends to be in a state of maximum network switching, enabling the exploration of a rich dynamical repertoire. The underlying idea is that during periods of fluctuation the system destabilizes, promoting the entrance of new information and consequently facilitating the exploration of a larger set of brain functional configurations. In line with our finding, a study by Demirtaş et al. (2016), using a similar analytical approach, reported increased temporal stability in patients with MDD. More generally, the link between reduced temporal variability and depression severity represented in our data is coherent with prior research implicating compromised metastability in other conditions including autism spectrum disorder (Fu et al., 2019, Watanabe and Rees, 2017) –with a high rate of depression comorbidity, schizophrenia (Kottaram et al., 2019, Koutsoukos and Angelopoulos, 2018, Lee et al., 2018), Alzheimer’s disease (Córdova-Palomera et al., 2017) and traumatic brain injury (Hellyer et al., 2015). In the last study, reduced levels of metastability after traumatic brain injury were associated with reduced cognitive flexibility and disrupted information processing. Interestingly, as proposed by other researchers in the field, higher levels of metastability might be acquired through mental training. A recent study (Escrichs et al., 2019) employing the intrinsic ignition framework demonstrated increased metastability in meditation experts while they were resting in the scanner. Similar results were reported in yoga practitioners (Gard et al., 2014) and high trait mindfulness individuals (Marusak et al., 2018). These findings are of interest, given the benefits of meditation for depression, anxiety and other negative emotional symptoms (Jain et al., 2015), which often increase the risk of major depression.

Contributions from specific functional networks

While it was not the main aim here, we explored whether whole-brain integration elicited by specific functional networks was associated with the severity of the depressive symptoms. Specifically, we found imbalances in the functional role of the salience network in dynamically integrating information over time to be linked to higher depressive ratings. This result seems consistent with the implications of the salience network in the orchestration of brain dynamics (Menon and Uddin, 2010, Uddin et al., 2011), and with studies reporting impairments of this network in MDD (Peterson et al., 2014, Sliz and Hayley, 2012), late-life depression (Yuen et al., 2014) and at-risk individuals (Kenett et al., 2018). In particular, similar to our findings, the study by Kenett et al. (2018) revealed a negative relationship between subclinical depression and integration abilities in the right anterior insula. The salience network is considered to mediate between-network interactions, constituting a key element for the integration of sensory, emotional, and cognitive information (Craig, 2009, Menon and Uddin, 2010, Uddin et al., 2011). Accordingly, altered functional integration of this network could in principle be linked to the inappropriate allocation of cognitive and attentional resources that are often linked to depressive symptoms (Wang et al., 2017a).

In addition to this finding, three other networks also showed attenuated temporal variability in association with the severity of depressive symptoms. These included the cingulo-opercular task control network, putatively involved in cognitive control (Cohen and D’Esposito, 2016, Sadaghiani and D’Esposito, 2015); the dorsal attention network, implicated in sustained attention and working memory (Corbetta and Shulman, 2002); and the memory retrieval network. Given that time-varying fluctuations are believed to play a central role in cognitive and behavioral functioning (Avena-Koenigsberger et al., 2017, Chang and Glover, 2010, Hellyer et al., 2015, Hutchison et al., 2013, Kringelbach et al., 2015), it is possible that the reduced temporal variability found in these networks relates to the difficulties in network switching seen in patients with depression, and might further account for the maladaptive processes of rumination and increased self-focus (Nolen-hoeksema et al., 1993). Finally, the fact that depression-related differences in these three networks only became apparent with respect to how much their ignition level fluctuated over time, but not with respect to the overall level of ignition, indicates that additional information may be learned from these non-static representations of otherwise dynamic brain activity.

Limitations, strengths and future work

Two parts of this study require further consideration: the population sample and the ignition approach. Concerning the population sample, this study moved beyond comparisons between healthy and patient populations to examine changes in brain activity along the continuum of depressive symptoms, as advocated by the research domain criteria approach. By using the interindividual response variability to a relationship breakup, it was possible to investigate depressive symptoms in participants without a clinical diagnosis. This strategy was based on the extensive evidence that after a relationship breakup most people respond with symptoms of depression (Field et al., 2009, Kendler et al., 1999, Najib et al., 2004, Verhallen et al., 2019). Yet, the inclusion of only individuals with a recent relationship breakup suggests the need for replication in a more heterogeneous sample. As far as the method is concerned, one of the advantages of the present study is the application of the intrinsic ignition analysis, which enabled the identification of global integration changes along the continuum of depressive symptoms. More importantly, this novel approach provided valuable information regarding the underlying spatiotemporal dynamics, encouraging us to acknowledge the presence and the relevance of these dynamics. Although we also explored the contributions of individual functional networks, our analysis does not explain how these networks form and dissolve over time to induce different levels of integration. Flexible cognition and behavior depend on the balance between integration and segregation (Cohen and D’Esposito, 2016, Deco et al., 2015). Thus, while this study focused only on integration, future research should explore the interplay between these two processes.

Conclusion

The goal of the present work was to investigate whether changes in the dynamical complexity were associated with the variability of depressive symptoms reported by individuals after a recent relationship breakup. The application of a novel intrinsic ignition analysis revealed that reduced global integration and reduced spatiotemporal dynamics, that is, less hierarchical organization and less functional variability over time, scales along the continuum of depressive symptoms even in participants without a clinical diagnosis. Given the growing evidence that demonstrates altered resting-state dynamics across neuropsychiatric disorders, our results in a nonclinical (yet vulnerable) population sample suggest the merit of investigating brain rigidity, understood as less complex brain dynamics, as a potential risk marker for mental health problems. Our approach may provide new opportunities for understanding depressive symptoms in the general population, offering dimensionality across health and disease.

Funding

This work was funded by a donation of Mr. Hazewinkel to Prof. dr. Gert J. ter Horst and the Research School of Behavioural and Cognitive Neurosciences, University of Groningen, the Netherlands. The funder had no role in study design, the collection, analysis and, interpretation of data, the decision to publish, or the writing of the manuscript.

CRediT authorship contribution statement

Sonsoles Alonso Martínez: Conceptualization, Formal analysis, Writing - original draft, Writing - review & editing. Jan-Bernard C. Marsman: Investigation, Resources, Data curation, Writing - review & editing. Morten L. Kringelbach: Methodology, Writing - review & editing. Gustavo Deco: Methodology, Writing - review & editing. Gert J. Horst: Conceptualization, Resources, Supervision, Funding acquisition.