Cortical thickness across the lifespan: Data from 17,075 healthy individuals aged 3–90 years

Abstract Delineating the association of age and cortical thickness in healthy individuals is critical given the association of cortical thickness with cognition and behavior. Previous research has shown that robust estimates of the association between age and brain morphometry require large‐scale studies. In response, we used cross‐sectional data from 17,075 individuals aged 3–90 years from the Enhancing Neuroimaging Genetics through Meta‐Analysis (ENIGMA) Consortium to infer age‐related changes in cortical thickness. We used fractional polynomial (FP) regression to quantify the association between age and cortical thickness, and we computed normalized growth centiles using the parametric Lambda, Mu, and Sigma method. Interindividual variability was estimated using meta‐analysis and one‐way analysis of variance. For most regions, their highest cortical thickness value was observed in childhood. Age and cortical thickness showed a negative association; the slope was steeper up to the third decade of life and more gradual thereafter; notable exceptions to this general pattern were entorhinal, temporopolar, and anterior cingulate cortices. Interindividual variability was largest in temporal and frontal regions across the lifespan. Age and its FP combinations explained up to 59% variance in cortical thickness. These results may form the basis of further investigation on normative deviation in cortical thickness and its significance for behavioral and cognitive outcomes.


Abstract
Delineating the association of age and cortical thickness in healthy individuals is critical given the association of cortical thickness with cognition and behavior. Previous research has shown that robust estimates of the association between age and brain morphometry require large-scale studies. In response, we used cross-sectional data from 17,075 individuals aged 3-90 years from the Enhancing Neuroimaging Genetics through Meta-Analysis (ENIGMA) Consortium to infer age-related changes in cortical thickness. We used fractional polynomial (FP) regression to quantify the association between age and cortical thickness, and we computed normalized growth centiles using the parametric Lambda, Mu, and Sigma method. Interindividual variability was estimated using meta-analysis and one-way analysis of variance. For most regions, their highest cortical thickness value was observed in childhood. Age and cortical thickness showed a negative association; the slope was steeper up to the third decade of life and more gradual thereafter; notable exceptions to this general pattern were entorhinal, temporopolar, and anterior cingulate cortices. Interindividual variability was largest in temporal and frontal regions across the lifespan. Age and its FP combinations explained up to 59% variance in cortical thickness. These results may form the basis of further investigation on normative deviation in cortical thickness and its significance for behavioral and cognitive outcomes.
Structural MRI is the most widely used neuroimaging method in research and clinical settings because of its excellent safety profile, ease of data acquisition and high patient acceptability. Thus, establishing the typical patterns of age-related changes in cortical thickness as reference data could be a significant first step in the translational application of neuroimaging. The value of reference data is firmly established in medicine where deviations from an expected range are used to trigger further investigations or interventions. A classic example is the body mass index (BMI) which has been instrumental in informing about risk for relating to cardio-metabolic outcomes (Aune et al., 2016).
The present study harnessed the power of the Enhancing Neuroimaging Genetics through Meta-Analysis (ENIGMA) Consortium, a multinational collaborative network of researchers organized into working groups, which conducts large-scale analyses integrating data from over 250 institutions (Thompson et al., 2017;Thompson et al., 2020).
Within ENIGMA, the focus of the Lifespan Working group is to delineate age-associations in brain morphometric measures extracted from MRI images using standardized protocols and unified quality control procedures harmonized and validated across all participating sites.
The ENIGMA Lifespan data set is the largest sample of healthy individuals available worldwide that offers the most comprehensive coverage of the human lifespan. This distinguishes the ENIGMA Lifespan data set from other imaging samples, such as the UK Biobank (http:// www.ukbiobank.ac.uk) which includes individuals over 40 years of age. In the present study, we used MRI data from 17,075 healthy participants aged 3-90 years to infer age-associated trajectories of cortical thickness. We also estimated regional interindividual variability in cortical thickness across the lifespan because it represents a major source of inter-study variation (Raz et al., 2010;Wierenga et al., 2020). Based on prior literature, our initial hypotheses were that in most regions the relationship between age and thickness will follow an inverse U-shape and will be influenced by sex.

| Study samples
De-identified demographic and cortical thickness data from 83 worldwide samples ( Figure 1) were pooled to create the data set analyzed in this study. For samples from longitudinal studies, only baseline MRI scans were considered. The pooled sample comprised 17,075 participants (52% female) aged 3-90 years; only participants with complete data were included (Table 1). All participants had been screened to exclude psychiatric disorders, medical and neurological morbidity and cognitive impairment. Information on the screening protocols and eligibility criteria is provided in Table S1.

| Image acquisition and processing
Prior to pooling the data used in this study, researchers at each participating institution (a) used the ENIGMA MRI analysis protocols, which are based on FreeSurfer (http://surfer.nmr.mgh.harvard.edu), to compute the cortical thickness of 68 regions from high-resolution T1-weighted MRI brain scans collected at their site; (b) inspected all images by overlaying the cortical parcellations on the participants' anatomical scans and excluded improperly segmented scans; (c) identified outliers using five median absolute deviations (MAD) of the median value (additional details in the supplement). Information on scanner vendor, magnetic field strength, FreeSurfer version and F I G U R E 1 ENIGMA Lifespan samples. Abbreviations are explained in Table 1 acquisition parameters for each sample as provided by the participating institutions is detailed in Table S1.

| Analysis of age-related changes in cortical thickness
We modeled the effect of age on regional cortical thickness using higher order fractional polynomial (FP) regression analyses (Royston & Altman, 1994;Sauerbrei, Meier-Hirmer, Benner, & Royston, 2006) implemented in STATA software version 14.0 (Stata Corp., College Station, TX). FP regression is one of the most flexible methods to study the effect of continuous variables on a response variable (Royston & Altman, 1994;Sauerbrei et al., 2006). FP allows for testing a broad family of shapes and multiple turning points while simultaneously providing a good fit at the extremes of the covariates (Royston & Altman, 1994). Prior to FP regression analysis, cortical thickness values were harmonized between sites using the ComBat method in R (Fortin et al., 2018). ComBat uses an empirical Bayes method to adjust for inter-scanner variability in the data while preserving biological variability. As the effect of scanner was adjusted using ComBat, we only included sex as a covariate in the regression models. Additionally, standard errors were adjusted for the effect of scanner in the FP regression. We centered the data from each brain region so that the intercept of an FP was zero for all covariates. We used a predefined set of power terms (−2, −1, −0.5, 0.5, 1, 2, 3) and the natural logarithm function, and up to four power combinations to identify the best fitting model. FP for age was written as age (p1, p2, … p6) 0 β where p in age (p1, p2, …p6) refers to regular powers except age (0) which refers to ln(age). Powers can be repeated in FP; each time a power s repeated, it is multiplied by another ln(age). As an example: age 0,1,1 ð Þ 0 β = β 0 + β 1 age 0 + β 2 age 1 + β 3 age 1 ln age ð Þ = β 0 + β 1 ln age ð Þ+ β 2 age + β 3 age ln age ð Þ 494 models were trained for each region. Model comparison was performed using a partial F-test and the lowest degree model with the smallest p-value was selected as the optimal model. Following permutation, critical alpha value was set at .01 to decrease the probability of overfitting. The age at maximum cortical thickness for each cortical region was the maximum fitted value of the corresponding optimal FP model.
Further, we divided the data set into three age-groups corresponding to early (3-29 years), middle (30-59 years) and late life (60-90 years). Within each age-group, we calculated Pearson's correlation coefficient between age and regional cortical thickness. Finally, we used the cocor package in R to obtain P-values for the differences in correlation coefficients between males and females in each age-group.

| Interindividual variation in cortical thickness
The residuals of the FP regression models for each cortical region were normally distributed. Using one-way analysis of variance we extracted the residual variance around the optimal fitted FP regression model so as to identify age-group differences in interindividual variation for each cortical region. Separately for each age-group (t), we calculated the mean age-related variance of each cortical region using P ffiffiffi ffi where e 2 denotes the squared residual variance of that region around the best fitting FP regression line for each individual (i) of that age-group, and n the number of observations in that age-group. Because the square root of the squared residuals was positively skewed, we applied a natural logarithm transformation to the calculated variance. To account for multiple comparisons (68 regions assessed in three age-groups), a Bonferroni adjusted p-value of 0.0007 as chosen as a cut-off for a significant F-test. To confirm that the scanner effect did not drive the interindividual variability analyses, we also conducted a meta-analysis of the SD of the regional cortical thickness in each age-group, following previously validated methodology (Senior, et al., 2016). To test whether interindividual variability is a function of surface area (and possibly measurement error by FreeSurfer) we plotted the SD values of each region against their corresponding average surface area.
LMS is considered a powerful method for estimating centile curves based on the distribution of a response variable at each covariate value (in this case age). GAMLSS uses a penalized maximum likelihood function to estimate parameters of smoothness (effective degrees of freedom) which are then used to estimate the λ, μ, and σ parameters.
The goodness of fit for these parameters in the GAMLSS algorithm is established by minimizing the Generalized Akaike Information Criterion (GAIC) index.

| RESULTS
3.1 | Association of age with cortical thickness Figure 2 shows the shape of the association of age with cortical thickness in each lobe, while the corresponding information on all cortical regions is provided in File S1. For most regions, the highest value for cortical thickness was observed in childhood; age and cortical thickness showed a negative linear correlation, with the slope being steep until the third decade of life (Table S2). By contrast, the entorhinal and temporopolar cortices showed an inverse U-shaped relation with age bilaterally while in the anterior cingulate cortex (ACC) showed an attenuated U-shape. In general, age and its FP combinations explained up to 59% of the variance in mean cortical thickness (Table S2). Age explained the smallest proportion of the variance for entorhinal (1-2%) and temporopolar (2-3%) cortices but the largest proportion of variance for the superior frontal and precuneus gyri (50-52%).
We observed significant sex differences in the slopes of age- Further, sex differences were also noted at the regional level in the early-and middle-life groups. In the early-life group, the slope of the association between age and cortical thickness was steeper in males than in females in the bilateral cuneus, lateral occipital, lingual, superior parietal, postcentral, and paracentral, precuneus, and pericalcarine gyri (all p < .0007). In middle-life age-group, the slope was steeper in males than in females in the bilateral pars orbitalis and pars triangularis as well as left isthmus of the cingulate, pars opercularis, precuneus, rostral middle frontal, and supramarginal, and right fusiform, inferior temporal, inferior parietal, lateral occipital, lateral orbitofrontal, rostral anterior cingulate, superior frontal, supramarginal regions, and the insula (all p < .0002) (Figures 3 and S1, Table S3).

| Interindividual variation in cortical thickness
Across age-groups (early, middle, and late life), interindividual variability in regional cortical thickness, as measured by pooled SD, was between 0.1 and 0.2 mm. Details are provided in Table S4, Figures 4 and S2. High interindividual variation was mainly confined bilaterally in the entorhinal, parahippocampal, transverse temporal, temporopolar, frontopolar, anterior cingulate and isthmus, and pars orbitalis regions. We confirmed the replicability of these findings in each agegroup by conducting meta-analysis following the procedures set-out by Senior et al. (2016).
Finally, we observed a nonlinear association between regional cortical surface area and interindividual variability with variability F I G U R E 2 Illustrative Fractional Polynomial Plots for the association of age and cortical thickness. We present exemplars from each lobe as derived from fractional polynomial analyses of the entire data set. Details regarding the association of age and thickness for all cortical regions (for the entire data set and separately for males and females) are given in the supplementary material being typically higher in regions with smaller surface areas ( Figure S3).

| Centile curves of cortical thickness
Representative centiles curves for each lobe are presented in Figure 5.
Centile values for the thickness of each cortical region, stratified by sex and hemisphere, are provided in Tables S5 to S7 and File S2.

| DISCUSSION
In the present study, we provide the most comprehensive characterization of the association between age and regional cortical thickness across the human lifespan based on multiple analytic methods (i.e., FP analysis, meta-analysis and centile calculations) and the largest data set of cortical thickness measures available from healthy individuals aged 3 to 90 years. In addition to sample size, the study benefited from the standardized and validated protocols for data extraction and quality control that are shared by all ENIGMA sites and have supported all published ENIGMA structural MRI studies (Thompson et al., 2020).
Most regional cortical thickness measures reached their maximum value between 3 and 10 years of age, showed a steep decrease during the second and third decades of life and an attenuated or plateaued slope thereafter. This pattern was independent of the hemisphere and sex. A recent review (Walhovd, Fjell, Giedd, Dale, & Brown, 2017) has highlighted contradictions between studies that report an increase in cortical thickness during early childhood and studies that report a decrease in cortical thickness during the same period. The results from the current study help reconcile previous findings as they show that the median age at maximum thickness for most cortical regions is in the lower bound of the age-range examined here. However, these findings must be considered in the context to the fewer data points available for those below the age of 10 years.
The general pattern of greater cortical thinning with advancing age was similar in both sexes. When participants were divided in early-, middle-and late-life groups, sex differences in the slope between age and cortical thickness was noted primarily for the midlife group. In this age-group, which included individuals aged 30-59 years, the slope was steeper in males than in females. This sexdifference has not been reported in other studies (Fjell et al., 2015;Raz et al., 2005;Raz et al., 2010;Storsve et al., 2014) which generally had smaller samples (<2000), shorter observation periods or examined age-related trajectories of cortical thickness after the effect of sex was regressed-out (e.g., Fjell et al., 2009). Although the sexdifferences reported here may be incidental, they resonate with findings of generally higher cognitive reserve in women as they enter later-life (Mauvais-Jarvis et al., 2020).
In the entorhinal and temporopolar cortex there were minimal age-related changes until the seventh to eighth decades of life; thereafter both regions showed age-related decrease in cortical thickness.
Although the FreeSurfer estimation of cortical thickness in these regions is often considered suboptimal (compared with the rest of the brain), we note that our findings are consistent with a prior multicenter study of 1,660 healthy individuals (Hasan et al., 2016). Further, the current study supports results from the National Institutes of Health MRI study of 384 individuals that found no significant change in the bilateral entorhinal and medial temporopolar cortex between the ages of 4-22 years (Ducharme et al., 2016). A further study of 207 healthy adults aged 23-87 years also showed no significant cortical thinning in the entorhinal cortex until the sixth decade of life (Storsve et al., 2014). These observations suggest that the cortex of the entorhinal and temporopolar regions is largely preserved across the lifespan in healthy individuals. Both these regions contribute to episodic memory while the temporopolar region is also involved in semantic memory (Rolls, 2018). Degenerative changes of the temporopolar cortex have been reliably associated with semantic dementia, which is characterized by loss of conceptual knowledge about realworld items (Hodges & Patterson, 2007 A consistently higher degree of interindividual variation was observed in the most rostral frontal regions (frontopolar cortex and pars orbitalis), in the ACC and in several temporal regions (entorhinal, parahippocampal, temporopolar, and transverse temporal cortex). To some degree, greater variability in several of these regions may reflect measurement challenges associated with their small size ( Figure S3).
Nevertheless, the pattern observed suggests that greater interindividual variability may be a feature of proisocortical and periallocortical regions (in the cingulate and temporal cortices) that are anatomically connected to prefrontal isocortical regions, and particularly the frontopolar cortex. This prefrontal isocortical region is considered evolutionarily important based on its connectivity and function in humans and nonhuman primates (Ongür, Ferry, & Price, 2003;Semendeferi et al., 2011). The frontopolar region has several microstructural characteristics, such as a higher number and greater width of minicolumns and greater interneuron space, which are conducive to facilitating neuronal connectivity (Semendeferi et al., 2011). According to the popular "gateway" hypothesis, the lateral frontopolar cortex implements processing of external information ("stimulus-oriented" processing) while the medial frontopolar cortex attends to self-generated or maintained representations ("stimulusindependent" processing) (Burgess, Dumontheil, & Gilbert, 2007).
Stimulus-oriented processing in the frontopolar cortex is focused on multitasking and goal-directed planning while stimulus-independent processing involves mainly metalizing and social cognition (Gilbert, Gonen-Yaacovi, Benoit, Volle, & Burgess, 2010). The other regions (entorhinal, parahippocampal, cingulate, and temporopolar) with high interindividual variation in cortical thickness are periallocortical and proisocortical regions that are functionally connected to the medial frontopolar cortex (Gilbert et al., 2010;Moayedi, Salomons, Dunlop, Downar, & Davis, 2015). Notably, the periallocortex and proisocortex are considered transitional zones between the phylogenetically older allocortex and the more evolved isocortex. Specifically, the entorhinal F I G U R E 5 Illustrative normative centile curves of cortical thickness. We present exemplar sets of centile curves for each lobe as derived from LMS of the entire data set. Normative centile curves for all cortical regions (for the entire data set and separately for males and females) are given in the supplementary material cortex is perialiocortical (Insausti, Muñoz-López, Insausti, & Artacho-Pérula, 2017), the cingulate and parahippocampal cortices are proisocortical and the cortex of the temporopolar region is mixed (Blaizot et al., 2010;Petrides, Tomaiuolo, Yeterian, & Pandya, 2012).
Considered together, these regions are core nodes of the default mode network (DMN; Raichle et al., 2001). At present, it is unclear whether this higher interindividual variation in the cortical thickness of the DMN nodes is associated with functional variation, but this is an important question for future studies.
The results presented here are based on the largest available brain MRI data set worldwide covering the human lifespan. However, none of the pooled samples in the current study was longitudinal. We fully appreciate that longitudinal studies are considered preferable to cross-sectional designs when aiming to define age-related brain morphometric trajectories. However, a longitudinal study of this size over nine decades of life is not feasible. In addition to problems with participant recruitment and retention, such a lengthy study would have involved changes in scanner types, magnetic field strengths, and acquisition protocols in line with necessary upgrades and technological advances. Nevertheless, it is possible to test the alignment between the results presented here and data from longitudinal cohorts, many of which are also available through the ENIGMA consortium. We consider this an important direction for follow-up studies. We took several steps to mitigate against site effects. First, we ensured that we used age-overlapping data sets throughout. Second, standardized analyses and quality control protocols were used to extract cortical thickness measures at all participating institutions.
Third, we estimated and controlled for the contribution of site and scanner using ComBat prior to conducting our analysis. The validity of the findings reported here is reinforced by their alignment with the results from short-term longitudinal studies of cortical thickness (Shaw et al., 2008;Storsve et al., 2014;Tamnes et al., 2010;Thambisetty et al., 2010;Wierenga et al., 2014). The generalizability of our findings for the older age-group is qualified by our selection of individuals who appear to be aging successfully in terms of cognitive function and absence of significant medical morbidity. Nevertheless, despite the efforts to include only healthy older individuals, the observed pattern of brain aging may still be influenced by subclinical mental or medical conditions. For example, vascular risk factors (e.g., hypertension) are prevalent in older individuals and have been associated with decline in the age-sensitive regions identified here (Raz et al., 2005). Thus, we cannot conclusively exclude the possibility that such factors may have contributed to our results. In addition, a wide range of factors have been associated with cortical morphology throughout the lifespan. Key among them are genetic factors (Grasby, 2020;Teeuw et al., 2019) and indices of socioeconomic status (Chan et al., 2018;Modabbernia et al., 2020;Ziegler et al., 2020) and possibly race (Zahodne et al., 2015). These factors were not modeled here as the relevant information was not collected in a systematic and harmonized fashion across contributing cohorts. It is therefore unclear to what extent they might have influenced the general pattern of age-related associations with cortical thickness reported in the current study; qualifying their possible effects is a priority for future investigations. Cellular studies show that the number of neurons, the extent of dendritic arborization, and amount of glial support explain most of the variability in cortical thickness (la Fougère et al., 2011;Pelvig, Pakkenberg, Stark, & Pakkenberg, 2008;Terry, DeTeresa, & Hansen, 1987). MRI lacks the resolution to assess microstructural tissue properties but provides an estimate of cortical thickness based on the MR signal. Nevertheless, there is remarkable similarity between MRI-derived thickness maps and postmortem data (Fischl & Dale, 2000). Finally, we present the centile curves to stimulate further research in developing normative reference values for neuroimaging phenotypes which should include investigation of measurement errors and reproducibility. In this context, the centile curves should not be used clinically or to make inferences about single individuals.
The findings of the current study suggest several avenues of further research. MRI-derived measures of cortical thickness do not provide information on the mechanisms that underlie the observed agerelated associations. However, the results provided here could be used to study further factors that may lead to deviations in cortical thickness way from the expected age-appropriate range. Additionally, the results of the current study provide a new avenue for investigating the functional correlates, either cognitive or behavioral, of agerelated changes and interindividual variation in regional cortical thickness.
In summary, using existing cross-sectional data from 17,075 individuals we performed a large-scale analysis to investigate the agerelated changes in cortical thickness. The size and age-coverage of the analysis sample has the potential to inform about developmental and aging changes in cortical morphology and provide a foundation the study of factors that may lead to deviations from normative patterns.

ACKNOWLEDGMENTS
This study presents independent research funded by multiple agencies. The funding sources had no role in the study design, data collection, analysis, and interpretation of the data. The views expressed in the manuscript are those of the authors and do not necessarily repre- Healthcare.

DATA AVAILABILITY STATEMENT
The ENIGMA Lifespan Working Group welcomes expression of interest from researchers in the field who wish to use the ENIGMA samples. Data sharing is possible subsequent to consent for the principal investigators of the contributing datasets. Requests should be directed to the corresponding authors.