Advanced diffusion imaging for assessing normal white matter development in neonates and characterizing aberrant development in congenital heart disease

Background Elucidating developmental trajectories of white matter (WM) microstructure is critically important for understanding normal development and regional vulnerabilities in several brain disorders. Diffusion Weighted Imaging (DWI) is currently the method of choice for in-vivo white matter assessment. A majority of neonatal studies use the standard Diffusion Tensor Imaging (DTI) model although more advanced models such as the Neurite Orientation Dispersion and Density Imaging (NODDI) model and the Gaussian Mixture Model (GMM) have been used in adult population. In this study, we compare the ability of these three diffusion models to detect regional white matter maturation in typically developing control (TDC) neonates and regional abnormalities in neonates with congenital heart disease (CHD). Methods Multiple b-value diffusion Magnetic Resonance Imaging (dMRI) data were acquired from TDC neonates (N = 16) at 38 to 47 gestational weeks (GW) and CHD neonates (N = 19) aged 37 weeks to 41 weeks. Measures calculated from the diffusion signal included not only Mean Diffusivity (MD) and Fractional Anisotropy (FA) derived from the standard DTI model, but also three advanced diffusion measures, namely, the fiber Orientation Dispersion Index (ODI), the isotropic volume fraction (Viso), and the intracellular volume fraction (Vic) derived from the NODDI model. Further, we used two novel measures from a non-parametric GMM, namely the Return-to-Origin Probability (RTOP) and Return-to-Axis Probability (RTAP), which are sensitive to axonal/cellular volume and density respectively. Using atlas-based registration, 22 white matter regions (6 projection, 4 association, and 1 callosal pathways bilaterally in each hemisphere) were selected and the mean value of all 7 measures were calculated in each region. These values were used as dependent variables, with GW as the independent variable in a linear regression model. Finally, we compared CHD and TDC groups on these measures in each ROI after removing age-related trends from both the groups. Results Linear analysis in the TDC population revealed significant correlations with GW (age) in 12 projection pathways for MD, Vic, RTAP, and 11 pathways for RTOP. Several association pathways were also significantly correlated with GW for MD, Vic, RTAP, and RTOP. The right callosal pathway was significantly correlated with GW for Vic. Consistent with the pathophysiology of altered development in CHD, diffusion measures demonstrated differences in the association pathways involved in language systems, namely the Uncinate Fasciculus (UF), the Inferior Fronto-occipital Fasciculus (IFOF), and the Superior Longitudinal Fasciculus (SLF). Overall, the group comparison between CHD and TDC revealed lower FA, Vic, RTAP, and RTOP for CHD bilaterally in the a) UF, b) Corpus Callosum (CC), and c) Superior Fronto-Occipital Fasciculus (SFOF). Moreover, FA was lower for CHD in the a) left SLF, b) bilateral Anterior Corona Radiata (ACR) and left Retrolenticular part of the Internal Capsule (RIC). Vic was also lower for CHD in the left Posterior Limb of the Internal Capsule (PLIC). ODI was higher for CHD in the left CC. RTAP was lower for CHD in the left IFOF, while RTOP was lower in CHD in the: a) left ACR, b) left IFOF and c) right Anterior Limb of the Internal Capsule (ALIC). Conclusion In this study, all three methods revealed the expected changes in the WM regions during the early postnatal weeks; however, GMM outperformed DTI and NODDI as it showed significantly larger effect sizes while detecting differences between the TDC and CHD neonates. Future studies based on a larger sample are needed to confirm these results and to explore clinical correlates.


