Automated brain morphometric biomarkers from MRI at term predict motor development in very preterm infants.

Very preterm infants are at high risk for motor impairments. Early interventions can improve outcomes in this cohort, but they would be most effective if clinicians could accurately identify the highest-risk infants early. A number of biomarkers for motor development exist, but currently none are sufficiently accurate for early risk-stratification. We prospectively enrolled very preterm (gestational age ≤31 weeks) infants from four level-III NICUs. Structural brain MRI was performed at term-equivalent age. We used a established pipeline to automatically derive brain volumetrics and cortical morphometrics - cortical surface area, sulcal depth, gyrification index, and inner cortical curvature - from structural MRI. We related these objective measures to Bayley-III motor scores (overall, gross, and fine) at two-years corrected age. Lasso regression identified the three best predictive biomarkers for each motor scale from our initial feature set. In multivariable regression, we assessed the independent value of these brain biomarkers, over-and-above known predictors of motor development, to predict motor scores. 75 very preterm infants had high-quality T2-weighted MRI and completed Bayley-III motor testing. All three motor scores were positively associated with regional cortical surface area and subcortical volumes and negatively associated with cortical curvature throughout the majority of brain regions. In multivariable regression modeling, thalamic volume, curvature of the temporal lobe, and curvature of the insula were significant predictors of overall motor development on the Bayley-III, independent of known predictors. Objective brain morphometric biomarkers at term show promise in predicting motor development in very preterm infants.

Introduction

Globally, slightly more than 10% of pregnancies result in preterm birth (Chawanpaiboon et al., 2019). Very preterm infants (gestational age (GA) 31 weeks) are at greatly enhanced risk for cerebral palsy (CP) and also for less severe motor impairments, which, depending on GA, affect around 40% of preterm-born children (Williams et al., 2010, Moore et al., 1995). Such impairments are often not diagnosed until at least two years of age, which means that clinicians and therapists miss a critical window for intervention, during which neuroplasticity is maximal. Given the proven efficacy of early intervention to enhance motor outcome (Spittle et al., 2015), a number of early motor biomarkers have been proposed (Parikh, 2016); but many require independent validation.

Of the possible structural brain biomarkers available at term-equivalent age (TEA), abnormality or injury on early cranial ultrasound or structural magnetic resonance imaging (MRI) have been widely used (van’t Hooft et al., 2015, Hintz et al., 2015, Woodward et al., 2006); however neither is sufficient for accurate prediction. As such, researchers have turned to other structural brain features beyond overt injury. Brain volumes are consistently decreased in extremely-low-birth-weight or preterm cohorts (Grunewaldt et al., 2014, Valkama et al., 2001, Lind et al., 2010, Lind et al., 2011, Gui et al., 2019, Zhang et al., 2015), and there is evidence that these volumetric changes prognosticate motor ability (Gui et al., 2019, Spittle et al., 2010, Limperopoulos et al., 2014, Keunen et al., 2016, Setänen et al., 2016, Bolk et al., 2018, Moeskops et al., 2017, Loh et al., 2020). Preterm infants also exhibit altered cortical morphometrics, including inner cortical curvature (Makropoulos et al., 2016, Kline et al., 2019); surface area (Zhang et al., 2015, Makropoulos et al., 2016, Kline et al., 2019, Kapellou et al., 2006, Rathbone et al., 2011, Sølsnes et al., 2015, Sripada et al., 2018, Engelhardt et al., 2015), gyrification (Zhang et al., 2015, Kline et al., 2019, Engelhardt et al., 2015, Shimony et al., 2016, Slaughter et al., 2016), and sulcal depth (Zhang et al., 2015, Kline et al., 2019, Engelhardt et al., 2015, Shimony et al., 2016) at term. We previously identified a strong independent association of inner cortical curvature and cortical surface area with two-year cognitive and language outcomes in our very preterm cohort (Kline et al., 2020), suggesting that these biomarkers may aid in the prediction of other early neurodevelopmental outcomes. An association between cortical surface area and motor ability has been previously established in preterm children and adults (Grunewaldt et al., 2014, Sripada et al., 2015). However, to our knowledge, no prior study has assessed the individual value of inner cortical curvature, sulcal depth, or gyrification index to prognosticate motor ability in preterm infants. A classifier trained on combined cortical morphometrics and tissue volumetrics from structural MRI at TEA had good success predicting abnormal two-year motor function in a preterm cohort (Moeskops et al., 2017). However, the authors did not explore the individual predictive power of specific cortical morphometrics (e.g. curvature) or regional brain volumes (e.g. thalamic volume).

