Automated processing pipeline for neonatal diffusion MRI in the developing Human Connectome Project

Introduction

A key challenge in neuroscience is to understand how the human brain grows and adapts during development (Kaiser, 2017). Magnetic Resonance Imaging (MRI) and the advent of connectome mapping methods have revolutionized the way we can approach this challenge by allowing non-invasive multi-modal explorations in living subjects (Miller et al., 2016; Van Essen and Ugurbil, 2012). Connectome studies require large datasets, and to date have been most successful in defining cerebral topography (Glasser et al., 2016) or linking connectivity to behaviour in adults (Smith et al., 2015). The developing Human Connectome Project (dHCP, http://www.developingconnectome.org) is set to create and make openly available a 4-dimensional map of early structural and functional developmental changes. A large dataset of multi-modal MRI scans from fetuses and newborn infants aged 20–44 weeks post-menstrual age is currently being acquired, alongside genetic, clinical and developmental information, with the aim of allowing researchers to explore genetic and environmental influences on brain development, and to further our understanding of the relation between cerebral connectivity and neurocognitive function.

In this work, we focus on newborn infants and on structural connectivity assessed using diffusion MRI (dMRI). This is a powerful technique to probe brain connections and microstructure during early development (Ball et al., 2013a, 2013b, 2015; Batalle et al., 2017b; Chung et al., 2016; Dubois et al., 2014a; Huppi and Dubois, 2006; Mishra et al., 2013; Pandit et al., 2013, 2014; Smyser et al., 2016; Yu et al., 2016). However, dMRI in newborn infants poses several challenges (Heemskerk et al., 2013; Pannek et al., 2012). Imaging newborn infants can be difficult and requires specialized medical facilities, operators, equipment and acquisition protocols (Hughes et al., 2017). Image analysis is challenging, and tools developed for adult brain data have to be refined or extended. Problems that are commonly overlooked or negligible when processing adult data, for example within-volume motion, need to be considered when analysing neonatal dMRI (Fig. 1A). Excessive motion can also cause severe signal dropouts that do not reflect the anatomically related anisotropic behaviour of water molecules; moreover, in combination with short repetition times it can cause significant intensity changes due to spin-history effects. The spin-history effects can be further aggravated by thick and overlapping slices due to the resulting reduced time between consecutive excitations of any given spin. Small and variable head sizes modulate the overall signal-to-noise ratio (SNR) and have to be accounted for (Fig. 1B). Cerebral free water content is higher in the infant brain, anisotropy is lower, and inhomogeneous maturation alters the signal profile significantly within a few weeks (Fig. 1C). Bespoke solutions are thus required for the successful connectomic analysis of neonatal dMRI.

In this work, we present a fully automated processing pipeline developed and optimised for neonatal dMRI data. The pipeline takes as input raw data and information about the scanning protocol, corrects for distortions, artefacts and subject motion, and super-resolves the originally acquired thick, but overlapping slices. The pipeline then estimates fibre orientations with methods that consider the inherently low anisotropy in white matter of the neonatal brain. These results are then used by an automated tractography framework to segment 16 (13 of which are bilateral) association, projection, limbic, cerebellar and commissural tracts in age-matched template space (Serag et al., 2012). This allows mapping tract-specific microstructural and morphological changes at the population level. Given the large number of subjects in the dHCP, an automated quality control (QC) framework has been implemented to assess quality and allow the detection of residual issues or inconsistencies with the data. The proposed data processing is applied to an exemplar dataset of infants imaged at approximately term-equivalent age.

Subjects and dMRI acquisition

MR images were acquired as a part of the dHCP; the study was approved by the National Research Ethics Committee and informed written consent given by the parents of all participants. Infants were fed and wrapped and encouraged to sleep naturally in the scanner; no sedation was administered. Images were acquired in the Evelina Newborn Imaging Centre, located within the Neonatal Intensive Care Unit of the Evelina London Children's Hospital. Pulse oximetry, temperature and heart rate were monitored throughout and ear protection was provided for each infant (President Putty, Coltene Whaledent, Mahwah, NJ; MiniMuffs, Natus Medical Inc., San Carlos, CA), and a paediatrician or neonatal nurse present throughout the scan. A dedicated patient handling system has been developed (Hughes et al., 2017) which encourages natural sleep and reduces the propensity of the infant to move during image acquisition.

In this exemplar study the images analysed were obtained from 140 infants (51 male, 89 female) born at median (range) 40 (36.6–42.1) and imaged at 40.9 (37.6–44.4) weeks post-menstrual age (Sup. Fig. 1). All images were reviewed by a paediatric neuroradiologist prior to inclusion in the dataset.

Preprocessing pipeline

Quality control

When dealing with large datasets such as the dHCP, automated quality control (QC) is of great importance for detecting data processing issues and inconsistencies. We have developed a dMRI QC framework that reads the outputs at different stages of the processing pipeline and generates summary data quality statistics (Bastiani et al., submitted). QC is performed at both the subject and the group level. Table 1 gives an overview of all the estimated QC indices.

At the data import stage the automated QC framework checks whether all 300 expected volumes have been acquired. The data are considered incomplete and not processed further when less than 11.3% (i.e., 34 volumes) of the total expected number of volumes have been acquired. This is the minimum number of volumes we need to be able to estimate an off-resonance field as two b0 volumes with opposite PE directions are acquired when reaching 34 volumes. If this threshold is reached, we are able to run through the distortion correction step and generate QC metrics. It does not imply optimality in any sense and the reduced number of volumes will be reflected in reduced values of all CNR metrics.

Absolute and relative motion are computed and reported. The former quantifies the displacement between the reference volume and each of the other volumes, the latter the displacement between two consecutive dMRI volumes. We found that for dHCP subjects born at term, the median across subjects of average absolute motion was 1.52 mm (std: 0.62, min: 0.56, max: 6.31 mm), i.e., slightly more than a voxel. As a comparison, we computed the average absolute within-session motion for a sub-set (n = 50) of subjects from the adult Human Connectome Project (HCP) and found a median of 0.83 mm (std: 0.29, min: 0.11, max: 3.21 mm).

Several outlier-related QC metrics are computed. Total number of outliers are split according to b-value shell, PE direction and single volume. The number of times a single slice has been marked as outlier is also computed to check for stationary spatial inconsistencies in the acquisition due to image reconstruction issues or hardware faults. When compared to HCP data, the median percentage of total outlier slices is almost 5 times higher (1.89% for dHCP, 0.39% for HCP).

To check the accuracy of susceptibility-induced distortion corrections and of the registration from dMRI space to subject native anatomical space, we also compute the correlation across brain masked voxels between the registered average b0 volume and the high-resolution T2-weigthed volume. The distribution of these correlation values helps in identifying problematic subjects.

Once registration to the subject's native anatomical space is complete, white and grey matter segmented binary masks are brought to dMRI data space. These are used to compute tissue-specific average contrast-to-noise ratios (CNRs) and SNRs. We use the data predictions obtained from the non-parametric Gaussian process embedded in the EDDY tool, evaluated using the final estimates of distortion fields, to compute CNR and SNR. Specifically, we obtain the voxel-wise CNR by dividing the standard deviation of the predicted signals by the standard deviation of the residuals (i.e., the difference between raw and predicted data) in undistorted space. We obtain the CNR separately for each acquired shell. Voxel-wise SNR is computed using only the b0 volumes as the ratio between the average predicted signal and the standard deviation of the residuals in undistorted space. The voxel-wise values are then averaged within WM/GM. Both CNR and SNR are useful to quantify the quality of a dataset (Fig. 7). For 140 dHCP subjects, the median CNR for the b = 1000 s/mm² shell in white matter is approximately 1.0 (std: 0.19, min: 0.12, max: 1.36), which compares to 0.45 (std: 0.2, min: 0.11, max: 0.73) for 50 HCP subjects. For 140 dHCP subjects, the median SNR in white matter is approximately 16 (std: 3.36, min: 6.65, max: 26.87), which compares to 10 (std: 1.2, min: 7.2, max: 11.3) for 50 HCP subjects. Differences in voxel size, though, play a role here as the acquired voxel volume for the dHCP is ∼3.5 times higher than for the HCP.

To investigate any potential differences or issues, we have also processed 25 extra preterm brains (age at scan range: 29–35 weeks). After running the full pipeline, we used the automated QC framework to assess differences in the distributions of quality metrics. No significant differences that would point to potential issues when processing such data were found (Sup. Fig. 4).

Microstructure modelling

The diffusion tensor (DT) model (Basser and Jones, 2002; Basser and Pierpaoli, 1996) is still the most widely used method to probe the brain's microstructure. However, the Gaussian displacement assumption underlying DTI breaks at high b-values. To overcome this limitation and account for deviations from Gaussian displacement, we use a diffusion kurtosis-augmented (DKI) tensor model (Jensen et al., 2005; Paydar et al., 2014). A single parameter is added to model mean kurtosis using a 2^nd order Taylor expansion of the log-attenuation in powers of b-values. This approach improves the overall fit to the multi-shell data and estimates a mean Kurtosis map (Fig. 8) that quantifies the deviation from Gaussianity of water molecules' displacement and can reflect different degrees of tissue heterogeneity (Steven et al., 2014). From the 6 independent components of the DT we derive the eigenvalues and eigenvectors and generate fractional anisotropy (FA), axial (AD), radial (RD) and mean diffusivity (MD) maps. Fig. 8 shows a longitudinal atlas of the estimated DT-derived microstructure parameters at 38–44 weeks post-menstrual age, with 20 subjects per age-bin (see Sup. Fig. 5 for single subjects' results).

We also apply the neurite orientation dispersion and density imaging (NODDI) model (Tariq et al., 2016; Zhang et al., 2012). NODDI describes the dMRI signal as a summation of three different contributing pools of water, i.e., restricted (intra-neurite), hindered (extra-neurite) and isotropic and estimates an orientation dispersion index to account for fibre fanning (Fig. 9 and Sup. Fig. 6 for single subjects' results). To model the distribution of fibre orientations it uses either Watson or Bingham distributions. We use the Bingham distribution as it allows extra flexibility by describing fibre dispersion along two orthogonal axes (Tariq et al., 2016). By introducing the orientation dispersion index and estimating the volume fractions of the intra- and extra-neurite compartments, NODDI can model tissue microstructure and geometry that may alter during development. Due to the very long computation times required to fit such a model, we use a GPU-based library for accelerated nonlinear optimisation (Hernandez Fernandez et al., 2017, www.fmrib.ox.ac.uk/∼moisesf/cudimot).

Fibre orientation modelling

In addition to microstructure modelling, we also estimate the fibre orientation densities (FODs) that will be used for tractography. We perform the estimation using a parametric framework, based on the assumption that the voxel-wise dMRI signal is the convolution of the FOD and a single-fibre response function (Behrens et al., 2003; Tournier et al., 2004).

Differences in microstructural properties across tracts, for instance between projection and association white matter tracts, together with axonal pruning and glia development (Dehaene-Lambertz and Spelke, 2015; Dubois et al., 2008, 2014a; Huang et al., 2006; Lebel and Deoni, 2018; Mukherjee et al., 2002; Oishi et al., 2013; Wozniak and Lim, 2006), make the use of a single common single fibre response function sub-optimal. We perform parametric deconvolution using a multi-compartment model that treats the FOD as a sum of delta functions (Behrens et al., 2007; Jbabdi et al., 2012; Sotiropoulos et al., 2016). Compared to previous parametric models that use an infinitely anisotropic (stick) single-fibre response model, we use a wider profile using the response from an axially symmetric tensor (i.e. with equal second and third eigenvalues). This extension captures better the inherently low anisotropy of neonatal dMRI data. We allow changes in the anisotropy of this single-fibre response kernel (Sotiropoulos et al., 2016), as its width (λ_R, the ratio of the minor and major axis of the tensor) is a model parameter which is estimated. We provide prior knowledge for this parameter, using a Gaussian prior for λ_R, with the prior mean and standard deviation being estimated from the data. To do that, we sample from a white matter region where one coherent fibre population is expected. Specifically, we used a standard mask of the corpus callosum and retained, for each subject, the 50% voxels with highest FA values. We then compute the mean and standard deviation of the voxel-wise ratio λ_R across subjects. In the exemplar dataset of 140 infants the average λ_R of 0.3 and an average standard deviation of 0.08; values respectively 2.3 and 2.7 times larger than those estimated for HCP subjects (Sotiropoulos et al., 2016), reflecting the lower anisotropy in white matter.

This parametric model allows voxel-specific single-fibre response kernels (all informed by the same prior) and can estimate fibre crossings in voxels within deep white matter regions that show low degrees of anisotropy. To estimate all the different parameters and their uncertainties, we use Bayesian inference and a Rician noise model. The estimation framework is capable of removing unnecessary fibre orientations by using automatic relevance determination (ARD) priors, as in (Behrens et al., 2007).

Fig. 10 shows averaged volume fraction maps of secondary and tertiary fibre orientations from 38 to 44 weeks post-menstrual age, where 2^nd and 3^rd fibre populations become more prominent over time, particularly in regions where adult brain shows multiple fibre populations crossing such as periventricular deep white matter regions. The figure shows the results obtained from a model selection procedure (ARD; see Behrens et al., 2007). This method automatically rejects extra fibre compartments that are not strongly supported by the data. Therefore, because those regions show a very low diffusion contrast-to-noise ratio (CNR) due to low fibre coherence and/or high water content, it is very challenging to reliably estimate more than one fibre compartment. The increased ability to detect crossings at older ages is indicative of higher diffusion contrast in the data, caused by higher structural coherence and myelination in white matter. This is probably reflecting the well-known increase of myelination that reduces the extra-axonal water fraction and the development of short cortico-cortical projections (Takahashi et al., 2012; Xu et al., 2014).

Tractography

Given the significant structural connectivity changes that occur in the neonatal brains, it is of interest to build an atlas of white matter connections, standardized across different developmental stages. We have developed standard-space tractography protocols that allow to automatically reconstruct 16 white matter tracts (see Table 2) across subjects and across different ages.

Specifically, we defined a set of standardized ROIs in template space (Serag et al., 2012) for week 44, i.e., the available time point closest to adult age (Fig. 11). Based on white matter atlases (Catani and Thiebaut de Schotten, 2008; Wakana et al., 2004) and using the segmentation results from the dHCP structural pipeline, we defined sets of “seed” (regions where tracts start), “waypoint” (regions through which tracts should go through) and “stop” (regions where tracts terminate) masks for each tract (de Groot et al., 2013). “Exclusion” (regions that tracts should avoid) masks were also used for all tracts, comprising at least the whole mid-sagittal plane, except for the commissural projections.

Prior to starting tracking, the ROIs of each selected tract are warped to the template's time point matched to the subject's age at scan (in weeks) using non-linear registration, FNIRT (Andersson et al., 2003). This step is done to ensure that all the tracking results are stored in age-matched template space for further processing. Probabilistic tractography is then performed in subject's native diffusion space and the results are resampled from diffusion space to age-matched template space using the non-linear warp fields. When the tracking is complete, we obtain tract probability maps by normalising the resulting volumes storing the number of streamlines passing through each voxel by the total number of successful streamlines, i.e., those that were not rejected by any exclusion criterion. Finally, we average all the age-matched tract probability maps to obtain group maps. Fig. 12 shows such average tract probability maps for all the bilateral tracts. The tool can easily be extended/modified and is computationally efficient as it uses GPUs (Hernandez Fernandez et al., 2016).

Defining subject-specific tracts using tractography can be very useful for assessing tract-specific microstructural changes in the neonatal brain (Dehaene-Lambertz and Spelke, 2015; Dubois et al., 2008, 2014a; Huang et al., 2006; Lebel and Deoni, 2018; Mukherjee et al., 2002; Oishi et al., 2013; Pecheva et al., 2017; Wozniak and Lim, 2006). Most studies investigating white matter microstructure in newborns use atlas-based ROIs (Kersbergen et al., 2014) to define the core of the tracts or manually segment them (Kunz et al., 2014; Paydar et al., 2014). As an example of the use of pipeline data we used the results from the automated framework to define age-specific white matter masks to extract tract-specific microstructural indices. Fig. 13 and supplementary figure 6 show the results of a linear regression analysis relating tract-specific microstructural indices obtained from DKI and NODDI in the exemplar dataset of 140 subjects (see supplementary table 1 for a full list of regression coefficients). We added total brain volume as an extra covariate to regress out partial volume effects and de-meaned both co-variates. Tractography results were averaged across individuals in age-matched template space, thresholded at 0.01 and binarised to obtain tract-specific masks. The differences in the magnitude of the age-related regression coefficients reflect the different developmental pathways of individual white matter tracts, which is well known to be asynchronous both from post-mortem (Kinney et al., 1988; Rados et al., 2006) and in vivo dMRI studies (Dehaene-Lambertz and Spelke, 2015; Dubois et al., 2008, 2014a; Huppi and Dubois, 2006; Kersbergen et al., 2014; Kunz et al., 2014; Mukherjee et al., 2002; Oishi et al., 2011; Paydar et al., 2014).

Relating the intercepts between different microstructural indices obtained from the same analysis technique also reveals clusters based on tract category. A general increase in FA and decrease in MD with age is matched by a significant increase of intra–neurite volume fraction and a decrease of isotropic volume fraction contamination. The MD and isotropic volume fraction changes point to the generalized reduction of water content occurring in the developing brain (Neil et al., 1998).

Discussion

In this work, we presented a comprehensive and automated pipeline to analyse neonatal dMRI data as acquired in the developing Human Connectome Project. Previous developmental studies have addressed single aspects of foetal (Oubel et al., 2012; Studholme, 2011) or neonatal data analysis, such as pre-processing (Dubois et al., 2014b; Holdsworth et al., 2012; Morris et al., 2011), microstructural modelling (Huppi et al., 1998; Kunz et al., 2014) and structural connectivity (Howell et al., 2018; Tymofiyeva et al., 2012). Our pipeline allows to analyse consistently a large number of subjects, despite the challenges that characterize neonatal dMRI data analysis (Fig. 1). We optimised existing methods and developed new ones to successfully pre-process the datasets. With a single tool, we can correct susceptibility and eddy-current-induced distortions, between (Andersson and Sotiropoulos, 2016) and within-volume (Andersson et al., 2017) motion artefacts and signal dropouts (Andersson et al., 2016). Moreover, we introduced an automated QC framework that allows to detect issues and data inconsistencies occurring at different stages of the pipeline and to quantify overall data quality. Such a framework is crucial when dealing with big cohorts of largely non-compliant subjects that make manual QC very time consuming and subjective.

The neonatal brain is not just a miniature adult brain (Batalle et al., 2017a), but poses several analysis challenges. The dMRI signal profile is very different when compared to the adult. Newborns are characterised by higher water content and lower myelination. Moreover, grey and white matter also show different microstructural features. Previous work argued that using a low b-value (∼700 s/mm²) might be optimal for neonatal dMRI due to the higher water content (Pannek et al., 2012). For the dHCP, we have used a data-driven approach to determine the optimal combination of diffusion encoding orientations and b-values (Tournier et al., 2015; Tournier et al., 2018, under revision). This allowed us to maximise the contrast between different tissue types and to use advanced analysis techniques that can estimate complex fibre geometries.

Furthermore, newborns typically show higher levels of motion when compared to adult subjects. Artefacts that are not so noticeable in adults become prominent in neonates, such as within-volume (sudden) motion, slice misalignment, and signal dropouts. The gyrification pattern, contrast and head sizes vary significantly throughout development. Therefore, we have developed new approaches that address such bespoke scenarios. Additionally, we used age-specific templates and transformations, as well as a specialized MR coil that improves SNR when imaging small brains.

Most previous studies have used DTI to assess microstructural changes and estimate structural connectivity profiles. The data quality and acquisition strategy followed in the dHCP allow the use of more advanced multi-shell analysis techniques, which in principle can estimate richer microstructural features and multiple fibre orientations within the same voxel. To realise these benefits in tractography, we used an approach for estimating fibre crossings, which accounts for the low anisotropy of the neonatal brain by assuming a non-perfectly anisotropic response function. This is crucial when estimating fibre orientations using a spherical deconvolution approach, which typically assume the same deconvolution kernel across the whole brain (Parker et al., 2013). This assumption does not hold in the neonatal brain, where localised changes in anisotropy occur throughout development due to changes in myelination and development of cortico-cortical projections. We therefore utilised an approach that relaxes this assumption and allows variations of the deconvolution kernel through space.

Using the voxel-wise estimated fibre orientations, we developed a neonatal-specific automated tractography tool. This allows to automatically segment white matter into a set of subject-specific long-range projections. These can be used to define regions of interests for further tract-specific microstructural or morphological analyses.

All the data used in this work will be shared with the scientific community. Raw and processed dMRI, structural and rfMRI datasets will be made available at http://www.developingconnectome.org as part of the official dHCP public data releases. This includes atlases and ROIs for, e.g., automated tractography. Processing pipelines will also be made available on the same platform. The current version of the dHCP neonatal dMRI pipeline can already be downloaded from: https://git.fmrib.ox.ac.uk/matteob/dHCP_neo_dMRI_pipeline_release.

In this work, we have described the current version of the dMRI neonatal processing pipeline. However, there are ongoing developments to address additional challenges in the data analysis stream. For instance, because of the high degree of motion present in neonatal data, the susceptibility field can change significantly during the acquisition. This means that a single fieldmap estimate is not sufficient and should be updated in the presence of big rotations and translations (Andersson et al., 2018). Moreover, registration to anatomical space can be improved. Diffusion MRI data provides additional information when compared to structural acquisitions (e.g., T1 or T2-weighted). This can be used to improve the registration to subject's anatomical space using, e.g., the diffusion tensor (Zhang et al., 2006). Tractography results can also be improved by using methods that allow the estimation of complex asymmetric sub-voxel fibre patterns (Bastiani et al., 2017). Both cortical microstructure and connectome reconstruction can benefit from improved surface registration, such as advanced techniques that allow to register surfaces based on multiple modalities (Robinson et al., 2013), together with templates that preserve cortical locations across ages.

Conclusions

We have presented an automated pipeline to pre-process neonatal dMRI data that is capable of effectively dealing with the main challenges presented by neonatal data and providing high-quality distortion correction, brain microstructure and connectivity results across early developmental ages, and provided exemplars of the output of the pipeline and the use of the data. The pipeline will allow researchers to use the dHCP datasets to explore important questions about the development of the human brain.

Appendix A. Supplementary data

The following is the supplementary data related to this article: