Associations Between Neonatal Brain Structure, the Home Environment, and Childhood Outcomes Following Very Preterm Birth

Background Very preterm birth is associated with an increased risk of childhood psychopathology and cognitive deficits. However, the extent to which these developmental problems associated with preterm birth are amenable to environmental factors or determined by neurobiology at birth remains unclear. Methods We derived neonatal brain structural covariance networks using non-negative matrix factorization in 384 very preterm infants (median gestational age [range], 30.29 [23.57–32.86] weeks) who underwent magnetic resonance imaging at term-equivalent age (median postmenstrual age, 42.57 [37.86–44.86] weeks). Principal component analysis was performed on 32 behavioral and cognitive measures assessed at preschool age (n = 206; median age, 4.65 [4.19–7.17] years) to identify components of childhood psychopathology and cognition. The Cognitively Stimulating Parenting Scale assessed the level of cognitively stimulating experiences available to the child at home. Results Cognitively stimulating parenting was associated with reduced expression of a component reflecting developmental psychopathology and executive dysfunction consistent with the preterm phenotype (inattention-hyperactivity, autism spectrum behaviors, and lower executive function scores). In contrast, a component reflecting better general cognitive abilities was associated with larger neonatal gray matter volume in regions centered on key nodes of the salience network, but not with cognitively stimulating parenting. Conclusions Our results suggest that while neonatal brain structure likely influences cognitive abilities in very preterm children, the severity of behavioral symptoms that are typically observed in these children is sensitive to a cognitively stimulating home environment. Very preterm children may derive meaningful mental health benefits from access to cognitively stimulating experiences during childhood.


MRI acquisition and processing
Infants underwent MRI at term-equivalent age on a Philips 3T system (Philips Medical Systems, Best, The Netherlands) using an eight-channel phased array head coil. High resolution anatomical T2-weighted turbo spin echo was acquired with TR = 8670 ms, TE = 160 ms, flip angle 90°, slice thickness of 2 mm with 1 mm overlap, and in-plane resolution 0.86 × 0.86mm.
Images from successful T2 scans (N=486) were visually inspected to detect and exclude images with artefacts (N=10). Infants with major lesions (N=40) or PMA at scan ≥ 45 weeks (N=50) were excluded in order to minimise registration failures. After these exclusions, a further 2 subjects were removed due to processing failures. The processing pipeline included field bias correction and tissue type segmentation using the neonatal processing pipeline described in Makropoulos et al. (2014) (1). A study-specific template was generated from a subset of 161 participants using Advanced Normalization Tools (ANTS) software (2). Images were then registered to the study-specific template using multimodal Symmetric Normalisation (SyN) (2). Both T2-weighted images and T2-based tissue type segmentation were used as input modalities in order to improve image registration (3).
Deformation tensor fields (i.e. warps) from the non-linear registration were used to obtain a logarithm transformation of Jacobian determinant maps, reflecting local expansion/shrinkage of each voxel with respect to the template (4). Log Jacobian maps were spatially smoothed with a Gaussian kernel of 4mm full width at half maximum. In order to reduce computational load, maps were down-sampled to 2mm isotropic resolution. To further reduce the number of included voxels and aid in localisation of effects, only brain tissue voxels defined within the neonatal version of the Automated Anatomical Labelling (AAL) atlas (5,6) were included in the analysis, resulting in an inclusion 37947 voxels. All voxel values were exponentiated before submitting to NNMF analysis in order to ensure non-negative input data.

Cognitively Stimulating Parenting Scale
Parents completed a questionnaire adapted from the Cognitively Stimulating Parenting Scale reported in Wolke et al. (7). It consists of 21 items included in the Home Observation for Measurement of the Environment (HOME) Inventory (8) and has been shown to have acceptable internal consistency (Cronbach α = 0.77) (20). Briefly, the Cognitively Stimulating Parenting Scale assesses the availability and variety of experiences that promote cognitive stimulation in the home. 16 binary items capture the child's access to stimulating objects such as educational toys (6 items), parental interactions such as teaching numbers or colours (7 items), and parental behaviours such as reading habits (3 items). One item ("access to cassette player") was adapted to reflect technological advances appropriate at time of testing ("access to cassette/CD/DVD player"). Five additional Likert-scale items assess the frequency of cognitively stimulating experiences such as family excursions (4 items) and number of books in the home (1 item). All questionnaire items can be found in Supplemental Table 1.
Yes/No items were scored as 1/0, and Likert-scale items were scored as 0 for answers 0-3 and as 1 for answers 4-6. A total sum score was calculated from all items (minimum 0, maximum 21).