We aimed to independently validate and extend current knowledge about the predictive value of structural MRI biomarkers for motor outcome in very preterm infants. We examined the relationships between whole-brain and regional cortical morphometrics and volumetrics and overall, gross, and fine motor aptitude at two years of age. We then used machine learning methods to identify the top predictive biomarkers for each motor outcome. We assessed their combined predictive ability in multivariable regression, corrected for known confounders of motor development. Our population-based cohort, ability to derive objective morphometric and volumetric features using an automated pipeline, and higher MRI field-strength (3 T) than previous studies should allow us to pinpoint the most predictive morphometric MRI biomarkers of motor development available at TEA, which will lead to improved predictive models that facilitate earlier intervention in this high-risk cohort.

Materials and methods

Participants

Between December 2014 and April 2016, we prospectively enrolled a consecutive cohort of 110 very preterm infants (GA 31 weeks) from four regional, level-III neonatal intensive care units (NICUs) – Nationwide Children’s Hospital, Ohio State University Medical Center, Riverside Methodist Hospital, and Mount Carmel St. Ann’s Hospital. Babies were excluded if they had cyanotic heart disease or congenital or chromosomal anomalies of their central nervous system (e.g., Dandy-Walker malformation, Chiari II malformation, trisomy 21). Infants that remained hospitalized in the NICUs of Ohio State, Riverside, or St. Ann’s beyond 44 weeks postmenstrual age (PMA) were also excluded (N = 7). This study was approved by Nationwide Children’s Hospital Institutional Review Board. A parent or guardian of each infant gave written informed consent before study enrollment.

MRI data acquisition

We acquired all MRI data at Nationwide Children’s Hospital using a 3 T Siemens Skyra scanner and a 32-channel phased array head coil. Subjects were imaged at a mean (SD) PMA of 40.4 (0.6) weeks. A skilled neonatal nurse and neonatologist attended all scans. Each infant was fed 30 min prior to the scan and was fitted with silicone earplugs (Instaputty, E.A.R. Inc, Boulder, CO) and swaddled in a blanket and a vacuum immobilization device (MedVac, CFI Medical Solutions, Fenton, MI) to promote natural sleep. We used the following acquisition parameters: axial T2-weighted MRI: echo time 147 ms, repetition time 9500 ms, flip angle 150°, voxel dimensions 0.93x0.93x1.0 mm3, scan time 4:09 min; 3-dimensional magnetization-prepared rapid gradient-echo T1-weighted MRI: echo time 2.9 ms, repetition time 2270 ms, echo spacing time 8.5 ms, inversion time 1750 ms, flip angle 13°, voxel dimensions 1.0x1.0x1.0 mm3, scan time 3:32 min; axial susceptibility-weighted MRI : echo time 20 ms, repetition time 27 ms, flip angle 15°, voxel dimensions 0.7x0.7x1.6 mm3, scan time 3:11 min.

Motor outcomes

Our primary outcome of interest was the overall motor score from the Bayley Scales of Infant and Toddler Development–Third Edition (Bayley-III) (Bayley, 2006) standardized examination. Our secondary outcomes of interest were the gross and fine motor subscale scores. The Bayley-III overall motor score is a scaled composite of the gross motor and fine motor subtest scores, ranging from 40 to 160, with a mean (SD) of 100 (15) based on a normative population (Piñon, 2010). The entire cohort was scheduled for Bayley-III testing between 22- and 26-months corrected age. Four infants did not return during their appointment window and were instead tested between 33- and 36-months corrected age. These four subjects were nonetheless included in the analysis, as the Bayley-III motor scores obtained in this later time window would be equally accurate or more accurate than earlier scores.

MRI data processing

