A Novel Method for High-Dimensional Anatomical Mapping of Extra-Axial Cerebrospinal Fluid: Application to the Infant Brain

Cerebrospinal fluid (CSF) plays an essential role in early postnatal brain development. Extra-axial CSF (EA-CSF) volume, which is characterized by CSF in the subarachnoid space surrounding the brain, is a promising marker in the early detection of young children at risk for neurodevelopmental disorders. Previous studies have focused on global EA-CSF volume across the entire dorsal extent of the brain, and not regionally-specific EA-CSF measurements, because no tools were previously available for extracting local EA-CSF measures suitable for localized cortical surface analysis. In this paper, we propose a novel framework for the localized, cortical surface-based analysis of EA-CSF. The proposed processing framework combines probabilistic brain tissue segmentation, cortical surface reconstruction, and streamline-based local EA-CSF quantification. The quantitative analysis of local EA-CSF was applied to a dataset of typically developing infants with longitudinal MRI scans from 6 to 24 months of age. There was a high degree of consistency in the spatial patterns of local EA-CSF across age using the proposed methods. Statistical analysis of local EA-CSF revealed several novel findings: several regions of the cerebral cortex showed reductions in EA-CSF from 6 to 24 months of age, and specific regions showed higher local EA-CSF in males compared to females. These age-, sex-, and anatomically-specific patterns of local EA-CSF would not have been observed if only a global EA-CSF measure were utilized. The proposed methods are integrated into a freely available, open-source, cross-platform, user-friendly software tool, allowing neuroimaging labs to quantify local extra-axial CSF in their neuroimaging studies to investigate its role in typical and atypical brain development.

1. Introduction

1.1. General Information on CSF

Cerebrospinal fluid (CSF) is a clear, colorless fluid that circulates in the brain, providing necessary mechanical and immunological protection to the brain. In addition to its protective purpose, recent findings have shown that CSF circulation plays a crucial role in brain development and function prenatally and during the lifespan (Jessen et al., 2015; Lun et al., 2015). Figure 1 illustrates CSF circulation in the subarachnoid space around the brain, spinal cord, and in the ventricles of the brain. Following CSF production by the choroid plexus in the ventricles, it circulates from the lateral, third, and fourth ventricles to the cisterns of the brain. CSF flow continues to the subarachnoid space, where it covers the cortical convexities of the brain. CSF then flows back into the parenchyma, where it interacts with the interstitial fluid within the perivascular spaces. Finally, CSF returns to the subarachnoid space, where it is absorbed through meningeal lymphatic vessels and arachnoid granulations.

FIGURE 1

Figure 1. Illustration of cerebrospinal fluid circulation in the ventricles, subarachnoid space surrounding the brain, interstitial space within the brain parenchyma, and draining into the meningeal lyphatic system and arachnoid granulations. Figure adapted with permission from Shen (2018).

It is now realized that healthy CSF circulation serves two essential functions to the brain. The first is a regulatory function by delivering growth factors and signaling molecules critical to brain development (Lun et al., 2015). Second, CSF circulation provides a filtration mechanism by removing neuroinflammatory proteins and metabolic waste byproducts of neuronal function, which would otherwise accumulate (Iliff et al., 2012). Disrupted CSF circulation has been shown to be involved in neurodegenerative conditions such as Alzheimer's disease, ischemic and traumatic brain injury, and neuroinflammatory conditions such as multiple sclerosis (Simon and Iliff, 2016). More recently, abnormalities in CSF were also linked to the onset of neurodevelopmental disorders (NDDs), including autism spectrum disorder (ASD) (Shen et al., 2013, 2017, 2018).

1.2. MRI-Based CSF Biomarkers

The volume of different CSF compartments can be accurately measured from in vivo structural magnetic resonance imaging (structural MRI or sMRI), which can serve as indirect markers of abnormal CSF production and absorption. Current findings in lateral ventricles (LV) volume related to ASD have indicated no consistent, significant group differences in children (Vidal et al., 2008) or adults (McAlonan et al., 2002). In contrast, there is evidence in ASD for increased volume of CSF located outside the ventricles (Hallahan et al., 2009), as well as increased volume of global CSF across the entire brain (McAlonan et al., 2004). Notably, two studies of infants at high familial risk for ASD have reported increased global volume of CSF in the subarachnoid space (Extra-Axial CSF or EA-CSF) (Shen et al., 2013, 2017). Increased EA-CSF volume at 6 months (see example in Figure 2) of age, prior to the defining behavior symptoms of ASD, was observed in infants who were later diagnosed with ASD (Shen et al., 2013). Further, EA-CSF remained abnormally elevated at 12 and 24 months of age (Shen et al., 2013). Moreover, greater EA-CSF volume at 6 months was also associated with more severe autism symptoms at the time of diagnosis at 2 years of age (Shen et al., 2013). Such EA-CSF findings were later confirmed through replication in a larger, independent cohort of infants (Shen et al., 2017). The previous EA-CSF studies relied on a novel method in infant MRIs to quantify the volume of EA-CSF in the dorsal subarachnoid space above the horizontal plane of the anterior-posterior commissure, thereby avoiding ventral regions that contain cisterns, sinuses, and vasculature that should not be classified as EA-CSF.

FIGURE 2

Figure 2. (A) T1-weighted MRI of a typically-developing infant with a normal MRI at 6 months of age. (B) T1-weighted MRI of an infant with enlarged EA-CSF volume at 6 months of age, who was diagnosed with ASD at 2 years of age. The dark regions between the brain folds and skull indicate increased volume of EA-CSF.