Introduction
In the young adult human brain, white matter (WM) contains 149,000 to 176,000 km of myelinated axons (Marner et al., 2003). WM is the wiring structure of the brain that allows effective communication between different cortical gray matter areas that have distinct functional characteristics. The organization and maturation of axonal pathways in regional WM structures follow heterogeneous trajectories beginning in the prenatal to early postnatal stages. For example, postmortem studies have shown that the Posterior Limb of Internal Capsule (PLIC) has partial myelination seen at birth (Brody et al., 1987). On the other hand, in the Corpus Callosum (CC), myelination progresses much later during the postnatal period (Kinney et al., 1988). Further, the Superior Longitudinal Fascicle (SLF), an association pathway, has a slow maturation period that extends into adulthood (Zhang et al., 2007). Origins of several neurodevelopmental and psychiatric disorders have been associated with the dysregulation in the organization and maturation of the WM regions (Dietrich et al., 1988;Tkachev et al., 2003;Pujol et al., 2004;Cascio et al., 2006;Courchesne et al., 2007;Brennand et al., 2011).
Similarly, individuals with Congenital Heart Disease (CHD) make up a population that has been associated with a history of neurologic and neuro-developmental impairments such as motor and visuospatial skills, as well as cognitive impairments such as working memory, attention, and language (Bellinger et al., 1999;Limperopoulos et al., 1999;Hovels-Gurich, 2002;Bellinger et al., 2003a;Bellinger et al., 2003b;Hövels-Gürich et al., 2007. Moreover, previous neonatal brain imaging studies have shown evidence of delayed development of pyramidal tracts, presence of white matter injury associated with hypoxia, lower gray and white matter volume, abnormal cerebral blood flow, and abnormal metabolism in CHD neonates (Mahle et al., 2002;Licht et al., 2004;Pas te et al., 2005;Partridge et al., 2006;McQuillen et al., 2007;Petit et al., 2009;Ortinau et al., 2012;Dimitropoulos et al., 2013;Dehaes et al., 2015). Given the importance of the WM during development and its association with neurodevelopmental disorders, it is crucial to accurately characterize the normal development in typically developing neonates and potentially to detect aberrant development of regional WM structures in neonates with CHD. Increasing our understanding about affected WM regions will allow for better characterization of WM abnormalities which are needed to guide the search for potential etiologies and to understand cognitive outcomes. Advanced MRI techniques such as diffusion Magnetic Resonance Imaging (dMRI) have been crucial in characterizing the early organization of regional WM in neonates, and are essential for investigating abnormalities in early WM development. Moreover, past studies have used measures such as Diffusion Tensor Imaging (DTI) and Neurite Orientation Dispersion and Density Imaging (NODDI) to study the normal development and to explore WM abnormalities in CHD (Kunz et al., 2014).
While DTI (Basser et al., 1994) is a classic method used to determine WM structural integrity, it has several limitations: 1) it assumes the existence of only a single fiber population oriented in one particular direction (Tuch et al., 2002), and 2) it assumes Gaussian diffusion with no restriction to motion of water molecules, which is an oversimplification of the underlying biological process. While FA derived from DTI has been reported to be sensitive to myelination of axons (Chenevert et al., 1990), it is not specific to one particular type of abnormality, i.e., several different biological processes can produce similar types of alterations in FA. Hence, it cannot provide specific information about the nature of the microstructural abnormalities (Westin et al., 2016).
To overcome the limitations of DTI, several compartmental models of diffusion have been proposed. For example, the NODDI (Zhang et al., 2012) model attempts to characterize neural tissue structure with three different compartments representing restricted, hindered, and isotropic diffusion. In addition, NODDI models the dispersion of axonal fibers with the use of an Orientation Dispersion Index (ODI). One key limitation of the NODDI model is that it requires an a-priori chosen value to be set for the diffusivity along the axons, which may not be accurate in all regions of the developing neonate brain and may be difficult to estimate when development or disease alters brain structure (Lampinen et al., 2017). Nevertheless, NODDI has been used to study the microstructure of the neonate brain (Kunz et al., 2014). However, the potential of this model to characterize better regional WM tract development in normal neonates and to detect changes in neonates with CHD remains unexplored.
An alternative approach to DTI and NODDI is the Gaussian Mixture Model (GMM). Unlike DTI, the GMM does not assume Gaussian diffusion and, unlike NODDI, the GMM is independent of tissue models, i.e., it does not make any assumption about the number of compartments or their diffusivities. Instead, the GMM estimates the Ensemble Average diffusion Propagator (EAP), which describes the probability distribution of the displacement of water molecules within an experimentally set diffusion time (Cheng, 2012;Özarslan et al., 2013;Ning et al., 2015). The EAP allows calculation of several tissue microstructure related scalar measures, such as the Return-to-Origin Probability (RTOP) and Return-to-Axis Probability (RTAP) (Ning et al., 2015). RTOP is the probability that a water molecule returns to its starting position (the origin) within the experimental diffusion time. Thus, in a highly restricted or hindered medium, we expect high RTOP values, whereas, in the case of free diffusion, the RTOP values will be low, as there is very little probability that the water molecules will return to their starting position. It has also has been reported in the physics literature that RTOP is inversely proportional to the volume filled by air in a porous material (Özarslan et al., 2013;Ning et al., 2015). In the case of complex biological brain tissue, it provides a measure of the total volume within which water diffusion occurs. This, in turn, is related to the cellular and axonal volume, size, and myelination. For instance, a densely packed set of myelinated axons will leave very little extra-cellular space within which water diffuses (since most of the space is occupied by myelin and cell membranes), leading to small diffusion volume but higher RTOP (since volume is inversely proportional). Additionally, one can compute RTAP, which measures the probability of water molecules returning back to the axis or line representing the principal diffusion direction. RTAP is inversely proportional to the transverse cross-sectional area of the fiber bundles and hence is related to axon diameter, packing and amount of myelination. Since the GMM assumes very little about the diffusion process, it can be very sensitive to abnormalities in the diffusion process but lacks the ability to provide model specific information such as fiber orientation dispersion. The potential for GMM to detect and to characterize white matter abnormalities in CHD is further explored in this work.
In this paper, we explore the potential of both NODDI and the GMM to characterize regional WM development in TDC neonates as well as to detect abnormalities in specific WM regions in neonates with CHD. We also compare these novel measures with the standard DTI based measure of fractional anisotropy (FA) and mean diffusivity (MD). Thus, the present study also aims to compare different dMRI models to determine: 1) the sensitivity of different diffusion measures in their ability to detect age-related (cross-sectional) regional WM changes in TDC neonates, and 2) which measures detect regional WM differences in CHD compared to TDC neonates.

