Normative spatiotemporal fetal brain maturation with satisfactory development at 2 years – Nature

Participants and fetal US scans

Three-dimensional US scans of the fetal head were acquired as part of the FGLS of the INTERGROWTH-21st Project²⁹ to explore their use for undertaking 2D measurements, and for further exploratory analyses^30,31. A detailed description of FGLS and its inclusion criteria has been published elsewhere².

In brief, INTERGROWTH-21st was a population-based project, conducted between 2009 and 2016 in eight delimited urban areas: Pelotas (Brazil), Turin (Italy), Muscat (Oman), Oxford (UK), Seattle (USA), Shunyi County in Beijing (China), the central area of Nagpur (India) and the Parklands suburb of Nairobi (Kenya). Participating women, who initiated antenatal care before 14 weeks’ gestation, were selected on the basis of WHO criteria for optimal health, nutrition, education and socioeconomic status needed to construct international growth standards³. Hence, they had low-risk pregnancies that fulfilled well-defined and strict inclusion criteria at both population and individual levels².

The last menstrual period (LMP) was used to calculate gestational age (±7 days) provided that: (1) the date was certain; (2) the woman had a usual 24–32 day menstrual cycle; (3) she had not been using hormonal contraception or breastfeeding in the preceding 2 months and (4) any discrepancy between the gestational ages on the basis of LMP and crown–rump length, measured by US at 9⁺⁰ to 13⁺⁶ weeks(^+days) from the LMP was fewer than or equal to 7 days, using the formula described by Robinson and Fleming⁴⁹. To ensure that crown–rump length measures were interpreted consistently, the Robinson–Fleming formula was loaded into all the study US machines; whenever another machine had to be used locally for crown–rump length measurement, a conversion table extracted from the same formula was provided. The crown–rump length technique was also standardized across sites and all ultrasonographers were trained uniformly.

US scans were then performed every 5 ± 1 weeks from 14⁺⁰ weeks’ gestation to delivery (that is, 14–18, 19–23, 24–28, 29–33, 34–38 and 39–42 weeks’ gestation). Dedicated sonographers performed the US scans using identical, commercially available, equipment (Philips HD-9, Philips Ultrasound), with a curvilinear abdominal 3D transducer (V7-3). A transvaginal probe was not used because it would have been culturally inappropriate in many settings. The US probe was positioned such that the central axial view was collected at the level of the thalami, and the angle of insonation was adjusted to include the entire skull (roughly 70°) for a typical volume acquisition time of 4 s. We conducted centralized hands-on training of sonographers, and the Oxford-based Ultrasound Quality Control Unit regularly carried out site-specific standardization procedures to ensure proper use of the US equipment and protocol adherence.

Extended Data Fig. 1a presents a flowchart of the data inclusion pipeline used to generate the 4D (3D + time) probabilistic atlas based on 899 healthy singleton fetuses in the FGLS database that were appropriately grown and born at term (451 (50.2%) females and 448 (49.8%) males). The maternal and perinatal characteristics of this subpopulation were similar to the total FGLS population (n = 4,321) (Supplementary Tables 3–5 and Extended Data Fig. 8).

Infant follow-up

Across all sites, standardized clinical care and feeding practices were implemented using the INTERGROWTH-21st Neonatal Group protocols (www.intergrowth21.org.uk). Exclusive breastfeeding up to 6 months and appropriate nutritional support for infants born preterm were promoted during and after pregnancy. Detailed information was obtained from the mother at the age of 1 and 2 years about the infant’s health, severe morbidities, hospitalizations, length of breastfeeding, timing of the introduction of solid food, age at weaning, feeding practices and food intake, using specially produced forms (www.intergrowth21.org.uk). The proportion of infants receiving breast milk, and vitamin and mineral supplements, and those following a special diet were estimated at the ages of 1 and 2. Similarly, at age 1 and 2, the infant’s weight, length and head circumference were measured following WHO protocols⁵⁰, and their age- and sex-specific z-scores and centiles were compared to the WHO Child Growth Standards³⁵. These anthropometric measures, as indicators of general nutrition at the age of 2, are strongly predictive of later attained height, development and human capital⁵¹.

