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

Abstract This study aimed at replicating a previously reported negative correlation between node flexibility and psychological resilience, that is, the ability to retain mental health in the face of stress and adversity. To this end, we used multiband resting-state BOLD fMRI (TR = .675 sec) from 52 participants who had filled out three psychological questionnaires assessing resilience. Time-resolved functional connectivity was calculated by performing a sliding window approach on averaged time series parcellated according to different established atlases. Multilayer modularity detection was performed to track network reconfigurations over time, and node flexibility was calculated as the number of times a node changes community assignment. In addition, node promiscuity (the fraction of communities a node participates in) and node degree (as proxy for time-varying connectivity) were calculated to extend previous work. We found no substantial correlations between resilience and node flexibility. We observed a small number of correlations between the two other brain measures and resilience scores that were, however, very inconsistently distributed across brain measures, differences in temporal sampling, and parcellation schemes. This heterogeneity calls into question the existence of previously postulated associations between resilience and brain network flexibility and highlights how results may be influenced by specific analysis choices.


Supplementary File 1 MRI data preprocessing
After reconstruction BOLD fMRI data were organized using the Brain Imaging Data Structure (BIDS; Gorgolewski et al., 2016) and preprocessed using fMRIPrep version 20.1.1 (Esteban et al., 2018a(Esteban et al., , 2018b; RRID:SCR_016216) which is based on Nipype 1.5.0 (Gorgolewski et al., 2011;Gorgolewski et al., 2018;RRID:SCR_002502). Many internal operations of fMRIPrep use Nilearn 0.6.2 (Abraham et al., 2014; RRID:SCR_001362) mostly within the functional processing workflow; all other software tools used are specified below. For more details of the pipeline see the section corresponding to workflows in fMRIPrep's documentation. In the following, preprocessing is described in more detail by text that was automatically generated by the fMRIPrep software and then manually edited with respect to readability.

Functional data preprocessing
For each subject's BOLD run the following preprocessing was performed: First, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep: Head-motion parameters with respect to the BOLD reference (transformation matrices and six corresponding rotation and translation parameters) were estimated before any spatiotemporal filtering using MCFLIRT (FSL 5.0.9;Jenkinson et al., 2002). A B0nonuniformity map (or fieldmap) was estimated based on a phase-difference map calculated with a dual-echo GRE (gradient-recall echo) sequence processed with a custom workflow of Susceptibility Distortion Correction workFlows (SDCFlows; Markiewicz et al., 2020) inspired by the epidewarp.fsl script and further improvements in Human Connectome Project Pipelines (Glasser et al., 2013). The fieldmap was then co-registered to the target EPI (echo-planar imaging) reference run and converted to a displacements field map (amenable to registration tools such as ANTs) with FSL's FUGUE and other SDCflows tools. Based on the estimated susceptibility distortion a corrected EPI (echo-planar imaging) reference was calculated for a more accurate co-registration with the anatomical reference. The BOLD reference was then coregistered to the T1w reference using bbregister (FreeSurfer) which implements boundarybased registration (Greve & Fischl, 2009). Co-registration was configured with six degrees of freedom. The BOLD time-series were resampled onto FreeSurfer's the fsnative and fsaverage surfaces. The BOLD time-series were resampled onto their original native space by applying a single composite transform to correct for head-motion and susceptibility distortions. These resampled BOLD time-series will be referred to as preprocessed BOLD in original space or just preprocessed BOLD.
The BOLD time-series were resampled into standard space generating a preprocessed BOLD run in MNI152NLin2009cAsym space. First a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep: Several confounding time-series were calculated based on the preprocessed BOLD, i.e., framewise displacement (FD), derivative of RMS variance over voxels (DVARS), and three region-wise global signals. FD was computed using two formulations following Power (absolute sum of relative motions; Power et al., 2014) and Jenkinson (relative root mean square displacement between affines; Jenkinson et al., 2002). FD and DVARS were calculated for the functional run using their implementations in Nipype (following the definitions by Power et al., 2014). The three global signals were extracted within the CSF, WM, and whole-brain masks. Additionally, a set of physiological regressors were extracted to allow for component-based noise correction (CompCor; Behzadi et al., 2007). Principal components are estimated after high-pass filtering the preprocessed BOLD time-series (using a discrete cosine filter with 128s cut-off) for the two CompCor variants: temporal (tCompCor) and anatomical (aCompCor). TCompCor components are then calculated from the top 5% variable voxels within a mask covering the subcortical regions. This subcortical mask is obtained by heavily eroding the brain mask, which ensures it does not include cortical GM regions. For aCompCor, components are calculated within the intersection of the aforementioned mask and the union of CSF and WM masks calculated in T1w space after their projection to the native space of the functional run (using the inverse BOLD-to-T1w transformation). Components are also calculated separately within the WM and CSF masks. For each CompCor decomposition, the k components with the largest singular values are retained such that the retained components' time series are sufficient to explain 50 percent of variance across the nuisance mask (CSF, WM, combined, or temporal).
The remaining components are dropped from consideration. The head-motion estimates 5 calculated in the correction step were also placed within the corresponding confounds file. The confound time series derived from head motion estimates and global signals were expanded with the inclusion of temporal derivatives and quadratic terms for each (Satterthwaite et al., 2013). All resamplings can be performed with a single interpolation step by composing all the pertinent transformations (i.e., head-motion transform matrices. susceptibility distortion correction, and co-registrations to anatomical and output spaces). Gridded (volumetric) resamplings were performed using antsApplyTransforms (ANTs) configured with Lanczos interpolation to minimize the smoothing effects of other kernels (Lanczos, 1964). Non-gridded (surface) resamplings were performed using mri_vol2surf (FreeSurfer).

Multilayer Modularity details
In general, community detection and node assignment can be solved algorithmically by maximizing the modularity quality function Q. Q allows to quantify the goodness of fit of a certain partition, with high values indicating high quality partitions (Furtunato, 2010).
Following the idea that random graphs do not show a cluster structure, the quality function compares the number of intracommunity edges to what would be expected at random (i.e., under a null model), to quantify communities (Fortunato, 2010;Mucha et al., 2010). The quality function Q can be formalized as follows: where !" and !" are the observed and expected weights of the connection between node i and node j. Here, !" is an element of the adjacency matrix and !" refers to a specific null model.
The Kronecker delta function % ! , " ( equals 1 if node i and j are assigned to the same module (M) and 0 if otherwise. Noteworthy, the modularity function suffers from a "resolution limit", meaning that communities that are rather small compared to the whole graph may be undetectable (Fortunato & Barthelemy, 2007). To overcome this caveat, one may introduce a structural resolution parameter that weights the null model. Optimizing ( ) with large values of leads to the detection of many small communities, whereas small values of result in few large communities. Here, we used an expanded modularity function that can be used with multilayer networks, formalized as follows: where !"# describes elements of of layer (time window) s and denotes the Newman-Girvan configuration null model (i.e., the expected connectivity strength at chance level, with k = node degree at layer s; m = sum of all edges at layer s) that is scaled by the resolution parameter # of layer s, % !# , "# ( and % !# , "$ ( = 1 if node i and node j are in the same module (M) or 0 if they do not belong to the same module (here !# depicts the community assignment of node i in layer s; "#$ is the interlayer coupling parameter connecting node j in layer s to itself in layer r, defining the similarity of communities across layers with large values of "#$ leading to increased homogeneity across layers, as nodes will prefer to stay in one community (Vaiana & Muldoon, 2018).

Correlation results on a nodal level without covariates.
Figures S1 and S2 depict the correlations between brain network measures on a nodal level and resilience questionnaires calculated without covariates. These figures can be compared to   8

Node -network assignment in the AAL, Power264 atlas
To additionally test whether results reported here might have been driven by the specific parcellation we used (i.e., Schaefer 100), we repeated the pre-processing and denoising pipeline with two additional atlases, i.e., the AAL parcellation with 90 regions (Tzourio-Mazoyer et al., 2002) and the atlas of (Power et al., 2011) with 264 regions, both of which had been used in the study by Long and colleagues (2019). For the AAL90 atlas, we referred to the node -RSN assignment provided by Long et al. (2019) in their Table S1. However, according to that Table, some nodes were assigned to multiple subnetworks (e.g., node middle frontal gyrus), which we avoided. Accordingly, our subnetwork assignment may vary in details from the one used in Long et al. (2019). Taking Long et al.'s (2019) Table S1 as a starting point, we assigned each node to the subnetwork which appears first in the respective column, leading to a total number of eight unique networks, i.e., 'Sensorimotor', 'Visual', 'Salience', 'Auditory', 'Subcortical', 'Frontoparietal', 'Default-mode', and 'Cingulo-opercular'. 12 nodes were not assigned to any sub-network (i.e.,'None'), including the amygdala and hippocampus. For the Power atlas (nodes 1,3,8,9,119,126,183,and 184) from further analyses, as some participants were lacking BOLD signal for these nodes. The remaining nodes were assigned to one of the following resting state networks (RSN): 'somatomotor hand', 'somatomotor mouth', 'cingulo-opercular', 'auditory', 'default mode', 'memory retrieval ', 'visual', frontoparietal', 'salience', 'subcortical', 'ventral attention', 'dorsal attention', or 'cerebellar'. Analogous to the AAL atlas, 21 nodes that could not be assigned to any of these categories were sorted into a 'None' category. For both parcellations, we again calculated partial correlations between dynamic brain measures and the respective resilience questionnaires on a nodal, RSN, and global scale.

Effect of different intra-and interlayer parameters on measures of brain dynamics
The community detection in multilayer networks is tuned by the intra-and inter-layer parameters and , respectively. Whereas the intralayer parameter determines how many communities will be detected (i.e., larger values leading to the detection of smaller communities) at a given timepoint, the interlayer resolution parameter scales the network dynamics over time, i.e., larger values leading to increased homogeneity over time and nodes consequently tend to switch communities less often (Vaiana & Muldoon, 2018). We here investigate the stability of our correlation results for neural flexibility and promiscuity measures, by systematically varying these parameters. In line with previous studies (e.g., Pedersen et al., 2018;Yin et al., 2020) we considered intralayer parameters ( ) from 0.8 to 1.4 and interlayer parameters ( ) from 0.4 to 1.0 in steps of 0.2, while fixing the respective other parameter to 1 (using MB data with TR = .675 and the Schaefer100 atlas). The detection of smaller communities (i.e., larger values of ) led to higher nodal flexibility (see Figure S3), whereas it led to smaller nodal promiscuity (i.e., nodes tend to change their community assignments more often over time, but these interactions are not distributed across all other modules but appear to be rather limited to a small number of communities. The spatial similarity (i.e., correlations) of nodal flexibility values (range r = .52 to .92, all p < .0001) and promiscuity (range r = .40 to .92, all p < .0001) across the range of parameters studied suggest a moderate to good correspondence across values (see Figure XX). Across layers, lower values of the interlayer parameter led to increased node flexibility, whereas node promiscuity showed the opposite pattern, i.e., decreasing values with lower omega (see Figure S3). However, despite numerical changes, the spatial similarity of node flexibility (r > .95, p < .0001) and promiscuity (r > .96, p < .0001) for all parameters of was high. Both patterns suggest that nodes switch between communities more often over time, but that this increase is driven by a constant switching between a small number of communities and not a broad interaction with various communities (i.e., simultaneous low promiscuity with high flexibility). More generally, these results demonstrate relative robustness of node flexibility and node promiscuity measures against variations in the parameters gamma and omega.

Effect of different intra-and interlayer parameters on the correlation between measures of brain dynamics and psychological resilience
While varying the interlayer parameter from 0.8 to 1.4 in steps of 0.2 and computing nodal flexibility and promiscuity based on the original data (TR= .675s; Schaefer100 atlas), we did observe one significant correlation between flexibility of a visual node (LH_Vis_7) and the BRS (r = .50, p = .02) when using a parameter value of = 1.4. For the remaining two levels of the intralayer parameter (note that we do not report = 1 here, as it refers to the initial analysis in the manuscript) we did not observe any significant correlations. The significant correlation indicates that higher flexibility is associated with higher self-reported resilience, which is directly opposed to what has been reported before (e.g., Long et al., 2019;Paban et al., 2019). Furthermore, across levels of , there were no significant correlations between node promiscuity and the resilience questionnaires (all p > .35). We did not observe any significant correlation between node flexibility (all p > .10) or node promiscuity (all p > .11) and resilience while setting to either 0.4, 0.6, or 0.8.