Material deprivation is associated with neural resilience for late-life depression

Depressive symptoms

A touchscreen questionnaire was implemented to collect information on sociodemographic characteristics and mental health. Recent depressive symptoms were measured using four items (RDS-4) that assess depressed mood, disinterest, restlessness, and tiredness in the past two weeks (Dutt et al., 2022). This continuous measure is highly comparable to several standardized self-report depression scales including PHQ-9, CES-D, and MASQ-30 (Dutt et al., 2022), and shows an area under the curve of 0.79 for its correlation with depression diagnosis (Khubchandani et al., 2016). In response to each of these four questions, participants indicated their experiences from “not at all” (scoring 1) to “nearly every day” (scoring 4) such that the total symptom score ranged from 4 to 16.

Imaging preprocessing and multimodal imaging features

All brain imaging data from the UKB were acquired using standard Siemens Skyra 3T scanners with a standard Siemens 32-channel RF receive head coil. Detailed information for UKB imaging acquisition can be found in the UK Biobank Imaging Documentation, hosted on the Oxford FMRIB UK Biobank resource page (https://www.fmrib.ox.ac.uk/ukbiobank/). A fully automated processing and QC pipeline was developed for UKB brain imaging data, which included T1, T2, FLAIR, susceptibility-weighted MRI, resting-state MRI, task-evoked MRI and diffusion MRI (Alfaro-Almagro et al., 2018). Additionally, this pipeline also generated a set of imaging-derived phenotypes (IDPs) such as cortical and subcortical structure volumes, microstructural measures in major tracts, and functional connectivity metrics. In this study, a total of 4,632 IDPs was included as multimodal imaging features in the prediction models for estimating brain-depression score (see section 3.1). Approximately three quarters of these features were derived from the resting-state and task-evoked fMRI data, and the remaining from structural and diffusion MRI data:

• – Resting-state fMRI features: Amplitude of ICA100 nodes and ICA25 nodes; Edges of full correlation matrix from ICA100 and ICA25; Edges of partial correlation matrix from ICA100

• – Task-evoked fMRI features: Median and 90th percentile BOLD signals for shapes and faces, as well as shape-face contrasts using a group mask and an amygdala mask respectively; Median and 90th percentile Z-statistics for shapes and faces, as well as shape-face contrasts using a group mask and an amygdala mask respectively

• – Structural MRI features: Volumes of cortical and subcortical (including sub-segments) structures; Total volume of white matter hyperintensities

• – Diffusion MRI features: Mean FA, MD, MO, L1-L3, ICVF, OD, ISOVF based on Standardized FA Skeleton; Weighted-mean FA, MD, MO, L1-L3, ICVF, OD, ISOVF in White-matter tracts

Late-life Depression Risk factors

A list of risk factors for late-life depression (LLD) was curated based on the literature that covers information on demographics, lifestyle, medical conditions, adverse experiences, and psychosocial factors. Measures of these variables were collected from participants either via a touchscreen questionnaire or a verbal medical history interview on the same day of imaging acquisition. Specifically, several variables such as self-reported health, long-standing illness and Townsend deprivation index, which is a census-based index of material deprivation calculated by the combination of four indicators of deprivation: non-home ownership, non-car ownership, unemployment and overcrowding (Townsend, Phillimore and Beattie, 1988) were available directly from the UKB, whereas aggregate measures for healthy lifestyle, sleep quality, vascular risk factors, adverse or traumatic experiences, social isolation and loneliness were derived using multiple items as follows:

• – A healthy lifestyle score was constructed based on smoking status, physical activity, diet, and alcohol consumption that are well documented as depression risk factors (Sarris et al., 2020; van Lee et al., 2020; Kang et al., 2021). Based on national recommendations, participants were given 1 for healthy, and 0 for unhealthy behaviors. Detailed coding information can be found in (Lourida et al., 2019).

• – A sleep (low) risk score was calculated using five sleep questions. Low-risk sleep factors were defined as a) having an early chronotype, b) sleeping 7–8 hours per day, c) never having or rarely having insomnia symptoms, d) not reporting snoring, and e) not reporting frequent daytime sleepiness (Fan et al., 2020). Participants received a score of 1 if their behaviors were classified as low risk for that factor and a sum score was calculated across five factors, where higher scores represent healthier sleep patterns or low risks of sleep issues (Hepsomali and Groeger, 2021).

• – An aggregate measure of vascular risk factors was calculated for each participant by counting instances of having a high BMI (>25), having a high waist-hip ratio (WHR>0.85 for females and WHR>0.90 for males), having ever smoked, and a self-reported diagnosis of hypertension, diabetes, or hypercholesterolemia (Cox et al., 2019). The resulting sum score of instances represent potential vascular risks, where higher scores indicating higher risks.

• – scores were calculated separately for childhood, adulthood, and lifetime experiences with dichotomization of responses to each individual question (Yapp et al., 2021). Specifically, responses of “never true” to negative experiences (e.g., hit hard) scored a 0 and responses of “rarely true” and more often (i.e., “sometimes true”, “often”, “very often true”) scored a 1, which was reversed for positive experiences (e.g., in a confiding relationship). Binary responses (“yes”, “no”) were scored 1 and 0 respectively. Separate sum scores were calculated to indicate the magnitude of traumatic experiences in different time periods and included as separate predictors in the models.

• – Psychosocial factors comprised two measures: loneliness and social isolation. Participants were classified as lonely if they reported feeling lonely often and if they could confide to someone close only occasionally (e.g., less than once every few months), and socially isolated if they met at least two criteria of 1) living alone, 2) visiting their family or friends less than once a month, and 3) participating in none of the listed leisure/social activities (Mutz, Roscoe and Lewis, 2021). These psychosocial factors were included in statistical models as separate predictors.