Neurodevelopment assessment

We assessed neurodevelopment at 2 years of age using the INTER-NDA (www.inter-nda.com), an international, psychometrically valid and reliable, standardized tool, targeted at children 22–30 months of age, which measures several dimensions of early development using a combination of directly administered, concurrently observed and caregiver reported items⁵². It was designed to be implemented by non-specialists across multinational settings, and includes a reduced number of culture-specific items comprising six domains measuring cognition, language, fine and gross motor skills, and positive and negative behaviour, in an assessment time of 15 min on average.

The INTER-NDA has been validated against the Bayley Scales of Infant Development III edition⁵³. On the basis of established guidelines⁵⁴, it showed ‘good’ agreement with interclass coefficient correlations across domains ranging between 0.75 and 0.88. Its norms are the first international standards of early child development, constructed according to the prescriptive WHO approach using data from five of the eight INTERGROWTH-21st study sites⁴⁰. So far, more than 18,000 children in 22 countries have been assessed using the INTER-NDA.

Attentional problems and emotional reactivity were measured on the respective subscales of the Preschool Child Behavior Checklist⁵⁵; responses were based on caregiver reports. Vision was assessed using the Cardiff Visual Acuity and Contrast Sensitivity tests for binocular vision⁵⁶. These are indicative of the integrity of the central visual pathway, and as directly observed neurodevelopmental markers, are unlikely to be affected by cultural influences and co-occurring disturbances in cognitive, hearing and language skills.

Motor development was assessed against four WHO milestones that are less likely to be affected by recall bias: sitting without support, hands knees crawling, standing alone and walking alone⁵⁷. Trained staff collected the data on a form with pictures of the relevant child positions and corresponding definitions. Parents were asked to report the age in months and weeks when they ‘first observed’ or ‘never observed’ the milestones. We assessed the age (in months) at which WHO gross motor milestones were first achieved.

All INTER-NDA assessors were trained and standardized centrally. All assessors were subject to a protocol adherence and reliability assessment following training; only those with protocol adherence scores in excess of 90% and inter-rater reliability of more than 0.8 conducted assessments. The administration of the above tests was supported by a tablet-based data collection and management system. Field staff were unaware of the INTER-NDA domain and total scores for individual children and sites. Data were uploaded onto secure servers as soon as each assessment was completed and compared to the international normative values established by the INTERGROWTH-21st Project⁴⁰.

Image curation

Volumetric images were selected from the INTERGROWTH-21st database on the basis of image quality criteria alone and AILN was blinded to the study site information during image selection. Briefly, inclusion criteria required the fetal head to occupy at least 30% of an artefact-free image (defined as free from the aberrations introduced by spontaneous fetal movement during scanning, acoustic shadows caused by the skull’s convex shape or reverberations from structures in the proximal hemisphere or maternal tissues), and clear visibility of structures in the distal hemisphere, for example, thalami, CoP, Sylvian fissure and ChP (details in Supplementary Table 6).

To detect artefact-free cases to achieve the best quality atlas, AILN manually assessed all 48,813 3D volumes; each volume took 1–2 min to assess depending on the overall quality, that is, a total of roughly 1,200 h. Artefacts included fetal motion, strong acoustic shadows and poor contrast, which limited visualization of brain structures in the distal hemisphere. The final selection of 899 cases depended on image quality, which was further assessed on the basis of the criteria provided in Supplementary Table 6 after the initial screening. Each assessment took on average 3 min (rough range 1–6 min), that is a total of roughly 210 h. We did not have access to a tool to automate the process.

We ultimately excluded more than 55,000 scans obtained between 14⁺⁰ and 30⁺⁶ weeks’ gestation on the basis of image quality (Supplementary Table 6 and Extended Data Fig. 1a). Extended Data Fig. 1c shows the distribution of US scans across these gestational ages for each study site.