Materials and methods
The cross-sectional study included 16 Typically Developing Controls (TDC) (mean age: 42.03 ± 2.28 weeks; Gender: 12 females, 4 males) and 19 with Congenital Heart Disease (CHD) (mean age: 39.54 ± 1.08 weeks; Gender: 3 females, 16 males). All the CHD and TDC neonates were screened at the Boston Children's Hospital and have birth gestational week (GW) greater than 36 weeks. The TDC cohort included typically developing neonates recruited from a well baby nursery with uneventful deliveries, including Apgar scores greater than 8 at 5 minutes and no clinical signs or symptoms concerning any brain disorder. The CHD cohort included three different subtypes of CHD: 6 neonates with transposition of great arteries (TGA), 10 single ventricle physiology, and 3 bi-ventricle physiology. The details of the specific cardiac anomalies are summarized in Table 1. In addition, injuries observed in T1, T2, diffusion and Susceptibility Weighted Images are summarized in Table 2. MR acquisition was performed on a 3 T TimTrio Siemens system using a 32-channel head coil and a simultaneous multislice acquisition sequence (Setsompop et al., 2012). A multi-b-value DWI protocol was used with a total of 81 axial slices acquired with four non-diffusion weighted images (at b = 0 s/mm 2 ). Acquisition parameters were: TR = 3700 ms, TE = 104 ms, flip angle = 90°, 2 mm isotropic spatial resolution, with two b-value shells at b = {1000, 2000} s/ mm 2 each shell having 30 gradient directions, with a total acquisition time of about 6 min.

Data preprocessing
All DWI images first went through a semi-automatic quality control pipeline (in-house MATLAB script) that detects signal intensity drop in each slice by measuring the divergence of signal intensity between adjacent slices. Any dMRI volume with significant motion artifact or signal loss was manually inspected and removed. Next, we performed head motion and eddy current correction for each data set using the FSL FLIRT software (Jenkinson and Smith, 2001;Jenkinson et al., 2002) where the gradient directions were appropriately corrected using the rotation parameters obtained from rigid registration to the b = 0 image.

