Moderating Effect of Cortical Thickness on BOLD Signal Variability Age-Related Changes

The time course of neuroanatomical structural and functional measures across the lifespan is commonly reported in association with aging. Blood oxygen-level dependent signal variability, estimated using the standard deviation of the signal, or “BOLDSD,” is an emerging metric of variability in neural processing, and has been shown to be positively correlated with cognitive flexibility. Generally, BOLDSD is reported to decrease with aging, and is thought to reflect age-related cognitive decline. Additionally, it is well established that normative aging is associated with structural changes in brain regions, and that these predict functional decline in various cognitive domains. Nevertheless, the interaction between alterations in cortical morphology and BOLDSD changes has not been modeled quantitatively. The objective of the current study was to investigate the influence of cortical morphology metrics [i.e., cortical thickness (CT), gray matter (GM) volume, and cortical area (CA)] on age-related BOLDSD changes by treating these cortical morphology metrics as possible physiological confounds using linear mixed models. We studied these metrics in 28 healthy older subjects scanned twice at approximately 2.5 years interval. Results show that BOLDSD is confounded by cortical morphology metrics. Respectively, changes in CT but not GM volume nor CA, show a significant interaction with BOLDSD alterations. Our study highlights that CT changes should be considered when evaluating BOLDSD alternations in the lifespan.


INTRODUCTION
Normal aging is associated with marked functional and structural neuroanatomical alterations in cortical thickness (CT), gyrification, cortical surface area (CA), gray (GM), and white matter volume (WM) (Salat et al., 2009;Thambisetty et al., 2010;McGinnis et al., 2011;Hogstrom et al., 2013). Importantly, magnetic resonance imaging studies (MRI) show that the magnitude and rate of change of these cortical morphometry metrics is not constant across the cortex but rather it varies with age and brain region (Raz, 2005;Jiang et al., 2014;Storsve et al., 2014;Dotson et al., 2016), and is reported to accelerate with increasing age (Driscoll et al., 2009;Jiang et al., 2014). For example, a longitudinal study of alterations in cortical morphometry in older adults found accelerated changes with increasing age in temporal and occipital cortices (Storsve et al., 2014). Furthermore, other studies report that the temporal lobes are most vulnerable to age-related morphometric changes, and that these changes reflect age-related cognitive impairment (Singh et al., 2006;Fjell et al., 2009;Pacheco et al., 2015). There is considerable evidence that neuroanatomical alterations reflect underlying functional alterations, especially in cognition (Fjell et al., 2006;Rossini et al., 2007;Ziegler et al., 2010;Achiron et al., 2013). In fact, functional magnetic resonance imaging (fMRI) studies, which rely on the blood oxygen level-dependent (BOLD) signal as a correlate of neuronal activity, report that changes in cortical morphology across adult lifespan impact the hemodynamic properties of the brain. For example, there are cortical laminar differences in BOLD signal (Tian et al., 2010;Huber et al., 2015), thicker cortical regions were reported to have a lower relative oxygen extraction fraction (Zhao et al., 2016). Therefore, since aging is associated with significant neuroanatomical alterations, these should be considered when assessing function (i.e., cognitive ability) using the BOLD signal.
The Standard Deviation of the BOLD signal can be used to estimate variability (hereafter, "BOLD SD "), and is believed to reflect the brain's dynamic ability to undergo fast moment-to-moment switching through network reconfigurations (Garrett et al., 2010Grady and Garrett, 2014). It is an emerging index of cognitive health in aging, with higher regional BOLD SD being associated with enhanced performance on certain cognitive tasks (i.e., task switching) but not on others (i.e., distractor inhibition) (Armbruster-Genc et al., 2016;Guitart-Masip et al., 2016). Generally, increased BOLD SD is associated with younger age, faster and more consistent performance on cognitive tasks, and cognitive flexibility (Garrett et al., 2013a;Armbruster-Genc et al., 2016). The physiological mechanisms underlying BOLD SD remain largely unknown. For instance, decreased dopaminergic transmission is proposed to be associated with decreased BOLD SD in subcortical areas in older adults compared to younger ones (Guitart-Masip et al., 2016). Another study reported that higher BOLD SD is associated with superior pain coping capabilities (Rogachov et al., 2016). However, there are few studies investigating the interaction between age-related alterations in cortical morphology and BOLD SD . One study reported that increased microstructural integrity of WM pathways is associated with greater BOLD SD (Burzynska et al., 2015). Given that neuroanatomical alterations associated with normative aging are known to influence cognitive performance, and that they influence the BOLD signal, their impact as physiological confounds to BOLD SD should be investigated. This is particularly relevant because there is a considerable degree of inconsistency of methods used and results across studies investigating BOLD SD . In fact, some studies report greater regional BOLD SD in older adults (Garrett et al., 2010(Garrett et al., , 2011Nomi et al., 2017), individuals with stroke (Kielar et al., 2016), multiple sclerosis (Petracca et al., 2017), Alzheimer disease (Makedonov et al., 2013(Makedonov et al., , 2016Scarapicchia et al., 2018) and other neurological disorders (Zöller et al., 2017).
The objective of this study is to investigate the contribution of cortical morphology (i.e., CT, CA, and GM volume) to age-related BOLD SD changes in older adults. A longitudinal framework, consisting of two scan points, should help reduce some of the inter-individual variance in neuroanatomy by accounting for external factors such as lifestyle, and various socio-economic and demographic factors. To this end, we investigated global and regional differences in BOLD SD in a group of older adults scanned twice at an interval of approximately 2.5 years, and regressed the effect of cortical morphology by introducing CT, CA, and GM as covariates. We hypothesized that cortical morphology metrics show an interaction with BOLD SD . We predict that age-related neuroanatomical alterations in CT, CA, GM are physiological confounds to BOLD SD measures, and that consequently adjusting for these metrics may help "unmask" the functional relevance of BOLD SD .