As a result of acoustic shadows and reverberation artefacts, only the distal cerebral hemisphere typically contains clearly discernible brain structures in US images. All analyses were, therefore, performed separately for images presenting the left and right hemispheres in the distal region of the US scan, resulting in two atlas (template) maps per gestational week. Each hemispheric map was constructed from a minimum of ten US scans per gestational week (Supplementary Table 7).

Image preprocessing

Each individual US image was processed following a series of manual and automated procedures. As summarized in Fig. 1a, these mainly included: (1) rigid alignment to a standardized coordinate space, (2) brain extraction and (3) structural enhancement. Before entering this pipeline, the 3D volumes were resampled to an isotropic voxel size of 0.6 × 0.6 × 0.6 mm³ using trilinear interpolation, from a median size of 0.32 × 0.51 × 0.85 mm³. To ensure that all fetal brains were included whole, the images were cropped to 160 × 160 × 160 voxels around the centre of the brain.

Rigid alignment

Fetal brains show inter-participant variability in shape, size and localization. In US images, the position of the US probe relative to the head directly affects the orientation of the imaged brain. Establishing the structural coordinate system is crucial as it reduces the degrees of freedom in the non-rigid transformations required for atlas construction. To compare anatomy across participants, we performed semi-automated brain alignment in two steps to bring the brains into structural correspondence. First, a deep learning-based, convolutional neural network was used to localize the brain within the 3D volume, exclude maternal and extracranial tissues, and linearly align each brain to a common 3D coordinate space⁸. A secondary manual correction step ensured inter-participant co-alignment across all images by following the convention of the stereotaxic space of the MNI-Colin27 template⁵⁸. The cropped brain images were rigidly aligned to a standardized 3D coordinate space, using a seven-parameter linear transformation (three translation terms, three rotation terms, one scale term to preserve the aspect ratio). This was achieved using MATLAB’s graphical user interface toolkit (MathWorks, 2021), which involved locating the midsagittal plane (that coincides with the longitudinal fissure separating the two cerebral hemispheres), and aligning key neuroanatomical landmarks therein (procedure shown in Supplementary Video 2). Once aligned, it was easy to determine the left and right hemispheres. We found that manually locating three landmark points on the corpus callosum and cavum septum pellucidum complex and rigidly registering these to a fixed point-based template achieved sufficient structural alignment across all images within a specific gestational age (Supplementary Video 2).

Fetal brain masking

Extracranial tissues (for example, eye orbits, scalp, maternal tissues) were removed using brain masks from the CRL MRI brain atlas constructed from 81 fetuses scanned between 19 and 39 weeks’ gestation⁹. We manually registered the CRL atlas template of the corresponding gestational week to each aligned US image, and propagated the brain extraction mask to each image. The CRL brain extraction mask included cortical grey matter, white matter, subcortical grey matter structures, cerebrospinal fluid, lateral ventricles and cerebellum. Care was taken to ensure that the mask outline was aligned to coincide with the skull’s inner boundary. It is worth noting that the CRL atlas contains templates from 21 to 37 weeks’ gestation, so template-matching was only possible within the gestational age range that overlapped with our atlas (21 to 31 weeks). For fetuses with scans collected between 14 and 20 weeks’ gestation, the brains were masked using an isotropically scaled version of the earliest CRL atlas brain mask (at 21 weeks’ gestation).

Image enhancement and/or hemisphere selection

All images were intensity normalized using histogram matching to a preselected, age-matched reference volume using the imhistmatchn function implemented in MATLAB. Speckle reduction and ridge enhancement (sulci and soft tissue boundaries) were achieved by filtering the images using the contrast-invariant monogenic signal constructed with a multiscale log-Gabor filter. This produced an edge map for each image, which provided a secondary image channel containing enhanced structural information for the atlas construction step (Extended Data Fig. 5b). The processed 3D volumes were classified as clearly capturing the left or right cerebral hemisphere in their distal portion. We reliably ascertained whether the right or left cerebral hemisphere was captured in the image by observing the presentation of the fetal head in the aligned brain volumes. As such, the proximal hemisphere was excluded by cropping the image 10 voxels to the side of the longitudinal fissure, thus removing most of the proximal hemisphere. This step was applied to the images and corresponding edge maps.