Diffusion Tensor Imaging Model
Diffusion tensors were estimated at each voxel in Slicer Version 4 (http://www.slicer.org), using weighted linear least squares fitting to obtain the eigenvectors and eigenvalues, which were then used to calculate Fractional Anisotropy (FA) and Mean Diffusivity (MD).

Neurite Orientation Dispersion and Density Imaging (NODDI)
The NODDI model decomposes the dMRI signal into three subdivisions or compartments: the intra-cellular signal, the extra-cellular signal, and the isotropic compartment (Zhang et al., 2012). The model incorporates the dispersion of fibers while modeling the intra-cellular compartment using a Watson distribution function, which can be used to estimate the fiber ODI (Zhang et al., 2012). The intra-cellular compartment (V ic ) of the NODDI model represents diffusion within the axons and cells. The isotropic compartment (V iso ) estimates the volume fraction of freely diffusing extracellular water such as cerebrospinal fluid. The normalized signal, A, in the NODDI model is given by the following equation: where, A ic is the signal contribution from the intra-cellular compartment and A ec is the signal due to diffusion in the extra-cellular space. We should note that the NODDI model assumes an a-priori fixed diffusion coefficient, 1.7 × 10 −3 mm 2 / s, along the principal diffusion direction for both the intra-and extra-cellular compartments. We used the default fixed parameters while estimating all the measures (Orientation Dispersion Index, Intra-cellular Volume Fraction, and Isotropic Volume Fraction) for the NODDI model using the MATLAB software that is publicly available (Zhang et al., 2012).

Gaussian Mixture Model
The GMM estimates the 3-dimensional probability distribution of the displacement of water molecules in a given experimental time.
Analytical expressions for computing measures such as RTOP and RTAP were proposed in Ning et al. (2015), which we used in this work. Unlike compartmental models, this method utilizes a directional radial basis function for reconstructing the diffusion signal as a linear combination of Gaussian basis functions centered across q-space, and not just at the origin. As noted previously, these measures are sensitive to the axonal and cellular size, density, and volume. All of these measures were computed using the publicly available software (https://github.com/ LipengNing/RBF-Propagator/).

Atlas registration
To analyze the microstructure in different brain regions, we used the John Hopkins Neonate Atlas (JHNA) (Oishi et al., 2011) to define the regions of interests (ROIs) in the neonate brains. This atlas contains 122 WM and gray matter ROIs, which have been used to study brain development in neonates with an age range of 37 to 53 GW (Oishi et al., 2011). All diffusion measures were computed in the native subject space, and were up-sampled linearly to the voxel size of 0.6 mm 3 to match the resolution of JHNA. Advanced Normalization Tools (ANTs) (Avants et al., 2008) were used to first affine transform the JHU-neonate FA atlas to each neonate subject FA space, followed by a non-rigid registration. The registration matrices obtained during these transformations were applied to the JHU-neonate label map to obtain subjectspecific definitions of several WM regions. In this work, we focused our analysis on deep WM regions because deformable atlas registration was robust and consistent in registering the JHU-neonate label primarily in the deep WM regions. Mean MD, FA, V ic , V iso, ODI, RTAP, RTOP were calculated from 22 ROIs (11 left hemispheres and 11 right hemispheres) defined by the JHU atlas for each individual subject.

Table 1
Details of cardiac anomaly including specific diagnosis, group and obstruction in congenital heart disease neonates (N = 19). Our first analysis focused on understanding the changes in diffusion measures in TDC neonates at different postnatal weeks (age) in all the 11 WM regions. We grouped the 11 ROIs into two groups: 1) projection fibers, and 2) association and callosal fibers (see Fig. 2). Regression analysis with respect to GW was performed for MD, FA, V ic , V iso, ODI, RTAP and RTOP for the TDCs using R software package. False discovery rate (FDR) correction was used to account for multiple comparisons. We particularly focused on three WM regions namely, PLIC, CC and SLF to perform qualitative comparison to determine which ROI has the highest and lowest values during the early postnatal developmental stage.

WM abnormalities in CHD
Since there were age differences between the CHD and TDC neonates, we first regressed out the effect of age from all TDC and CHD subjects. Subsequently, a t-test was used to estimate statistical differences between the two groups for each of the 7 measures and the 22 WM regions (a total of 154 t-tests). FDR correction was used to account for multiple hypothesis testing.  Table 3. Below, we provide more details about the results for each diffusion model. We also provide mean and standard deviation values for both TDC and CHD neonates for each of the diffusion measures in (Table 4 for TDC, and Table 5 for CHD).   1-A.7), we present the results from Fig. 2 in a different format to appreciate the slope (increase or decrease) of diffusion measures in different WM regions between TDC neonates at different ages. In particular, we ordered the slope of each region (for each measure) in decreasing order to fully understand the differential WM maturation between TDCs at different ages.

