The midpoint of cortical thinning between late childhood and early adulthood differs between individuals and brain regions: Evidence from longitudinal modelling in a 12-wave neuroimaging sample

Charting human brain maturation between childhood and adulthood is a fundamental prerequisite for understanding the rapid biological and psychological changes during human development. Two barriers have precluded the quantification of maturational trajectories: demands on data and demands on estimation. Using high-temporal resolution neuroimaging data of up to 12-waves in the HUBU cohort (N = 90, aged 7-21 years) we investigate changes in apparent cortical thickness across childhood and adolescence. Fitting a four-parameter logistic nonlinear random effects mixed model, we quantified the characteristic, s-shaped, trajectory of cortical thinning in adolescence. This approach yields biologically meaningful parameters, including the midpoint of cortical thinning (MCT), which corresponds to the age at which the cortex shows most rapid thinning - in our sample occurring, on average, at 14 years of age. These results show that, given suitable data and models, cortical maturation can be quantified with precision for each individual and brain region.


Introduction 23
The human cortex undergoes protracted microscopic and macroscopic structural changes 24 between childhood and adulthood 1 . Individual differences in cortical structure have been 25 associated with a range of phenotypic differences, including physical and mental health 2-5 , 26 neurodevelopmental disorders as well as cognitive performance in childhood and 27 adolescence 2,6 . Although other measures of brain structural development such as brain 28 volume or white matter connectivity provide complementary insights into brain maturation, 29 cortical thinning is one of the most widely used proxies of brain maturation 7-9 . We here use 30 the term apparent cortical thickness throughout to highlight that Magnetic Resonance 31 Imaging (MRI) studies measure a proxy of cortical thickness, that follows a similar spatial 32 patten to that observed in histological studies 10 , but has absolute values that may be 33 influenced by signal intensities and contrasts 11 . 34 Previous work has generally observed developmental decreases in cortical thickness 35 from childhood (after 2-3 years of age) to early adulthood 8 , with longitudinal studies showing 36 that the rate of cortical thinning increases in adolescence 9,12 . There is an emerging consensus 37 on a characteristic s-shaped, non-linear trajectory at the population level ( Figure 1) 9,12 . This 38 process of cortical thinning is thought to reflect a range of underlying biological processes 13 , 39 including increasing myelination of the deeper cortical layers 14 and decreasing synaptic 40 Hypotheses concerning individual differences in the process of cortical thinning and 45 their link to psychological development feature prominently in neurodevelopmental 46 theories. For instance, it has been hypothesized that a life history of adversity may lead to 47 accelerated 16,17 or delayed 18,19 cortical maturation. The developmental mismatch theory 48 suggests that a mismatch in maturity between subcortical and cortical brain regions (with 49 frontal regions commonly thought to thin last) may help explain the prevalence of risk-taking 50 behaviour during adolescence when this maturation disparity is thought to be maximal 20,21 . 51 At the group comparison level, hypotheses posit that girls demonstrate earlier cortical 52 maturation than boys 22-24 , with these differences hypothesized to underlie developmental 53 differences in behavioural and psychopathological phenotypes. Similarly, Nunes et al. (2020) 54 hypothesized that children with autism spectrum disorder are characterized by accelerated 55 brain maturation 25 . Overall, hypothesized differences in cortical maturation are central to 56 some of the most influential neurodevelopmental theories. However, the data and analytic 57 methods currently used to capture maturation are no match for ambitions in understanding 58 and applying the construct of cortical maturation. 59 Longitudinal data is costly and time-consuming to procure, therefore, most 60 neuroimaging studies to-date rely on cross-sectional data. Cross-sectional data precludes 61 the investigation of developmental changes and cross-sectional measures necessarily 62 conflates distinct sources of cortical thickness differences (baseline thickness, the onset of 63 maturation, as well as speed and total amount of thinning). Only under extremely restrictive 64 assumptions (e.g., identical brain thickness in early childhood and late adolescence, identical 65 Most longitudinal studies investigating individual differences to date have used linear 75 modelling approaches, such as linear mixed effects models, using the percentage of, or 76 absolute change, in cortical thickness between two ages as an indicator of maturation of a 77 given brain region, with larger changes commonly equated to more protracted 78 maturation 29,30 . However, these estimates are typically confounded by the initial thickness 79 of a region -thicker regions can thin more than thinner regions. Even when such confounds 80 are controlled for, the absolute change in thickness remains dependent on the precise age 81 range studied. Alternative, nonlinear approaches, such as Generalized Additive Mixed 82 Models (GAMMs) 31 , can capture complex nonlinear relationships and are excellent tools for 83 predictive purposes. Nonlinear mixed models offer a similarly flexible approach for capturing 84 complex nonlinear relationships to GAMMs. A particular advantage of nonlinear mixed 85 models is that they yield readily interpretable parameters, informative, of, e.g., ages of rapid 86 development ( Figure 1), making them an attractive, but currently underused, tool for 87 developmental neuroscientists. 88 To address the challenge of quantifying cortical maturation at the individual level, we 89 here leverage a unique dataset (with up to 12 longitudinal measurements between late 90 childhood and early adulthood) and a quantitative framework currently underutilized in 91 cognitive neuroscience (non-linear random effects modelling) to demonstrate that cortical  To check the quality of the FreeSurfer outputs, we followed the Enigma protocol 156 Finally, the Euler number, which indicates the topological complexity of the 168 reconstructed cortical surface, was extracted from FreeSurfer for each scan and used as a 169 proxy for in-scanner motion and a quantitative measure of image data quality in our 170 statistical analyses, to account for potential systematical bias in data quality across age. The 171 Euler number has been correlated with visual image quality inspection scores as well as with 172 cortical thickness in a regionally heterogeneous pattern across datasets 43 . The Euler number 173 was extracted for each hemisphere and summed to produce one value per scan. 174 175

Statistical analyses 176
We here modelled cortical thinning between childhood and adulthood using nonlinear In a second step, we assessed differences in MCT estimates across different brain 212 regions. We assessed whether MCTs across brain regions could be constrained to equality 213 using Confirmatory Factor Analysis (CFA), as implemented in lavaan 52 for R (version 0.6-9). 214 We used full information maximum likelihood with robust standard errors to account for 215 missingness and nonnormality. We estimated a one-factor CFA model in which factor 216 loadings were freely estimated for each brain region. We then compared this model to one 217 where the factor loadings were still freely estimated, but the intercepts were constrained to 218 equality using the likelihood ratio test. A significant likelihood ratio test here indicates a loss 219 in model fit for a constrained model, providing evidence that brain regions differ in their 220

MCTs. 221
Third, we used Exploratory Factor Analysis to assess the dimensionality of MCTs across 222 regions and to identify maturational factors capturing developmental trends across regions. 223 This was implemented through parallel analysis via the psych 53 package for R (version 2.1.6), 224 using an oblique oblimin transformation. 225 Finally, we tested the popular hypothesis of a posterior-anterior gradient of 226 development 54 by correlating MCTs for each region with their y-coordinates in MNI space. 227 Coordinates were obtained from the brainGraph 55 package in R (version 3.0.0). In an 228 additional exploratory analysis we also correlated MCTs for each region with their z-229 coordinates in MNI space to test for a potential dorsal-ventral gradient. 230 231

232
Cortical thickness showed nonlinear changes between childhood and early adulthood across 233 cortical regions (Figure 3). The characteristic sigmoid, or s-shaped, curve, of cortical thinning, 234 was apparent across most brain regions. We captured this shape using the four-parameter 235 logistic function: 236

237
The function yields four biologically meaningful parameters: The upper asymptote (AUpper, 238 maximal apparent thickness in mm), lower asymptote (ALower, minimal apparent thickness in 239 mm), and Hill, the slope of change, and the Inflection. We were particularly interested in the 240 latter parameter, which corresponds to the MCT, and was used here as an index of cortical 241 maturation (see Figure 1). The MCT can be compared across individuals, sexes, and brain 242 regions. To this end, we first modelled cortical thickness averaged across the cortex, 243 identified the average pattern of maturation, and assessed sex differences. Next, we 244 modelled cortical thickness for different brain regions to investigate patterns of maturation 245 across the cortex. Sex and in-scanner motion, as a proxy for image quality, was controlled for 246 in all analyses. 247 Because the sigmoid is asymptotic, there is no age at which the brain is mature. 248 Instead, the brain develops throughout the age range investigated here (7 -21 years).    Table 2). 280 281