Estimating the delta of brain-depression score

Our approach was inspired by the estimation of brain-age delta from neuroimaging features that is defined as the difference between the estimated brain age and chronological age in a given individual, which has been used to indicate underlying problems in outwardly healthy people and related to the risk of cognitive ageing or age-associated brain disease (Franke et al., 2010; Cole and Franke, 2017; Baecker et al., 2021). In this study, we employed Multivariate Adaptive Regression Splines models (Friedman and Roosen, 1995; Friedman, 2007) to estimate the brain-depression score (BDS) as informed by multimodal imaging features and quantified the difference between the estimated and the actual depression symptom scores (i.e., ΔBDS), with a positive ΔBDS indicating neural vulnerability to depression (i.e., actual reported symptom score higher than BDS) and negative ΔBDS indicating neural resilience to depression (i.e., actual symptom score lower than BDS).

Multivariate Adaptive Regression Splines (MARS) is a flexible regression technique that can capture the intrinsic nonlinear and multidimensional relationship of variables with an ensemble of linear functions joined together by one or more spline basis functions, where the number of basis functions and the parameters associated with each function (e.g., product degree and knot locations) are determined by the data (Friedman and Roosen, 1995; Friedman, 2007). Specifically, MARS builds a model in two phases: the forward pass and the backward pruning, similar to growing and pruning of tree models. In the forward pass, MARS starts with a model consisting of just the intercept term (i.e., the mean of the response values), followed by the assessment of every single predictor to find a basis function pair that produces the maximum improvement in the model error. This process iterates until either the model reaches a predefined limit number of terms, or the error improvement reaches a predefined limit. The result of the forward pass is a MARS basis matrix with rows of observations and columns of basis functions (e.g., hinge functions). To avoid overfitting by the full terms in the basis matrix from the forward pass, the backward pass is to find the subset of these terms that gives the best generalized cross validation (GCV) via a stepwise term deletion procedure. This backward pruning process continues until only one term remains (the intercept term) and the final model with the best GCV is selected (Friedman and Roosen, 1995; Friedman, 2007).

In this study, IDPs from multimodal imaging data and depressive symptom sum scores were included in MARS models as independent (X) and dependent (Y) variables. Additionally, age, sex, head size, head motion (i.e., the averaged head motion across space and time points) during fMRI acquisition (i.e., for both the resting-state and task-evoked fMRI), scanner site, scanner table position and data acquisition dates were included in all MARS models as confounding variables. Per partition of the full data in the nested cross-validations (see details below), principal component analysis (PCA) was employed to decompose high-dimensional X and components collectively explaining more than 50% variance were retained. As the outcome measure Y (i.e., symptom sum score) in our study is highly skewed with a long tail, transformation of Y was performed per data partition before modeling to enforce Gaussianity. This was realized by using a data-drive approach that finds the optimal normalization method from a suite of possible transformation options including the Box-Cox transformation, Yeo-Johnson transformation, the ordered quantile technique, Arcsinh transformation, exponential transformation, square root transformation and the Lambert W x F transformation (Peterson, 2021). We also applied dummy coding for all factorial variables such as sex and scanner site before data partitioning to ensure the same number of predictors across models.

Nested cross-validations were applied to increase robustness and generalizability of our MARS estimations, with 3 iterations in the outer loop and 10 folds per iteration in the inner loop. As mentioned earlier, data processing including PCA (on X) and transformation (on Y) was performed within each iteration of the outer loop and parameters obtained from the training and validating set were applied to the held-out testing set. Per training fold, a grid search was performed to identify the optimal combination of two parameters: the maximum degree of interactions among terms (ndegree) and the number of terms retained in the final model (nprune). Given the available data points and the empirical evidence that 3^rd-degree interactions never benefited model fit using a subset of data, expansion of interaction degree was restricted to 2 (e.g., testing ndegree = 1 and = 2) and nprune up to 5 times of the predictor numbers (i.e., varied between 690 and 695 due to changing numbers of principal components in each fold). Per prediction model on the testing set, we further partitioned data to estimate prediction uncertainty as the noise and model variance of the out-of-fold predicted values over 5 iterations of 50 cross-validations, using variance models of MARS (Milborrow, 2015), and an averaged r-squared of 0.19 in contrast to 0.09 from the training models.

For each participant, we quantified ΔBDS as the difference between the predicted BDS (ŷ) and the reported symptom sum score (Y) while accounting for the prediction uncertainty including irreducible “aleatoric” or noise variance and model variance of prediction or “epistemic” uncertainty such that . By definition, the resulting ΔBDS for each participant indexed less depressed or resilient brain patterns if the quantity was negative, and more depressed or vulnerable brain patterns if positive. Here, we focused on individuals with high neural resilience (i.e., top 10 percentile of the negative ΔBDS) and individuals with high neural vulnerability (i.e., top 10 percentile of the positive ΔBDS). This was determined to drive a good balance between sufficient statistical power (i.e., n= 144 with 11 regressors) and high clinical relevance, and to reduce dependency of BDS as a function of original symptom score (i.e., minimize high correlation between ΔBDS and original symptom score).

Testing risk factor effects

All analyses were conducted for individuals with high neural resilience and vulnerability separately. We first conducted partial rank correlation analyses to test the associations between each risk factor and the ΔBDS while controlling for the reported symptom sum score (i.e., resulting in effects not influenced by the sum score). This is a common practice to deconfound the associative effects driven by chronological age rather than the brain-age delta in the brain age literature (Smith et al., 2019). False Discover Rate (FDR) corrections were applied to adjust for multiple testing. We further conducted quantile regression analyses to identify unique contributions of individual risk factors, while accounting for the effects of reported symptom sum score.