White matter regional abnormalities in CHD
We compared each of the diffusion measures between the TDCs and CHD subjects. An overall table of p-values for all measures is given in (Table 6), where we have reported the adjusted p-values after FDRbased multiple comparisons.
MD did not show any significant difference between TDC and neonates with CHD ( Fig. 2.a). FA was statistically lower in CHD compared to TDC neonates in bilateral CC, ACR, UF, left SLF and left RIC ( Fig. 2.b). V ic was lower in the CHD subjects in the bilateral SFOF and left PLIC (Fig. 2.c). Orientation Dispersion Index (ODI) was statistically higher in CHD neonates in the left CC (Fig. 2.d), while V iso showed no statistically significant differences for any WM structure (Fig. 2.e). RTAP was significantly lower in CHD in bilateral CC, UF, SFOF, and left IFOF (Fig. 2.f), while RTOP was statistically lower in CHD in several regions: bilateral CC, UF, SFOF, left IFOF, right ALIC and left ACR ( Fig. 2.g). Table 6 shows that measures derived from GMM detected more regional differences in white matter structures between CHD and TDC groups with more regions found to be statistically different using this model compared to the NODDI model. Further, the Cohen's d reported  Table 7, shows that the magnitude of effect sizes is larger for the GMM compared to the NODDI and DTI models.

Discussion
One of the key goals of our study was to characterize regional WM microstructure in typically developing neonatal brains using new advancements in dMRI. Thus far, past studies had used DTI, CHARMEDlight and NODDI measures (Kunz et al., 2014;Qiu et al., 2015). In this study, we not only characterized regional WM structures with DTI measures (FA, MD) and NODDI measures (Vic, Viso, ODI), but we also used tissue model-free GMM measures (RTAP, RTOP). This allowed us to compare the sensitivity of different diffusion methods to study specific microstructural changes in the early postnatal weeks.

Projection fibers
The results obtained from DTI measures demonstrate changes in the WM regions between TDC neonates at different postnatal weeks. More specifically, an increase in FA and a decrease in MD are observed across all WM regions with an increase of age between TDC neonates. Furthermore, the PLIC, where partially mature myelin can be observed at birth (Brody et al., 1987), shows the highest value for FA and V ic versus the lowest value for MD compared to other projection, association, and callosal fibers throughout the age range in our study. The results from RTOP and RTAP showed similar developmental trends in PLIC as the FA derived from DTI. The presence of myelination in PLIC decreases the extra-cellular space, making diffusion of water more restricted and hindered. Restriction of water with developing myelin increases the anisotropy of WM regions (Neil et al., 1998), leading to higher FA, V ic , RTOP and RTAP values in PLIC. Known trends in maturation from dorsal to ventral (Kinney et al., 1988) were also observed in MD, FA, V ic , RTAP, and RTOP.

Corpus callosum
The results from our study show that CC, unmyelinated at birth (Kinney et al., 1988), has lower FA and V ic values than PLIC but higher values than the other WM regions. Similarly, we observed the lowest ODI values in the CC compared to all other WM regions, which reflect high compactness and less dispersion of fibers. While the CC has high directional coherence of fibers, the partial volume from the ventricles might contribute to the lower FA compared to PLIC. We observed the highest isotropic volume fraction, V iso , in the CC, indicating the existence of partial volume effects.
Moreover, the results from RTOP and RTAP showed much lower values in the CC compared to PLIC. We note that both these measures are sensitive to the overall diffusional volume. Thus, due to the presence of myelinated axons in PLIC, the extracellular volume is lower and the contribution of restricted diffusional fraction higher, leading to higher RTOP/RTAP values. On the other hand, in CC, which has mostly unmyelinated axons at birth, the extracellular volume fraction is slightly higher, leading to lower RTOP and RTAP values. Another reason for lower values in CC might be due to partial volume with ventricles, which are in close proximity.

Association fibers
We also observed differences in the pattern of maturation across different measures within association fibers of the ventral language pathway (UF and IFOF) showing greater maturation than those of the dorsal language pathway (SLF). Similarly, a past DTI study has shown that the SLF has a slow maturation rate extending into childhood (Zhang et al., 2007). This can also be observed in our results during the early postnatal weeks, as SLF has the lowest FA, V ic , RTOP, and RTAP values compared to other association fibers. Thus, our results show a heterogeneous maturation of different WM regions across the TDC population.