1.3. EA-CSF Quantification

The earlier results indicate that the quantification of global EA-CSF could be important in understanding the nature of brain and CSF pathology and its relation to ASD symptoms. However, the global EA-CSF measure does not provide an anatomical localization of the effect. A localized EA-CSF extraction would provide measurements suitable for localized group analysis or localized discriminative analysis. Moreover, the ability to obtain anatomically precise measures would provide greater insight into the underlying physiological and anatomical mechanisms, as well as a more fine-grained ability to detect abnormalities in NDDs. One older method to extract local measurements of EA-CSF is through voxel-based morphometry (VBM), which utilizes a statistical approach of parametric mapping (Ashburner and Friston, 2000). VBM methods are computationally efficient as they do not involve surface reconstruction of complex cortical surfaces, but the accuracy and precision of localized EA-CSF measurements from VBM are severely limited by the voxel resolution and are sensitive to volumetric registration errors, which are known to be abundantly present in most cortical areas due to the inherent cortical folding variability. Moreover, such voxel-based EA-CSF measurements cannot be easily correlated with other cortical surface-based measurements (e.g., cortical thickness and surface area) that have been shown to hold value as early biomarkers for NDDs such as ASD (Hazlett et al., 2017). Hence, the ability to extract high-dimensional surface-based local EA-CSF measurements would allow for a better understanding of how such biomarkers are related to each other, leading to optimal combinations for accurate early prediction of NDDs using deep learning techniques of multiple measures.

1.4. Proposed Method of Local EA-CSF

To overcome the limitations mentioned above, as shown in Figure 3, this paper presents a novel framework for extracting surface-based local EA-CSF measurements from sMRI. The proposed framework first computes a probabilistic tissue segmentation of white matter (WM), gray matter (GM), and CSF. A hard segmentation is obtained from the tissue probability maps. These hard segmentations are then used to reconstruct polyhedral models of the outer CSF hull surface as well as the WM and GM surfaces. A Laplacian partial differential equation (PDE) is solved between a defined inner surface and the CSF hull surfaces to generate a vector field that is used to create streamlines connecting the surfaces. Along these streamlines, the CSF space is sampled, and CSF probability values are integrated to generate local EA-CSF measures at each point cortical surface. To the best of our knowledge, the proposed framework is the first to address the problem of extracting local EA-CSF measurements in a way that is suitable for localized surface-based analysis.

FIGURE 3

Figure 3. The proposed framework for the extraction of local EA-CSF from structural MRI. The proposed processing pipeline combines probabilistic brain tissue segmentation, cortical surface reconstruction, as well as streamline-based quantification to produce accurate and reliable local EA-CSF measurements.

2. Methods

2.1. Participants and MRI Acquisition