We used the automated, open-source pipeline from the developing Human Connectome Project (dHCP) (Makropoulos et al., 2018) to process all T2-weighted MRI scans. The dHCP pipeline performs tissue segmentation and volume estimation for each main tissue type and all cortical and subcortical structures. It also calculates the following metrics for each cortical region of the dHCP’s neonatal atlas: surface area, sulcal depth, gyrification index, inner cortical curvature (equivalent to curvature of the white matter surface), and cortical thickness, 510 metrics in all. The dHCP utilizes the 50-region Gousias neonatal atlas and further subdivides some structures, yielding 87 total regions (Gousias et al., 2012, Gousias et al., 2013). Briefly, the pipeline bias corrects and brain extracts the T2 images, segments tissue types and brain regions using the Draw EM algorithm of Makropolous et al. (Makropoulos et al., 2014), extracts white and pial surfaces, performs surface inflation and spherical projection via a reimplementation of a model used by Freesurfer (Fischl et al., 1999), and generates metrics of cortical geometry from the extraction and inflation process. We visually inspected each segmented scan and rejected any scans with tissue identification errors, improper segmentation, or artifacts. Very preterm infants with moderate or severe ventriculomegaly were also excluded, as large ventricles were often inaccurately identified as white matter by the pipeline.

Data analysis

Our main analysis identified the most promising motor biomarkers in four steps. Because the dHCP generates a very large number of regional and merged morphometrics and volumetrics, variable reduction in several stages was necessary to find the most highly-correlated predictors and reduce the likelihood of type 1 errors. Step 1) Application of a priori knowledge: as our initial cortical morphometric feature set, we chose four features that are known to be altered in prematurity (Kline et al., 2019); surface area, gyrification index, sulcal depth, and inner cortical curvature, for the whole brain and bilaterally for all four brain lobes. We also included curvature, sulcal depth, and surface area for the insula bilaterally. Previously, we had found cortical thickness from the dHCP to not be correlated with outcomes from our cohort (Kline et al., 2020); so we did not include thickness in the current analysis. As our initial set of volumetrics, we chose overall brain volume, bilateral volume of all four brain lobes and the insula, and volume of the main tissue types: cortical grey matter, white matter, deep grey matter (consisting of the caudate, thalamus, subthalamic nucleus, lentiform nucleus, and the internal capsule bilaterally), corpus callosum, CSF, and ventricles. We also included volumes for the main subcortical structures: the amygdala, hippocampus, caudate nucleus, lentiform nucleus, and subthalamic nucleus bilaterally and the corpus collosum, thalamus, brainstem and cerebellum unilaterally. Step 2) Correlation with outcome: From the 72 initial features described in step 1, we performed partial-correlation analysis, corrected for PMA at MRI scan, to identify variables that were significantly associated with each of the Bayley-III motor scales. We also accounted for false discovery rate in this step, using Benjamini Hochberg correction with a significance level of p < 0.05. Step 3) Iterative lasso regression for feature selection: For each motor scale, we used least absolute shrinkage and selection operator (Lasso) regression (Tibshirani, 1996, Tibshirani et al., 2012) to select the top features from the subset identified in step 2. Lasso regression is a form of linear regression that applies a penalty (alpha) to the model coefficients in the regression cost function. Increasing alpha drives more coefficients toward zero, encouraging sparse solutions. To find the optimal penalty that balanced error due to overfitting and error due to bias in our models, we first used cross-validated lasso model fitting with an 80–20 train-test split. After identifying the penalty that minimized the mean squared error, we fit the lasso model a final time for each motor scale, using the optimized alpha. From this model, we identified the top 3 predictors with the largest magnitude coefficients, denoting the largest contribution to the model. Step 4) Linear regression and bootstrapping: We fit a linear regression model for each motor scale using the top 3 coefficients (limited to 3 to prevent overfitting) from lasso regression and adjusting for known predictors and confounders. Given their established association with motor outcome, we retained male sex, gestational age at birth, PMA at MRI scan, birth hospital, and global brain abnormality score in all final models irrespective of their significance, to ensure the independent predictive power of our morphometric biomarkers. Briefly, global brain abnormality score sums abnormalities in the cortical grey matter, cerebral white matter, deep grey matter, and cerebellum, and higher scores indicate greater injury (Kidokoro et al., 2013, Tamm et al., 2020). We verified that the final models met all assumptions of linear regression, including no multicollinearity, normality, additivity, no autocorrelation, and homoscedasticity. We used bootstrap resampling with 10,000 repetitions to ensure that the confidence intervals for our final models were stable.

We also performed a secondary analysis to contrast the method outlined above with another popular data reduction method, principle component analysis (PCA). Starting from the 51 variables that were significantly correlated with Bayley-III motor scores (p < 0.05) in step 2 of the main analysis, we first normalized the data using python’s standard scaler. We then applied PCA to reduce the dimensionality of the dataset while preserving as much variance as possible. In PCA analysis, principal components are the eigenvectors of the data’s covariance matrix. The first principal component lies in the direction that maximizes the variance of the projected data in lower dimensional space. Each subsequent component captures as much of the remaining variance as possible, with the constraint that it must be orthogonal to all previous principal components (PCs). We continued to extract PCs until we had accounted for 95% of the variance in the intial morphometric data. We then used linear regression analysis to examine the predictive power of any PCs that accounted for at least 1% of the data variance for prognostication of motor scores. We contrasted these results with the results of our main analysis.