Model comparison for development
We observed the expected pattern of change with age across the three models used. However, when comparing PLIC with CC using FA, we see much fewer differences between the values between these two WM regions as compared to the greater difference seen with GMM. This increase in FA could be either due to increase in myelination from the PLIC or coherence from the CC. The values obtained from V ic were the closest between PLIC and CC, indicating that there might be a possibility of partial volume effects from the lateral ventricles in the CC affecting the DTI and GMM measures, and removing such effects in the NODDI model brought the values closer between these two regions (by modeling an isotropic diffusion fraction V iso ). However, it is also important to note that estimated parameter values for the NODDI model might be suboptimal as the diffusivity values that are by default fixed into the NODDI model may not be accurate for neonate's brains. For example, our estimate of the axial diffusivity in the TDC neonates in the CC is 0.0021 mm 2 /s, whereas the default value is 0.0017 mm 2 /s. Further, as shown in Jelescu et al. (2015), the estimation of the NODDI parameters might have several confounds, which could lead to over-or under-estimation of the volume fractions in neonatal data. We also note that fixing the diffusivity to the one found in the CC of TDC neonates also may not be optimal, as the axial diffusivity varies across different regions. In our study, the average axial diffusivity for SFOF was 0.0015 mm 2 /s, much lower than the axial diffusivity of 0.0021 mm 2 /s observed in the CC. This demonstrates that fixing axial diffusivity may be an inherent limitation of the NODDI model.
In summary, we found that all measures were sensitive to age-related changes in TDC neonates in the regional WM structure. However, DTI and NODDI models simplify the modeling of biological tissue thereby limiting its sensitivity to capture all the variations during the developmental period of the brain. Since GMM does not assume any tissue model, it allows more freedom to capture the variability of each WM region.

Differential development of WM in CHD
In this study, we identified region-specific WM abnormalities related to CHD. The differences between CHD and TDC were most prominent in the CC as well as the association fibers of the ventral language pathway (UF, SFOF, left IFOF) and dorsal language pathway (SLF). Furthermore, assessing WM regions with seven different dMRI measures, provided more information as to the nature of WM structural abnormalities in CHD neonates.

Corpus callosum
The DTI results demonstrate the lower integrity of WM regions in CHD neonates. A previous DTI study in neonates with CHD demonstrated slower maturation of WM regions including the posterior part of the CC (Mulkey et al., 2014). Similarly, our result demonstrates lower FA in the CC among CHD neonates. Although FA is not specific to any particular type of biological change, lower FA might suggest that the compactness of CC is lower in CHD neonates. This also reflects the results obtained from ODI, as the ODI values in left CC for CHD neonates are significantly higher than in TDC neonates, suggesting greater dispersion of fiber bundles. Furthermore, we observe significantly lower RTOP and RTAP in the CC of CHD neonates. Since the corpus callosum in neonates primarily contains densely packed unmyelinated axons, lower RTOP and RTAP could be due to larger extra-axonal and extracellular volume, which implies fewer axons and cells in the CC of neonates with CHD. On the other hand, increased ODI implies incoherent axonal layout. These results demonstrate the possible factors contributing to the abnormalities seen in CHD, given the relationship between CC and several developmental disorders including developmental language delay (Fabbro et al., 2002;Paul, 2010).

Ventral language pathway
In addition, the results from DTI demonstrate lower FA in the UF of CHD neonates implying delayed or abnormal maturity of this WM region. RTAP and RTOP were also significantly lower in this region. A previous DTI study in young adults with CHD showed lower FA in the left UF compared to controls, and the lower FA was positively correlated with verbal memory (Brewster et al., 2015). In our study, we were able to detect differences in FA, RTOP, and RTAP measures in the UF at birth among neonates with CHD. These results may indicate abnormal maturation pattern of UF during early postnatal weeks in CHD. Similarly, we observe statistically significant and lower RTOP and RTAP values in the left IFOF of neonates with CHD. We note that both the DTI and NODDI measures failed to detect these differences in the IFOF. Furthermore, UF and IFOF are believed to be part of the ventral language pathway involved in speech recognition, representation of lexical concepts and semantic processing of language (Chang et al., 2015). The ventral pathway is believed to be bilaterally organized within each hemisphere performing computationally different roles (Hickok and Poeppel, 2007). Although these roles are unclear, some studies have shown the dominance of right hemisphere for processing slow temporal features, whereas the left hemisphere supports rapid temporal features (Abrams et al., 2008). With increasing evidence of individuals with CHD facing neuro-developmental disorders including language dysfunctions, these results observed in the ventral pathway of language, with strong evidence of abnormalities in the left hemisphere, could be among the possible factors contributing to these impairments.

