Probing the association between resting-state brain network dynamics and psychological resilience

Abundant literature in human clinical neuroscience has established a link between changes in intrinsic functional connectivity of large-scale brain networks and psychological disorders such as depression (Mulders et al., 2015) or schizophrenia (Dong et al., 2018; Li et al., 2019; see also Menon, 2011, for a general overview). As a consequence, brain bases of preserving mental health in the face of stress and adversity (resilience) have also become a focus of interest (e.g., Southwick & Charney, 2012). Two recent studies (Long et al., 2019; Paban et al., 2019) have reported associations between psychological resilience and brain network dynamics by using multilayer modularity, a relatively new tool from the evolving field of network neuroscience that integrates spatial and temporal information in a graph theoretical framework (Muldoon & Bassett, 2016). In both studies, resilience was assessed with a frequently used questionnaire, the Connor–Davidson Resilience Scale (CD-RISC; Connor & Davidson, 2003), and network dynamics under task-free (resting-state) conditions were examined by detecting time-evolving patterns of nonoverlapping and coherent modules and by quantifying the frequency with which brain nodes switched between modules (node flexibility; Bassett et al., 2011).

The link between resilience and brain network dynamics is motivated by (a) an assumed relationship between neuronal and cognitive flexibility (e.g., Braun et al., 2015; Chen et al., 2016) and by (b) associations of resilience with higher cognitive flexibility, both in theoretical models (Southwick & Charney, 2012) and in empirical work (e.g., Genet & Siemer, 2011). While this would predict greater resilience in more (cognitively or neuronally) flexible persons, it has also been argued that cognitive flexibility might not be universally adaptive and that resilience may depend on an interplay between flexibility and stability (Parsons et al., 2016). This proposal receives support from observations of changes along flexibility-stability dimensions in psychiatric conditions, both behaviorally (e.g., rigid behaviors like rumination in depression; Nolen-Hoeksema et al., 2008) and neuronally (e.g., increased network flexibility, variability in patients with autism or schizophrenia; Gifford et al., 2020; Harlalka et al., 2019). These and similar results suggest that the extreme tails of the ‘flexibility-stability’ distribution (Kashdan & Rottenberg, 2010) may indeed be related to psychopathology and that, accordingly, adaptive behavior in the ‘normal’ (i.e., unaffected) range depends on a balance between flexibility and stability. Whereas the above-cited work by Long et al. (2019) and Paban et al. (2019) took a primarily neuroscientific perspective as a starting point, the psychological perspective offered here is not fully consistent with their empirical results. However, it may offer a valuable conceptual framework for a neurocognitive model of resilience and therefore calls for further empirical research to clarify the current inconsistencies.

Paban et al. (2019) measured task-free EEG and conducted network analyses in source space. Negative correlations between node flexibility and resilience were observed in the alpha, beta, and delta frequency bands, including superior parietal cortex, medial orbitofrontal cortex, and cuneus. A subsequent resting-state functional MRI (rs-fMRI) study (Long et al., 2019) also reports lower node flexibility in more resilient persons, primarily in visual cortices and the left medial-orbital superior frontal gyrus. Except for a partial overlap in visual regions (lingual gyrus), localization results differed between studies and Long et al. (2019) did not replicate correlations in ‘higher order’ (superior parietal, inferior frontal) areas reported by Paban et al. (2019). Whereas some of these inconsistencies may result from inherent differences between methods (e.g., differences in temporal resolution between fMRI and EEG), others may reflect specific methodological choices by the authors. For example, Long et al. (2019) studied BOLD ‘dynamics’ using only 12 nonoverlapping time windows (each of 20 sample points length, derived from 250 measurements of TR = 2 sec).

Even though gold standards for analyses of dynamic multilayer modularity are yet to be established, recent studies suggest that more data are required for sensitive and reliable estimation of network dynamics from BOLD fMRI (Hindriks et al., 2016; Yang et al., 2021). To ameliorate such methodological shortcomings and to further our understanding of the relationship between resilience and brain network dynamics, we here replicate and extend these results by using temporally highly resolved fMRI (multiband/MB sequence; MB factor 4; TR = .675 sec) from 52 healthy young adults who completed three resilience questionnaires. We first replicated Long et al.’s (2019) analysis pipeline as closely as possible, by downsampling data to a TR of 2.025 sec and using the same analysis parameters. Following this, we explore effects of specific analysis choices (like different windowing schemes) on network flexibility. Lastly, we repeat correlation analyses with optimized denoising and the full temporal resolution of the MB data, resulting in a time series of 753 overlapping windows (Figure 1). Node flexibility was calculated as the number of times a node changes its community assignment between windows, divided by the total possible number of changes (Bassett et al., 2011). Given the results of Long et al. (2019), we expected that resilience and node flexibility should be negatively correlated.