Atlas construction

We created a digital 4D spatiotemporal atlas to characterize fetal brain maturation, encoding the structural variability expected of a healthy population at different gestational ages. We used a diffeomorphic Demons-based approach^59,60 to estimate the mapping that brought each brain image’s anatomy into spatial correspondence with all other brain images in the population set. As a consequence of the rapid structural evolution of the brain during intrauterine life, we opted for multi-channel groupwise registration^61,62: a template-free approach to reduce bias introduced by selecting a fixed reference brain image, and the ability to supplement the US intensity images with (several) extra noise-suppressed edge images. That is, no initial template was selected as input to the non-rigid registration; instead, the atlas was completely derived from the set of population images following an iterative optimization procedure. For each gestational age, the algorithm simultaneously generated a deformation map for each population image, and an estimated representation of the average brain template, improving the structural alignment (at the pixel level) with each iteration. The two hemispheres were brought together in the final templates without any anatomically overlapping regions (Fig. 1b and Extended Data Fig. 5c). The ten voxels of the proximal hemisphere that had been retained in the image preparation step were removed to enable the hemispheric templates to be joined along the longitudinal fissure.

Specifically, for each week and each cerebral hemisphere, the brain template was constructed so that each individual brain was minimally deformed, while the structures were maximally aligned. The algorithm took as input a set \({M}_{{\rm{i}}}^{{\rm{a}}}=\,{\{({{\rm{I}}}_{{\rm{i}}},{{\rm{H}}}_{{\rm{i}}})\,\}}_{\{\,{\rm{i}}=1\}}^{{\rm{N}}}\) comprising N pairs of input images and corresponding edge maps, where a is the gestational week at which each scan was acquired. Each input image I_i is a 3D B-mode US scan defined in Euclidean space, and its corresponding structural edge map H_i, was constructed using multiscale feature asymmetry with a log-Gabor filter with kernels of scales \(\lambda =\{0.075,0.125,0.175\}\) to enhance fissures while suppressing speckle. The objective was to find a set of non-rigid transformations \({D}_{{\rm{i}}}^{{\rm{a}}}\), each of which mapped the individual images I_i to a group-representative average image \({\hat{I}}^{{\rm{a}}}\) by simultaneously minimizing the intensity distance across the images at each voxel location (Extended Data Figs. 5 and 10), with additional regularization to enforce a diffeomorphic mapping between the atlas and each individual image⁵⁹ (Extended Data Fig. 3c).

To obtain the age group’s average image \({\hat{I}}^{{\rm{a}}}\), the deformation maps (D_i) were applied to their corresponding images, and the images were combined using a voxel-wise mean (Extended Data Fig. 10). Knowing that the fetal brain evolves over weeks, and that some tissues are transient and may not be consistently observed across the entire gestational age range, we opted to generate a separate atlas for each gestational week. This resulted in 17 brain atlas template maps, spanning 14 to 31 weeks’ gestation (Fig. 1b and Supplementary Video 1).

Structural atlas labels

Age-specific label maps were generated by manual segmentation of the volumetric atlas images. These were conducted by four authors (A.I.L.N. and F.A.M. for TBV; L.S.H. for subcortical structures; and M.K.W. for CoP), in consultation with histology-based atlases of human fetal brains^63,64, and were independently verified by two co-authors (A.I.L.N. and W.S.). The atlases were first manually labelled in the axial plane, and iteratively modified in the coronal and sagittal planes by L.S.H. and M.K.W. The following structures were labelled: CoP, CB, ChP and intracranial space (TBV).

The atlas labels were used to train convolutional neural network models to segment the 3D brain mask⁶, subcortical⁷ and cortical²⁷ structures automatically. The trained models were then applied to an out-of-sample data set of 3D brain volumes to extract volumetric measures for evaluation of the reliability of the normative patterns of fetal brain growth. Intra-rater variability (between 85 and 91% agreement) was computed for the manual segmentation work related to structural atlas labelling⁷.

