Neurodevelopmental Trajectories in Children With Internalizing, Externalizing and Emotion Dysregulation Symptoms

Introduction Childhood and adolescence are crucial periods for brain and behavioral development. However, it is not yet clear how and when deviations from typical brain development are related to broad domains of psychopathology. Methods Using three waves of neuroimaging data within the population-based Generation R Study sample, spanning a total age range of 6–16 years, we applied normative modeling to establish typical development curves for (sub-)cortical volume in 37 brain regions, and cortical thickness in 32 brain regions. Z-scores representing deviations from typical development were extracted and related to internalizing, externalizing and dysregulation profile (DP) symptoms. Results Normative modeling showed regional differences in developmental trajectories. Psychopathology symptoms were related to negative deviations from typical development for cortical volume in widespread regions of the cortex and subcortex, and to positive deviations from typical development for cortical thickness in the orbitofrontal, frontal pole, pericalcarine and posterior cingulate regions of the cortex. Discussion Taken together, this study charts developmental curves across the cerebrum for (sub-)cortical volume and cortical thickness. Our findings show that psychopathology symptoms, are associated with widespread differences in brain development, in which those with DP symptoms are most heavily affected.


INTRODUCTION
Over the course of childhood and adolescence, both the brain and behavior undergo tremendous development. Regarding the relationship between the developing brain and atypical behavior, a body of evidence has associated differences in brain morphology to multiple domains of psychopathology (1)(2)(3)(4)(5). These studies have assessed multiple measures of brain morphology, including cortical volume and cortical thickness. However, the brain regions that have been identified are widespread and vary substantially across studies (5,6). Additionally, the direction of effect also differs across studies, meaning that some studies find positive relationships between cortical thickness/volume and psychopathology, whereas others find negative associations (1,(7)(8)(9)(10).
The age at which children are assessed may potentially be a crucial factor to unravel why effects across studies differ in both location and direction. Non-linear patterns in brain development across age may partially underlie differences in the direction of observed effects. Total brain volume, for example, increases until adolescence, where it reaches a plateau and starts to decline (11), whereas gray matter volume reaches this peak in early childhood (12). Additionally, evidence suggests that distinct brain lobes and regions within lobes, develop at their own pace (11,13). The asynchronous development of regions and lobes may be an explanation for the effect differences observed in the brain regions involved in these studies. Recent work has therefore used data-driven normative modeling, a technique that can be used to derive typical development curves for brain morphology (14)(15)(16). Emerging evidence suggests that deviations from typical brain development, estimated using these normative models, improves prediction of psychopathology over predictions based on raw brain morphology measures (17).
The direction of the effect in earlier work on internalizing symptoms, seems to be dependent on the age range that is used in studies. Namely, studies including younger age ranges generally observed positive associations (1,7,19), whereas in older age ranges negative associations are observed (8,9). In contrast, for externalizing symptoms, the majority of studies, including a meta-analysis for cortical and subcortical gray matter volume, point toward lower volume and thickness in children with externalizing disorders (5,10,22,(24)(25)(26)(27)(28)(29), while a few report higher cortical volume or thickness (10,31) and one reports a non-linear relationship (23).
The aims of this study were (i) to establish normative development curves for cortical thickness and (sub-)cortical volumes, covering the gray matter of the cerebrum, and (ii) to study to what extent deviations from typical development are related to psychopathology symptoms in a large populationbased cohort of children and adolescents. We hypothesized that all three domains of psychopathology would be related to deviations from normative development of brain morphology. Specifically, we hypothesized that for internalizing symptoms, alterations would be most prominent in the rostral middle frontal cortex, OFC, amygdala and hippocampus; for externalizing in the ACC, OFC, amygdala, hippocampus and striatum; and for DP symptoms in the OFC, amygdala and striatum. Further, we hypothesized that the direction of these deviations varies with age for internalizing symptoms, with positive deviations at younger and negative deviations at older ages. For the externalizing domain we hypothesized that, in line with most prior work, higher symptoms are related to negative deviations from normative development at all ages.

Participants
This study is embedded in the neuroimaging component of the Generation R Study, a large, longitudinal, population-based cohort with an observational design. The recruitment strategy has been described elsewhere (32)(33)(34). In brief, women living within specific zip codes of Rotterdam with a delivery date between April 2002 and January 2006 were invited to participate. The families are still being followed. When children were 6-10 years old (T1), 8-12 years old (T2), and 13 to 16 years old (T3), neuroimaging and behavioral data were collected. Children were included in the current study if they had good quality neuroimaging and behavioral data available in at least one wave of data collection. Neuroimaging data were excluded if any of the following conditions were present: dental braces, incidental findings that significantly alter brain morphology, or poor image quality. At T1, a total of 842 children were included, at T2, 2,708 children were included and at T3, 1,904 were included, resulting in a total sample of 5,454 scans from 4,415 children. A flowchart of the study sample is provided in Figure 1. The Generation R Study was approved by the medical ethics committee at the Erasmus MC and conducted according to the Declaration of Helsinki. Written informed consent, and when applicable assent, was obtained from the caregivers and their children.

Behavioral Assessment
Child behavior was assessed using the Child Behavior Checklist (CBCL). At T1, the CBCL version for children aged 1.5-5 years was used (35), and at T2 and T3, the CBCL version for children aged 6-18 years was used (36). Both versions are reliable and valid questionnaires to assess child behavior (35,36). The CBCL v1.5-5 has 99 items, and v6-18 has 112 items that are scored on a threepoint Likert scale (0 = not true, 1 = somewhat true, 2 = very true). From the CBCL v1.5-5, seven empirically derived syndrome scales were calculated. From the CBCL v6-18, eight syndrome scales were obtained. These syndrome scales were summed into three broad domains of psychopathology [internalizing, externalizing and dysregulation profile (DP) symptoms]. In the CBCL v1.5-5, the internalizing scale includes the emotionally reactive, anxious/depressed, withdrawn and somatic complaints syndrome scales, and in the CBCL v6-18 it is a sum-score of the anxious/depressed, withdrawn/depressed and somatic complaints syndrome scales. Externalizing symptoms were assessed with the attention problems and aggressive behavior syndrome scales in CBCL v1.5-5, and with the rule-breaking behavior and aggressive behavior syndrome scales in v6-18. Lastly, the DP is a comorbid profile, which is the summed score of the anxious/depressed, attention problems and aggressive behavior syndrome scales in both CBCL versions (35,36).

MRI Processing
Image processing for data from T1 to T3 was performed using FreeSurfer analysis suite v6.0.0 (http://surfer.nmr.mgh.harvard.edu/). All images were processed individually using FreeSurfer. FreeSurfer processing steps have been described in detail previously (37). Briefly, the analysis stream includes converting raw DICOM data to "MGZ-files, " skull stripping, intensity normalization, and voxel segmentation of gray matter, white matter and cerebrospinal fluid. Labeling of the gray matter regions was performed using the Desikan-Killiany atlas (38).

MRI Quality Assurance
MRI quality assurance has been described previously (33,34,39). To summarize, FreeSurfer image reconstructions were visually inspected by at least one rater. Based on how well FreeSurfer delineated the gray-white matter and the outer gray matter boundaries, each scan was rated on a Likert scale. Raters included master students, PhD students and postdoctoral researchers, who were all trained extensively, which was completed after correctly rating 30 scans of which quality was determined previously. At T1 and T2, scans were rated on a five point Likert scale (unusable, poor, sufficient, good, excellent). At T3, scans were rated on a three point Likert scale (poor, questionable, good). All scans that were unusable or of poor quality were excluded from the analyses. Quality assessment based on visual inspection was also compared to an automated quality assessment, which has been described previously for T1 and T2 data (40). Visual ratings were also  compared to this automated quality assessment at T3, as well as to the Euler number which can be extracted after FreeSurfer reconstruction (41), this comparison is depicted in Supplementary Figure S1.

Covariates
Multiple covariates were included in the analyses. Sex was derived from medical records at birth. Handedness was measured at each data collection wave, with the Edinburgh Handedness  Inventory (EHI) (42), from which a laterality quotient was obtained ranging from -1 (fully left-handed) to +1 (fully righthanded). Maternal education, household income and child national origin were assessed using a questionnaire. Maternal education and household income were assessed at T1 and used as proxies for socioeconomic status (SES). Maternal education was divided into three categories: low (no education/primary school), middle (high school/vocational training), and high (higher vocational training/university) and household income into two categories: below €2000,-per month and above €2000,per month. Child national origin was assessed at baseline, based on the birth country of the parents, it was categorized as Dutch and non-Dutch (African, American western, American non-western, Asian western, Asian non-western, Cape Verdean, Dutch Antilles, European, Indonesian, Moroccan, Oceania, Surinamese and Turkish).

Statistical Analyses
Our primary analyses assessed the relationship between deviations from typical development in cortical and subcortical regions of interest (ROIs) and multiple domains of psychopathology, using normative modeling. Specifically we included the following ROIs: the ACC (sum of rostral and caudal ACC), OFC (sum of lateral and medial OFC), rostral middle frontal cortex, the amygdala, hippocampus and the striatum (sum of putamen, caudate and nucleus accumbens). For cortical ROIs, we included measures of cortical thickness as well as cortical gray matter volume, for subcortical ROIs, gray matter volumes were included. In our secondary analyses, we explored the remaining (sub-)cortical regions labeled within FreeSurfer (38,43), following the same procedure as for our primary analyses. To reduce the total number of tests, brain measures were averaged across both hemispheres.
The analyses consisted of five steps. First, we residualized brain morphology measures for possible covariate effects using two different models. In model 1 the effects of sex and handedness were regressed out of brain and CBCL measures, in model 2 the effects of SES and child national origin were additionally regressed out. Second, these residualized brain morphology measures were used to fit our normative model. A common way to fit a normative model is to use Gaussian process regression with age, and have the model predict the brain measure from those inputs. Generally, the subjects used in these analyses are considered to have typical development and the model is then validated using a held out subset of typically developing subjects (14)(15)(16). The Generation R study, however, is a population-based sample that is not enriched for children with psychopathology. Thus, we fit the model on all participants, which has been described as a viable option previously (15). Importantly, the current sample did not merely include cross-sectional data, but also longitudinal data for many of the participants. We leveraged the longitudinal data by bootstrapping multiple unique combinations of subsets of scans as training sets. In the individual training sets, all participants were only included once, to prevent overfitting on a single person's development. To reach an approximately even distribution of participants across ages in each training set, approximately 50% of the scans acquired at T2 and T3 were not included in each individual training set. In each test set, 10% of the participants from T1, T2, and T3 were included.
We used Gaussian process regression (GPR) to fit (non)-linear normative trajectories to each brain measure across age. GPRs fit a Gaussian process to the given data points, such that any age, along a continuum (x-axis), is associated with a normal distribution for each brain region (y-axis). This approach is especially well-suited for normally distributed data. Given that we used a population sample, we assume that points for each brain region are normally distributed. Each brain measure was thus fit using a separate Gaussian process with GPytorch's (44) exact Gaussian processes module. The Gaussian process is continuous and can interpolate and extrapolate from the age range of a given set of points. An added benefit of associating a normal distribution with each brain measure given a certain age is that we can use the standard deviation of the distribution to calculate how confident the Gaussian process is in its prediction. It also allows us to interpret the distributions over time as each brain measure's normative trajectory. An important hyperparameter for Gaussian process regression is its kernel, which determines the shape of the line that the normal distributions are centered on. We empirically evaluate a variety of typical kernels for each brain region to limit assumptions about their normative trajectory. Before fitting the GPR, all brain measures were standardized to a mean of 0 and a standard deviation of 1, to accommodate direct comparisons of effect-size estimates across brain regions and measures. The age was rescaled between 0 and 1 based on the minimum and maximum age in the dataset. To assure we obtain the best fit for the trajectory of a given region, we evaluated multiple types of kernels on an unseen validation set. This validation set was a small (10% of each wave in the training fold) subsample of the training set. The kernels we used included a linear, Matern, radial basis function (RBF), and a rational quadratic kernel. We averaged the performance of each kernel over the validation sets to select the best kernel. The best kernel and complete training set, including the validation set, were used to train the final model. The final models for each training set were then used to predict the mean and standard deviance at each age in the test fold.
Third, the difference between these predicted mean and standard deviations for each morphological value and the true morphological values were used to calculate the z-scores for each participant. The formula for this calculation is shown in Where, y are the true values for the brain measures at each age,ŷ are the mean predicted brain measures at each age, and σ are the predicted variances at each age. Note that because the normative model predicts a normal distribution at each age, average predicted brain measure at each age is the most likely value of that brain measure. Fourth, the association between the deviations of each individual from normative development and psychopathology was tested, using separate linear mixed model analyses. Internalizing, externalizing and DP symptoms were entered as dependent variables, z-scores for all brain measures were entered as independent variables, and a random effect was applied for participant ID. Finally, these analyses were repeated with an interaction term for age, to assess whether differences in the slope of deviations from typical development were age-dependent.
To assess the robustness of the findings, two sensitivity analyses were performed. First, normative development curves for (sub-)cortical volume in each region were fit with the effect of total intracranial volume (ICV) regressed out of individual volumes. Z-scores obtained from this model were subsequently related to psychopathology symptoms, to assess whether the effects observed were global or specific to regions. Second, we assessed whether deviations from typical development are specific to psychopathology domains. Therefore, analyses were repeated for brain regions that showed a significant relationship with two or more individual psychopathology domains, in which all significant psychopathology domains were entered in the model simultaneously.
Lastly, as post-hoc analyses, both normative developmental trajectories and deviations due to psychopathology were established for surface area in all hypothesisdriven and exploratory regions of interest, after which deviations from typical development were related to all psychopathology domains.
Bootstrapping and regression analyses were performed in R version 3.6.3 (45), normative modeling was performed in Python version 3.9.0 (46). Missing data in the covariates were imputed 30 times with 30 iterations using multiple imputation through chained equations with the mice package (47). The false discovery rate (FDR) was controlled using the Benjamini Hochberg procedure (48). Primary analyses were corrected for a total of 27 tests, at q-value = 0.05. Exploratory analyses were separately corrected for a total of 180 tests (q-value = 0.05). Analyses using an interaction term were corrected for multiple testing following the same procedure for a total of 207 tests. Hypotheses and analyses for this project were publicly preregistered, a time-stamped version of this preregistration is available via: www.osf.io/aqc4s. Slight deviations from our initial preregistration are described in the Supplementary Material. Analysis scripts are publicly available via https://github.com/ eloygeenjaar/normative-smri-psychopathology.

Sample Characteristics
Sample characteristics are described in Table 1

Normative Development of Brain Morphology
The linear and/or non-linear development curves were fit for each residualized region (model 1 and model 2) included in the Desikan-Killiany atlas (38). Examples of the most common patterns that we observed are depicted in

Cortical and Subcortical Volume
Normative development curves between age 6-16 revealed an increasing slope for (sub-)cortical volume in the entorhinal, inferior temporal, middle temporal, temporal pole, hippocampus, pallidum and thalamic regions; a decreasing slope for the cuneus, frontal pole, isthmus cingulate, lateral occipital, lingual, paracentral, pars opercularis, pars triangularis, pericalcarine, precuneus, post central, supramarginal and transverse temporal regions; and an inverted U-shaped curve for the anterior cingulate, banks of superior temporal sulcus, caudal middle frontal, inferior parietal, orbitofrontal, pars orbitalis, precentral, posterior cingulate, rostral middle frontal, superior frontal, superior parietal, superior temporal, amygdala and striatal regions. Flat trajectories were observed in the fusiform, insula and parahippocampal regions. These patterns were consistent across model 1 and 2. Given that normative development curves were fit on 12 bootstrapped folds of the dataset, the optimal fit differed slightly between individual folds, however, patterns described were consistent across all folds.

Cortical Thickness
The normative models fit on cortical thickness data showed a decreasing slope from early to later neurodevelopment in the majority of regions (see Supplementary Figure S1). The steepest slope was primarily seen between 6 and 12 years of age. Noteworthy exceptions with a fairly flat slope across neurodevelopment were the entorhinal and temporal pole regions. These patterns were consistent across models and folds.

Cortical and Subcortical Volume
All a priori selected, hypothesis driven, regions of interest for cortical and subcortical volume showed a negative relationship with some psychopathology domains, meaning that psychopathology symptoms were related to negative deviations from typical brain development. After correction for multiple testing, negative deviations from typical development in the ACC were related to all psychopathology domains; negative deviations in the OFC, the rostral middle frontal cortex, and the amygdala were related to externalizing and DP symptoms. Lastly, negative deviations from typical development in hippocampal and striatal volume were related to DP symptoms. Full results are shown in Figure 5 (model 2), Table 2 and Supplementary Table S1 (model 1).

Cortical Thickness
Positive associations were observed between deviations from typical cortical thickness development in the OFC, and externalizing and DP symptoms after correction for multiple testing. No associations were observed between deviations from normative development in the ACC and rostral middle frontal cortex, and each of the psychopathology domains. Full results are shown in Figure 6 (model 2), Table 3 and Supplementary Table S2 (model 1).

Cortical and Subcortical Volume
Exploratory analyses revealed significant negative associations between deviations from typical development in several (sub-)cortical volume regions and psychopathology domains. After correction for multiple testing, all psychopathology domains were related to negative deviations from typical development in the precuneus. Negative deviations from typical development in the cuneus, fusiform, inferior parietal, inferior temporal, isthmus cingulate, lateral occipital, lingual, middle temporal, posterior cingulate, precentral, superior parietal and thalamus were observed for externalizing and DP symptoms. Lastly, negative deviations in the parahippocampal region were specific to externalizing symptoms and negative deviations in the insula, pars triangularis, pericalcarine, postcentral and superior temporal region were specific to DP symptoms. Full results are shown in Figure 5 (model 2), Table 4 and Supplementary Table S3 (model 1).

Cortical Thickness
Positive deviations were observed in the pericalcarine region in relation to internalizing symptoms. Additionally, positive deviations from typical development in the frontal pole and the posterior cingulate region were related to higher externalizing and DP symptoms. Deviations from typical development in cortical thickness for most regions were not related to psychopathology symptoms. Full results are shown in Figure 6 (model 2), Table 5 and Supplementary Table S4 (model 1).

Interaction Effect Age
After correction for multiple testing, a significant positive interaction effect between age and the deviations from typical development in cortical volume was observed in the fusiform and parahippocampal region in relation to externalizing symptoms. This indicates that with increasing age, deviations from typical development become larger in those with externalizing symptoms. Full results are shown in Supplementary Tables S5-S8.

Sensitivity Analyses
In the first sensitivity analysis we repeated our analyses on cortical and subcortical volume while correcting for ICV, to assess whether the observed relationships were global or region specific. While the relationships attenuated in many regions, some region specific deviations were identified. In our hypothesis driven regions, the relationship between the ACC development and all psychopathology domains, as well as the relationship between the amygdala and externalizing symptoms remained statistically significant after adjustment for ICV. Additionally, exploratory analyses indicated a significant relationship between the lingual region and all psychopathology domains; between the cuneus, inferior parietal, lateral occipital and precuneus, and externalizing and DP symptoms; and between the isthmus cingulate, pars triangularis, pericalcarine and posterior cingulate, and DP symptoms. Full results are shown in Supplementary Tables S9, S10.
In the second sensitivity analysis we assessed whether, for those regions that showed a significant relationship with multiple

Post-hoc Analyses
As post-hoc analyses, normative trajectories were established for cortical surface area. Normative trajectories showed an inverted U-shaped relationship for the majority of regions. Notable exceptions were positive trajectories in the anterior cingulate, entorhinal, fusiform, inferior temporal, insula, middle temporal, orbitofrontal, parahippocampal, pars opercularis, pars orbitalis, pars triangularis and temporal pole; and a negative trajectory for the transverse temporal region. Examples of the most common patterns are shown in Figure 7 and full results from the normative model are shown in Supplementary Figure S4. Relationships between deviations from typical development of surface area and psychopathology largely mirrored the associations between cortical volume and psychopathology, in which most regions showed negative deviations from typical development. Full results are presented in Figure 8 and Tables 6, 7.

DISCUSSION
This study aimed to assess the relationship between deviations from typical brain development and psychopathology symptoms using three waves of neuroimaging and behavioral data from a large population-based cohort. We applied normative modeling to derive typical development curves, which showed regional differences in the development of subcortical and cortical volume, as well as cortical thickness and surface area. Psychopathology symptoms were related to deviations from typical development in subcortical volumes, and widespread regions across the cortex for cortical volume and surface area, and some regions for cortical thickness.
Our hypothesis driven and exploratory analyses together revealed that deviations from normative development related to psychopathology symptoms are not restricted to the a priori defined regions of interest, but rather that these deviations were present in regions across the entire cortex, as well as subcortical gray matter volumes. This raises the idea that the observed associations might not be region specific, but rather represent a global effect on brain development. Indeed, the majority of these findings did not remain after additional correction for ICV, indicating that the effects obtained are mostly global. Some areas, however, show a significant relationship on top of this global effect. Taken together, the observed pattern suggests that, while some regions may be particularly important for the emergence of psychopathology symptoms or affected by downstream effects of psychopathology (39), associations between cortical volume and psychopathology are not necessarily restricted to these regions. Given that emotion and behavior require integration of information that involves many brain regions, these small, but widespread differences may together lead to differences in psychopathology symptoms. Thus, our findings bolster the importance of analyzing the entire cortex and subcortex when assessing the relationship between brain morphology and psychopathology in youth, as opposed to restricting analyses to a priori defined regions of interest.
In line with our hypotheses, children with externalizing symptoms deviated from typical development for (sub-)cortical volumes in the ACC, OFC and amygdala, and children with DP symptoms deviated from typical development in the OFC, amygdala and striatum. As opposed to our initial hypothesis, we did not find evidence for a relationship between internalizing symptoms and the hypothesized regions of interest. Additionally, the relationship between deviations from typical development in the hippocampus and the striatum, and externalizing symptoms did not reach statistical significance. In line with the widespread alterations in development of subcortical and cortical volume, we additionally observed associations between the development of the ACC, and internalizing and DP symptoms, the rostral middle frontal cortex and externalizing and DP symptoms, and the hippocampus and DP symptoms. Earlier research on brain morphology and DP symptoms is relatively scarce, which is likely to explain why fewer regions had been reported to be related to DP symptoms previously. Our findings showed that associations between brain morphology and the DP were even more widespread than those for internalizing and externalizing symptoms. After correction for ICV, our results indicated region specific deviations from typical development for cortical volume in the ACC in relation to all psychopathology domains. Indeed, the ventral part of the ACC has been shown to have an important role in emotion regulation, including contextual fear generalization and the top-down regulation of aggressive impulses (49)(50)(51)(52), and the dorsal ACC is crucial for cognitive control (53). Altered development of amygdala volume was only associated to externalizing and DP symptoms in this study, which is surprising given the role of the amygdala in processing of emotional cues that are important for all psychopathology domains (54). We had predicted that deviations in amygdala volume would also be associated with internalizing symptoms, although an earlier meta-analysis indicated that results on amygdala volume are somewhat inconsistent across studies, resulting in an absence of an effect in this meta-analysis (6). Regarding cortical thickness, deviations from typical development in the OFC were related to externalizing and DP symptoms, which is in line with the results obtained for cortical volume. However, in contrast with the findings for cortical volume, no associations between deviations from typical development in the ACC and any type of psychopathology were observed. Lastly, we hypothesized deviations from typical development in cortical thickness in the rostral middle frontal cortex to be related to internalizing symptoms, however, we did not find evidence for the presence of this relationship. Deviations from typical development in surface area in children with symptoms of psychopathology showed remarkably similar results as those obtained for cortical volume. The overlap is likely partially explained by the high correlations between cortical volume and surface area, which in our sample ranged between 0.26 and 0.94. Given the low correlation between cortical thickness and surface area, these findings point toward surface area as an important brain morphology measure to study in relation to psychopathology. Although surface area is studied less extensively than cortical volume and cortical thickness in relation to psychopathology, recent work also showed similar results between psychopathology, and both cortical volume and surface area, but no relationship with cortical thickness (55). Contrary to these findings and the current findings, other work has observed alterations in cortical thickness in relation to multiple domains of psychopathology, whereas alterations in surface area were specific to externalizing disorders (56). A critical difference is, however, that the latter study evaluated the association between childhood psychopathology and brain morphology in mid-adulthood, and thus may not generalize to developmental populations. It will be important for future work to extend the current findings by assessing the relationship between brain morphology and psychopathology at multiple ages across the lifespan. In both our hypothesis driven and exploratory analyses we showed that deviations from normative development were associated with psychopathology. Although the overlap in confidence intervals in the majority of regions do not suggest significant differences in effect sizes between psychopathology domains, a consistent pattern is observed with the largest effect sizes for DP symptoms and the smallest effect sizes for internalizing symptoms. This pattern was also observed in earlier work using the first and third wave of the current sample, for the relationship between cognitive performance, and internalizing, externalizing and DP symptoms (57,58). These findings align closely with recent evidence that some alterations in brain structure and function are shared across many psychiatric disorders (59,60). It is likely that the overlap in involved brain regions can partially be attributed to the high correlation among psychopathology domains. For example, internalizing and externalizing symptoms generally correlate with a coefficient of around 0.5 (61). Achenbach et al. (61) recommended adjustment for externalizing symptoms when internalizing symptoms are assessed and vice versa. Following this recommendation we performed sensitivity analyses for those regions in which deviations from typical development were related to multiple psychopathology domains. Our findings indicated that only DP symptoms were related to regional deviations from typical development above and beyond internalizing and externalizing symptoms. Thus, our findings add to the current knowledge that those with DP symptoms are most heavily affected in terms of symptomatology (62) and cognitive performance (57,58), that they are also most heavily affected in terms of deviations from brain development. A promising line of research that has emerged in recent years, aims to unpack the shared variance between individual dimensions of psychopathology into a general psychopathology factor (63,64). The variance in symptomatology that remains after extraction of the general psychopathology factor can then be viewed as more specific internalizing and externalizing symptoms. Indeed, recent work has used a similar approach to normative modeling as used in the current study, and showed that general psychopathology was associated with deviations from typical development for gray matter volume in widespread regions across the cortex, and additionally identified regions that associated with specific psychopathology dimensions (17). This study focused on gray matter volume and used clinical cases to assess associations with psychopathology. Thus, a promising extension of our study would be to assess whether deviations from typical development in relation to the general psychopathology factor are also more pronounced in cortical volume than cortical thickness.
Regarding the direction of effect, negative deviations from typical development were observed for cortical volume and surface area, and psychopathology, whereas positive deviations were observed in the relationship between cortical thickness and psychopathology. These directions were consistent across all psychopathology domains, in which very few interaction effects for age were observed. The absence of this interaction effect suggests that the direction of effect between psychopathology symptoms and deviations from typical development were largely stable between 6 and 16 years of age. Thus, while we had hypothesized an age-dependent effect on the deviations of typical development related to internalizing symptoms, we only provide evidence for age related effects in the development of the fusiform and parahippocampal volume in relation to externalizing symptoms. To provide more context to the direction of effect, both in cortical volume and thickness, it would be beneficial to not only establish typical development curves using cross-sectional data, but also using longitudinal MRI data. Indeed, some evidence suggests that longitudinal trajectories of cortical development are different in those with internalizing (1,8) and externalizing symptoms (24). Extending findings of the current study, by studying temporal changes in the deviations from typical development in relation to temporal changes in psychopathology, could provide unique insights in the bidirectional relationship between behavioral and brain development.
Findings from our normative model are in line with contemporary work that finds the average growth trajectory for gray matter volume to peak at 6-10 years of age, after which it declines (12,(65)(66)(67). This pattern is also observed for total brain volume, although earlier work did not model a decline in total brain volume after this peak (68) . The parietal and occipital lobe mirror this inverted U-shaped pattern in our normative model. Regional differences in developmental trajectories have been reported previously (11,13,69) and indeed we extend prior knowledge by showing that, similar to the earlier work (69), across a span of 6-16 years of age, many regions in the frontal and temporal lobe reach their peak volume after 6 years of age. Thus, our findings indicate that each brain region develops at its own pace, with brain regions located in the same lobe showing similar developmental trajectories. We also observed almost identical patterns in subcortical gray matter volume as previous work (12,66), which showed inverted U-shaped trajectories that reached a peak around 15 years of age. Our findings indicate either an inverted U-shaped or positive developmental trajectory for individual subcortical regions. Further, our model showed that cortical thickness has a negative developmental curve for the vast majority of regions, with only some regions showing a fairly flat trajectory. This is consistent with earlier findings showing a decrease in cortical thickness across this age range (70) and largely consistent with other work, although the peak of cortical thickness was either reached at earlier (66) ages for average cortical thickness or later ages for average as well as regional cortical thickness (65,71). Further, our findings showed that the rate of change differs across brain regions. For example, regions in the occipital lobe have a steep decline in cortical thickness between 6 and 12 years of age, after which they reach a plateau, whereas regions in the parietal and cingulate cortex have a more linear decline over time, of which most do not yet reach a plateau at 16 years of age. Regarding surface area, earlier work indicated that across the entire cortex, development follows an inverted Ushaped trajectory across childhood and adolescence (65,66). In line with these findings, our results indicate similar trajectories for most regions. However, our findings also indicate that for regions mostly in the frontal and temporal lobe, surface area continues to increase until after 10-12 years of age, mirroring the regional differences in development observed for cortical volume and thickness. Our study has several strengths and limitations. First, a key strength of the current study is that we used normative modeling to establish z-scores representing deviations from normative development, which were subsequently related to psychopathology symptoms. Earlier work has shown that these regional deviations form typical development provide a greater prediction accuracy for psychopathology symptoms than using raw measures of brain morphology (17). Second, this study included a large age range, spanning early childhood to midadolescence, which allowed us to extend contemporary findings (12,70) in an age range that has not been studied extensively. Third, we were able to adjust our analyses for many potentially confounding factors. Although we interpreted the second model, some associations observed in model 1, corrected for biological sex and handedness, attenuated upon adjustment for SES and child national origin. These differences in the results obtained from the first and second model indicate the importance of adjustment for potentially confounding factors. Fourth, the population-based setting of the current study is both a strength and limitation. To derive typical developmental trajectories of brain morphology, a population-based sample is ideal. However, the majority of participants in this study exhibit relatively low levels of psychopathology, limiting the power to detect associations that might exist in those with clinical psychopathology levels. Fifth, while the current sample covers a large and important age range for developmental studies, the data acquired was not equally distributed across all ages. This can potentially influence the results around ages where fewer data was available. Sixth, although the current findings were derived from one of the largest population-based samples covering childhood and adolescence, our findings warrant replication in other comparable cohorts. Seventh, the data at T1 was obtained on a different MRI scanner than the data at T2 and T3, which may have influenced the results. Finally, a limitation of the current study is the amount of inter-individual variability in brain morphology measures and subsequently in the obtained z-scores. This results in small effect sizes in the associations we observed. However, small effect sizes are consistently reported for studies on brain morphology and psychopathology, and in those that have obtained large effect sizes, the effect size is often inflated by the small sample size (72). In summary, this study charted regional typical development of subcortical and cortical volume, surface area and cortical thickness. Findings showed that deviations from this typical development curve were related to psychopathology symptoms in widespread regions of the cerebral (sub-)cortex. DP symptoms were related to regional deviations from typical brain development above and beyond internalizing and externalizing symptoms in cortical volume. Our findings underline the evidence that assessing deviations from typical development in terms of (sub-)cortical volume and thickness can provide insights in the coupling between brain and behavioral development.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because Generation R data is available to researchers upon reasonable request. Requests should be directed to the management team of the Generation R Study. Individual level data are not publicly available for privacy and ethical restrictions. Requests to access the datasets should be directed to generationr@erasmusmc.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by METC Erasmus MC. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.