Participants
All data obtained for the present study were obtained from the longitudinal Geneva Aging Study, after approval by the ethics committee of the Faculty of Psychology and Educational Sciences of the University of Geneva and the Swiss Ethic committee. Older subjects were initially recruited either from the University of the Third Age of Geneva or through newspaper and association advertisements for pensioners, as part of a larger longitudinal study. All participants gave written informed consent and older adults received a small amount of money as a compensation for their transportation fees.
Our initial sample consisted of 31 older adults scanned twice (mean age at first scan = 71.65 ± 6.03 years, mean age at second scan = 74.06 ± 5.99 years; 9 males). These subjects were chosen, within our pool, because they were the only ones that had undergone two T1 structural images and task fMRI scans, as well as other cognitive tests. Participants were screened for health problems with a questionnaire. The structural MRIs were inspected to rule out severe abnormalities (white matter changes, ventricular enlargement, tumors etc.). Three of the participants showed signs of Parkinson or lesion on their anatomical MRI, so they were excluded from the final sample (n = 28). All models and results in this current paper thus reflect 28 older adults (mean age at first scan = 71.61 ± 6.21 years, mean age at second scan = 74.07 ± 6.15 years; 7 males). The scans were 2.46 ± 0.69 years apart.

MRI Data Acquisition
Participants were scanned in a Siemens Trio 3T magnet. A BOLD fMRI task-rest sequence was administered using a reaction time paradigm, where the participant had to indicate on which side a cross was changing into a square, as fast as possible (Mella et al., 2013). The task consisted of eight experimental blocks, interspersed with eight resting blocks (respectively, 52-20 s). The BOLD activity was obtained using an echo planar imaging acquisition (echo time, TE = 30 ms, time repetition, TR = 2100 ms, flip angle = 80 • , field of view, FOV = 205 mm). Then, a structural T1-weighted MRI was acquired (TE = 2.27 ms, TR = 1900 ms, FOV = 256 mm, voxel size 1.0 mm × 1.0 mm × 1.0 mm).

Structural Data Processing
Structural T1-weighted MR images were analyzed using Freesurfer version 6.0, a widely used and freely available automated processing pipeline 1 , which allows surface-based three dimensional reconstruction and quantification of cortical morphology. The standard steps for analysis were implemented (using "recon-all" pipeline with the default set of parameters). Regional measures of GM, CT, and CA for each hemisphere were obtained using the automated anatomic parcellation procedure. Technical details are found in prior publications (Dale et al., 1999;Fischl et al., 1999aFischl et al., ,b, 2001Fischl and Dale, 2000;Zhang et al., 2001;Salat et al., 2004;Winkler et al., 2012). In brief, T1-weighted images underwent preprocessing steps including motion correction, brain extraction, intensity normalization, and Talairach transformation (Sled et al., 1998;Smith, 2002;Ségonne et al., 2004;Reuter et al., 2010). GM and WM surface boundaries were reconstructed to estimate the distance between them across the cortex (Fischl and Dale, 2000). The generated cortical models were inflated into spheres to be registered to a spherical atlas and parcellated into regions of interest using Destrieux atlas (Fischl and Dale, 2000;Segonne et al., 2007;Destrieux et al., 2010).