Data and code availability statement

Code used in this analysis and derived data that support the conclusions of this study are available upon direct request to the corresponding author.

Results

Subjects

Of the initial 110 eligible VPT infants, 17 were excluded due to structural brain injury or issues with their T2 MRI scan: ten exhibited ventriculomegaly, one exhibited posterior plagiocephaly (flat head syndrome); which interfered with dHCP volume estimation; five had significant artifacts or missing cortical boundaries on their T2 images, and one was missing a T2 image altogether. Of the remaining 93 infants, 75 (81%) returned for follow-up testing at a mean (SD) age of 24.7 (2.8) months corrected age and were included in the analysis. Seven infants (9.3%) had moderate injury on structural MRI based on their global brain abnormality scores (Makropoulos et al., 2018). The mean (SD) Bayley-III motor score for the final cohort was 94.7 (13.2), with a range from 46 to 122. Twelve infants (16%) scored more than 1 SD below the mean (<85). The mean (SD) raw fine and gross motor scores were 9.8 (2.5) and 8.3 (2.2), respectively. Additional baseline characteristics for the final cohort and for the infants lost to follow up are shown in Table 1. The two cohorts did not differ significantly on any baseline characteristics.

Table 1

Baseline characteristics of the final very preterm cohort with follow-up data included in the Bayley-III motor regression model (n = 75).

Characteristics	Final Cohort with Bayley Motor Follow-up (N = 75)	Lost to Follow-Up (N = 18)	p
Chorioamnionitis, n (%)	10 (13.3)	4 (22.2)	0.461
Antenatal steroids, complete course, n (%)*	40 (54.8%)	6 (33.3%)	0.121
Gestational age, weeks, mean (SD)	28.1 (2.4)	28.6 (3.0)	0.150
Birth weight, grams, mean (SD)	1105.7 (370.1)	1120.8 (385.1)	0.438
Male, n (%)	44 (58.7)	9 (50.0)	0.599
Severe retinopathy of prematurity, n (%)	8 (10.7)	1 (5.6)	1.000
Bronchopulmonary dysplasia, n (%)	38 (50.7)	9 (50.0)	1.000
Late-onset sepsis, n (%)	7 (9.3)	5 (27.8)	0.051
Postnatal steroids, n (%)	6 (8.0)	0 (0.0)	0.592
Low socioeconomic status, n (%)	41 (54.7)	14 (77.8)	0.109
Global brain abnormality score on structural MRI, mean (SD)	3.2 (2.7)	3.1 (3.5)	0.264

* Two infants in the final cohort did not have antenatal steroid information.

Cortical morphometrics

Total brain surface area was positively correlated with each Bayley motor scale, with PMA-corrected correlation coefficients of 0.52 for overall motor, 0.45 for gross motor, and 0.49 for fine motor. Total brain curvature was negatively correlated with each Bayley motor scale: −0.37 for overall motor, −0.34 for fine motor, and −0.33 for gross motor. Total brain sulcal depth was positively correlated with overall motor score (0.29) and gross motor score (0.28), but not with fine motor score (0.25). Total brain gyrification index was also significantly correlated with overall motor score (0.27) and gross motor score (0.29), but not with fine motor score (0.22).

Fig. 1 illustrates the strength of the correlation between each lobar cortical morphometric and all three scales, adjusted for PMA at MRI scan. Cortical surface area was positively correlated with every motor scale for every lobe of the preterm brain, with correlation coefficients ranging from 0.30 to 0.57. Cortical curvature was negatively correlated with each motor scale in several lobes of the brain, with correlation coefficients from −0.26 to −0.43. For overall motor score, correlations existed in the bilateral frontal and parietal lobes, the left temporal lobe, and the left insula. For fine motor score, correlations existed in the left frontal lobe, the left insula, the left temporal lobe, and the right parietal lobe. For gross motor score, correlations existed in the bilateral frontal and parietal lobes, the left insula, and the left temporal lobe. Gyrification index was correlated with overall motor score in the right occipital lobe and with gross motor score in the left parietal lobe and the right occipital lobe, with correlation coefficients from 0.26 to 0.29. Sulcal depth was correlated with overall motor score in the right parietal lobe and with gross motor score in the left frontal lobe, with correlation coefficients from 0.26 to 0.28.