Significant, but noisy, sex differences in trajectories 282
There were significant sex differences: Females started with thicker cortices than males (b = 283 -0.06, p = .001), had thicker corteces in adulthood (b = -0.07, p = .035) and showed earlier MCT (b = 1.92, p < .001, Figure 4). However, the CVs for these estimates were all well above 20% 285 (Table 1), indicating that these parameters were estimated with low precision. The low 286 precision was likely due to the small number of males in the sample (NMales = 37). were all excluded from subsequent analyses (see Supplementary Tables 3 and 4 for all 301 estimates). For these regions, the logistic may not be a good approximation of the functional 302 relationship between age and cortical thickness. For the pericalcarine, for instance, we find 303 that a linear model fits marginally better than the four-parameter logistic (ΔAIC = 0.50), 304 potentially reflecting the early maturation for this region. The enthorinal is known to be more 305 likely to be affected by eye-movement artifacts and may therefore show a poor model fit.  Table 5). 318 Next, we implemented a formal quantitative test to examine whether regions differed 334 in their MCTs. To do so, we compared a confirmatory factor model with intercepts for MCTs 335 across regions constrained to equality (reflecting the hypothesis that regions mature at the 336 same approximate age) to a model where they are estimated freely, reflecting the 337 hypothesis that regions mature at distinct ages). Despite the considerable added complexity, 338 we found that allowing MCTs to differ between regions, substantially improved model fit 339 (ΔΧ 2 (29) = 930.69, p < .001) indicating pronounced differences in maturational timing 340 between regions independent of their overall thickness.
In addition to MCT-differences across regions globally, we also examined whether 342 MCTs are linked between some or all regions: In other words, if a person is mature in one 343 region, are they then also more mature in all other regions, or are there clusters of brain 344 regions that covary in their relative maturity? To examine this question, we first examined 345 the fit of a one-factor confirmatory model, testing the hypothesis that a single factor could 346 capture MCTs across the brain. This model fit poorly (X 2 (377) = 1236.68, p < .001, CFI = 0.787, 347 Finally, we explored whether there is a regional ordering in the timing of maturation in 355 line with hypotheses of a posterior-to-anterior gradient across the cortex 7,56 . We correlated 356 inflection ages of each region with the region's average y-coordinate in MNI space, as 357 contained in mni.y in brainGraph 55 . We found no evidence for a significant correlation 358 between the spatial location and the MCT (t(27) = 0.12, p = .903, r = 0.02, Bayes Factor = 0.35). 359 In an additional exploratory analysis, we also analyzed whether there is a regional ordering 360 in the timing of maturation in line with hypotheses of a dorsal-ventral gradient across the 361 cortex 57 . We correlated inflection ages of each region with the region's average z-coordinate 362 in MNI space, as contained in mni.z in brainGraph 55 . We found no evidence for a significant 363 correlation between the spatial location and the inflection ages (t(19) = -0.65, p = .522, r = -364 0.15, Bayes Factor = 0.45). 365 Together, these analyses offer new insight into cortical maturation. We demonstrate, 366 for the first time, that it is possible to estimate non-linear maturation independent of overall 367 cortical thickness. Maturational trajectories differed between individuals and cortical 368 regions. The ability to estimate these differences offers a new window into elucidating long-

372
Using longitudinal data of up to 12-waves imaged between late childhood and early 373 adulthood and flexible nonlinear mixed models, we here show that cortical maturity as 374 indexed by the MCT can be separated from other cortical thickness parameters (i.e., 375 asymptotes and slope -hill) and estimated precisely and reliably for each individual and most 376 brain regions. We identified a characteristic, s-shaped, trajectory of cortical thinning: 377 Cortical thickness was high in childhood, followed by decreases in early adolescence, Our finding of an early MCT in several frontal areas is surprising, given previous, mostly 393 cross-sectional studies. There are several possible explanations for this result, which will 394 need to be tested in future research: First, while most frontal regions are usually imaged well, 395 the frontal pole is known to be difficult to image. We therefore advise to interpret the finding 396 for this region with caution. Second, we cannot rule out that data quality issues or model 397 misfit may have affected findings, although quality control procedures, statistical control for 398 in-scanner motion, diagnostic plots and precision estimates do not suggest that this was 399 likely to be the sole explanation. Third, while the scans were evenly distributed across most 400 ages, the interval between waves 11 and 12 was longer, and data was sparser at older ages.
to affect estimation of late maturing, rather than early maturing regions. Fourth, the late 403 maturation suggested by the, mostly, cross-sectional analyses to date, do not replicate in 404 appropriate longitudinal analyses like ours, because of the inherent limitation of cross-405 sectional data for the examination of developmental patterns. Finally, it may be that cross-406 sectional analyses have captured aspects of maturation that are independent of the MCT, 407 the time point of fastest thinning. It may be, for instance, that frontal regions asymptote 408 later than other region, even though the peak of thinning happens quite early. Altogether 409 this finding, along with recent work identifying structural brain networks and hubs [60][61][62][63] , 410 supports a complex systems account, in which the brain matures in a distributed pattern. 411 Future studies could use complex systems approaches like network models to identify how 412 maturation across regions produces changes in cognition and mental states across 413

development. 414
We found evidence for pronounced individual differences, with MCTs differing by 415 several years between individuals. This supports the notion that brain maturation is highly 416 variable and invites questions of potential predictors and outcomes. We here investigated 417 sex as a potential predictor of individual differences. There was some evidence for sex 418 differences, with females showing earlier maturation than males, by about 2 years. This 419 supports initial evidence for earlier maturation in girls from older longitudinal studies 46,47 and 420 may be linked to earlier puberty and different socialization in girls 64,65 . Our estimates of sex 421 differences were relatively noisy, however. Future studies could use nonlinear mixed models 422 in larger cohorts, such as ABCD 28 , to investigate the robustness of sex differences in 423 maturation. Future studies could also investigate other candidate predictors such as 424 environmental influences (e.g., adversity and education) to identify whether these 425 accelerate or decelerate maturation -a yet unsolved conundrum in developmental 426 science 18,19,66 . Investigations of potential outcomes of maturational differences, e.g., 427 cognitive performance or psychiatric diagnosis, would be similarly fascinating, and could 428 eventually include distal effects, such as cognitive and brain changes during ageing 1,67 . 429 It is worth noting that our analytical approach, using the four-parameter logistic 430 function, depends on the nature of cortical thinning: Other measures of brain maturation 431 (e.g., brain v0lume) likely show different developmental shapes and different maturational decreases, reported for volume changes in some cortical regions, would not allow for a 434 quantitative 'midpoint.' However, nonlinear mixed models are extremely versatile and can 435 be fit, in principle, to any functional relationship. This also includes more complex 436 relationships than that captured by the four-parameter logistic fit here, allowing, e.g., for 437 asymmetries in development before and after adolescence. In the future, these nonlinear 438 mixed models can be used to model white matter trajectories and other morphometric 439 measures of cortical and subcortical development to better understand similarities and 440 differences across brain tissues. This will yield a more precise understanding of how changes 441 in grey and white matter work in concert to produce functional changes in the brain and 442

behaviour. 443
Readers may be interested in understanding the power prerequisites for using 444 nonlinear mixed models. Power in nonlinear mixed models depends on the interaction 445 between several factors, including the number of participants, number of timepoints, 446 spacing between time points, missingness and model complexity 68,69 . While past studies of 447 optimal design in pharmacology indicate that at least three time points may be sufficient to 448 estimate simple nonlinear mixed models 70 , future studies will need to evaluate power for 449 plausible developmental functions and typical designs of neuroimaging studies in detail. 450 In conclusion, this study shows that apparent cortical thinning in adolescence is s-