Functional Data Processing
Processing of the functional data were performed using FSL version 5.0 (Analysis Group, FMRIB, Oxford, United Kingdom . Standard preprocessing were followed using FSL's FEAT and FSL's Melodic for functional data (Jenkinson et al., 2012). Briefly, for each participant preprocessing steps included motion correction, slice timing, spatial normalization, highpass temporal filtering (100 s), smoothing (kernel 5 mm FWHM), and linear registration (12 degrees of freedom) of the functional data to the high-resolution T1 structural image, and from T1 to 1 mm standard space (MNI 152). Additionally, FSL's Melodic was used to correct each functional image for artifacts using automatic dimensionality estimation via independent component analysis (ICA) . Change in cortical morphology between the two scanning sessions was determined as cortical morphology metrics GM (delta.volume), CT (delta.thickness), CA (delta.area) at timepoint 2-at timepoint 1.

BOLD Signal Variability Calculation
As part of the Geneva Dataset, the subjects were performing different cognitive tasks, and for the current study, only the fixation/rest blocks from the block design task fMRI were selected to calculate BOLD SD, as previously described (Garrett et al., 2010). BOLD SD analysis was restricted to the GM using participant specific GM mask obtained from FSL's FAST. First, fixation blocks were normalized so that the overall four-dimensional mean across brain and block was 100. Next, for each voxel, the block mean was subtracted to remove block-wise drift, followed by concatenation of all blocks. The standard deviation of the normalized mean of the concatenated fixation blocks was used to obtain BOLD SD values for each brain region (n = 148) of each subject, as 1 http://surfer.nmr.mgh.harvard.edu/ defined by the Destrieux Atlas (Destrieux et al., 2010), using in-house MATLAB code. Change (delta.variability) in variability between the two scanning sessions was calculated as variability at timepoint 2-variability at timepoint 1. BOLD SD encompassing both timepoints was introduced as "variability" (see section "Regional Model").

Statistical Analysis
Linear mixed effects models (LMMs) were used to investigate the potential confounding effect of cortical morphology, GM, CT, CA on BOLD SD (Bates, 2005;Zuur et al., 2011).The LMMs allow estimation of the effects of explanatory variables ("fixed effects") and their interactions on the dependent variable (i.e., BOLD SD ), while statistically controlling for the effects of randomly selected participants ("random effects") on the dependent variable (BOLD SD ). Multiple models were run and the likelihood-ratio test was used to (1) investigate if introducing subjects as random effects improves the fit of the model (2) to select the optimal combination of fixed effects fitted with maximum likelihood, while keeping the random effects structure the same. Therefore, the likelihood-ratio test via ANOVA was used to compare the goodness of fit of different models. R statistical software package (R Core Team, 2013) 2 was used for all statistical analyses. Correlation between cortical morphology measures were computed using cor function in R. All models were fitted using the "lmer" or "lm" function in R. "lmerTest" R package was used to obtain summary table and p-values for linear mixed models via Satterthwaite's degrees of freedom method (Kuznetsova et al., 2017). A spatiotemporal approach to LMMs allowed characterization of regionally specific variation across the brain (Bernal-Rusiel et al., 2013). This approach was implemented to investigate if there is a significant change in BOLD variability across time (1) all cortical regions (2) region specific. Random effect structure with subjects varying in their "baseline" BOLD variability was retained (random intercept, 1| ID). Additionally, to model a different rate of change in the expected response levels, time varying predictors were introduced by random slopes (i.e., thickness, time) ("Regional Model"). To reduce spatial correlation issues an LMM model with the described structure was applied at each spatial location (i.e., region of interest) independently. Each model produced a parameter which quantifies the mean change in BOLD variability (delta.variability) for that region. P-values were corrected for multiple comparisons using false discovery rate (FDR) at q = 0.05 (Benjamini and Yekutieli, 2001

Relations Between Cortical Thickness and BOLD Signal Variability Age-Related Changes
The functional and structural data of 28 participants was assessed. As expected, our study showed that CA and GM are highly correlated r = 0.904, p < 0.001 (95% Cl 0.901-0.908) and consequently collinearity was suspected. Multiple LMMs were utilized to assess the effect of neuroanatomical metrics CT, CA, and GM on BOLD SD age-related changes. Results from the linear mixed effect models run with likelihood-ratio test via ANOVA, are presented in Table 1. (A) indicates that subject intercept should be introduced as a random effect, while (B-E) show the steps that have led to the final model. Specifically , Table 1. (B,C) indicate that introducing CT or CA as covariates, separately, each significantly improve the fit of the model p < 0.0001, p < 0.05, respectively. However, from (D) it is apparent that adding CA to a model that already has CT as a covariate does not improve the fit of the model, meaning that CT only should be included in the final model (see section "Final Model"). (E) indicates that introducing GM as a covariate does not improve the fit of the model. Neither CT, GM, nor CA mean changes were significant as tested with LMMs. Regional Model LMM indicated that neither overall nor regionally specific mean BOLD SD change was significant after FDR correction.

Cortical Thickness and Its Association With BOLD SD
In this study we aimed at determining the contribution of cortical morphology to BOLD SD in order to better understand the physiological nature of brain function capture by BOLD SD . BOLD SD changes across the lifespan have been shown to be robust to certain vascular factor such as cerebral blood flow, BOLD cerebrovascular reactivity, maximal BOLD signal change (Garrett et al., 2017), as well as GM volume changes (Nomi et al., 2017). However, using LMMs, we found that functional alterations in aging as captured by BOLD SD are confounded by the structural metric CT. This is not surprising, considering that it is well known that BOLD signal change/activation is dependent on the laminar organization of the cortex, and consequently is influenced by its depth or thickness (Koopmans et al., 2010;Tian et al., 2010;Bandettini, 2012). In fact, neurovascular coupling is reported to vary by cortical depth and layer (Goense et al., 2012). In a longitudinal study investigating CT, GM volume, CA across the lifespan, Storsve et al. (2014) reported cortical morphology metric specific and region specific rates of mean annual percentage change (APC) in healthy adults aged 23-87 years. In most regions, GM volume has a mean APC of -0.51, CT of -0.35, and CA of -0.19. Other longitudinal studies in older adults, report similar reductions: in GM volume, mean APC ranging from -0.5 to -2.1 ± 1.6% (Tang et al., 2001;Fjell et al., 2009), in CT, mean APC -0.3% (Shaw et al., 2016), in WM tracts, mean APC ranging from -0.20 to -0.65 depending on the tract location (Storsve et al., 2016). In our sample, the change in cortical morphology from scan 1 to scan 2 (2.5 years apart) was not statistically significant. However, the findings reported by these longitudinal studies suggest that while the neuroanatomical alterations may be too subtle to reach statistically significance in the investigated timespan of 2.5 years, they may still contribute to BOLD SD age-related findings. Particularly, accounting for CT age-related changes may help "unmask" the functional value of BOLD SD.
While the relationship between cortical morphology metrics in aging is dynamic, GM volume changes were largely accounted by changes in CT rather than CA, highlighting the importance of tracking changes in CT (Storsve et al., 2014). In fact, studies show that CT and CA are genetically independent (Panizzon et al., 2009). Their neurodevelopment in the lifespan is largely independent of each other suggesting that they should be considered as separate metrics with different contributions to cortical volume (Im et al., 2008;Winkler et al., 2010;Eyler et al., 2011;Lemaitre et al., 2012;Storsve et al., 2014). Furthermore, studies show that that some cortices experience decrease in cortical morphology metrics at a higher rate than others (Storsve et al., 2014).
The hemodynamic properties of each brain region are highly correlated with its cortical structure (Ulrich and Yablonskiy, 2016;Wen et al., 2018). Studies based on quantitative gradient recalled echo (qGRE) analysis use the GRE signal decay rate parameter (R2t * ) to gather information on tissue-specific contributions to the fMRI signal. Regional R2t * metric variations are associated with variations in neuronal density and myelination (positive correlation), as well as glial cells and synapses concentration (i.e., negative correlation) (Wen et al., 2018). In normal aging R2t * increases in all cortical regions (Zhao et al., 2016). Thicker regions show the opposite pattern, respectively, decreased neuronal density, and higher concentration of glial cells and synapses relative to neurons. Thicker areas also tend to extract less oxygen from the blood, as measured by oxygen extraction ratios (i.e., expressed as local-to-global ratio). These findings indicate that laminar differences in cellular content impact neurovascular coupling mechanisms, which in turn may compromise the power of BOLD SD measurements to detect "real" changes in neuronal variability processing (Harris et al., 2011). Although, laminar differences in BOLD SD remain rather elusive, our study suggests the that CT should be considered in BOLD variability studies.

Age-Related Changes in BOLD SD
In our study, we did not find a significant change in BOLD SD , likely due to the low timespan between scans (2.5 years), and relatively low sample (n = 28). Overall, most studies reporting a decrease in BOLD SD suggest that this finding may indicate structural reductions in synaptic complexity and integrity, as well as functional decline in neural optimization and flexibility (Garrett et al., 2010). Burzynska et al. (2015) reported a positive correlation between increased microstructural integrity of WM and BOLD SD in healthy adults, consequently warranting the consideration of structural alterations in variability studies. Nomi et al. (2017) found increases in variability in bilateral anterior insula, anterior cingulate and ventral temporal cortex, and decreases in BOLD SD in sensory brain regions and other subcortical areas, in participants with ages ranging from 6 to 86 years. In their study they accounted for GM volume but not CT or CA. Two seminal studies investigated BOLD SD differences in aging alone (Garrett et al., 2010) and in relation to performance on cognitive tasks (Garrett et al., 2011) between a young group of participants (20-30 years) and an older one (56-85 years) and found both increases and decreases in regional variability with younger age alone, and younger age and better performance (Garrett et al., 2011) making it rather difficult to isolate specific key contributing regions. Most regions in these two studies showed the same trend but there were some inconsistencies. For example, superior frontal gyrus was reported to show greater variability with younger age and better performance (Garrett et al., 2011) in one study, while another one reported that its variability increases with age (Garrett et al., 2010). CT was not considered in any of these studies. The relationship between cognitive performance and/or flexibility and BOLD SD is definitely complex and task dependent. Behavioral studies indicate age-related differences in intra-individual variability on cognitive performance tests involving reaction time, and working memory. Older adults showing higher intra-individual variability on reaction time tests than younger adults, while the opposite is observed on working memory tests (Fagot et al., 2018). Additionally, the relationship between BOLD SD change, aging and cognitive health has not been a consistent finding. Some studies have found increased BOLD in healthy older adults (Baum and Beauchamp, 2014), Alzheimer's (Scarapicchia et al., 2018), attention deficit hyperactivity disorder (ADHD) (Nomi et al., 2018), and other neurological disorders (Zöller et al., 2017).

Other Possible Confounds in BOLD SD Studies
Most studies investigating BOLD variability utilized a cross-sectional design while we used a within-subject design. This design is a highlight of our study because it allows for the investigation of normative age-related differences, specifically intra-individual effects of the processing of aging on cortical morphology and BOLD SD , rather than process of aging simply age differences across groups. Importantly, the within-subject design helps account for numerous sources of individual differences that may affect BOLD variability such as differences in dopaminergic neurotransmission (Guitart-Masip et al., 2016;Alavash et al., 2018), significant differences in vascular features, pain sensitivity and coping (Rogachov et al., 2016), clinical symptomology (Ke et al., 2015;Nomi et al., 2018;Zöller et al., 2018), and even differences in tendency for financial risk taking (Samanez-Larkin et al., 2010). Given that in our study the participants were scanned at an interval of approximately 2.5 years, we were able to at least partially account for such factors. Furthermore, it is clear from the literature that there are numerous environmental factors [i.e., socioeconomic status, education (Chan et al., 2018)] that may contribute to individual neuroanatomical alterations, which in turn may confound BOLD variability studies. This further supports the findings of our study, specifically that CT is a neuroanatomical metric that should be accounted for.
Lastly, there may be inconsistency in results between studies due to: (1) using task-fMRI and resting-state fMRI, (2) calculation of BOLD signal variability using standard deviation vs. mean-square successive difference) (for review, Garrett et al., 2013b), (3) type of statistical analysis performed (e.g., partial least squares method, LMM, general linear model), and (4) not accounting for contribution of CT and the other mentioned sources of individual differences.