Principal component analysis (PCA)
Significance of the principal components (PCs) was assessed using permutation testing. For each of 5000 permutations, the subject order for each column was randomly shuffled and PCA recomputed, saving the proportion of variance explained by each permutation PC, thereby establishing the null model. Original PCs were deemed significant if the proportion of variance they explained exceeded chance (p < .05 across 5000 permutations).
In order to assess which of the significant PCs was consistent across independent subsamples, we conducted a split-half replication. For each of 1000 iterations, we randomly split the sample into two non-overlapping sets of subjects and re-computed the PCA in each subsample. Correlations between the PC loadings from the two PCA solutions were evaluated for reproducibility, and significance of PCs for each sub-sample was again assessed using permutation testing as described above.
To assess meaningful contribution of individual variables to the identified PCs, we defined the loading threshold as � 1 = 0.18 , equivalent to the loading value if all 32 variables included in the PCA loaded equally onto the same component (since the sum of squares of all loadings for one PC must sum to 1). Loadings greater than this value were considered to be meaningful. We also inspected correlations between PC loadings and original input variables (outcome assessments) to aid interpretation of the components.

Non-negative Matrix Factorisation
Non-negative matrix factorisation (NNMF) is an unsupervised multivariate analysis technique which aims to obtain a low-rank approximation of high-dimensional data. NNMF decomposes a non-negative input matrix into the product of two low-rank factor matrices (W and H). Due to non-negativity constraints imposed on each element of the resulting matrices, this procedure results in a parts-based representation of the data, whereby the parts are additively combined to reproduce the original input. This is a particularly advantageous property of NNMF as it increases interpretability and specificity of the resulting components compared to other dimension reduction techniques. It is especially suited to identifying components of covariance within data that is inherently non-negative, such as brain structural measures (9).
Here, we aimed to identify structural networks in the neonatal preterm brain in which regional brain volumes co-vary across individuals in a consistent way (N=384). We thus used a matrix containing voxelwise jacobian values (dimensions: 37947 voxels × 384 subjects) as input to the NNMF algorithm, resulting in a factorisation into a matrix W (dimensions: 37947 voxels × k networks), representing the individual contribution of each voxel to each of k networks, and a matrix H (dimensions: k networks × 384 subjects), containing subject-specific scores for each of k networks. The procedure to estimate the optimal rank k is detailed in the following.

Rank estimation procedure for Non-negative Matrix Factorisation (NNMF)
First, we ran NNMF on the original data for ranks 2 to 20, in steps of 1. NNMF was run 50 times for each rank, each time randomly masking 20% of values within the input matrix. Next, we randomly permuted the elements within each column of the input matrix, and repeated the NNMF estimation of this permuted dataset for ranks 2 to 20 a total of 50 times each, again masking a random subset of 20% of elements each time. For each NNMF estimation (both on the original and permuted datasets), we calculated the reconstruction error for each solution, defined as the Frobenius norm between the input data matrix (original or permuted) and the resulting NNMF approximation (W × H). Additionally, we calculated the root mean square error (RMSE) for the 20% held-out datapoints as well as the included datapoints for each estimation.
The reconstruction error is expected to decrease as a function of increasing rank k, both for the original data and randomly permuted data. We therefore aimed to ascertain the optimal rank k at which the decrease in reconstruction error for the original data no longer exceeds the decrease in reconstruction error for the permuted data (i.e., the point at which the decrease in error is no more substantial than the decrease observed in random data). For this, we performed a t-test at each rank comparing the gradient of the reconstruction error (capturing the decrease in error) for the original data and the permuted data (each consisting of 50 runs at that rank), and chose the highest rank at which this difference between original and permuted data was significant.
NNMF was run in R using the NNLM and NMF packages. NNLM was used (with initialisation based on random seeding) during the rank optimisation procedure, as it allows for individual matrix elements to be masked during approximation. NMF was used for the final estimation of structural covariance networks, using the previously ascertained optimal rank (with initialisation based on non-negative double singular value decomposition [NNDSVD] and no data held out).