To extend previous work, we also assessed node promiscuity and node degree, and their relation to resilience. Node promiscuity, the fraction of communities a node participates in at least once, and node flexibility are complementary measures that inform us whether brain dynamics per se (i.e., assessed via node flexibility) or the diversity of brain systems with which a node interacts (assessed via its promiscuity) may be associated with resilience. Node degree, in turn, is a proxy for dynamic connectivity that does not rely on modularity detection algorithms (and is thus invariant to potential algorithmic idiosyncrasy or parameter choices). However, previous work demonstrates that dynamic connectivity is related to various cognitive and behavioral traits, as well as clinical conditions like posttraumatic stress disorder (PTSD, often used as proxy for studies of resilience; Bolsinger et al., 2018; Jin et al., 2017; Lurie et al., 2020), which makes it a candidate marker for analyzing the neurobiology of resilience. To the best of our knowledge, relationships between resilience and promiscuity as well as degree have so far not been studied. Accordingly, we here aim at a more complete characterization of putative relationships between brain network dynamics and resilience.

Descriptive statistics for resilience scales are listed in Table 1. Psychometric properties were satisfactory, with internal consistencies (Cronbach’s alpha) between .7 and .85 (Table 1). Resilience scores in our sample were comparable to the respective original publications for the German versions (CD-RISC: 30.6, Sarubin et al., 2015; RS-13: 70, Leppert et al., 2008; Brief Resilience Scale (BRS): 3.58 and 3.37, Chmitorz et al., 2018). Resilience ratings were significantly intercorrelated: CD-RISC − RS, r = .60; CD-RISC − BRS, r = .61; BRS − RS, r = .45, all p < .001 (FDR-corrected for multiple comparisons).

We replicated the analysis pipeline of Long et al. (2019) as closely as possible, involving downsampling of multiband data by a factor of 3 to a virtual TR of 2.025 sec, use of the two parcellation schemes used by Long et al. (2019), that is, AAL90 and Power264 atlases, and the same denoising strategy (with 26 parameters) as used in that study. Age, gender, and framewise displacement were included as covariates of no interest (see Methods). Given that recent evidence implies that functionally derived parcellations outperform anatomical atlases like the automated anatomical atlas (AAL) (Dadi et al., 2019), we also conducted all analyses with a further atlas, the functionally defined parcellation scheme of Schaefer (100 nodes), which approximately matches the number of nodes in the AAL. Following Long et al. (2019), the downsampled time series was segmented into nonoverlapping windows of 20 TRs length, resulting in 15 time windows (as opposed to 12 in the original study). We did not replicate the correlation between global flexibility and CD-RISC reported by Long et al. (2019; all p > = .33), but observe a borderline significant correlation between global flexibility (AAL atlas) and the BRS resilience score (r = −.33, p = .05). Across the three parcellation schemes and resilience questionnaires, we did not observe any further significant effects, neither at global, subnetwork, or node level (Table 2).

Segmenting the time series into nonoverlapping windows results in a low number of sample points, generally not considered sufficient for analyzing network dynamics. Many studies today segment BOLD time series into overlapping windows, as this allows for an estimation of network dynamics with greater sensitivity and reliability (Hindriks et al., 2016; Lurie et al., 2020; Yang et al., 2021). When analyzing our downsampled data with a high number of overlapping windows (i.e., 251; see Methods), nodal flexibility was scaled by a factor of ∼1/10; cp. y-axis scaling of Figures 2A and 2B), as expected given the higher overlap between consecutive windows. However, high spatial similarity was preserved between node flexibility derived from downsampled data with overlapping versus nonoverlapping windows (r = .90, p < .0001; Figure 2C and 2D), indicating that there are no qualitative differences between windowing schemes. We thus use overlapping windows for all further analyses.

Analyses of network dynamics may be limited by lower temporal resolutions (which was here simulated by downsampling to allow for a direct replication of the results by Long et al., 2019). As compared to the original data (TR = .675 sec), downsampling shifted node flexibility (calculated using overlapping windows; see section Replication of Long et al. (2019)) toward higher values. However, high spatial similarity (i.e., correlation) exists between the distribution of nodal flexibility in original versus downsampled data (r = .95, p < .0001). As outlined above, we will in the following also investigate node promiscuity and degree as further measures of brain network flexibility. We thus also explore how these measures are affected by downsampling of the BOLD data: as for flexibility, downsampling led to higher promiscuity values, but a high spatial similarity was preserved (r = .93, p < .0001). We observed a significantly lower nodal degree in the downsampled data (0.06–0.61% significant connections; M = 0.20) compared to the original data (0.12–0.75%; M = 0.33), t(99) = 13.93, p < .0001, which is consistent with a recent report by Pedersen et al. (2018) that fewer data result in lower nodal degree values.

Given that no previous data exist to derive specific hypotheses concerning the relationship between resilience and node promiscuity as well as degree, all subsequently reported associations were tested two-sided using Spearman correlations. Correlations were calculated for all three parcellations, for original (TR = .675 sec) and downsampled (TR = 2.025) data, using an optimized denoising strategy (36 parameters; see Methods), overlapping time windows, and the same covariates as during the replication attempt reported above.