STUDY LIMITATIONS
The main limitations, as discussed above and further addressed by Garrett et al. (2013b); Scarapicchia et al. (2018), are the lack of standardization in acquiring and analyzing the fMRI data. Additionally, the sample used in our study consists of a relatively small sample of older adults and rather short scan-rescan time, consequently this limits our ability to make strong generalization to other age groups or longer time spans, respectively. Nevertheless, since studies show that the brain undergoes extensive annual structural alterations across the lifespan at different rates, our finding that CT contributes to BOLD SD alterations, remains an important consideration.

CONCLUSION
In conclusion, contrary to a view that BOLD variability is just "noise, " we consider it to be emerging as an important metric of normal aging. Keeping in mind that across the lifespan there are considerable cortical morphometry alterations and that cortical depth affects the BOLD signal, we reported that CT contributes to BOL DSD changes, in an older sample of health adults. An increasing number of studies are considering healthy regional BOLD SD changes in the context of functional networks (Rogachov et al., 2016;Nomi et al., 2017). Pursuing this direction while accounting for CT and other possible confounding factors such as dopaminergic neurotransmission, socioeconomic background etc., should reveal new insights into the mechanisms behind age-related neural processes.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript and/or the supplementary files.

AUTHOR CONTRIBUTIONS
DP, RE, and SR conceived and designed this specific project on the data collected by NM and AR. DP analyzed the data. All co-authors participated in writing of the manuscript.