Principal components of childhood outcomes
PCA on outcome variables with permutation testing identified four significant PCs (PC1-PC4) in the full sample, explaining a cumulative 64% of total variance. Further permutation testing showed PC1 and PC2 to be consistently significant in both halves across 1000 split-half analyses (100% of split-halves). PC3 was significant in both halves in 62.8% of split-half analyses, and significant in only one half in a further 37.1%. PC4 was inconsistently significant (both halves: 0.1%; one half: 45.3% of analyses). The mean correlation between loadings in the two halves were high for PC1 (R = .98) and PC2 (R = .85), medium for PC3 (R = .46) and low for PC4 (R = -.03). We therefore retained PC1, PC2, and PC3 for further analysis, jointly explaining a cumulative proportion of 59% of total variance. Individual behavioural loadings (thresholded at 0.18) as well as significant correlations with outcome variables for PC1-PC3 are depicted in Figure 1. Numerical values of all loadings can be found in Supplementary   Table S1.

Results of rank selection for NNMF
As expected, reconstruction error decreased as a function of rank both for the original and permuted data ( Figure S1A). Figure S1B shows the superimposed gradients of the reconstruction errors for original and permuted data, each estimated 50 times under 20% Wold holdouts. The difference between mean reconstruction error gradients was significant for ranks 1-15 and non-significant thereafter, indicating that at ranks 16 and above the information gained by adding one component no longer exceeds that gained in random data. Inspection of the gradient of RMSE for included (training) and held-out (test) datapoints also showed a levelling-out at around 15 networks ( Figure S2). We therefore used the 15-network solution for all further analyses.

Removal of twins and triplets
We repeated all analyses after removing possible effects of the high proportion of twin and triplet sets in our data. Where a complete twin set, or at least two triplets from the same triplet set were present, we removed one twin or triplet (or two triplets, in the case of complete triplet sets) at random, and repeated all analyses. Before testing the effect of SCN volumes on behavioural outcome PCs, we removed 13 individual twins and 3 individual triplets, resulting in a total sample of N=140 (down from N=157) individuals. Controlling for multiple comparisons (15 models), SCN 12 was, again, the only SCN showing a significant multivariate effect on behavioural outcomes (p < .001), which was driven by its univariate effect on PC2 ("cognitive" component) (β = 5.58, p = .007) (see Table S3).
Before testing the effect of cognitively stimulating parenting on behavioural outcomes, we removed 23 individual twins and 4 individual triplets, repeating the analysis on N=175 (down from N=206) individuals. Multivariate and univariate effects are found in Table S4, showing that a multivariate effect of cognitively stimulating parenting was driven by its effect on PC1 ("preterm phenotype"), in line with the main analysis.

Controlling for lesions
We repeated the analysis controlling for severity of brain lesions (major / minor / absent). As infants with major brain lesions were removed from analyses assessing brain-behaviour relationships, lesion severity was coded as minor vs absent in these analyses.
After multiple comparison correction (for 15 models, one for each SCN), SCN 12 was, again, the only SCN showing a significant multivariate effect on behavioural outcomes (p < .001), which was driven by its univariate effect on PC2 ("cognitive" component) (β = 5.18, p = .011) (see Table S3). Multivariate and univariate effects of cognitively stimulating parenting on behavioural outcomes are found in Table S4, showing that a multivariate effect of cognitively stimulating parenting (p = .002) was driven by its effect on PC1 ("preterm phenotype"), (β = -0.34, p < .001), in line with the main analysis.
In addition, we re-ran the analysis testing an effect of cognitively stimulating parenting on behavioural outcomes after excluding those participants with major lesions (i.e., the same subjects removed from all MRI analyses). In line with the main analysis, there was a significant multivariate effect of cognitively stimulating parenting (p = .008), which was driven by its effect 8 on PC1 ("preterm phenotype"), (β = -0.31 p = .003), suggesting that the effect is not biased by those individuals who suffered major lesions at birth.

Controlling for maternal age and education
Maternal education, rather than IMD, could be studied as a proxy for socio-economic status.
We therefore repeated our analysis controlling for maternal education instead of IMD. We also controlled for maternal age in this model. Maternal education was captured in terms of the age After multiple comparison correction (for 15 models, one for each SCN), SCN 12 was, again, the only SCN showing a significant multivariate effect on behavioural outcomes (p = .003), which was driven by its univariate effect on PC2 ("cognitive" component) (β = 6.04, p = .005) (see Table S3). Multivariate and univariate effects of cognitively stimulating parenting on behavioural outcomes are found in Table S4, showing that a multivariate effect of cognitively stimulating parenting (p = .013) was driven by its effect on PC1 ("preterm phenotype"), (β = -0.26, p < .010), in line with the main analysis.