Fig. 1

Whole-lobe partial correlations (PMA-adjusted) between cortical morphometrics and Bayley-III motor outcomes. Row A displays correlations between Bayley-III overall motor score and cortical surface area (column 1), inner cortical curvature (column 2), gyrification index (column 3), and sulcal depth (column 4). Rows B and C display the same correlations for the Bayley-III fine and gross motor subscale scores, respectively.

Volumetrics

Total brain volume was positively correlated with each motor scale: 0.43 for overall motor, 0.37 for gross motor, and 0.41 for fine motor. In terms of overall tissue volumes, white matter volume (0.37), cortical grey matter volume (0.42), and deep grey matter volume (0.52) were all positively correlated with overall motor score. White matter volume (0.33), cortical grey matter volume (0.35), and deep grey matter volume (0.50) were positively correlated with gross motor score. White matter volume (0.34) and deep grey matter volume (0.46) only were positively correlated with fine motor score. Volume of the ventricles (R = -0.10 for main motor score) and volume of overall CSF (R = 0.02 for main motor score) were not significantly correlated with any motor scale.

Fig. 2 illustrates the correlation strength between each lobar and subcortical volume and the Bayley-III motor scales, adjusted for PMA at MRI scan. All five bilateral lobar volumes were positively correlated with all motor scales, with the exception of the temporal lobe bilaterally for the gross motor scale. Lobar volume correlations ranged from 0.29 to 0.47 for overall motor, 0.33 to 0.41 for gross motor, and 0.28 to 0.45 for fine motor. Nearly all subcortical volumes tested were positively correlated with all motor scales. Overall motor score was positively correlated with volume: bilaterally in the amygdala, the caudate, the lentiform nucleus, and the subthalamic nucleus; in the left hippocampus; and in the brainstem, the cerebellum, the corpus collosum, and the thalamus, with correlations ranging from 0.29 to 0.55. Gross motor score was positively correlated with volumes: bilaterally in the caudate, the lentiform nucleus, and the subthalamic nucleus; in the right amygdala; and in the brainstem, cerebellum, and thalamus, with correlations ranging from 0.32 to 0.52. Fine motor score was positively correlated with volume: bilaterally in the amygdala, the lentiform nucleus, and the subthalamic nucleus; in the left caudate and the left hippocampus; and in the brainstem, cerebellum, corpus callosum, and thalamus, with correlations ranging from 0.27 to 0.50. Thalamic volume was the most highly-correlated subcortical volume for all motor scales.

Fig. 2

Partial correlations (PMA-adjusted) between brain volumes and Bayley-III motor outcomes. From left to right, row A displays correlations between Bayley-III overall motor score and 1) all brain volumes shown in 2-D from sagittal, coronal, and axial perspectives, 2) all brain volumes shown in 3-D from right and left perspectives, and 3) subcortical structures only shown in 3-D from right and left perspectives. Rows B and C display the same correlations for the Bayley-III fine motor subscale score and the Bayley-III gross motor subscale score, respectively.

Regression results

Our automated feature-selection method using lasso regression with an optimized penalty selected the following best three predictors for the overall motor scale: thalamic volume, curvature of the left temporal lobe, and curvature of the left insula. For the fine motor scale, the top three predictors were: surface area of the right occipital lobe, curvature of the left temporal lobe, and curvature of the left insula. For the gross motor scale, the top three predictors were: thalamic volume, curvature of the left temporal lobe, and gyrification index of the parietal lobe. Table 2 displays the final regression models for all three Bayley-III motor scales, with and without known predictors added. All cortical morphometrics and volumetrics remained significantly correlated with motor outcomes when known predictors were introduced into the models, illustrating the independent predictive power of these biomarkers. For each motor scale, the combination of just three automatically-extracted and automatically-selected cortical morphometric and volumetric features explained over one-third of the variance in motor scores. In the case of overall motor score, the three features explained nearly one-half of the variance. Global brain abnormality score on structural MRI and PMA at MRI scan were both negatively correlated with the overall and gross motor scales. The negative correlation with PMA at MRI scan is likely due to the sicker babies in our cohort being scanned at slightly later time points, when they were healthy enough to safely undergo MRI. Male sex was negatively correlated with gross motor score. The bootstrap bias-corrected confidence intervals for all three final models were comparable (Table 2), which supports their internal validity.