In the original data, we observed two significant positive correlations at the nodal level of the Power264 atlas between flexibility of a cingulo-opercular node (#54) and the CD-RISC (r = .52, p = .02), and between promiscuity of a visual node (#165) and the BRS resilience score (r = .57, p = .004). After downsampling (but now exploiting higher temporal resolution due to the overlapping windowing scheme), the AAL-based analysis yielded a significant association between the promiscuity of the left pallidum (node level, AAL) and BRS resilience (r = −.49, p = .02). In the Schaefer100 analysis, significant negative correlations emerged between the degree of a visual node (left Vis_7) and the RS score (r = −.48, p = .03), as well as between global degree and BRS resilience (r = −.35, p = .03). In the Power264 atlas we observed two borderline significant results, between the BRS resilience score and the degree of the subcortical resting-state network (RSN) (r = −.37, p = .05) and the cerebellar RSN (r = −.37, p = .05). No correlation was observed for any additional tested combination (see Table 3). Figure 3, Figure 4, and Figure 5 visualize nodal and RSN results in an exemplary manner for the Schaefer100 parcellation.

To further investigate putative factors that may have influenced reported results, we conducted several post hoc analyses on data derived from the Schaefer100 parcellation.

Recalculation of all node-level correlations without covariates yielded results comparable to the above-reported findings (cf. Supporting Information File 1, Figures S1 and S2). Correlations for original (MB) resolution were −.41 < r < .44, all p > .11. For downsampled data, the above-reported significant correlation for one visual node disappeared (−.43 < r < .30, all p > .14).

To explore how node flexibility changes as a function of window length, node flexibility was calculated for window lengths between 40 and 100 sec in steps of 20 (original/MB data). Similar to the effect of downsampling (see above), we observed an increase in node flexibility with decreasing window size (Figure 2E), however, with high spatial similarity of flexibility values across window sizes, all r > .87, p < .0001 (Figure 2F). We did not observe any significant correlations between node flexibility and resilience when varying the window size between 40 and 80 sec (all p > .13).

While repeating the correlation analyses (Schaefer100 parcellation, overlapping windows) after denoising with a 26-parameter regression model (following Long et al., 2019) as opposed to the 36-parameter model used above, no correlations with resilience were found at global (multiband: all p > .54; downsampled: all p > .74) or nodal level (multiband: all p > .15; downsampled: all p > .53). In the RSN analysis, significant negative correlations were found between degree of the limbic RSN and all resilience scales (multiband; CD-RISC: r = −.39, p = .03; BRS: r = −.44, p = .008; RS: r = −.40, p = .03; downsampled: all p > .16).

As the Schaefer atlas is available in multiple parcellations, it is best suited to investigate the influence of spatial resolution. Using 200 (vs. 100) nodes, we did not observe any correlation on the global level, neither for original (all p > .50), nor downsampled (all p > .13) data. While there were also no significant correlations on the RSN level for original data (all p > .29), degree of the limbic RSN was negatively correlated to both CD-RISC (r = −.57, p < .001) and RS (r = −.38, p = .04) in downsampled data. No significant effects emerged at nodal level (original: all p > .36; downsampled: all p > .10).

For the 26-parameter denoised data, we observed a significant correlation between global flexibility and mean framewise displacement (FD) (r = .38, p = .02, FDR-corrected) derived from the original (multiband) data, whereas no other brain measure was correlated with FD for original or downsampled data (all p > .11) and no correlations emerged after optimized (36-parameter plus despiking) denoising (original: all p > .63; downsampled: all p > .05).

The present study investigated associations between intrinsic functional brain connectivity dynamics and psychological resilience. More specifically, we explored whether a previously reported negative correlation between resilience and node flexibility as determined using resting-state BOLD fMRI (Long et al., 2019) can be replicated in multiband BOLD data with higher temporal resolution (TR = .675 sec). Additionally, we aimed at extending previous findings by including two further metrics of network dynamics, as well as two further questionnaires measuring slightly different theoretical conceptions of resilience. Correlation analyses were performed at different topological scales, that is, for whole brain, (functional) subnetworks, and nodes. Results did not support the previously reported negative association between node flexibility and resilience, neither in a direct replication matching as closely as possible methodological parameters of the original study, nor when extending the methodology in terms of higher temporal resolution, finer grained sliding window approach, or further network measures. In these analyses, we found a small number of significant correlations that were distributed rather inconsistently across different measures of brain dynamics and resilience, across atlases, and across sampling schemes: for multiband data (original time resolution; Power264 parcellation), flexibility of a cingulo-opercular node correlated with resilience measured using the CD-RISC questionnaire and promiscuity of a visual node correlated with BRS resilience, however, both in the opposite direction relative to Long et al. (2019). For downsampled data, we observed negative correlations between global degree and BRS and between the degree of a visual node and the RS score (Schaefer100 parcellation). After increasing the spatial resolution to 200 nodes (follow-up analysis), these effects disappeared, and negative resilience correlations emerged for the limbic network (CD-RSIC, RS). When using less stringent denoising (as in the original study), we also observed negative correlations with limbic network degree (all three resilience scales), for multiband but not downsampled data. When applying the AAL atlas, a negative correlation emerged with promiscuity for a subcortical node (pallidum; BRS). Of note, we also observed a borderline significant (p = .05) negative correlation between global flexibility and resilience (BRS) in the direct replication. We do not consider this strong support in favor of a replication since we found no similar effects using the Power264 or Schaefer100 atlases or with the other two questionnaires. In sum, these heterogeneous results cannot be considered a successful replication. The fact that the limbic network was involved in these correlations multiple times suggests that its role for resilience should be investigated further. On the other hand, these effects were not consistent, so that they should be treated with caution. Lastly, post hoc analyses conducted with multiband data and the Schaefer100 parcellation suggest invariance of results against changes in windowing schemes (size, overlap of windows). In the following, we will discuss potential limitations of the present study, which may be important when considering reasons for differences between present and previous results. We then discuss in more depth factors that may account for the failed replication of the results of Long et al. (2019).

Potential limitations of investigating network dynamics in BOLD data with low temporal resolution, as in the original study, have already been discussed in the Introduction. When analyzing temporally highly resolved multiband data, we chose the window size to be equal to 1/f_min (Leonardi & Van De Ville, 2015). We cannot rule out that such long window sizes might hamper the sensitivity to small changes, as calculating connectivity matrices over longer periods serves as a smoothing kernel minimizing the ability to resolve time-varying functional connectivity (TVFC) (Vergara et al., 2019). However, this disadvantage of MB-echo planar imaging (EPI) is more than compensated by the ability to acquire more data in the same time, the possibility to investigate faster dynamics, and to thus better disentangle neural and physiological signals (Yang & Lewis, 2021). Furthermore, we showed empirically that node flexibility was largely robust against changes in window size.

Potential limitations may also stem from the questionnaires used: short versions of two scales (CD-RISC: 10 items; RS: 13 items) were used due to timing constraints, whereas the long (25-items) CD-RISC was used by Long et al. (2019). However, we deem this unproblematic as validity studies show that the short versions can be considered valid and reliable instruments for measuring trait resilience (Chmitorz et al., 2018; Leppert et al., 2008; Sarubin et al., 2015). Our data confirm this, as all questionnaires showed satisfying levels of internal consistency.

Even though the overall direction of correlations between node flexibility and resilience was the same in two previous publications using fMRI (Long et al., 2019) and EEG (Paban et al., 2019), neural localizations partly differed. Adding to this heterogeneity, correlation results in the present study were highly inconsistent (see above) and results thus provide no strong evidence for an association between network dynamics and resilience. This demonstrates the importance of replication attempts also in network neuroscience studies of brain-behavior correlations. To support such work in the future, we have compiled a list of methodological features that we have found important in the present replication study (Table 4).

When applying less stringent denoising, a positive correlation between global flexibility and mean FD emerged (in multiband time resolution, overlapping windows), suggesting that increased motion inside the scanner may artificially alter functional connectivity and network flexibility. We did not observe any other correlation between dynamic brain measures and motion, suggesting that denoising approaches and their effectiveness may interact with temporal sampling schemes (i.e., MB vs. conventional EPI). This underscores the importance of incorporating recent empirical insights regarding denoising strategies for module detection and analyses of network reconfigurations (Lydon-Staley et al., 2019).

Both windowing schemes used in our study (overlapping vs. nonoverlapping) are among the most used approaches within the sliding window framework (Iraji et al., 2020). For window-based TVFC analyses, window length is a critical aspect, as it has to balance the ability of robustly estimating TVFC (which benefits from higher numbers of timepoints) against susceptibility to noise (Damaraju et al., 2020; Iraji et al., 2020). Extensive studies have evaluated the impact of different window types, lengths, and overlap on TVFC, however, without yet identifying a gold standard (e.g., Abrol et al., 2017; Shakil et al., 2016). At present, lengths between 30 sec and 100 sec are commonly used (see Figure S1 in Preti et al., 2017). Within this range, one would, however, assume that strong effects should be robust against specific analysis choices (compare, e.g., Damaraju et al., 2020, vs. Haimovici et al., 2017). This is supported by our follow-up analyses showing that results did not vary between different window schemes or sizes.

Recent research suggests that multilayer modularity in particular (Yang et al., 2021) and TVFC in general (Hindriks et al., 2016) are a function of scanning duration and that rs-fMRI sequences around 10 min with standard TRs (as used by Long et al., 2019) are not sufficient for reliable parameter estimation. As reliability strongly influences which correlations are detectable (Hedge et al., 2018), short measurements may be prone to false positive correlations. Yang et al. (2021) propose that at least 20 min of rs-fMRI is needed for reliable multilayer modularity analysis—in that study based on a TR of 1.45 sec and nonoverlapping windows of size 100 sec. Here, we measured 10 minutes of rs-fMRI with a TR of about half that length and used a substantially larger number of windows, intentionally selected to allow for reliable multilayer analyses. Yang et al. (2021) also suggest that the standard settings for intra- and interlayer parameters (γ = ω = 1) in dynamic network measures may not necessarily turn out to be optimal across datasets and across different network measures in terms of reliability. However, when varying the intra- and interlayer parameters beyond the standard setting (in the original/MB data, Schaefer100 atlas; see Supporting Information), we did not observe any significant correlations between node flexibility or promiscuity and resilience scores, except for a single correlation between the flexibility of a visual node and BRS when setting γ = 1.4 (r = .50, p = .02). This correlation, however, points in the opposite direction of what had been reported before (Long et al., 2019). We can thus conclude that the specific choice of intra- and interlayer parameters has no strong effect on detecting resilience-flexibility correlations.

It is also important that the behavioral variables of interest show satisfying psychometric properties. Internal consistencies for the resilience questionnaires were satisfactory, ranging from .70 to .85. When investigating individual differences, samples should not be too homogeneous on the target measures, as interperson variability is a precondition for detecting correlations. Even though CD-RISC data appear to be in a somewhat narrower range and more left-skewed than in Long et al. (2019), both the RS-13 and the BRS show sufficient variability (Table 1), which increases confidence that our results are not driven by a lack of heterogeneity. Moreover, resilience scores in our sample were distributed similar to the respective original publications of these measures. It is, however, difficult to compare our behavioral results directly to the studies by Paban et al. (2019) and Long et al. (2019), as the former did not report descriptive values (but high internal consistency: Cronbachs α = .90), whereas the latter did report descriptive statistics but no psychometrics.

To summarize, our results do not provide support for the previously reported negative association between node flexibility and psychological resilience. We extended previous studies by including additional measures of functional connectivity dynamics and two further resilience questionnaires and found only weak and inconsistent evidence of associations between network dynamics and resilience, mostly for node degree, a proxy for TVFC, but not for node flexibility or promiscuity. Our study highlights how specific degrees of freedom in the analysis of functional connectivity may influence the presence or absence of effects of interest. This underscores the need for testing the robustness and generalizability of proposed effects via replication.

Code and data have been deposited at https://doi.org/10.5281/zenodo.5761102 and will be made publicly available after an embargo period that ends December 31, 2024.

In total, N = 69 right-handed university students were enrolled in the study, of whom N = 60 completed the study protocol and were included for further analyses. All participants were native speakers of German, right-handed, and between 18 and 35 years old; absence of current psychiatric episode was assured with a structured interview (MINI; Sheehan et al., 1998). This sample size was based on considerations of statistical power for investigating across-participant relationships between BOLD activation and cognitive flexibility (see preregistration: https://osf.io/a64jn); the task-free (resting state) fMRI data reported here were acquired together with the preregistered task-based fMRI experiments. During image quality checks and preprocessing, eight datasets were excluded from further analyses (four due to low quality of questionnaire data and four due to motion artifact; see below for details of exclusion criteria), so that the reported analyses are based on data from 52 participants (26 females, 25 males, 1 diverse; 18–34 years; mean age 24.0 ± 3.7). Using the pwr package in R (Champely et al., 2017), we estimated that our final sample size of N = 52 has a power of > .99 to detect correlational effects of the size reported by Long et al. (2019), that is, of around r = .55. All participants provided written informed consent, and all procedures were approved by the Ethics Committee of the Department of Psychology of Goethe University Frankfurt, Germany.

Resting-state functional MRI (rs-fMRI) data were collected prior to task-based fMRI (which is not part of the present report; see previous paragraph) on a 3-T Siemens Prisma MR-Scanner equipped with a 32-channel head coil, using a multiband (MB factor = 4) EPI sequence with the following parameters: 900 volumes (10:14 min), TR = 675 ms, TE = 30 ms, voxel size = 3 mm³, flip angle = 60°, FoV = 222 mm, acquisition matrix = 74 × 74, 40 slices. During the rs-fMRI measurement, participants were asked to keep their eyes open and gaze at a white fixation cross, located at the center of a screen (Nordic Neuro Lab, 40″, 1,920 × 1,080, 60 Hz), to stay relaxed and not to think about anything specific. In a separate session, a T1 weighted (T1w) 3D structural MR scan was acquired with a MPRAGE sequence (4:26 min, voxel size = 1 mm³, TR = 1,900 ms, TE = 2.52 ms, acquisition matrix = 256 × 256, 192 slices) for purposes of coregistration between functional and structural data.

The psychological construct ‘resilience’ was quantified using German versions of three self-report questionnaires, the 10-item version of Connor–Davidson Resilience Scale (CD-RISC – 10 items; Connor & Davidson, 2003; German version by Sarubin et al., 2015), the Brief Resilience Scale BRS – 6 items (Smith et al., 2008; German version by Chmitorz et al., 2018), and the Resilience Scale (RS – 13 items; Wagnild & Young, 1993; German version by Leppert et al., 2008). Note that short versions of the questionnaires were selected, as the present study was part of a more extensive study protocol and the short versions have been shown to allow for a time-efficient data collection with comparable validity and reliability (Chmitorz et al., 2018; Leppert et al., 2008; Sarubin et al., 2015). We included the CD-RISC and BRS because of their good ratings in an evaluation study by Windle et al. (2011). The RS was additionally chosen to increase comparability with other research, as this questionnaire is frequently used in resilience research. In addition, the selected questionnaires also differ in that CD-RISC and RS define resilience as a personality trait, whereas resilience is understood as an outcome in the BRS (see original publications for more details). Resilience questionnaires were filled out during an informational preparation session for the fMRI measurements and administered online using Unipark software (EFS Survey, Questback GmbH). As completion time of online questionnaires has been identified as the most reliable indicator of data being meaningful or meaningless (Leiner, 2013), we evaluated the quality of questionnaire data by using a quality index provided by the Unipark system that compares the completion time of each participant with the average completion time of our sample. As preregistered (https://osf.io/c94y8) participants with an Unipark quality index of .20 or lower were excluded, resulting in exclusion of five participants (see also Leiner, 2013, for a similar criterion). Resilience scores were calculated according to the respective manuals. For each resilience scale, we also assessed internal consistency by calculating Cronbach’s alpha and its 95% confidence interval.

Quality of imaging data was assessed using both fMRIPREP’s visual reports as well as MRIQC 0.15.2rc1 (Esteban et al., 2017, 2019). T1w and functional images for each participant were visually checked for signal artifacts, whole-brain coverage, and correct alignment between structural and functional data. Following a procedure proposed by Faskowitz et al. (2019), functional data were excluded if marked as an outlier (i.e., exceeding ×1.5 the interquartile-range either from Q1 or Q3) in more than 50% of the MRIQC quality metrics: dvars, tsnr, snr, efc, aor, aqi (see the MRIQC documentation for more information about these metrics; Esteban et al., 2017). Given the sensitivity of resting-state analyses to movements and given that some of the aforementioned metrics are influenced by motion, we additionally included FD as a metric for quantifying motion artifacts (Maknojia et al., 2019). Due to the higher sampling rate of multiband EPI, motion parameters exhibit a high-frequency (HF) component resulting from head motion due to respiration (Power et al., 2019) that is usually not observable with standard single-band fMRI (Williams & Van Snellenberg, 2019). With this ‘spurious’ HF motion component, the head appears to be in constant motion and summary measures (such as FD) might be contaminated, leading to exaggerated flagging of ‘bad’ volumes. Note that this applies less to subsequent functional connectivity measures, as those are routinely band-pass filtered, whereas summary measures (e.g., FD) are calculated on raw motion parameters (Gratton et al., 2020; Power et al., 2019; Williams & Van Snellenberg, 2019). We therefore calculated filtered FD (FD_filt) from Butterworth-filtered raw head movement traces to better isolate ‘true’ head movements (Power et al., 2019). Following Lydon-Staley et al. (2019), we excluded subjects with mean FD_filt >0.2 mm or >20 volumes with FD_filt >0.25 mm to ensure high data quality. We excluded four subjects based on their FD_filt, while no subject was excluded based on MRIQC’s metrics.

T1w and rfMRI images were preprocessed using fMRIPREP 20.1.1 (Esteban et al., 2019), which is based on Nipype 1.5.0 (Gorgolewski et al., 2011). A boilerplate text released under a CC0 license describing preprocessing details can be found in the Supporting Information. For further pipeline information, see fMRIPREP’s documentation (Esteban et al., 2019). Due to the use of MB data and the high sampling rate, no slice time correction was applied (see, e.g., Glasser et al., 2013). Distortion corrected functional images in T1w space were further denoised using the XPC Engine 1.2.1 (Ciric et al., 2017). We implemented a denoising strategy that has been shown to be relatively effective in mitigating motion artifacts in the study of dynamic functional connectivity, (multilayer) subnetwork detection, and measures of module reconfiguration (Lydon-Staley et al., 2019). BOLD data was first demeaned, detrended, and despiked on a voxelwise basis, instead of using more aggressive censoring methods that may result in varying window lengths across participants (Hutchison et al., 2013), and then temporally filtered with a first-order Butterworth filter by using a passband of 0.01–0.08 Hz. These operations were followed by a confound regression, that included (a) six motion estimates derived from fMRIPREPs realignment, (b) mean signals from white matter (WM) and cerebrospinal fluid (CSF), (c) mean global signal, (d) temporal derivatives of these nine regressors, as well as (e) the quadratic terms of all 18 parameters, resulting in a 36-parameter model to obtain residual BOLD time series. All regressors were also band-pass filtered to avoid reintroducing noise caused by a frequency-dependent mismatch (Hallquist et al., 2013). For the direct replication of Long et al. (2019) and to directly test the influence of different denoising schemes, we also implemented a more lenient 26-paramter denoising approach that included 24 motion parameters following Friston et al. (1996) together with mean signals from WM and CSF. Within the XCP engine, the Schaefer (resolution: 100, 200; Schaefer et al., 2018), AAL90 (Tzourio-Mazoyer et al., 2002), and Power264 (Power et al., 2011) atlases were transformed to native T1 space and resampled to match the BOLD images (see Supporting Information for details). For each subject the whole brain was parcellated into distinct regions and a functional time series was extracted for each region, corresponding to the average across all voxels within that region.

All subsequent TVFC analyses were performed with the numpy (1.18.5; Harris et al., 2020) and scipy (1.5.0; SciPy 1.0 Contributors et al., 2020) packages in python 3.8.3.

To directly replicate findings described in Long et al. (2019) and to investigate if different sampling rates (i.e., MB-EPI vs. conventional EPI) yield different results, we downsampled our original resting-state BOLD time series by a factor of three to resemble a TR of 2 as closely as possible (3 × .675 sec = 2.025 sec). Raw motion traces were also downsampled accordingly, and FD was calculated without filtering. For the direct replication, we segmented the time series into 15 nonoverlapping windows with length of 20 TR (40.5 sec), analogous to the windowing scheme used in Long et al. (2019). As this low number of windows may hamper the ability to investigate brain network dynamics (Hindriks et al., 2016; Lurie et al., 2020; Yang et al., 2021), we additionally used overlapping windows that were shifted by a single time point (TR). All windows were tapered with a Hamming filter to reduce potential edge artifacts and to suppress spurious correlations (Shakil et al., 2016; Zalesky et al., 2014), and pairwise Pearson correlation coefficients between all nodes were calculated within each window.

For analyses of the ‘original’ MB data with high temporal resolution, we used a fixed window length of 148 time points (100 sec), which satisfies the frequency criterion that the length should be at least be equal to 1/f_min (Leonardi & Van De Ville, 2015) and allows for a full oscillation of the slowest frequency in the range of 0.01–0.08Hz. After downsampling, sliding window correlations with an adjusted window length of 50 time points (101.25 sec) and subsequent TVFC and modularity analyses were calculated as described above. To investigate whether the calculated network measures differ between original and downsampled data and to rule out any idiosyncratic algorithmic behavior (i.e., from multilayer modularity), we compared all three metrics at the nodal level, as these were the starting points for all subsequent measures. For node flexibility and node promiscuity we calculated the spatial similarity (i.e., Spearman correlation coefficient) between mean values per node (average over participants) in the original and downsampled data to explore whether nodes behave in comparable manner in both time series. As the downsampling creates ‘sparser’ data and disrupts the smoothness of the original data (which might result in more abrupt changes in community assignments), we anticipated a tendency toward higher values for both metrics in the downsampled data. For node degree, we expected lower values in the downsampled data, as a recent study by Pedersen et al. (2018) showed that nodal degree decreased with less available data and thus tested this hypothesis using a paired t test.

We calculated the standard deviation (SD) of each node × node correlation (i.e., node-node connection or edge) over time as a proxy for TVFC (often also referred to as ‘dynamic’ connectivity; Lurie et al., 2020). To test whether these SDs likely reflect ‘true’ dynamics, we benchmarked these estimates against phase randomized surrogate (null) data that preserved autocorrelation, power spectral density, and stationary cross correlation of the observed data (Lurie et al., 2020; Prichard & Theiler, 1994; Savva et al., 2019). More precisely, we created 500 surrogates for each subject by phase randomizing the empirical time series obtained from the 100 nodes. To preserve the correlative structure between node time series, all signals were multiplied by the same uniformly random phase and the SD for each edge was calculated, respectively (see Savva et al., 2019, for a detailed description of phase randomization). The empirically observed SDs were ranked against this null distribution and p values were obtained by dividing the number of times SD_surr >= SD_real by the number of surrogates. Given that each subject has 4,950 unique connections, FDR correction (p < 0.05) was applied to reduce type I errors.

To describe ‘dynamics’ in nodal space, we calculated the nodal degree metric as the binary sum of significant (i.e., ‘dynamic’) connections for each node (Pedersen et al., 2018). Of note, when testing TVFC with surrogate null data, one needs to be cautious as the absence of significance (relative to the null model) does not necessarily indicate the absence of dynamic connectivity. This interpretation depends heavily on the null model applied, and the data may contain meaningful fluctuations relative to other definitions of the null model (see Lurie et al., 2020, for a comprehensive overview).

For multilayer modularity, negative correlations in each node × node matrix (i.e., per subject and time window) were set to zero and correlations were z-transformed as done in previous studies (e.g., Finc et al., 2020; Pedersen et al., 2018). We constructed an ordinal multilayer network in which each layer (i.e., time window) represents a weighted adjacency matrix. To assess the spatiotemporal community structure and to track network reconfigurations, each node was linked to itself across layers. To detect communities (i.e., groups of nodes that are more densely connected to one another than to the rest of the network), we used the multilayer counterpart (Mucha et al., 2010) of the modularity function proposed by Newman and Girvan (2004). To optimize the multilayer modularity function, we used an iterative and generative Louvain like algorithm (implemented with code from Jeub et al., 2011–2019; see Supporting Information for details). The tunable parameters γ and ω were held constant across layers and set to unity, as had been done in previous studies (e.g., Braun et al., 2015; Yin et al., 2020). As the modularity approach is not deterministic, we repeated it 100 times for each participant and chose the run yielding the highest modularity value as, for example, done in Finc et al. (2020). We also tested the influence of varying the intra- and interlayer parameters beyond the standard setting in the original data (see Supporting Information) but found no substantial influence on the reported results.

Within the time-varying community framework, we assessed node flexibility as the number of times a node changes its community assignment between adjacent layers, normalized by the total number of possible changes. Node flexibility can be interpreted as a metric that allows one to quantify reconfigurations of functional connectivity patterns that a brain region undergoes over time (Braun et al., 2015). To enrich the spatiotemporal description of a node, we further calculated node promiscuity, defined as the fraction of communities a node participates in at least once, across all layers. This metric allows one to quantify the distribution of a node’s connections over time, for example, whether high flexibility stems from a switching between two communities or an evenly distributed allegiance to a larger number of different modules (Garcia et al., 2018; Papadopoulos et al., 2016). The higher the nodal promiscuity, the more modules a node participates in at least once across time. For each of the three node-specific measures of dynamic network reconfiguration, that is, flexibility, promiscuity, and degree, we additionally calculated the respective global measure as the average across all nodes, as well as network-specific measures by averaging across all nodes belonging to the respective RSNs, which varied across atlases.

All correlational analyses were performed with the pingouin statistics package (0.3.11; Vallat, 2018) in python 3.8.3. To test for associations between resilience and time-varying brain network measures, we calculated partial (Spearman) correlations with age, gender, and FD included as covariates of no interest. Note that unlike Long et al. (2019), years of education was not included as covariate, as all participants were students and years of education were thus not acquired in our study. Correlations were performed at the global, nodal, and network-specific levels for all pairwise combinations between resilience questionnaires (CD-RISC, BRS, RS) and network measures (flexibility, promiscuity, degree). Results were corrected for multiple statistical comparisons by using the false discovery rate (p < .05; Benjamini & Hochberg, 1995).

Lastly, to more directly assess the effect of specific analysis choices on observed correlations between network measures and resilience, we varied a number of such factors systematically, using as ‘standard’ for comparison the analysis pipeline based on MB data parcellated using the Schaefer100 atlas, 36-parameter denoising, and overlapping sliding window scheme. We investigated the effects of including versus excluding covariates, varying sliding window size (in a frequently used range between 40 and 100 sec in steps of 20 sec; Preti et al., 2017), and how the two different denoising pipelines affect resilience correlations. We also explored within the same parcellation scheme (the Schaefer atlas) the effects of increasing the spatial resolution (i.e., from 100 to 200 nodes).

We thank Dr. Cindy Eckart for providing valuable comments on an earlier version of this manuscript and our research assistants for their help in data collection.

Supporting information for this article is available at https://www.doi.org/10.1162/netn_a_00216.

Dominik Kraft: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Software; Visualization; Writing – original draft; Writing – review & editing. Christian J. Fiebach: Conceptualization; Project administration; Writing – original draft; Writing – review & editing.

This research was funded by the German Research Foundation (CRC 1193 grant no. INST 247/859-1 awarded to C.J.F.). D.K. is supported by the German Academic Scholarship Foundation.

Functional connectivity:

Statistical dependencies among (fMRI) time series from different regions, nodes, or networks.

Graph theoretical framework:

Mathematical tool to model and analyze pairwise relations between different objects (here: time series from brain regions and their functional connectivity).

Downsampling:

Creating a ‘slower’, i.e., less densely sampled time series by artificially selecting every xth time point.

Multilayer modularity detection:

Detection of highly connected nodes (i.e., modules) in temporal, ordinal networks.

Time-varying functional connectivity:

Functional connectivity that varies over time, often labeled as ‘dynamic connectivity’.

Multiband EPI:

Fast echo planar imaging (EPI) sequence which is characterized by the simultaneous acquisition of multiple slices. Multiband EPI shortens the acquisition time of a single brain volume and thus allows to acquire more data at a given time.

Sliding window correlation:

FC calculation over a defined window, which is shorter than the whole time series. Window type may vary in its appearance. Sliding the window over the time series can produce time resolved functional connectivity (one matrix per window).