We analyzed 153 structural MR brain images, which were obtained at 3 time points (6, 12, 24 months of age) from 51 typically-developing infants, who were at low risk for ASD (i.e., no family history of ASD, intellectual disability, or major psychiatric disorder, and who had an older sibling with typical development). These infants were scanned longitudinally at 6, 12, and 24 months of age as part of a National Institutes of Health-funded, multi-site, Autism Centers of Excellence (ACE) Network study: the Infant Brain Imaging Study (IBIS). Data collection sites had study protocols approval from their Institutional Review Boards (IRB), and all enrolled subjects had informed consent provided by parent/guardian. The MRI scans were acquired at 4 different sites (University of North Carolina at Chapel Hill, University of Washington in Seattle, Washington University in St. Louis, and Children's Hospital of Philadelphia), each equipped with 3T Siemens Tim Trio scanners (Wolff et al., 2012). The scan sessions included T1-weighted (T1w) (160 sagittal slices with TR = 2,400 ms, TE = 3.16 ms, flip angle = 8°, field of view 224 × 256) and T2-weighted (T2w) (160 sagittal slices with TR = 3,200 ms, TE = 499 ms, flip angle = 120°, field of view 256 × 256) MRI scans. All datasets had the same spatial resolution of 1 × 1 × 1 mm³. Only participants with scans from all three time points were included to enable an accurate longitudinal study of extracted local EA-CSF trajectories. The total sample of N = 51 infants (153 scans) included N = 32 males (96 scans) and 19 females (57 scans); see Table 1 for complete demographic information.

TABLE 1

Table 1. Breakdown of ages, sex, and number of scans in the current dataset of typically developing infants from the IBIS sample.

Multiple procedures for quality control were employed to assess scanner stability and reliability across sites, times, and scanner upgrades. A Lego (Lego Group, Billund, Denmark) brick-based phantom (Fonov et al., 2011) was scanned monthly at each study site and analyzed to assess image quality and quantitatively address site-specific local distortions. Also, two adult living phantoms were scanned once per year at each scanner and after any significant scanner update. The data for these phantoms were evaluated for scanner stability across sites and time (Gouttard et al., 2008) and are also to assess stability for the local EA-CSF measure.

2.2. Image Processing and Surface Generation

2.2.1. Initial Preprocessing

The raw T1w and T2w brain images were corrected for intensity non-uniformity using the N4 algorithm (Tustison et al., 2010) (Figures 4A,B). Correction of geometric distortions was also applied for the optimal processing of multi-site longitudinal data (Fonov et al., 2010). T1w and T2w images were rigidly transformed to a prior pediatric 1-year-old atlas in stereotaxic space. A prior intensity growth map was applied to the 12-month T1w and T2w scans to improve the poor contrast of the WM/GM boundary from under-myelination (Kim et al., 2013).

FIGURE 4

Figure 4. An example of the probabilistic tissue segmentation obtained for an input infant sMRI scan. (A) T1-weighted scan, (B) T2-weighted scan, (C) WM, GM, and CSF label map, (D) CSF probability map, (E) CSF probability map with ventricular CSF space removed.

2.2.2. Skull Stripping

The brain mask necessary to perform skull stripping was performed using a multi-atlas approach that combines multiple candidate brain masks obtained via deformable registration of a prior set of atlases (each consisting of a T1w, T2w and brain mask label image). Five brain masks were utilized, namely FSL-BET (Smith, 2002), two in-house prior atlases, and two atlases of the CIVET pipeline (Kim et al., 2005). The deformable registration was computed via the ANTs registration toolkit (Avants et al., 2011) using both T1w and T2w data. The candidate brain masks were then combined, weighting each equally, to result in a majority vote fusion.

2.2.3. Tissue Segmentation

Different brain tissues were then segmented using a framework of atlas-moderated expectation-maximization (EM) implemented in the AutoSeg toolkit (Wang et al., 2014). Particularly, a deformable registration was applied to propagate a prior template and prior tissue probability maps for WM, GM, and CSF from MNI space into individual T1w data (Fonov et al., 2011). Then, an EM based tissue segmentation was performed (Van Leemput et al., 1999) to obtain a label map with segmentations for WM, GM, and CSF (Figure 4C). Ventricular CSF space (lateral ventricles, third and fourth ventricles) was then removed from the CSF posterior map (Figure 4D) by deformably co-registering a single prior template with an existing ventricular area mask and using the registered mask to remove the ventricle (Figure 4E).

2.2.4. Surface Reconstruction

Cortical surfaces were reconstructed with an adapted version of the CIVET workflow (Kim et al., 2005). The cortical surface model consisted of high-resolution triangle meshes (81,920 triangles and 40,962 vertices) in each hemisphere, and cortical surface correspondence among subjects was established via spherical registration to an average surface template (Robbins et al., 2003). CIVET was applied to 12 and 24 months old data following tissue segmentation with AutoSeg to construct WM and GM surfaces (Figures 5A,B). However, tissue segmentation for 6-month-old subjects did not yield reliable WM vs. GM segmentation because the WM and GM have almost the same intensity level in both T1w and T2w scans of isointense-phase infants (around 6–8 months of age). Cortical surfaces at 6 months were determined longitudinally via the corresponding 12 months visit to solve this problem. Using ANTs (Avants et al., 2011) deformable registration with normalized cross-correlation (metric radius 2 mm, Gaussian smoothing of 3 mm of the deformation map) of joint T1w and T2w data (both image sources were equally weighted), the pre-processed, brain masked MRI data of 12-month-old subjects was registered to data from the same subject at age 6 months. This registration was applied to the cortical surfaces of the 12-month-old subjects to propagate them into the 6 months space (see Figure 6). In addition to the WM and GM surfaces, a smoothed middle surface (Figure 5C) was obtained by averaging WM and GM surfaces and then two iterations of averaging based surface smoothing. Moreover, the outer CSF hull surface (Figure 5D) was generated by first dilating the intracranial mask, followed by a surface reconstruction using standard marching cubes algorithm (Lorensen and Cline, 1987) and a subsequent Laplacian surface smoothing. Finally, all the reconstructed surfaces were visually QC'ed with a surface cut overlay on the MR images by a single rater (MM).

FIGURE 5

Figure 5. An example of the cortical surfaces reconstructed from the input infant structural MRI. (A) WM Surface, (B) GM Surface, (C) Middle Surface, and (D) CSF Hull Surface.

FIGURE 6

Figure 6. The process of generating cortical surfaces for the 6-month-old structural MRI scans.

2.3. Extraction of Local Extra-Axial Cerebrospinal Fluid

2.3.1. Solving Laplace PDE

Following the reconstruction of the cortical surfaces, the next step is to solve a Laplace's equation between the inner surface (S_inner) and a corresponding outer surface (S_outer). While the CSF hull surface (Figure 5D) is defined as S_outer, in this work, we utilize the middle surface (Figure 5C) as S_inner instead of the GM surface to accommodate for potential surface reconstruction errors. Both S_inner and S_outer are assumed to have spherical topology, i.e., can be stretched and warped without breaking to form a sphere surface. Laplace's equation is a second-order PDE solved for a scalar field u(x) that is enclosed between boundaries S_inner and S_outer. The Laplace PDE takes the following form

$\begin{array}{ll}△u={\nabla }^{2}u(x)=0,& (1)\end{array}$

where u(x) = u_L for x ∈ S_inner and u(x) = u_H for x ∈ S_outer. To correctly measure local EA-CSF, a boundary condition map that defines the Laplace PDE boundary condition must be defined in an anatomically consistent manner. In the solution of the Laplace PDE in the proposed framework, the solution domain is bounded by the Dirichlet condition and the Neumann condition. The Dirichlet condition specifies values of the solution itself on the boundary, while the Neumann boundary condition defines values for the first-order derivative of the solution. The interface with the Dirichlet condition defines S_inner and S_outer where streamlines start and arrive, and the Neumann condition defines an open boundary that is parallel to the streamlines (see Figure 7). A consistent boundary map generation is ensured using surface-based pre-processing steps that were applied to create a boundary label map in the image domain (Lee et al., 2016). The Laplace PDE is iteratively solved in the created image voxel grid using the Jacobi method (Causon and Mingham, 2010).

FIGURE 7

Figure 7. A solution of the (A) Laplace equation with two different boundary conditions. As observed in the color map (B) and the contour image (C), the solution isolines are parallel to the Dirichlet boundary and perpendicular to the Neumann boundary. Figure from Lee et al. (2016).

2.3.2. Streamline-Based Local EA-CSF

After obtaining the Laplace PDE solution, the next step is the computation of the local EA-CSF between S_inner and S_outer, which we define as CSF accumulated along the lines connecting the two surfaces. Such lines need to be orthogonal to the PDE solution isolines at each point to obtain a biologically plausible path. To achieve that, streamlines that are tangent to the normalized gradient field of the PDE solution are utilized to provide an analogy to cortical columns and to establish a one-to-one correspondence between S_inner and S_outer. Such streamlines are then constructed explicitly by the integration of the Lagrangian vector field. A fourth-order Runge-Kutta (RK4) integration method (Yaakub and Evans, 1999) is used in generating the streamlines to minimize local truncation error and provide faster convergence.

To achieve sub-voxel accuracy, we process the starting and ending segments of the streamlines to fit them perfectly within the boundary of the defined inner and outer surfaces. Finally, a local EA-CSF measure is computed for each vertex v by accumulating CSF probability P at each point k on the streamline l_v associated with v. A linear approximation is utilized to account for the streamlines non-uniformity

$\begin{array}{ll}{\text{EA-CSF}}_{\text{v}}={\sum }_{\text{k}\in {\text{l}}_{\text{v}}}\left(\frac{\text{P}(\text{k})+\text{P}(\text{k}+1)}{2}\right)\times {\text{Δ}}_{\text{k}},& (2)\end{array}$

where Δ_k representing the Euclidean distance between point k and the successive point k+1. Figure 8 provides an example for the computed streamlines and sampled CSF probability map.

FIGURE 8

Figure 8. (A) Streamlines generated using a fourth-order Runge-Kutta (RK4) integration method. (B) Local EA-CSF measure is computed by accumulating CSF along the generated streamlines.

As CSF values are accumulated along the streamlines, and locations/voxels with no CSF contents do not contribute to the accumulate EA-CSF measure, our proposed streamline approach is rather robust to the exact location of the inner surface as long as all CSF regions of interest lay between the inner and outer surfaces. For that purpose we chose the middle cortical surface as the inner surface (see section 2.2). Choosing the gray matter surface as inner surface would likely exclude minor EA-CSF sections due to segmentations errors of the gray matter boundary.

2.3.3. Statistical Analysis

The extracted raw local EA-CSF maps were first mapped to the common MNI surface template space (Lyttelton et al., 2007) for additional processing and analysis. As an initial standard processing step in the study of cortical surface measurements, a geodesic heat kernel-based smoothing (FWHM of 20 mm) was applied to the EA-CSF maps (Chung et al., 2005). The effects of sex and age on local EA-CSF were tested using a longitudinal mixed-effects model in SurfStat, which is a toolbox for statistical analysis of cortical surface measurements applying random field theory for statistical inference (Chung et al., 2010). The longitudinal linear mixed model included a subject-specific random intercept to induce equal correlations between observations on the same subject. Slope terms were also added to model the fixed effects of sex, of age, as well as of sex and age interactions. In particular, with the local EA-CSF as the dependent variable Y, the following linear mixed model was fitted for each subject i:

$\begin{array}{ll}{Y}_i={\beta }_0+{\beta }_{\text{1}}{\text{Sex}}_i+{\beta }_{\text{2}}{\text{Age}}_i+{\beta }_{\text{3}}{\text{Sex}}_i{\text{Age}}_i+U_i+{\epsilon }_i,& (3)\end{array}$

where U_i captures estimates for the subject-specific random effect and ε_i is the independent noise term in every observation. Standard false discovery rate (FDR) (Benjamini and Hochberg, 1995) correction was applied to correct for the multiple comparisons in the model in Equation ( 3). Figure 9 illustrates the within-subject correlation of local EA-CSF across age, as revealed by the linear mixed model. High correlations are shown across most of the brain regions, particularly in the frontal, parietal, and temporal lobes. In the presence of highly correlated areas, subject-specific random effects need to be incorporated in the linear mixed model.

FIGURE 9

Figure 9. The correlation of local EA-CSF within-subject across ages 6–24 months. High correlations in the frontal, parietal, and temporal lobes were observed.

3. Experimental Results

3.1. Local EA-CSF Reproducibility

The stability and reliability of the proposed local EA-CSF measure were tested using a dataset with a large set of scan/rescan MRIs. Two human phantoms (young male subjects, age 26, and 27), were scanned with the same pulse sequence at four different sites using Siemens 3T Tim Trio scanners, at irregular intervals over 2.5 years. This resulted in 35 MRI scans for subject I and 31 MRI scans for subject II. The tissue segmentation, brain surface reconstruction, and computation of local EA-CSF maps were performed independently. Local EA-CSF maps were then analyzed using the local coefficient of variation (CV) as a measure of stability. The CV for a vertex v is defined as the ratio between the standard deviation and mean of the extracted local EA-CSF across diffident scans of the same subject:

$\begin{array}{ll}C{V}_v=\frac{{\sigma }_v}{{\mu }_v}\times 100.& (4)\end{array}$

It is noteworthy that the EA-CSF spaces in these two adult subjects were visually smaller than those we observed in our infant dataset. In general, in our experience infants have visually larger EA-CSF spaces than adults. It is thus likely our results yield a conservative estimate of the expected reproducibility in the infant settings.

The CV analysis showed excellent stability with mean across-site CV of 1.15 and 1.56% for all cortical regions in Case I and Case II respectively (Figure 10). Higher CVs were observed in few regions, including left supramarginal gyrus, left postcentral gyrus, left gyrus rectus, right postcentral gyrus, right superior temporal gyrus, and right precentral gyrus. Local EA-CSF variability in these regions was mainly linked to imperfect CSF tissue segmentation. It is worthwhile to mention that the global ICV measures showed CV values around 1% (Bryson et al., 2008; Hazlett et al., 2012) in the same human datasets. Hence, the proposed local EA-CSF extraction framework provided a stable local measure as compared to global ICV in adult brains.

FIGURE 10

Figure 10. Coefficients of variation (σ_v/μ_v× 100) for local EA-CSF maps of two sets of adult brains. Mean coefficient of variation of 1.15 and 1.56% were observed for (A) Case I and (B) Case II, respectively. Regions with high coefficient of variation are linked to CSF segmentation issues.

3.2. Anatomical Mapping of Local EA-CSF Across Age

Figure 11 shows the mean and standard deviation of local EA-CSF maps at 6, 12, and 24 months. At all ages, local EA-CSF is most abundant over the central and precentral sulci. This observation is consistent with visual inspection of hundreds of infant MRI brain images (Shen et al., 2013, 2017). We also observed finer grained findings that are consistent at all the studied ages, such as increased EA-CSF along the superior frontal gyrus, fusiform gyrus, and the calcarine sulcus, as well as decreased EA-CSF in the inferior parietal lobule and the middle temporal gyrus. These observations are novel. As a group, this sample of typically developing infants showed anatomical consistency of the average, local EA-CSF patterns across time (left column in Figure 11). Yet, differences across individuals at each age are quite large as evident from the mean coefficients of variation of 44.6, 42.2, and 42.9% at 6, 12, and 24 months, respectively.

FIGURE 11

Figure 11. The mean and standard deviation of local EA-CSF measured at 6 months (first row), 12 months (second row), and 24 months (third row). At all ages, local EA-CSF is most abundant over the central sulcus, precentral sulcus and longitudinal fissure. The highest inter-subject variability of local EA-CSF was found in the medial and ventral temporal areas.

Consistent with our previous global EA-CSF report (Shen et al., 2017), we observed a significant decrease in local EA-CSF over time. Figures 12A,D shows the linear change rates of local EA-CSF per month between 6 and 24 months of age. The mean negative change rate of all regions that showed statistically significant differences across age was –0.011 per month. Multiple regions of the frontal lobe showed significant negative change rates of EA-CSF, while the region with the greatest negative change rate (–0.046 per month) was observed in the left superior temporal gyrus. Figures 12B,C shows the FDR-corrected T scores (FDR threshold Q < 0.01) for the age effect resulting from the longitudinal mixed-effects model. Overall 16.9% of the brain showed a significant decrease in local EA-CSF in the first 2 years of infant's life. Local EA-CSF showed a highly negative correlation with age in most of the frontal lobe areas, including bilateral middle frontal gyrus and bilateral superior frontal gyrus. Negative correlation with age was also observed in some temporal regions, such as the bilateral superior temporal gyrus. No cortical areas showed a significantly positive association of age and the local EA-CSF.

FIGURE 12

Figure 12. The longitudinal age effect on local EA-CSF. (A) Frontal regions showed a mean negative change rate of local EA-CSF of –0.011. The greatest negative change rate of –0.046 per month was observed in the left superior temporal gyrus. The corresponding (B) raw T-scores (red: negative age effect; blue: positive age effect) and (C) FDR corrected T-scores in regions with statistically significant age effect (FDR Q < 0.01, only regions with a negative age effect survived FDR correction); (D) EA-CSF plot across age at the surface location with maximum T score, (T_max = 8.57, located in left frontal superior gyrus).

3.3. Sex Differences in the Anatomical Mapping of Local EA-CSF

Several brain regions showed a statistically significant main effect of sex across age (FDR threshold Q < 0.01), with higher local EA-CSF in males compared to females in regions such as: right middle frontal gyrus, inferior and medial orbital gyri, bilateral inferior frontal gyri, right insula, and right superior temporal gyrus, right supplementary motor cortex, and right superior frontal gyrus. 6.2% of the brain showed more local EA-CSF in males compared to females. There were no regions where females showed larger EA-CSF (See Figures 13A,B). It is noteworthy that our longitudinal mixed model did not covary for head size, for example, via intracranial cavity volume.

FIGURE 13

Figure 13. Regions with sex differences in local EA-CSF. (A) FDR corrected T-scores in regions with statistically significant sex effect (FDR Q < 0.01, red: male > female; blue: male < female). (B) Sex differences at location with maximum T score, (T_max = 6.03, located in the right superior frontal gyrus). (C) Regions with statistically significant sex by age interaction (FDR Q < 0.01). (D) Local EA-CSF profile at location with maximum sex by age effect (located in the left superior frontal gyrus).

In addition to the main sex effect, brain regions showed a significant sex by age interaction, whereby local EA-CSF decreased at a faster rate in males compared to females, including: bilateral calcarine fissure, bilateral parahippocampal gyri, bilateral middle temporal gyri, right superior parietal gyrus, right precuneus. Figures 13C,D shows the results of the longitudinal analysis of sex and age interaction.

4. Summary and Concluding Discussion

In this paper we propose a novel framework for extracting surface-based local EA-CSF measurements from MR brain images. The proposed framework is the first to address the problem of obtaining local EA-CSF measurements in a way that is suitable for localized surface-based analysis. The proposed processing relies on a probabilistic tissue segmentation approach to generate a CSF probability map that is used to reconstruct the outer CSF hull surface. A Laplacian partial differential equation is solved between the inner cortical surface and the CSF hull surfaces to generate a vector field that is used to create streamlines connecting the surfaces at sub-voxel accuracy via a fourth-order Runge-Kutta approach. The starting and ending segments of the streamlines are then processed to fit them correctly within the boundary of the defined inner and outer surfaces. Along these streamlines, the CSF probability values are accumulated to quantify EA-CSF measures at each vertex on the cortical surface mesh. The proposed local EA-CSF extraction tool was used to study early postnatal brain development in typically developing infants from the IBIS dataset. Due to a high within-subject correlation, a longitudinal linear mixed model is proposed to incorporate fixed effects of age and sex as well as sex and age interactions.

The experimental results obtained from the scan/rescan human dataset show that the proposed local EA-CSF measure is reliable in a scan/rescan setting and produces reasonably stable results. The stability of the proposed method is further confirmed by the consistency in local EA-CSF patterns across the 3-time points used in studying local EA-CSF trajectories in the first 2 years of infancy. The experimental results reveals several findings through the proposed processing and analysis pipeline. First, local EA-CSF in several cortical regions, mainly in the frontal lobe areas, shows a statistically significant negative correlation with age. The longitudinal analysis also reveals several cortical regions with statistically significant higher local EA-CSF in males compared to females. Most of these regions also show a more substantial decrease in local EA-CSF across age. However, few cortical areas show higher negative local EA-CSF change rate in male subjects compared to females. Such localized findings confirm that the proposed local EA-CSF extraction pipeline reveals specific regions of significant change that would not be possible to be observed using the previous global EA-CSF approach.

The quantification of local EA-CSF measurements relies on streamlines that are generated based on solving an isotropic Laplace PDE between inner and outer surfaces on a created voxel grid. However, such isotropic PDE solution is purely based on the boundaries implied by the inner and outer surfaces. Hence, a limitation of the proposed framework is that the generated streamlines are not constrained to areas containing CSF. The method also fails to account for partial volume effects that lie between the two boundaries. In the future, we plan to improve the proposed local EA-CSF quantification by replacing the isotropic Laplace PDE with an anisotropic version (Joshi et al., 2018) in which the diffusion coefficient varies spatially in proportion to the fraction of CSF in each voxel. This allows for generating streamlines that follow a realistic path through areas containing CSF, leading to improvements in the accuracy and reliability of extracted local EA-CSF measurements.

Data Availability Statement

The source code and scripts of our method will be released as open source with the next version of the Auto-EACSF tool on github (ongoing development). All tools developed at the UNC Neuro Image Research and Analysis Laboratory are open source. All raw MRI datasets and associated demographic information employed in this study is available at: NIH/NDA: https://nda.nih.gov/edit_collection.html?id=19.

Ethics Statement

The studies involving human participants were reviewed and approved by the Institutional Review Boards (IRB) of all participating institutions. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author Contributions

Method development and computational analysis was mainly performed by MM under guidance by MAS and MDS, with additional guidance and methodological input by GG and SP. Parts of processing was performed by SK and JG. Study design, data acquisition was performed by AE, SD, AE. RM, KB, RS, HH, and JP. Interpretation of results was performed by MM, MAS, and MDS. Document preparation was performed by all co-authors.

Funding

This study was supported by grants from the National Institutes of Health (R01-HD055741, T32-HD040127, U54-HD079124, U54-HD086984, and R01-EB021391), Autism Speaks, and the Simons Foundation (#140209). MDS is supported by a U.S. National Institutes of Health (NIH) career development award (K12-HD001441). The sponsors had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are sincerely grateful to all the families and children who have participated in the Infant Brain Imaging Study (IBIS). The Infant Brain Imaging Study (IBIS) Network is an NIH funded Autism Centers of Excellence project and consists of a consortium of 9 universities in the U.S. and Canada. Members and components of the IBIS Network include: J. Piven (IBIS Network PI), Clinical Sites: University of North Carolina: HH, C. Chappell, M.D. Shen, M. Swanson; University of Washington: SD, AE, D. Shaw, T. St. John; Washington University: KB, J. Constantino; Children's Hospital of Philadelphia: RS, J. Pandey. Behavior Core: University of Washington: AE; University of Alberta: L. Zwaigenbaum; University of Minnesota: J. Elison, J. Wolff. Imaging Core: University of North Carolina: MS; New York University: GG; Washington University in St. Louis: RM, J. Pruett. Data Coordinating Center: Montreal Neurological Institute: AE, D.L. Collins, V. Fonov, L. MacIntyre; S. Das. Statistical Analysis Core: K. Truong. Environmental risk core: John Hopkins University: H. Volk. Genetics Core: John Hopkins University: D. Fallin; University of North Carolina: MS.

References

Ashburner, J., and Friston, K. J. (2000). Voxel-based morphometry–the methods. Neuroimage 11, 805–821. doi: 10.1006/nimg.2000.0582

PubMed Abstract | CrossRef Full Text | Google Scholar

Avants, B. B., Tustison, N. J., Song, G., Cook, P. A., Klein, A., and Gee, J. C. (2011). A reproducible evaluation of ants similarity metric performance in brain image registration. Neuroimage 54, 2033–2044. doi: 10.1016/j.neuroimage.2010.09.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B Methodol. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Bryson, S. E., Zwaigenbaum, L., McDermott, C., Rombough, V., and Brian, J. (2008). The autism observation scale for infants: scale development and reliability data. J. Autism Dev. Disord. 38, 731–738. doi: 10.1007/s10803-007-0440-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Causon, D., and Mingham, C. (2010). Introductory Finite Difference Methods for PDEs. London, UK: Bookboon.

Google Scholar

Chung, M. K., Robbins, S. M., Dalton, K. M., Davidson, R. J., Alexander, A. L., and Evans, A. C. (2005). Cortical thickness analysis in autism with heat kernel smoothing. Neuroimage 25, 1256–1265. doi: 10.1016/j.neuroimage.2004.12.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, M. K., Worsley, K. J., Nacewicz, B. M., Dalton, K. M., and Davidson, R. J. (2010). General multivariate linear modeling of surface shapes using surfstat. Neuroimage 53, 491–505. doi: 10.1016/j.neuroimage.2010.06.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonov, V., Evans, A. C., Botteron, K., Almli, C. R., McKinstry, R. C., Collins, D. L., et al. (2011). Unbiased average age-appropriate atlases for pediatric studies. Neuroimage 54, 313–327. doi: 10.1016/j.neuroimage.2010.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonov, V. S., Janke, A., Caramanos, Z., Arnold, D. L., Narayanan, S., Pike, G. B., et al. (2010). “Improved precision in the measurement of longitudinal global and regional volumetric changes via a novel MRI gradient distortion characterization and correction technique,” in International Workshop on Medical Imaging and Virtual Reality (Beijing: Springer), 324–333.

Google Scholar

Gouttard, S., Styner, M., Prastawa, M., Piven, J., and Gerig, G. (2008). “Assessment of reliability of multi-site neuroimaging via traveling phantom study,” in International Conference on Medical Image Computing and Computer-Assisted Intervention (New York, NY: Springer), 263–270.

PubMed Abstract | Google Scholar

Hallahan, B., Daly, E., McAlonan, G., Loth, E., Toal, F., O'brien, F., et al. (2009). Brain morphometry volume in autistic spectrum disorder: a magnetic resonance imaging study of adults. Psychol. Med. 39, 337–346. doi: 10.1017/S0033291708003383

PubMed Abstract | CrossRef Full Text | Google Scholar

Hazlett, H. C., Gu, H., McKinstry, R. C., Shaw, D. W., Botteron, K. N., Dager, S. R., et al. (2012). Brain volume findings in 6-month-old infants at high familial risk for autism. Am. J. Psychiatry 169, 601–608. doi: 10.1176/appi.ajp.2012.11091425

PubMed Abstract | CrossRef Full Text | Google Scholar

Hazlett, H. C., Gu, H., Munsell, B. C., Kim, S. H., Styner, M., Wolff, J. J., et al. (2017). Early brain development in infants at high risk for autism spectrum disorder. Nature 542:348. doi: 10.1038/nature21369

PubMed Abstract | CrossRef Full Text | Google Scholar

Iliff, J. J., Wang, M., Liao, Y., Plogg, B. A., Peng, W., Gundersen, G. A., et al. (2012). A paravascular pathway facilitates csf flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid β. Sci. Transl. Med. 4:147ra111. doi: 10.1126/scitranslmed.3003748

PubMed Abstract | CrossRef Full Text | Google Scholar

Jessen, N. A., Munk, A. S. F., Lundgaard, I., and Nedergaard, M. (2015). The glymphatic system: a beginner's guide. Neurochem. Res. 40, 2583–2599. doi: 10.1007/s11064-015-1581-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Joshi, A. A., Bhushan, C., Salloum, R., Wisnowski, J. L., Shattuck, D. W., and Leahy, R. M. (2018). “Using the anisotropic laplace equation to compute cortical thickness,” in International Conference on Medical Image Computing and Computer-Assisted Intervention (Granada: Springer), 549–556.

PubMed Abstract | Google Scholar

Kim, J. S., Singh, V., Lee, J. K., Lerch, J., Ad-Dab'bagh, Y., MacDonald, D., et al. (2005). Automated 3-d extraction and evaluation of the inner and outer cortical surfaces using a laplacian map and partial volume effect classification. Neuroimage 27, 210–221. doi: 10.1016/j.neuroimage.2005.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. H., Fonov, V. S., Dietrich, C., Vachet, C., Hazlett, H. C., Smith, R. G., et al. (2013). Adaptive prior probability and spatial temporal intensity change estimation for segmentation of the one-year-old human brain. J. Neurosci. Methods 212, 43–55. doi: 10.1016/j.jneumeth.2012.09.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J., Kim, S. H., Oguz, I., and Styner, M. (2016). “Enhanced cortical thickness measurements for rodent brains via lagrangian-based rk4 streamline computation,” in Proceedings of SPIE–the International Society for Optical Engineering, Vol. 9784 (San Diego, CA: NIH Public Access).

PubMed Abstract | Google Scholar

Lorensen, W. E., and Cline, H. E. (1987). “Marching cubes: a high resolution 3d surface construction algorithm,” in ACM Siggraph Computer Graphics, Vol. 21 (Anaheim, CA: ACM), 163–169.

PubMed Abstract | Google Scholar

Lun, M. P., Monuki, E. S., and Lehtinen, M. K. (2015). Development and functions of the choroid plexus–cerebrospinal fluid system. Nat. Rev. Neurosci. 16, 445–457. doi: 10.1038/nrn3921

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyttelton, O., Boucher, M., Robbins, S., and Evans, A. (2007). An unbiased iterative group registration template for cortical surface analysis. Neuroimage 34, 1535–1544. doi: 10.1016/j.neuroimage.2006.10.041

PubMed Abstract | CrossRef Full Text | Google Scholar

McAlonan, G. M., Cheung, V., Cheung, C., Suckling, J., Lam, G. Y., Tai, K., et al. (2004). Mapping the brain in autism. a voxel-based MRI study of volumetric differences and intercorrelations in autism. Brain 128, 268–276. doi: 10.1093/brain/awh332

PubMed Abstract | CrossRef Full Text | Google Scholar

McAlonan, G. M., Daly, E., Kumari, V., Critchley, H. D., Amelsvoort, T. V., Suckling, J., et al. (2002). Brain anatomy and sensorimotor gating in asperger's syndrome. Brain 125, 1594–1606. doi: 10.1093/brain/awf150

PubMed Abstract | CrossRef Full Text | Google Scholar

Robbins, S., Evans, A. C., Collins, D. L., and Whitesides, S. (2003). “Tuning and comparing spatial normalization methods,” in International Conference on Medical Image Computing and Computer-Assisted Intervention (Montreal, QC: Springer), 910–917.

PubMed Abstract | Google Scholar

Shen, M. D. (2018). Cerebrospinal fluid and the early brain development of autism. J. Neurodev. Disord. 10:39. doi: 10.1186/s11689-018-9256-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, M. D., Kim, S. H., McKinstry, R. C., Gu, H., Hazlett, H. C., Nordahl, C. W., et al. (2017). Increased extra-axial cerebrospinal fluid in high-risk infants who later develop autism. Biol. Psychiatry 82, 186–193. doi: 10.1016/j.biopsych.2017.02.1095

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, M. D., Nordahl, C. W., Li, D. D., Lee, A., Angkustsiri, K., Emerson, R. W., et al. (2018). Extra-axial cerebrospinal fluid in high-risk and normal-risk children with autism aged 2–4 years: a case-control study. Lancet Psychiatry 5, 895–904. doi: 10.1016/S2215-0366(18)30294-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, M. D., Nordahl, C. W., Young, G. S., Wootton-Gorges, S. L., Lee, A., Liston, S. E., et al. (2013). Early brain enlargement and elevated extra-axial fluid in infants who develop autism spectrum disorder. Brain 136, 2825–2835. doi: 10.1093/brain/awt166

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, M. J., and Iliff, J. J. (2016). Regulation of cerebrospinal fluid (csf) flow in neurodegenerative, neurovascular and neuroinflammatory disease. Biochim. Biophys. Acta 1862, 442–451. doi: 10.1016/j.bbadis.2015.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S. M. (2002). Fast robust automated brain extraction. Hum. Brain Mapp. 17, 143–155. doi: 10.1002/hbm.10062

PubMed Abstract | CrossRef Full Text | Google Scholar

Tustison, N. J., Avants, B. B., Cook, P. A., Zheng, Y., Egan, A., Yushkevich, P. A., et al. (2010). N4itk: improved n3 bias correction. IEEE Trans. Med. Imaging 29, 1310–1320. doi: 10.1109/TMI.2010.2046908

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Leemput, K., Maes, F., Vandermeulen, D., and Suetens, P. (1999). Automated model-based tissue classification of mr images of the brain. IEEE Trans. Med. Imaging 18, 897–908. doi: 10.1109/42.811270

CrossRef Full Text | Google Scholar

Vidal, C. N., Nicolson, R., Boire, J.-Y., Barra, V., DeVito, T. J., Hayashi, K. M., et al. (2008). Three-dimensional mapping of the lateral ventricles in autism. Psychiatry Res. Neuroimaging 163, 106–115. doi: 10.1016/j.pscychresns.2007.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Vachet, C., Rumple, A., Gouttard, S., Ouziel, C., Perrot, E., et al. (2014). Multi-atlas segmentation of subcortical brain structures via the autoseg software pipeline. Front. Neuroinformatics 8:7. doi: 10.3389/fninf.2014.00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolff, J. J., Gu, H., Gerig, G., Elison, J. T., Styner, M., Gouttard, S., et al. (2012). Differences in white matter fiber tract development present from 6 to 24 months in infants with autism. Am. J. Psychiatry 169, 589–600. doi: 10.1176/appi.ajp.2011.11091447

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaakub, A., and Evans, D. J. (1999). A fourth order runge–kutta rk (4, 4) method with error control. Int. J. Comput. Math. 71, 383–411. doi: 10.1080/00207169908804817

CrossRef Full Text | Google Scholar