Table 2

Morphometric biomarkers and prediction of 1) Bayley-III motor scores with and without adjustment for known predictors and confounders (global injury score on structural MRI, gestational age, and male sex).

Predictors	Models with morphometrics only		Models with morphometrics and covariates
	Coefficients (95% CI)	P	Coefficients (95% CI)	Bootstrap 95% CI	P
Bayley-III overall motor (n = 75)	Adjusted R2: 0.47		Adjusted R2: 0.51
Thalamic volume (cm3)	5.97 (3.67, 8.28)	<0.001	5.53 (2.54, 8.52)	(2.58, 8.48)	<0.001
Curvature of the left temporal lobe	−17.34 (-24.63, −10.05)	<0.001	−17.00 (-24.32, −9.67)	(-24.75, −9.24)	<0.001
Curvature of the left insula	−8.32 (-14.12, −2.52)	0.006	−6.61 (-12.30, −0.91)	(-12.29, −0.92)	0.024
Global abnormality score			−1.24 (-2.27, −0.20)	(-2.43, −0.04)	0.021
Gestational age			−0.32 (-1.44, 0.80)	(-1.52, 0.88)	0.569
Male sex			−0.72 (-5.26, 3.83)	(-5.08, 3.65)	0.755
Postmenstrual age at MRI scan			−4.46 (-8.55, −0.38)	(-8.47, −0.46)	0.032
Hospital of birth			0.26 (-1.93, 2.46)	(-1.87, 2.40)	0.810
	Coefficients (95% CI)	P	Coefficients (95% CI)	Bootstrap 95% CI	P
Bayley-III fine motor (n = 75)	Adjusted R2: 0.36		Adjusted R2: 0.41
Surface area of the right occipital lobe (cm2)	0.11 (0.05, 0.16)	<0.001	0.08 (0.02, 0.14)	(0.02, 0.14)	0.008
Curvature of the left temporal lobe	−3.01 (-4.56, −1.47)	<0.001	−2.98 (-4.52, −1.44)	(-4.70, −1.26)	<0.001
Curvature of the left insula	−1.78 (-2.98, −0.57)	0.004	−1.40 (-2.60, −0.21)	(-2.65, −0.15)	0.022
Global abnormality score			−0.16 (-0.38, 0.06)	(-0.39, 0.07)	0.146
Gestational age			0.07 (-0.16, 0.29)	(-0.15, 0.28)	0.563
Male sex			0.33 (-0.60, 1.26)	(-0.61, 1.27)	0.481
Postmenstrual age at MRI scan			−0.70 (-1.52, 0.12)	(-1.51, 0.10)	0.091
Hospital of birth			0.43 (-0.03, 0.89)	(-0.02, 0.87)	0.067
	Coefficients (95% CI)	P	Coefficients (95% CI)	Bootstrap 95% CI	P
Bayley-III gross motor (n = 75)	Adjusted R2: 0.42		Adjusted R2: 0.52
Thalamic volume (cm3)	1.11 (0.71, 1.51)	<0.001	1.15 (0.64, 1.65)	(0.59, 1.70)	<0.001
Curvature of the temporal lobe (L)	−2.32 (-3.62, −1.02)	0.001	−2.04 (-3.27, −0.80)	(-3.20, −0.87)	0.002
GI of the parietal lobe (L)	1.95 (0.48, 3.43)	0.010	2.65 (1.23, 4.07)	(1.08, 4.23)	<0.001
Global abnormality score			−0.24 (-0.41, −0.06)	(-0.42, −0.06)	0.008
Gestational age			−0.06 (-0.25, 0.13)	(-0.29, 0.17)	0.539
Male sex			−0.82 (-1.61, −0.03)	(-1.59, −0.05)	0.043
Postmenstrual age at MRI scan			−0.75 (-1.45, −0.05)	(-1.52, 0.01)	0.035
Hospital of birth			−0.23 (-0.60, 0.14)	(-0.61, 0.15)	0.223

*Curvature is a unitless metric that ranges from + 1 to −1.

PCA analysis