Dorsal language pathway
Our results also demonstrate significantly lower FA in the left SLF of neonates with CHD. SLF is part of the dorsal language pathway which is responsible for the phonological processing of language (Chang et al., 2015). Unlike the ventral pathway, the dorsal pathway is strongly left hemisphere dominant (Hickok and Poeppel, 2007). The result in our study implies lower integrity of the left SLF in neonates with CHD, which could also contribute to language dysfunctions faced later in life. Although, we observed lower RTOP and RTAP in the SLF, there was no statistically significant difference after correction for multiple comparisons. It is important to note that the SLF is among the least mature WM regions observed during the first postnatal weeks, and has slower maturation and myelination rate than the regions of the ventral Fig. 1. a. Regression analysis of MD with postnatal age in TDC neonates shows decreasing MD for bilateral ALIC, PLIC, RIC, ACR, SCR, PCR, SLF, SFOF, IFOF and right UF (FDR corrected p < 0.01). Other decreases were not significant. MD measures the amount of water diffusion in the brain tissue. b. Regression analysis of FA with postnatal age in TDC neonates shows increasing FA for bilateral ALIC (FDR corrected p < 0.01). Other increases were not significant. FA measures the relative diffusivity of water molecules along the WM fibers compared to the perpendicular or radial direction. c. Regression analysis of V ic with postnatal age in TDC neonates shows increasing V ic for bilateral ALIC, PLIC, RIC, ACR, SCR, PCR, SFOF, IFOF, right CC, left SLF and left UF (FDR corrected p < 0.01). Other increases were not significant. V ic is derived from the NODDI model, and measures the intracellular volume fraction. d. Regression analysis of V iso with postnatal age in TDC neonates shows no significant changes with age. V iso is derived from the NODDI model, and measures the isotropic volume fraction. e. Regression analysis of ODI with postnatal age in TDC neonates shows increasing ODI for left SLF (FDR corrected p < 0.01). The remaining changes were not significant. ODI is derived from the NODDI model, and measures the fiber orientation dispersion index. f. Regression analysis of RTAP with postnatal age in TDC neonates shows increasing RTAP. Statistical significances were seen bilaterally in ALIC, PLIC, RIC, ACR, SCR, PCR, SFOF, IFOF, and left SLF (FDR corrected p < 0.01). The remaining increases were not significant. RTAP is derived from the Gaussian Mixture Model, and measures the return to axis probability. g. Regression analysis of RTOP with postnatal age in TDC neonates shows increasing RTOP. Statistically significant increases with age were seen bilaterally for ALIC, PLIC, RIC, SCR, PCR, SFOF and on the right for ACR and left for SLF (FDR corrected p < 0.01). The remaining increases were not significant. RTOP is derived from the Gaussian Mixture Model, and measures the return to origin probability. pathway in postnatal years (Zhang et al., 2007;Brauer et al., 2013). Given the slow maturation rate of SLF, the interpretation of abnormalities in early postnatal weeks for neonates with CHD is difficult.

Model comparison
Measures derived from the NODDI model showed far fewer statistical differences than those from DTI and GMM measures, indicating a lower sensitivity of the NODDI measures in capturing differences between CHD and TDC. Further, GMM-derived measures showed the highest sensitivity in terms of the effect size differences observed between TDC and CHD population. By carefully changing the default parameters, the sensitivity of NODDI might improve, yet it still might be suboptimal for several regions as the axial diffusivity is regionally very different (as discussed earlier). Further, axial diffusivity also might differ substantially in developing population and those affected by diseases. These reasons make it difficult to optimize this parameter in the NODDI model. Hence, in our study, we used default parameters setting in the NODDI mode. Overall, our study leveraged advanced diffusion imaging to obtain several parameters which allowed us to observe region-specific WM maturation in TDC and detect abnormalities in neonates with CHD in early post-natal weeks.

