Abstract
Stress following preterm birth can disrupt the emerging foundation of the neonatal brain. The current study examined how structural brain development is affected by a stressful early environment and whether changes in topological architecture at term-equivalent age could explain the increased vulnerability for behavioral symptoms during early childhood. Longitudinal changes in structural brain connectivity were quantified using diffusion-weighted imaging (DWI) and tractography in preterm born infants (gestational age <28 weeks), imaged at 30 and/or 40 weeks of gestation (N = 145, 43.5% female). A global index of postnatal stress was determined based on the number of invasive procedures during hospitalization (e.g., heel lance). Higher stress levels impaired structural connectivity growth in a subnetwork of 48 connections (p = 0.003), including the amygdala, insula, hippocampus, and posterior cingulate cortex. Findings were replicated in an independent validation sample (N = 123, 39.8% female, n = 91 with follow-up). Classifying infants into vulnerable and resilient based on having more or less internalizing symptoms at two to five years of age (n = 71) revealed lower connectivity in the hippocampus and amygdala for vulnerable relative to resilient infants (p < 0.001). Our findings suggest that higher stress exposure during hospital admission is associated with slower growth of structural connectivity. The preservation of global connectivity of the amygdala and hippocampus might reflect a stress-buffering or resilience-enhancing factor against a stressful early environment and early-childhood internalizing symptoms.
SIGNIFICANCE STATEMENT The preterm brain is exposed to various external stimuli following birth. The effects of early chronic stress on neonatal brain networks and the remarkable degree of resilience are not well understood. The current study aims to provide an increased understanding of the impact of postnatal stress on third-trimester brain development and describe the topological architecture of a resilient brain. We observed a sparser neonatal brain network in infants exposed to higher postnatal stress. Limbic regulatory regions, including the hippocampus and amygdala, may play a key role as crucial convergence sites of protective factors. Understanding how stress-induced alterations in early brain development might lead to brain (re)organization may provide essential insights into resilient functioning.
Introduction
Preterm birth has long-lasting adverse effects on brain development and can increase the risk for psychiatric symptoms later in life (Spittle et al., 2009; Eikenes et al., 2011; Loe et al., 2013; Fischi-Gómez et al., 2015). The development of the preterm brain is contingent on several (clinical) factors (Doesburg et al., 2013; Ranger and Grunau, 2014; Ranger et al., 2015; Chau et al., 2019). Few longitudinal studies have explored the complex interaction between postnatal stress, brain development, and behavioral functions following preterm birth. We examined the impact of postnatal stress following preterm birth on brain development and how alterations in topological architecture influence vulnerability to behavioral symptoms during early childhood (two to five years).
Early-life adversities may alter brain maturation trajectories during a critical developmental period (Nelson and Gabard-Durnam, 2020; Lammertink et al., 2021). Cross-sectional studies investigating the effects of preterm birth on brain structure and function have shown lower white matter integrity in association tracts and projection fibers (Zwicker et al., 2013; Menegaux et al., 2017; Vollmer et al., 2017; Duerden et al., 2018). Preterm birth has also been related to an upregulation of functional connectivity between stress-related and stress-vulnerable regions. Specifically, the temporal cortex, thalamus, anterior cingulate gyrus, hippocampus, and amygdala showed higher functional connectivity following preterm birth (Papini et al., 2016; Johns et al., 2019; De Asis-Cruz et al., 2020).
Advances in graph theory enable researchers to reveal meaningful information about the topological architecture of the neonatal brain. Studies showed, for instance, that the properties of fundamental community structures (i.e., groups of densely connected regions reflecting subsystems or “building blocks” of a network) of a preterm born infant seem to be similar to typically developing fetuses and neonates (Song et al., 2017; Turk et al., 2019). Connectomic studies also highlight a more segregated and less integrated network organization in preterm-born infants (Ball et al., 2013a,b; Groppo et al., 2014; Sa de Almeida et al., 2021) and children (Fischi-Gomez et al., 2016; de Kieviet et al., 2021). The alterations in brain connectivity following preterm birth may play a significant role in developing future psychopathology (Van Essen and Barch, 2015; Kaufmann et al., 2017; Gilchrist et al., 2022).
In this study, we first examined the influence of stress on the development of premature brain connectivity and, second, whether alterations in macroscale network architecture at term-equivalent age may predict vulnerability for anxiety-related symptoms during early childhood (two to five years of age). We examined diffusion imaging and tractography from preterm infants and data on postnatal stress related to their hospitalization. We aim to identify specific differences in resilient and vulnerable infants that may enable resilient individuals to maintain relative mental wellbeing during early childhood.
Materials and Methods
Subjects
A total of 268 preterm infants were included and scanned between 28 and 32 and/or 39 and 42 postmenstrual age. Data collection was part of standard clinical care, with permission obtained to use these data for clinical research from the medical ethical review committee of the University Medical Center Utrecht (METC Utrecht). Preterm infants with chromosomal and/or congenital anomalies were excluded. Details and demographics of the main and validation datasets are outlined in Table 1.
Main dataset
Data of N = 145 preterm infants born infants clinically diagnosed as “extremely preterm” with a gestational age <28 weeks were included in our study, admitted to the Neonatal Intensive Care Unit (NICU) between 2013 and 2019 at the Wilhelmina Children's Hospital Utrecht, The Netherlands. Infants were scanned using 45 diffusion-encoding gradient directions.
Validation dataset
A replication sample containing N = 123 preterm infants born infants with a gestational age <28 weeks was included to assess the robustness of our results. Infants were admitted to the NICU between 2008 and 2013 and were scanned using 32 diffusion-encoding gradient directions.
Magnetic resonance imaging (MRI)
MRI data included the examination of 3T structural anatomic T2-weighted imaging and diffusion-tensor imaging (main dataset: dMRI, n = 45 directions; validation dataset, n = 32 directions; 3T Achieva MR scanner). Images were obtained as part of a 35-min scanning session.
T2 data were acquired using a Turbo Spin Echo (TSE) sequence, using parameters: TR = 6112 ms, TE = 120 ms, voxel resolution in millimeters 0.53 × 0.64 × 2 for 30 weeks and TR = 4851 ms, TE = 150 ms, voxel resolution in millimeters 0.78 × 0.89 × 1.2 for 40 weeks. dMRI data were acquired at 2 mm isotropic resolution and SENSE factor of 2 in 2 shells; 45 noncollinear directions for the main dataset, with a b-value of 800 s/mm2 and one nondiffusion-weighted image (non-DWI; b = 0) with TR 6500 ms and TE 80 ms; and 32 noncollinear directions for the validation dataset, with a b-value of 800 s/mm2 and one non-DWI (b = 0) with TR 5685 ms and TE 70 ms.
Infants were immobilized by wrapping them into a vacuum cushion. MiniMuffs (Natus Europe) and earmuffs (EM's kids Everton Park) were used to reduce noise and the infant's propensity to move during image acquisition. Before scanning, preterm infants scanned at 30 weeks were either sedated with 30 mg/kg oral chloral hydrate or not sedated at all, whereas infants scanned at 40 weeks were all sedated with 50–60 mg/kg oral chloral hydrate. Scanning was halted if the infant woke up, and clinicians attempted to re-settle the infant without taking them out of the patient immobilization system. A neonatologist or physician assistant was present at all times during the examination.
Data processing
Structural images
Volumetric tissue segmentation of gray and white matter, and labeling of subcortical and cortical areas, were performed on the T2 image (voxel resolution in millimeters 0.53 × 0.64 × 2 for 30 weeks and 0.78 × 0.89 × 1.2 for 40 weeks) using the structural pipeline from the developmental human connectome project (dHCP; http://www.developingconnectome.org/). The dHCP pipeline utilizes an “expectation-maximization” scheme that combines structure priors and an intensity model of the images (Makropoulos et al., 2018). A total of 47 (sub)cortical gray matter labels were automatically generated during segmentation (Fig. 1).
A total of 47 gray matter regions are segmented by the structural pipeline of the developmental Human Connectome Project (dHCP). HPL; hippocampus left, HPR; hippocampus right, AML; amygdala left, AMR; amygdala right, ATLML; anterior temporal lobe medial part left, ATLMR; anterior temporal lobe medial part right, ATLLL; anterior temporal lobe lateral part left; ATLLR; anterior temporal lobe lateral part right, GPAL; gyri parahippocampalis et ambiens anterior part left, GPAR; gyri parahippocampalis et ambiens anterior part right, STGL; superior temporal gyrus middle part left, STGR; superior temporal gyrus middle part right, MITGAL; medial and inferior temporal gyri anterior part left, MITGAR; medial and inferior temporal gyri anterior part right, LOGAL; lateral occipitotemporal gyrus/anterior fusiform left, LOGAR; lateral occipitotemporal gyrus/anterior fusiform right, CBL; cerebellum left, CBR; cerebellum right, BRS; brainstem, INL; insula left, INR; insula right, OLL; occipital lobe left, OLR; occipital lobe right, GPPL; gyri parahippocampalis et ambiens posterior part left; GPPR; gyri parahippocampalis et ambiens posterior right, LOGPL; lateral occipitotermporal gyrus/posterior fusiform part left, LOGPR; lateral occipitotermporal gyrus/posterior fusiform part right, MITGPL; medial and inferior temporal gyri posterior part left, MITGPR; medial and inferior temporal gyri posterior part right, STGPL; superior temporal gyrus posterior part left, STGPR; superior temporal gyrus posterior part right, CGAL; cingulate gyrus anterior part left, CGAR; cingulate gyrus anterior part right, CGPL; cingulate gyrus posterior part left, CGPR; cingulate gyrus posterior part right, FLL; frontal lobe left, FLR; frontal lobe right, PLL; parietal lobe left, PLR; parietal lobe right, CNL; caudate nucleus left, CNR; caudate nucleus right, THL; thalamus left, THR; thalamus right, SNL; subthalamic nucleus left, SNR; subthalamic nucleus right, LNL; lentiform nucleus left, LNR; lentiform nucleus right.
DWI tractography
DWIs were corrected for eddy current distortions, motion-induced signal drop-out, and head motion using a nonparametric approach using FSL (FSL EDDY; Andersson and Stamatios, 2016). The b0 image (voxel-size 2 × 2 × 2 for the main dataset, voxel-size 1.41 × 1.41 × 2.00 for validation dataset, b = 0 s/mm2) was registered to the T2-weighted image for anatomic alignment of the DWI images using FLIRT with a boundary-based-registration (BBR) cost function (Greve and Fischl, 2009). The linear transformation matrix was combined with a nonlinear warp registration using FSL FNIRT (Andersson et al., 2007) to map the diffusion space to an age-matched template. A single tensor model was used to estimate the main diffusion direction in each voxel (Basser et al., 1994) based on the 45 DWIs (b = 800 s/mm2; 32 directions for the validation dataset). An FA and MD whole-brain map was created based on the fitted tensors. White matter pathways were reconstructed using FACT (fiber assignment by continuous tracking; Mori and Van Zijl, 2002). Tractography involved starting eight streamline seeds in each white matter voxel, with fiber tracking, continued along the main diffusion direction of each voxel until a streamline showed high curvature (>65°), exited the brain mask, and/or when a streamline entered a voxel with low FA (<0.05). The mean FA value of a streamline was computed as the weighted average FA value, including all voxels that a streamline passed. Individual brain networks consisting of 47 gray matter regions and their interconnecting pathways were created by combining the subcortical and cortical segmentation map with reconstructed white matter tractography streamlines. For all combinations of regions, their interconnecting streamlines were mapped, with the weight of each region-to-region connection taken as the non-zero mean FA of the selected streamlines. Connections with a low connectivity strength (lowest 5%) were taken as potential false-positive reconstructions and set to 0. A group-based threshold was applied, retaining connections present in at least 50% of the participants, balancing the number of false-positive and false-negative structural connections (de Reus and van den Heuvel, 2013). Results were validated using different levels of group-based consensus thresholds (50–90%, steps of 5%).
Three summary measures were used to detect outliers among connectivity matrices: the presence of odd connections, the absence of common connections, and the average fractional anisotropy. We calculated the interquartile range (IQR) for each group separately by subtracting the 25th percentile from the 75th percentile (i.e., IQR = Q3–Q1). Participants with a score below Q1 – Q2 × IQR or above Q2 + Q2 × IQR for any of the three measures were considered outliers. This quantification led to the removal of seven outliers at 30 weeks of gestation and 19 outliers at 40 weeks of gestation.
Behavioral measures
Postnatal stress
NICU-related stress was quantified according to prior studies (Grunau et al., 2009; Duerden et al., 2018) and has been associated with alterations in stress-regulatory systems, including the hypothalamic-pituitary-adrenal axis and the autonomic nervous system (Grunau et al., 2004, 2006; Grunau, 2013; Provenzi et al., 2016; Lavanga et al., 2021), as well as childhood behavioral problems (Ranger et al., 2014). Data on invasive and stressful procedures were automatically extracted from the digital medical system. We generated a proxy for stress exposure by computing a global index of NICU-related stress using a principal component analysis on six parameters: skin-breaking procedures (i.e., heel lance, arterial and venous punctures, peripheral venous line insertion), total days of invasive mechanical ventilation, and suctioning of the nose and mouth. Each row (i.e., subject) was weighted on the total days of NICU admission. The main component explained 72.5% of the variance, with factor loading ranging from 0.74 to 0.91.
Residualized approach to postnatal stress
Participants were invited for standard clinical follow-up at 2.5 and/or 5.75 years of age. The main and validation datasets had behavioral follow-up data, resulting in a total of respectively 71 and 91 infants with an MRI at term equivalent age and data on behavioral symptoms (for an overview, see Table 2). During the clinical follow-up, parents reported the level of internalizing symptoms of their child, such as depression and anxiety, using the Child Behavior Checklist (CBCL; Achenbach and Rescorla, 2001). The CBCL is a parent-report questionnaire used to assess the frequency of dysfunctional behavior exhibited by the child in the past six months. Caregivers rate their children's behavior by answering questions about their child on a three-point scale (0–2), zero being “not true” one being “somewhat or sometimes true,” and two being “very true or often true.” If children did not have a behavioral symptom assessment at 5.75 years of age, we used the 2.5 years assessment (moderate correlation between the two time points; r = 0.45, p < 0.001; Fig. 2A). The follow-up also included other assessments not part of the current study, such as motor development and intelligence.
A, Significant and positive association between internalizing symptoms assessed at two and five years of age. B, Residualization approach; orange observations are categorized as stress-overreactve (vulnerable), and green observations are characterized as stress-underreactive (resilient).
Resilience was quantified as a metric of mental health by indexing the internalizing symptoms subscale of the CBCL, taking into account the degree of NICU-related stressor exposure using simple linear regression. We observed a significant positive association between postnatal stress and early childhood internalizing symptoms (t(11,151) = 4.08, p < 0.001). The fitted regression line (Fig. 2B) reflected the normative level, with participants positioned above the linear line (i.e., positive residual) expressing an over-reactivity of behavioral symptoms to stressor exposure in the neonatal period and data points below the linear line (i.e., negative residual) representing individuals with under-reactivity to stressor exposure (Amstadter et al., 2014; Van Harmelen et al., 2017).
Preterm-born individuals were classified accordingly: the resilient group showed fewer behavioral symptoms than expected, and a vulnerable group showed more behavioral symptoms problems than expected.
Experimental design and statistical analyses
Analyses (connectome development and group differences, see below) were corrected for confounding factors, including gender, birthweight (z scores [Hoftiezer et al., 2019]), mean FA, gestational age, age at scan, degree of brain injury (i.e., intraventricular hemorrhage), neonatal surgeries, administration of prenatal and postnatal corticosteroids (i.e., accelerates lung maturation), and days of morphine. There was no joint dependency among covariates (variance inflation factor ≤1.5).
Stress and connectome development
Longitudinal changes in whole-brain structural connectivity between 30 and 40 weeks of gestation were examined using a time × postnatal stress interaction model using network-based statistic (NBS), a permutation-based method specifically designed to statistically assess network differences (Zalesky et al., 2010). We created an NBS linear-mixed model (R package NBR version 0.1.4; Gracia-Tabuenca and Alcauter, 2020) that took into account the repeated-measures nature of the structural connectivity data by including a per-participant random intercept and by modeling all within-subject covariates [gender, gestational age, age at scan, degree of brain injury (i.e., intraventricular hemorrhage), surgeries, administration of prenatal and postnatal corticosteroids (i.e., accelerates lung maturation), and administration of morphine in days] and predictors (i.e., postnatal stress, time) as fixed effects. The model was applied to all non-zero Ni × Nj connections of the individual networks (lower triangle; consensus-based threshold). The N × N matrix of F-statistics and matching p-values associated with the interaction effect was thresholded at a p-value of p < 0.05. NBS defines the largest connected component, and the size of the largest component is tested against a null-model of permuting subject labels 10,000 times. The subsequent null distribution was used to calculate an FWE-corrected p-value for the largest identified component. We used the main sample and validated the findings in a separate, independent population (Table 1).
Group differences between resilient and vulnerable individuals
Differences in network organization between resilient and vulnerable individuals were assessed by examining global and local network metrics from the individual structural matrices at term-equivalent age (R packages igraph, braingraph; R Core Team, 2021). A separate general linear model was specified to test for significant group difference for each network metric and was compared with permuted data (four graph-level and four vertex-level; 47 regions), building a null distribution. Graph-level analyses were permuted 10,000 times, and vertex-level measures were permuted 5000 times. The contrast was thresholded on false discovery rate (FDR) adjusted p < 0.001 to correct for multiple comparisons.
Local graph parameters, including clustering coefficient, nodal efficiency, betweenness centrality, and communicability (Rubinov and Sporns, 2010), were calculated. Global measures included (mean) clustering coefficient, modularity, strength, and global efficiency. Clustering coefficient describes the tendency of regions to cluster together in triangles and is computed by the ratio between the number of connections between region i and its neighbor regions and the total number of possible connections with neighbors (Rubinov and Sporns, 2010). The global measure is computed by taking the mean clustering coefficient of all individual regions in the network. Nodal efficiency describes for every region in the network the inverse harmonic mean of the minimum path length between a given region i and all other regions j (Achard and Bullmore, 2007). Higher nodal efficiency is indicative of a higher capability of information integration, and these regions can also be categorized as a hub. The global measure is computed by taking the mean of nodal efficiency of all individual regions in the network. Betweenness centrality describes the influence of a region in the communication between pairs of regions and is measured by the frequency with which a region falls between pairs of other regions on their shortest interconnecting path (Rubinov and Sporns, 2010). This measure reflects the potential influence of controlling global connectivity between nondirectly connected regions. Communicability describes how well a region communicates with every other region in the network and is computed by the weighted sum of all paths and walks between regions i and j (Estrada and Hatano, 2008). High communicability indicates multiple and strong alternative paths connecting the region with other regions. Modularity describes the degree to which a network can be organized into modules of densely interconnected regions but sparsely connected between modules and is computed by the difference between the number of edges within a community and a random network of the same degree sequence (Rubinov and Sporns, 2010). High modularity reflects a highly segregated network. Strength describes the total sum of the weights of all individual nodal connections in the network. Together, these provide a good understanding of the connectivity and influence of a particular region on the network.
Predictive classification
Random-forest regression with conditional inference trees (RFR-CIT) was used to assess the level of predictive power of node-wise centrality measures in correctly differentiating resilient (stress-underreactive) and vulnerable (stress- overreactive) individuals. RFR-CIT was used as it has been shown to have advantages in working with relatively small sample sizes, to be relative insensitive to multicollinearity, and to be less sensitive to overfitting (Breiman, 2001; Liaw and Wiener, 2002). The optimal performance of the model (i.e., feature selection and parameter tuning) was assessed using ten times repeated 10-fold cross-validation [65% of the dataset (n = 105), R package caret version 6.0–92; Kuhn, 2008)]. We divided the training dataset randomly into 10 folds, training the model on 90% of the data and validating on the rest. In each fold, a different one-tenth of the data were selected for validation. The training data sample is shuffled before each repetition, and results are aggregated across all folds from all runs (totaling 100 different models). The optimal model is selected with the smallest root mean squared error. The performance of the trained model was finally tested using a holdout dataset (35% of subjects; n = 57).
The first predictive classification model combined all four local centrality measures (i.e., betweenness centrality, communicability, nodal efficiency and clustering coefficient) of 47 gray-matter nodes. Analyses were repeated for each node-wise centrality measure separately. Differences in features because of technical variability in acquisition protocol were removed while preserving biological variability using ComBat (R packages neuroCombat; R Core Team, 2021) before model fitting using empirical Bayes and parametric adjustment (Johnson et al., 2007; Fortin et al., 2017, 2018). We estimated ComBat parameters separately for the training and test dataset to prevent data leakage. The statistical significance of a difference in accuracy between the two models was assessed using the McNemar test.
Results
The sample consisted of a main (N = 145, Mage = 26.53, Msd = 0.97, 43.5% female) and validation (N = 123, Mage = 26.54, Msd = 1.00, 39.8% female) dataset of preterm born individuals. Both the main (n = 71) and validation (n = 91) dataset have follow-up data on parent-reported internalizing symptoms. Key demographics of the two samples are presented in Table 1.
Sample demographic and neonatal clinical details of participants (N = 268)
The effects of postnatal stress on the development of whole-brain structural connectivity
We performed NBS (for details, see Materials and Methods) to identify subnetworks of edge-wise effects that showed significant alterations in growth depending on the degree of postnatal stress exposure. NBS analysis revealed one significant cluster of connections, involving 48 connections, with slower growth in connectivity strength from 30 to 40 weeks of gestation for individuals exposed to higher stress (p = 0.003, consensus-based threshold; Figs. 3A, 4A). The cluster spanned both hemispheres, involving 20 brain regions such as the amygdala, thalamus, caudate nucleus, and cortical regions such as the insula, fusiform, parahippocampal gyrus, anterior/posterior cingulate cortex, parietal lobe, and frontal lobe. Figure 3 provides a matrix of the vertices and edges involved. The subnetwork reduced in size but remained significant across prevalence thresholds (Fig. 4E). Also, postnatal stress significantly affected white-matter connectivity at term-equivalent age, with higher stress resulting in lower structural connectivity in a subnetwork of 49 connections (Fig. 3B, 4C, p = 0.014).
Matrix of largest significant subnetwork (NBS threshold of 50%) identified in the main dataset. Edges with a cross are part of the subnetwork showing a significant time × stress effect (matrix shows the δ in mean connectivity between 30 and 40 weeks of gestation; A) or a significant main effect of stress (B) at term-equivalent age. An overview of abbreviations can be found in Figure 1. FA = fractional anisotropy.
A, Schematic representation of the longitudinal time × stress effects, orange representing high stress (highest 25%) and green representing low stress (lowest 25%). B, Replication of time × stress effects in an independent sample (32 directions diffusion protocol). C, Negative effect of postnatal stress on structural connectivity at term-equivalent age (included 49 connections; 45 directions diffusion protocol). D, Replication of stress effects in an independent sample (32 directions diffusion protocol). E, Robustness of NBS findings across a range of prevalence thresholds (50%: p < 0.05; prevalence threshold of 60%: p < 0.05, prevalence threshold of 70%: p < 0.05, two-sided permutation testing, 10,000 permutations).
Replication
The NBS findings were replicated in an independent sample, providing robust evidence for the effects of postnatal stress on the growth of white-matter connectivity. We masked the connectivity matrix such that connections were retained if they were part of the subnetwork identified in the main sample. Then, we calculated a non-zero mean of connectivity strength and tested the effects of postnatal stress on changes in connectivity strength between 30 and 40 weeks of gestation. We observed a significant stress × time interaction such that higher levels of postnatal stress were associated with slower growth in connectivity strength (estimate = −0.007(0.003), F(1,37) = 4.79, p = 0.035, 95% CI [−0.014, −0.001]; Fig. 4B). Higher stress was associated with significantly lower levels of white-matter connectivity at term-equivalent age (t(13,96) = −2.44, p = 0.016; Fig. 4D).
Network architecture at term-equivalent age reveals differences between resilient and vulnerable individuals
A total of 41 and 42 neonates were classified as stress under-reactive (now being referred to as resilient), and 30 and 49 infants were classified as stress over-reactive (now being referred to as vulnerable) based on the normative levels of stress-reactivity (the relationship between postnatal NICU-related stress and long-term behavioral symptoms; see above, Residualized approach to postnatal stress). There were no group differences in birth weight, age at birth, age at scan, corticosteroids, days of morphine administration, and mean FA (Table 2). There was, however, a slight difference in gender in the main dataset (included as a covariate). The reported findings below were thresholded on 75% prevalence, i.e., connections were included if they were reported in at least 75% of the participants. The reported results are based on structural connectivity at term-equivalent age.
Sample demographic and neonatal clinical details of resilient and vulnerable infants
Global graph theoretical measures
Analyses revealed no significant group effects in measures of global network architecture.
Local graph theoretical measures
We observed significant group effects on local network measures. Group differences were region-specific such that both reduced and increased centrality were observed in vulnerable relative to resilient individuals (Table 3).
Group-difference on nodal centrality measures for contrast vulnerable > resilient
We first examined the contribution of regions in local network organization as measured by “nodal clustering.” Vulnerable infants, relative to resilient, showed a lower clustering of several cortical brain regions overall, including the posterior cingulate cortex (t(69) = −5.48, p < 0.001), parahippocampal gyrus (t(69) = −5.25, p < 0.001), frontal lobe (t(69) = −6.29, p < 0.001), and parietal lobe (t(69) = −7.45, p < 0.001). In contrast, higher clustering was observed in the hippocampus (t(69) = 7.19, p < 0.001), amygdala (t(69) = 4.8, p < 0.001), and medial anterior temporal lobe (t(69) = 6.3, p < 0.001). It is important to note that only differences in the posterior cingulate cortex and parietal lobe were successfully replicated in the validation sample. Statistical details of group differences found in the main and validation dataset can be found in Table 3.
We assessed the contribution of regions in global communication across the brain through “betweenness centrality.” On average, vulnerable infants showed a lower centrality of the hippocampus (t(69) = −9.5, p < 0.001) and the anterior fusiform (t(69) = −7.45, p < 0.001), whereas a higher centrality was observed in the brain stem (t(69) = 3.76, p < 0.001), posterior cingulate cortex (t(69) = 5.72, p < 0.001), and parietal lobe (t(69) = 6.11, p < 0.001; Table 3). These results suggest differential susceptibility in connections central to global brain communication.
We further examined global network connectivity through “communicability.” Vulnerable individuals showed, on average, lower communicability of the hippocampus (t(69) = −16.03, p < 0.001), amygdala (t(69) = −3.74, p < 0.001), medial anterior temporal lobe (t(69) = −5.73, p < 0.001), and subthalamic nucleus (t(69) = −11.44, p < 0.001; Fig. 5 (created with Simple Brain Plot; Scholtens et al., 2021); Table 3). A higher global connectivity was observed in the posterior fusiform (t(69) = 9.65, p < 0.001) and parietal lobe (t(69) = 3.73, p < 0.001). Bivariate analyses, excluding covariates, showed similar results (e.g., hippocampus (t(69) = −18.05, p < 0.001), amygdala (t(69) = −5.49, p < 0.001), parietal lobe (t(69) = −11.44, p < 0.001).
Distribution of group differences in communicability values between vulnerable and resilient infants (left) and regions colored according to T-value [right; vulnerable < resilient (green); Simple Brain Plot, Scholtens et al., 2021].
Resilient and vulnerable infants did not differ on measures of nodal efficiency.
Classification model
RFR-CIT was used to investigate the potential predictive power from local network metrics. The combined set of local network measures (e.g., communicability) of the 47 (sub)cortical gray matter regions was able to correctly classify vulnerable and resilient individuals with an accuracy of 83.9% (p < 10−5, κ = 0.679, area under the curve (AUC) = 0.919; Table 4). For model classification and calibration, see Figure 6. Individual centrality measures showed a distinct pattern with higher accuracy for measurements of global connectivity (e.g., communicability) than local connectivity measures (e.g., nodal efficiency; for an overview, see Table 4).
Classification of vulnerable versus resilient infants using the (35%) holdout dataset
A, Difference between the mean of predicted probabilities of the vulnerable and resilient group. B, True frequency of the positive label against its predicted probability, the x-axis represents the average predicted probability in each bin and the y-axis is the proportion of sample whose class is the vulnerable class (fraction of positives). Error bars depict 95% confidence intervals.
Discussion
Preterm-born infants have a life-long increased risk for stress-related psychopathology characterized by anxiety and socio-emotional problems (Arpi and Ferrari, 2013; Upadhyaya et al., 2021). Our findings suggest that higher stress exposure during NICU admission is associated with slower growth of structural connectivity in regions such as the amygdala, hippocampus, insula, and posterior cingulate cortex. However, some infants seem to be less affected by postnatal stress exposure; resilient infants at term-equivalent age show higher global connectivity of regions central for bottom-up emotion regulation. The classification model showed an excellent predictive accuracy of group membership, early childhood resilience, or vulnerability, using local network measures at term-equivalent age shortly following exposure. To summarize, the postnatal, stressful environment contributes to significant alterations in brain development. However, only a proportion of infants show a higher susceptibility to behavioral problems in early childhood.
Our findings underscore the impact of postnatal stress on the growth of structural brain connections in corticolimbic pathways across both hemispheres. The affected brain regions are involved in (bottom-up) emotion regulation and processing, such as the amygdala, insula, hippocampus, parahippocampal gyrus, and posterior cingulate cortex. These findings align with evidence from other neuroimaging studies showing a delayed development in white matter pathways following preterm birth relative to full-term controls (Dodson et al., 2017; Bouyssi-Kobar et al., 2018; Duerden et al., 2018). We show evidence that in addition to the effects of prematurity, stressful early exposure significantly contributed to a more pronounced impact on delayed development in a subnetwork of connections.
Our findings suggest that higher global connectivity characterizes resilience in several brain regions, including the hippocampus. Measures of global connectivity estimate communication efficiency among all nodes in a network, enabling the integration and distribution of neural information between spatially distant brain regions (Sporns, 2018). Vulnerable individuals showed lower global connectivity of the hippocampus. This region is a crucial regulator of the hypothalamic-pituitary-adrenal axis activation and plays a critical role in storing and retrieving emotional memories (Chan et al., 2014; Duval et al., 2015). Prior studies on (early-life) trauma indicated that the hippocampus is particularly vulnerable to chronic pain and stress, with lower volumes and a hypoconnectivity following early-life trauma (Andersen et al., 2008; Shin and Liberzon, 2010). We observed a similar pattern for the amygdala and subthalamic nucleus. The amygdala is part of the (medial) temporal lobe and densely connected with the prefrontal cortex (PFC) and has extensive anatomic connections with the paraventricular thalamus and hippocampus. This region plays a critical role in perception, regulation, and plasticity of emotion (Davis and Whalen, 2001; Yang et al., 2017). A less interconnected amygdala in vulnerable infants might seem contradictory, as it does not agree with studies showing evidence of lower amygdala connectivity in resilient trauma-exposed adults (Roeckner et al., 2021). However, a less interconnected amygdala might also be evidence of a decreased inhibitory control of more segregated cortical regions, including the ventromedial PFC (vmPFC; Johnstone et al., 2007; Rogers et al., 2017; Andrewes and Jenkins, 2019). The lower centrality of the frontal lobe in vulnerable infants substantiates this interpretation. Hence, the increased global connectivity of the hippocampus and amygdala might be a key system in a healthy adaptation with resilience following early disturbances of preterm birth.
Previous studies in children and adults with depression and anhedonia consistently reported a lower capacity for global connectivity (Cullen et al., 2014; Yang et al., 2017). The subthalamic nucleus interconnects with the amygdala and hippocampus, receives convergent cortical and pallidal projections (Accolla et al., 2016), and plays a role in threat appraisal (Serranová et al., 2011). Although the subthalamic nucleus received attention concerning Parkinson's disease, the increased social and affective alterations following deep-brain stimulation have been implicated in the emergence of enhanced affective processing and decreased depressive symptoms (Schneider et al., 2003; Smeding et al., 2006). Higher global connectivity of these regions might be beneficial in retaining short-term mental wellbeing following preterm birth. Although the current study contributes to the resilience literature following preterm birth, our findings cannot be generalized to developmental stages other than early childhood. Future studies investigating the (neuroprotective) mechanisms by which global connectivity is higher in resilient infants, children, and adolescents are warranted.
There is no direct inverse relationship between resilience and vulnerability (Bolsinger et al., 2018). Because of the lack of resilience studies following preterm birth, vulnerability studies seem to present the best available approximation for the concept of resilience in preterm-born individuals. In line with studies on trauma exposure, preterm-born individuals with more problem behavior seem to show reduced hippocampal connectivity and lower volumes (Aanes et al., 2015; Rogers et al., 2018) and a lower interconnected amygdala (Rogers et al., 2017). Further, preliminary interventional studies focusing on neuroprotection, such as music and massage therapy, reported significantly improved white matter maturation of the uncinate fasciculus (Sa de Almeida et al., 2020). Hence, these results implicate that higher global connectivity of the amygdala and hippocampus may lead to symptom attenuation and is consistent with our observation that infants differ in a small number of regions or pathways that may facilitate resilience.
The topological measures provide a wealth of information and have the prominent capability of characterizing the behavior and organization of complex brain systems (Farahani et al., 2019). The differences in neural representations between resilient and vulnerable infants enable the accurate classification of group membership. The current study shows that the connectome has distinguishable features in topological architecture at term-equivalent age. Specifically, measures of global connectivity (e.g., communicability) seem to have higher predictive accuracy than local connectivity measures (e.g., nodal efficiency). This feature suggests that especially edges involved in communication among spatially distant regions are central to resilience's (short-term) etiology. Demonstrating these patterns highlights that resilience and vulnerability occur in the context of unique neurobiological differentiability. However, these findings are not evidence of a neuroimaging-based biomarker of resilience. There is no convergence on the proposed resilience mechanisms, and it is unclear whether alterations are a hallmark of general pathology following prematurity.
The structural alterations in emotion regulation and processing are not limited to postnatal stress. Studies have consistently reported alterations in limbic and prefrontal microstructure in neonates and infants exposed to third-trimester prenatal stress (e.g., maternal anxiety; Rifkin-Graboi et al., 2013, 2015; Posner et al., 2016; Dennis et al., 2019; Lautarescu et al., 2020). Unfortunately, the current study does not include information on prenatal stress exposure. Although prenatal (psychosocial-induced) and postnatal (pain-induced) stress are different dimensions of distress, the commonalities in alterations in brain development following exposure are not surprising. Both types of stressors share significant conceptual and physiological overlaps. Both phenomena challenge the body's homeostasis (Mörelius et al., 2016; Jahnke et al., 2021; Lammertink et al., 2021) and can be conceptualized as a stress overload, termed allostatic overload, with the brain's limbic system standing as a central mediator. This aligns with the striking overlap in neuroanatomical networks involved with chronic pain and stress modulation, specifically the corticolimbic system, including the hippocampus, thalamus, amygdala, and vmPFC (Timmers et al., 2019). Hence, these commonalities suggest a significant role of stress physiology in the onset and persistence of chronic pain (for extensive review, see Abdallah and Geha, 2017).
Several methodological issues should be taken into consideration when interpreting our findings. First, the preterm-born population is increasingly susceptible to (chronic) disease. The degree of illness and stress exposure is often tightly linked. Although we adjusted for clinical confounders of prematurity, it remains a challenge to disentangle the effects of illness and stress. Several uncontrolled confounders (including socioeconomic status, race, and ethnicity) might hinder the generalizability of the current results. Noteworthy, we did not find an association between parental educational level and child internalizing symptoms (maternal education: χ2(29,147) = 30.24, p = 0.402; paternal education: χ2(29,144) = 35.18, p = 0.199; for educational levels see Table 2 (Steenis et al., 2015)). We were unable to control for prenatal stress exposure (other than prenatal glucocorticoid administration), such as maternal anxiety. It is unclear whether the effects of postnatal stress can be accounted for by the continuation of stress from the prenatal period. Second, although our unique dataset enables us to investigate individual differences in longitudinal white-matter development, we could only model linear change. Studies involving three or more time points (Remer et al., 2017) can fit several growth slopes, including quadratic, facilitating a more nuanced understanding of how postnatal stress affects brain development (King et al., 2018). For instance, a quadratic growth pattern would mean that the effects of postnatal stress emerge during a specific developmental period and then decline or disappear during a particular period and then reappear later. Despite this methodological limitation, our results provide convincing evidence that between 30 and 40 weeks of gestation, postnatal stress significantly reduces linear growth in a subnetwork of connections. Third, resilient infants might have been healthier than vulnerable infants. Although the stress-residualization approach controls for postnatal stress exposure, resilient infants could still have experienced fewer complications and clinical procedures than vulnerable infants. Bivariate analyses showed that infants did not differ on a large set of clinical parameters (Table 2). Fourth, the small sample size and lack of a diverse geographical population might hinder the generalizability and accuracy of the classification model. Another important factor is the choice of validation method. Although we used a widely accepted train-validate-test approach, the initial train/test split might induce biased model estimates and explain the higher reported classification accuracy. Despite these valid limitations, a study compared the degree of overfitting in small sample sizes using several validation methods [e.g., (nested/partially nested) cross-validation, train/test split], and they reported that both nested and train/test split validation approaches results in robust and unbiased performance estimates (Vabalas et al., 2019).
Our longitudinal findings suggest that postnatal stress leads to sparser brain connectivity after preterm birth. Alterations in specific brain areas impacting bottom-up emotion regulation might render preterm infants resilient to internalizing symptoms in early childhood. These findings show the detrimental impact of postnatal stress, but we also recognize the relative plasticity of the preterm brain in structural changes that indicate resilience. The current results suggest that resilience might appertain to a preserved global connectivity, which is a first step in identifying the biological basis of mental wellbeing following preterm birth.
Footnotes
F.L. was supported by the Wilhelmina Children's Hospital Grant D-17-010007. M.P.v.d.H. was supported by the Netherlands Organization for Scientific Research (NWO) VIDI Grant 452-16-015 and the European Research Council Grant ERC-2020-CoG 101001062. E.J.H. was supported by the European Research Council Grant ERC-2015-CoG 682591.
The authors declare no competing financial interests.
- Correspondence should be addressed to Femke Lammertink at femke.lammertink{at}gmail.com