19 principle components captured greater than 95% of the variance of our morphometric data. PC1 explained 56% of the variance in the total dataset, PC2 explained 9%, and PC3 explained 5%. Only the first 14 PCs accounted for 1% or greater of the dataset variance. Of these, only PCs 1 and 2 were significantly correlated with Bayley-III motor score in univariate (PMA-corrected) linear regression. In multivariate linear regression corrected for PMA, the combination of these 2 significant PCs explained 39% of the variance in two year motor scores (Table 3). In order of significance, the top variable contributions to PC1 were: overall brain volume, volume of the cortical gray matter, volume of the right parietal lobe, overall surface area, and volume of the deep gray matter. The magnitude of the contribution of all 5 top variables to PC1 were nearly identical and were all positive. The rest of the volume and surface area metrics all contributed strongly to this PC, including the volume of the thalamus. For PC2, the top variable contributions came from overall curvature, curvature of the left frontal lobe, curvature of the left parietal lobe, curvature of the right parietal lobe, and curvature of the right frontal lobe. All curvature variable contributions to this PC were positive.

Table 3

Principal components of the morphometric dataset and prediction of Bayley-III motor scores with and without PMA adjustment.

	Coefficients (95% CI)	P	Coefficients (95% CI)	P
Bayley-III gross motor (n = 75)	Adjusted R2: 0.37		Adjusted R2: 0.39
PC1	1.19 (0.74, 1.65)	<0.001	1.31 (0.85, 1.78)	<0.001
PC2	−2.44 (-3.60, −1.28)	<0.001	−2.53 (-3.68, −1.39)	<0.001
PMA			−3.94 (-8.39, 0.50)	0.08

Discussion

In this study, we illustrate that motor ability at two-years corrected age in a preterm cohort can be predicted with brain morphometrics and volumetrics automatically extracted from structural MRI at term. This is especially true for thalamic volume, cortical surface area, and cortical curvature. Greater global and lobar cortical surface area, greater overall volume of the whole brain and of most tissue types, and greater subcortical volumes were positively correlated with all scales of two-year motor ability. Greater cortical curvature was negatively correlated with all scales, both globally and in the majority of lobes. Greater gyrification and sulcal depth were positively correlated with overall motor score and gross motor score, globally and in one to two brain lobes each. In multivariable regression, just three automatically-extracted and automatically-selected cortical morphometric and/or volumetric biomarkers accounted for between a third and a half of the variance in fine, gross, and composite scales of motor development at 2-years corrected age.

Inner cortical curvature is a novel prognostic biomarker for motor ability. Makropoulos et al. previously reported increased inner cortical curvature throughout the preterm brain compared to term controls (Makropoulos et al., 2016), and we corroborated this finding in our very preterm cohort (Kline et al., 2019). Curvature of the left insula and curvature of the left temporal lobe were both retained in overall motor and fine motor models. Curvature is likely such a valuable metric for prediction because it generalizes some of the major brain maturational abnormalities associated with prematurity. Premature infants consistently exhibit shallower sulci and decreased gyrification compared to their full-term peers. Shallower sulci are missing sulcal wall area (low curvature), but have the same amount of sulcal trough area (high curvature), which results in higher average curvature (Shimony et al., 2016), although the individual anatomy of sulci and gyri can vary considerably. The insula is involved in some aspects of motor and somatosensory control (Uddin et al., 2017), while regions of the temporal lobe help to mediate visuomotor coordination (Tankus and Fried, 2012). Taken together, the curvature of the insula and the curvature of the temporal lobe explained 28% of the variance in the overall motor model and 22% of the variance in the gross motor model.