Limitation
Our study is limited in that 3 subjects were imaged postoperatively and 4 of 19 subjects (two imaged preoperatively and two imaged postoperatively) showed a recent injury on diffusion imaging. However, these acute lesions were small and, since recent, only two lesions were associated with decreased FA. Therefore, these lesions are likely to have little effect on our findings. Our study is also limited in assessing functional impairments associated with language development in the Fig. 2. a. Group comparison for MD between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF), Uncinate Fasciculus (UF). No statistical significant differences are observed in any of the tracts. Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum. b. Group comparison for FA between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF), Uncinate Fasciculus (UF). FA is significantly lower for CHD in bilateral CC, UF and left SLF (FDR corrected p < 0.01). Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum. c. Group comparison for V ic between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF), Uncinate Fasciculus (UF). V ic is significantly lower for CHD in bilateral SFOF (FDR corrected p < 0.01) Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum. d. Group comparison for V iso between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF), Uncinate Fasciculus (UF). No statistical significance is observed in any region. Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum. e. Group comparison for ODI between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF), Uncinate Fasciculus (UF). ODI is significantly higher for CHD in the left CC (FDR corrected p < 0.01). Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum. f. Group comparison for RTAP between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF) and Uncinate Fasciculus (UF). RTAP is significantly lower for CHD in the bilateral CC, UF, SFOF and left IFOF (FDR corrected p < 0.01). Note the improved separation of CHD from TDC with this measure. Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum. g. Group comparison for RTOP between CHD and TDC neonates in Corpus Callosum (CC), Inferior Fronto-occipital Fasciculus (IFOF), Superior Fronto-occipital Fasciculus (SFOF), Superior Longitudinal Fasciculus (SLF), Uncinate Fasciculus (UF). RTOP is significantly lower for CHD in the bilateral CC, UF, SFOF and left IFOF (FDR corrected p < 0.01). Graph displaying statistical data based on the minimum, first quartile, median, third quartile, and maximum.
(caption on next page) early postnatal weeks because speech production specific to native language only progresses between 6 months and the first year of the postnatal year (Dubois et al., 2016). Therefore, a future goal would be to perform a longitudinal study among CHD neonates, which in turn might help us understand whether the abnormalities associated with WM regions seen in early postnatal weeks might be associated with neurodevelopmental deficits including language development later during childhood and early adulthood. This might allow early interventions among neonates with CHD to help improve the cognitive outcomes associated with these WM regions because plasticity of the brain is very high during the early phase of development. Another limitation is that our groups were not gender-balanced, which may also contribute to differences between TDC and CHD populations.

Conclusion
Our results demonstrate the feasibility of using advanced multiple bvalue diffusion imaging to characterize regional white matter maturation in neonatal brains and to detect micro-structural abnormalities associated with CHD. Our analysis of 22 different WM regions (both hemispheres) in TDC neonates showed significant changes in all diffusion metrics but MD and Viso in 12 different regions (i.e., FA, Vic, ODI, RTOP and RTAP changed with age). Thus, in this cross-sectional study, these 5 measures were sensitive to microstructural changes with age. Our analysis also showed statistically significant group differences between CHD neonates and TDC, with RTAP and RTOP providing the strongest statistical significance in terms of effect sizes and the number of abnormal regions. Furthermore, a major finding from our study indicates the presence of WM differences (lower axonal and cellular packing density and volume) in bilateral CC and UF and left IFOF and left SLF in CHD. To the best of our knowledge, this is the first study to use advanced diffusion measures to better characterize regional WM differences in congenital heart disease in a systems based approach. The left lateralization of these differences in CHD suggests that the delayed language development often observed in CHD may be connected to the abnormalities seen in our study in the early postnatal weeks. Further, this study also compared different models used in the diffusion MRI field and found that the GMM and DTI measures are more sensitive at detecting differences in tissue structure than the NODDI model, as measured in terms of the number of regions found to be abnormal. Long-term follow-up of the CHD neonates is planned to assess their cognitive outcomes especially in the domain of language to determine the predictive power of these differences. Moreover, our work may aid in the monitoring of future early interventions targeting functions related to these specific WM regions.

Acknowledgements
This project was funded by the National Institute of Health (NIH) from grants R01MH097979 (PI Rathi), R01HD076258 (PI Grant) and R21HD058725 (PI Grant) as well as partially funded by Abbott Nutrition through the Center for Nutrition, Learning, and Memory at  Table 7 Cohen's d score reflecting the effect size magnitude. |d| < 0.2 "negligible", |d| < 0.5 "small", |d| < 0.8 "medium", |d| > 0.8 "large". Blue color indicates effect size for statistically significant ROIs. the University of Illinois (the trial being registered at clinicaltrials.gov as NCT02058225).Informed consent This study was approved by the Boston Children's Hospital IRB and informed consent was taken from all subjects' parents or legal guardians participating in this study.