The CBV, CoPV and CoPA values reported only cover the period between 18 and 27 weeks’ gestation. Before 18 weeks’ gestation no MRI reference exists and the voxel spacing (0.6 mm), made it difficult to delineate the CB and CoP. After 27 weeks’ gestation, the increased ossification of the skull (for example, the petrous part of the temporal bone for visualizing the CB) reduced image quality.

Hemispheric asymmetry

We characterized the emergence of local developmental asymmetries between the cerebral hemispheres. For each gestational age, the affinely registered brain images consisting of a clearly visualized right hemisphere were flipped across the midsagittal plane (longitudinal fissure) creating mirror images that matched the left hemisphere maps spatially \(({I}_{{\rm{R}}}\to {I}_{{\rm{R}}}^{{\prime} }\approx {I}_{{\rm{L}}})\) (Extended Data Fig. 5d). The analyses were conducted separately for each gestational week to show the timing and regions showing cerebral lateralization.

Left and right data distributions

The final selection sample of 899 cases showed a bias towards left visible hemispheres (Extended Data Fig. 1e), which is in keeping with US data, although mostly near term, showing that about two thirds of fetuses are in the left occiput position in utero^65,66,67. To verify the effect sizes between scan data sampled from the two hemispheres (for each gestational age), we applied a two-sample Kolmogorov–Smirnov test. Rejection of the null hypothesis confirmed sufficient similarity between the gestational age and brain volumes from each hemisphere. The null hypothesis was rejected for all gestational ages, suggesting that the investigation of asymmetry was valid.

Anthropometric asymmetry

The distribution of all volumetric measures for each hemisphere was assessed for normality, conditional on gestational age. Gaussian additive models were fitted to each of the five volumetric measures, separately for each cerebral hemisphere. We then tested whether to reject the null hypothesis of equality (that is, absence of asymmetry) between the left and right hemispheres by computing the Cohen’s d estimates for each brain region with the statsmodels Python package (v.0.13.2). The null hypothesis at a nominal 5% level of significance was tested. Significance was found only for ChPV, which produced an F statistic of 27.424 on 4 and 847 degrees of freedom, with P < 0.001.

Spatial asymmetry

To identify spatial patterns of cerebral asymmetry, we performed tensor-based morphometry, which can reveal the local volumetric change (expansion or contraction) between a target and source image. Tensor-based spatial statistics of asymmetry were calculated by applying the logarithm to the diffeomorphic Jacobian determinant maps that resulted from the non-rigid registration step (Atlas construction above). Application of the logarithm encourages the distribution of the deformation fields to be zero-mean and symmetric, which enables ease in interpreting relative tissue growth and/or loss⁶⁸.

Voxel-wise permutation tests were performed on all scans at each gestational age to show which hemisphere had the greater US signal across all scans collected at a given gestational age. We used non-parametric ‘Monte Carlo’ permutation testing as implemented in the FSL RANDOMISE method⁶⁹, and applied threshold-free cluster enhancement to the statistical maps⁷⁰, to enhance the brain areas that showed spatial contiguity. This approach is appropriate when the null distribution is not known a priori, and has been shown to handle noise and spatially correlated signals⁷⁰. We tested for main group effects (left hemisphere, +1, right hemisphere, −1), while including residualized gestational age as a covariate in the general linear model. Five thousand permutations were performed for each contrast, and the regional clusters surviving a conservative family-wise error rate correction threshold of P < 0.05 (two-tailed and permutation-based) were deemed as sites of significant asymmetry at the given age. All permutation testing was conducted within the mask of the left cerebral hemisphere (Extended Data Fig. 5d). This tensor-based morphometry approach normalizes for differences in brain volume, and so any detected regions of statistical significance show local morphological, rather than size, differences. In the generated statistical maps (Fig. 5a), positive values shown in the left or right hemisphere indicated either a leftwards (L > R) or rightwards asymmetry (R > L).

We performed secondary analysis of the cerebral subregions to explore the spatial and temporal patterns of asymmetry. Age-matched structural parcellation templates from the CRL atlas were rigidly aligned to our US-based atlas using a shape-preserving similarity transform (scaling, translation, rotation). For each gestational age group, we performed cluster-level analysis to identify the cortical regions showing significant asymmetry. Cluster tables summarize the percentage of cluster-specific voxels contained within each region in the parcellation map (Supplementary Table 8). A surface-based representation of this result was achieved by labelling each cortical region with the percentage of significant voxels within the overlapping cluster. To facilitate visualization of the largest asymmetric regions, only the parcellated regions with at least 10% of voxels having survived permutation testing are shown on the surface. Figure 5d shows the timing of fetal brain lateralization during the second trimester, and the regions of significant difference between the two hemispheres.

Temporal patterns of maturation

To examine the emergence and evolution of internal brain structures, voxel-wise statistical analyses were performed on pairs of weekly atlases, each separated by 2 weeks, by gathering the scans from the two timepoints (a and b = a + 2) and generating a single groupwise atlas template from all scans (Extended Data Figs. 3 and 4). The atlases were spatially normalized using a global affine transformation to remove size effects. The groupwise registration step yielded voxel-level deformation maps that would map each scan to the central (median) age. That is, the structures in each scan were reconfigured such that the anatomies in the earlier and later gestational timepoints were deformed to the same template, representing the brain at a gestational age between a and b. Discovery of age-group-specific differences was achieved in two ways. First, we computed the log-Jacobian maps (\({J}_{i}^{t}\), which show regions of local structural changes associated with growth or shrinkage) for each deformation field map, and subtracted the mean maps of the two groups:

Positive regions indicate structural expansion from age a to b, and negative values correspond to regions undergoing shrinkage (Fig. 4a and Extended Data Figs. 3 and 4).

Second, we performed voxel-based morphometry with permutation testing implemented in FSL RANDOMISE⁶⁹ with threshold-free cluster enhancement⁷⁰ to highlight regions that were statistically significantly different between the two gestational timepoints. We tested for main group effects by constructing a design matrix with residualized gestational age. Again, 5,000 permutations were performed for each contrast, and only the voxels surviving a conservative threshold of P < 0.05 were considered as significantly evolving between the two timepoints (Extended Data Figs. 11 and 12).

Within-population structural variance

For a given gestational age, each fetus has a spatial map of the amount of voxel-level deformation required for their brain to match that of the population average. By aggregating these deformation fields, we performed voxel-level PCA to determine the breadth of healthy phenotypic structural presentation.

We visually detected a progressive increase in inter-participant variability across the gestational period (Fig. 3e) and for all volumetric measures (Fig. 6), which was confirmed by White’s Lagrangian test for heteroscedasticity⁷¹.

Ethics

The INTERGROWTH-21st Project and its ancillary studies were approved by the Oxfordshire Research Ethics Committee ‘C’ (reference no. 08/H0606/139), the research ethics committees of the individual participating institutions, as well as the corresponding regional health authorities where the project was implemented. All mothers provided written informed consent for the use of their clinical data. The sponsors had no role in the study design, data collection, analysis, interpretation of the data, or writing of the paper. The following authors had access to the full raw data set: R.B.G., S.H.K., A.I.L.N., A.P., and J.V. The corresponding author had full access to all the data and final responsibility for submitting the paper.

Data analysis software statement

Statistical analysis was carried out with the Python statsmodel package (v.0.13.2) and the FSL RANDOMISE tool (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Randomise; v.6.0.5). The deep learning models used to perform whole-brain extraction and alignment are available on Github (https://github.com/felipemoser/kelluwen), as is the model used to segment the subcortical structures (https://github.com/lindehesse/FetalSubcortSegm_Code). The atlas was constructed using a script written in MATLAB (v.R2022a), adapted from an implementation of diffeomorphic log-demons image registration (https://www.mathworks.com/matlabcentral/fileexchange/39194-diffeomorphic-log-demons-image-registration). All data analysis scripts were written in Python (v.3.9.6). Plots were generated using the Python seaborn package (v.0.12.1), and cortical surface maps were created using the Python-based ggseg package (v.0.1).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.