The association between brain volumes and motor ability in preterm cohorts has been well-established. Some studies have identified this relationship at the same developmental timepoint, rather than using volumes to prognosticate later motor outcomes (Limperopoulos et al., 2014, Loh et al., 2017). At seven years, basal ganglia and thalamic volumes were associated with motor function (Loh et al., 2017); and at three years, volumes of the left and right sensorimotor cortices were associated with gross motor score (Limperopoulos et al., 2014). Still other prognostic studies have examined associations between only a few select brain volumes at TEA and motor outcome, for instance reduced cerebellar volume and abnormal general movements at three months (Spittle et al., 2010) or unmyelinated white matter volume and two-year motor outcomes (Keunen et al., 2016). Gui et al. showed that a classifier trained using TEA brain volumes was better able to predict psychomotor developmental index score at 18–24 months (Gui et al., 2019) than a classifier trained on social and demographic variables alone, but they did not examine which brain volumes were the most predictive. Specific brain volumes from TEA have been shown to prognosticate motor ability into late childhood. Volume of the cerebellum and brainstem predicted fine motor skills at age six (Bolk et al., 2018), and volume of the frontal lobes, basal ganglia, thalamus, cerebellum, and brainstem were correlated with motor scores at age 11 (Setänen et al., 2016). Our results corroborate and extend prior studies by showing positive correlation between numerous subcortical volumes at TEA — the volume of the amygdala, brainstem, caudate, cerebellum, hippocampus, lentiform nuclei, subthalami nuclei, and thalamus—and motor outcomes at two-years corrected age. Like Loh et al. (Loh et al., 2017), we identified the thalamus as a top predictive volumetric biomarkers for motor development. Thalamic volume was a top-three predictive biomarkers for overall motor and gross motor ability. In bivariable regression, thalamic volume alone accounted for 27% of the variance in overall motor score and 26% of the variance in gross motor score. The thalamus is a crucial motor hub that relays signals from the cerebellum and basal ganglia to the cerebral cortex, which may explain why it is such a prominent motor biomarker.

We found that surface area of the very preterm brain was positively correlated with all Bayley-III motor scales globally and bilaterally in all five lobes, making it the most universally-correlated metric tested. Comparing our results with those in the literature, the robust relationship we identified between cortical surface area and motor ability appears to persist into childhood and early adulthood. Grunewalt et al. identified reduced cortical surface area and poorer motor skills in a cohort of extremely-low-birth-weight children at age ten (Grunewaldt et al., 2014). In preterm-born young adults (18–22 years of age), Sripada et al. demonstrated a positive relationship between measures of visual motor integration and cortical surface area in all lobes, with associations in the temporal lobe, the insula, and the occipital lobe being the most consistently significant (Sripada et al., 2015). In the present study, the surface area of the right occipital lobe explained 18% of the variance in the fine motor subscale, presumably reflecting the high level of visual input and feedback needed to perform a fine motor task. The fact that surface area was not retained in the overall or gross motor models was likely due to the high degree of collinearity between thalamic volume and surface area. In our final cohort, overall surface area was significantly correlated with thalamic volume (r = 0.79, p < 0.0001). Conversely, overall curvature was not (r = -0.10, p = 0.41), potentially explaining why curvature biomarkers were retained in the final models with thalamic volume.

The results of our secondary PCA analysis bolster our main analysis, by reinforcing that surface area and volume biomarkers are of primary importance but that they capture similar (collinear) features in this morphometric dataset. The 39 significant volume and surface area metrics from step 2 of the main analysis all contributed strongly to PC1, which in turn explained 56% of the morphometric dataset variance. Furthermore, the PC results corroborate the secondary but still important contribution of curvature to motor prediction, with curvature metrics being the primary contributors to PC2, which explained 9% of the morphometric dataset variance. However, PCA inherently sacrifices some accuracy in order to simplify and streamline the data. This is likely why the top two principal components explained less of the variance in two year motor scores than our carefully selected individual biomarkers. The relative ease of interpretability of individual regional biomarkers compared to PCA-transformed data also supports the validity of our main approach for clinical model building.

Our study was limited by a modest sample size and 19% loss to follow-up for our primary analysis. However, the infants lost to follow-up did not differ from included infants on any baseline characteristics, which should minimize attrition bias. Given that only seven subjects in our final cohort developed cerebral palsy, we lacked sufficient study power to examine prediction of this important outcome.

Conclusions

In conclusion, of the numerous biomarkers available on structural MRI at TEA, thalamic volume and/or whole lobar surface area combined with whole lobar measures of curvature appear to be the best prognostic biomarkers for early motor development in very preterm infants. Gyrification index of the left parietal lobe also improves prediction of gross motor ability. Using just three automatically extracted metrics, we were able to account for nearly one-half of the variance in overall Bayley-III motor score at two-years corrected age. Overall, cortical morphometrics and volumetrics automatically-derived from structural MRI at term show promise to prognosticate later motor ability in very preterm infants. Larger cohort studies are needed to validate these biomarkers and to develop additional advanced MRI biomarkers to enhance the models’ predictive power and facilitate aggressive early intervention for the most at-risk infants.

CRediT authorship contribution statement

Julia E. Kline: Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing - original draft, Writing - review & editing, Resources, Software. Venkata Sita Priyanka Illapani: Resources, Validation, Writing - review & editing. Lili He: Data curation, Writing - review & editing. Nehal A. Parikh: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.