Modeling healthy male white matter and myelin development: 3 through 60months of age.

An emerging hypothesis in developmental and behavioral disorders is that they arise from disorganized brain messaging or reduced connectivity. Given the importance of myelin to efficient brain communication, characterization of myelin development in infancy and childhood may provide salient information related to early connectivity deficits. In this work, we investigate regional and whole brain growth trajectories of the myelin water fraction, a quantitative magnetic resonance imaging measure sensitive and specific to myelin content, in data acquired from 122 healthy male children from 3 to 60months of age. We examine common growth functions to find the most representative model of myelin maturation and subsequently use the best of these models to develop a continuous population-averaged, four-dimensional model of normative myelination. Through comparisons with an independent sample of 63 male children across the same age span, we show that the developed model is representative of this population. This work contributes to understanding the trajectory of myelination in healthy infants and toddlers, furthering our knowledge of early brain development, and provides a model that may be useful for identifying developmental abnormalities.

Introduction

The elaboration of the myelin sheath around neuronal axons, and the associated white matter maturation, is a cornerstone of human neurodevelopment. Myelinated white matter forms efficient communication pathways that shape the integrated neural systems responsible for higher order functioning (Fornari et al., 2007; Johnson and Munakata, 2005). Given myelin's critical role in brain communication, processes that disrupt its development may result in reduced brain connectivity and inefficient interneuronal communication. In turn, this may lead to altered neuronal functioning, and may contribute to some neurodevelopmental and psychiatric disorders, including autism and attention deficit and hyperactivity disorder (Courchesne, 2004; Haroutunian and Davis, 2007; Konrad and Eickhoff, 2010).

Myelination during the first five years of life is a rapid and dynamic process. Prior histological studies have established that myelination begins in the cerebellum and brainstem in utero (Barkovich et al., 1988; Flechsig, 1901; Paus et al., 2001; Yakovlev and Lecours, 1967). Following birth, myelination proceeds caudocranially from the splenium of the corpus callosum, optic radiations and internal capsule by 3–4 months; occipital and parietal lobes by 5–6 months; temporal and frontal lobes by 9–11 months (Flechsig, 1901; Yakovlev and Lecours, 1967); and continues into the second decade of life (Barnea-Goraly et al., 2005; Bartzokis et al., 2010). However, while retrospective histological studies provide the most faithful characterization of myelin development, they suffer significant limitations. They are i) inherently cross-sectional; ii) difficult to combine into a single temporal timeline, owing to differences in staining techniques and inconsistent brain coverage; iii) difficult to obtain from large specimen pools spanning the infant age-range; iv) preclude investigation of underlying structure–function relationships; and v) labor intensive. Further, they may not necessarily reflect healthy development as these studies are conducted post-mortem.

Recent in vivo magnetic resonance imaging (MRI) techniques, including conventional T1- and T2-weighted structural imaging, diffusion tensor (DT)-MRI, and magnetization transfer imaging (MTI), have become popular for investigating early brain development (Giedd and Rapoport, 2010; Giedd et al., 1999; Knickmeyer et al., 2008) and, especially, white matter maturation (Geng et al., 2012; Lebel et al., 2012; van Buchem et al., 2001). These non-invasive techniques provide detailed anatomical tissue contrast and micro-structural insight that affords a more sensitive and direct means of examining white matter development. However, these methods also have their disadvantages. While conventional MRI (T1- and T2-weighting) have shown alterations in the gray/white matter contrast (Huang et al., 2006; Paus et al., 2001) temporally mirroring myelination, these qualitative observations are influenced by a variety of micro-structural and biochemical elements (MacKay et al., 2009). DT-MRI offers a quantitative approach, with metrics including fractional anisotropy (FA), mean diffusivity (MD), and axonal and radial diffusivity (AD and RD, respectively). Changes in these metrics during development have often been associated with myelination, however these measures are also associated with changes in the local tissue architecture (Jones et al., 2013; Mädler et al., 2008). Many of these measures are also derived directly from the tensor model of diffusion that does not apply to all brain voxels, making the interpretation difficult (Wheeler-Kingshott and Cercignani, 2009). Similarly, while the magnetization transfer ratio (MTR) has been shown to correlate strongly with myelin content (Moll et al., 2011; van Buchem et al., 2001; Zaaraoui et al., 2008), the MTR is also influenced by other processes including edema and/or inflammation (Gareau et al., 2000; Vavasour et al., 2011).

Multi-component analysis of relaxation time data, also termed multi-component relaxometry (MCR), may provide a more sensitive measure of myelin content. MCR decomposes the measured MRI signal into contributions from distinct micro-structural water compartments. Prior MCR studies have consistently reported at least two water compartments: a fast-relaxing water pool attributed to water trapped between the myelin-lipid bilayers; and a slower-relaxing water pool attributed to intra-/extra-cellular water (MacKay et al., 2006; Whittall et al., 1997). Quantification of the signal from the myelin-bound water, termed the myelin water volume fraction (MWF), has been shown to strongly correlate with histological assessments of myelin content (Gareau et al., 2000; Laule et al., 2006, 2008; Odrobina et al., 2005; Webb et al., 2003) and provide improved myelin specificity compared to DT-MRI measures or MTR (Gareau et al., 2000; Mädler et al., 2008; Vavasour et al., 2005).

While MCR has traditionally been performed using multi-echo T2 decay data, a more recent approach, termed mcDESPOT (multi-component driven equilibrium single pulse observation of T1 and T2), has been proposed (Deoni et al., 2008). Though at the expense of a more complicated signal model that must include the effects of water exchange, mcDESPOT offers the potential advantages of improved SNR, reduced acquisition times, and increased spatial resolution and volumetric coverage compared to the established T2 approach. While comparison of multi-echo T2 and mcDESPOT MWF values present a known discrepancy, with mcDESPOT values being consistently larger (Zhang et al., 2013), they do, however, correspond qualitatively with histological myelin content measures in a Shaking Pup model of dysmyelination (Hurley et al., 2010), and have been used to investigate structure–function impairment in MS (Kitzler et al., 2012; Kolind et al., 2012) and other demyelinating disorders (Kolind et al., 2013). More recently, the mcDESPOT has been applied to the study of white matter maturation and healthy infant neurodevelopment (Deoni et al., 2011, 2012), revealing a strong consistency with the known spatial–temporal pattern of myelination.

A continuous and probabilistic model of myelination could alleviate these concerns. Derivation of an appropriate growth model would allow estimation of the typical mean MWF, and variance, for any age. This work sought to develop such a model by comparing common growth functions fit to measured MWF data. The most appropriate model was then used to generate a continuous, four-dimensional “atlas” of healthy MWF development, allowing calculation of the ‘typical’ average and standard deviation MWF maps at any desired age. As proof-of-concept, individual and group-averaged MWF maps were statistically compared to the growth model derived MWF maps, with no significant differences found. The developed atlas, therefore, represents the first continuous model of myelin maturation in healthy male infants; provides an important step for understanding the typical myelination trajectory; and provides a framework from which to identify the earliest of white matter changes.

Materials and methods

Subjects

MRI data analyzed in this work are part of an ongoing longitudinal study investigating white matter maturation in healthy, typically developing children and its relationship to behavioral development (Deoni et al., 2012). Informed parental consent was obtained in accordance to ethics approval from the Institutional Review Board of the host institution. Enrolled infants met the following inclusion criteria: uncomplicated single birth between 37 and 42 weeks; no exposure to alcohol or illicit drugs during pregnancy; no familial history of major psychiatric or depressive illness; no diagnosis of major psychiatric, depressive or learning disorder in participant; and no pre-existing neurological conditions or major head trauma. In total, 122 healthy male infants and toddlers between 70 and 1809 days of age (mean = 690.14 days, corrected for a 40-week gestation) were analyzed. Table 1 provides an age-group break down of these participants.

Table 1

Age information of the 122 scanned participants. Ages given are given in days and corrected for gestational period.

	3 months	6 months	9 months	12 months	15 months	18 months	21 months	24 months	30 months	36 months	42 months	48 months	54 months	60 months
Number of subjects	21	15	11	6	5	5	6	2	15	9	7	9	6	5
Age (range)	70–124	140–211	240–313	316–383	413–487	497–573	593–677	720–724	772–987	1030–1728	1205–1328	1368–1526	1541–1683	1728–1809
Age (mean)	101.67	178.53	271.91	346.17	466.4	522.6	634	722.00	887.87	1158.56	1284.71	1430.22	1620.67	1758.00
Gestational period (mean)	39.31	39.58	39.51	39.95	39.8	39.54	40.17	37.14	39.47	39.02	39.84	38.90	39.29	39.51

Measuring MWF in infants

Whole-brain MWF maps were acquired using the rapid mcDESPOT (Deoni et al., 2008) imaging technique. Imaging protocols for this age-range, using acoustically-muffled sequences, have been presented previously (Deoni et al., 2011, 2012), and comprise 8 T1-weighted spoiled gradient echo images (SPGR or spoiled FLASH), 2 inversion-prepared (IR)-SPGR images and 16 T1/T2-weighted steady-state free precession (SSFP or TrueFISP) images.

MWF maps were calculated from these data using a three-pool signal model estimating intra/extra-axonal water; myelin-associated water; and a non-exchanging free water pool (Deoni et al., 2013). Corrections for flip angle (B1) and off-resonance (B0) inhomogeneities were also performed (Deoni, 2010). Total imaging times ranged from 19 min for the youngest infants to 24 min for older and larger children.

Children under 4 years of age were scanned during natural, non-sedated, sleep; while children over this age were able to watch a favorite TV show or movie. All data was acquired on a 3 T Siemens Tim Trio scanner equipped with a 12 channel head RF array. To minimize intra-scan motion, children were swaddled with an appropriately sized pediatric MedVac vacuum immobilization bag (CFI Medical Solutions, USA) and foam cushions were placed around their head. Scanner noise was reduced by limiting the peak gradient amplitudes and slew-rates to 25 mT/m/s. A noise-insulating insert (Quiet Barrier HD Composite, UltraBarrier, USA) was also fitted to the inside of the scanner bore. MiniMuff pediatric ear covers and electrodynamic headphones (MR Confon, Germany) were used for all scanned children, while a pediatric pulse-oximetry system and infrared camera were used to continuously monitor the sleeping infants during scanning. After acquisition, image data was assessed for motion artifacts (blurring, ghosting, etc).

MR analysis and myelin trajectory modeling

Following calculation of the 122 MWF maps, all maps were non-linearly aligned to a study specific template (Deoni et al., 2012) using the Advanced Normalization Tools software package (Avants et al., 2008) and smoothed with a 3 mm Gaussian kernel. Non-brain parenchyma was removed using FSL's brain extraction tool (BET) (Smith, 2002). Regional masks for bilateral frontal, temporal, parietal, and occipital lobes, cingulum, and cerebellar white matter, as well as the genu, splenium and body of the corpus callosum were derived from the MNI adult template (Mazziotta et al., 2001), co-registered to the study template, and superimposed upon each infant's MWF maps (Deoni et al., 2012). Mean MWF values for each mask were obtained for each infant and plotted against the infant's gestational-corrected age. To examine whether the nonlinear transformation affected the quantitative values, mean MWF values were also extracted from native space after applying the inverse warp transformation to the regional masks. These native-space values were then compared to the standard space values.

To the plotted data, the growth models of Gompertz (1825), Stannard et al. (1985), Richards (1959), Bleasdale and Nelder (1960), as well as more generic functions such as the logistic model, and hyperbolic tangent were fit. In addition to these, a modified form of the Gompertz model was also included in the group of growth models. Although diverse, each of these sigmoidal functions shares similar characteristics: 1) a lag or period of slow growth; 2) a period of rapid exponential growth; and 3) a period of reduced growth-rate (Fig. 1A). From these functions, biologically relevant metrics, such as lag period, growth rate, and maximum size, may be derived (Fig. 1A).

Fig. 1

(A) Properties of a sigmoid function. Although many functional forms exist, each shares three similar characteristics: 1) an initial lag or period of slow growth, 2) a period of rapid exponential growth, and 3) a period of reduced growth-rate. These properties and shape of the overall curve are governed by the free parameters of the model. (B) Modified Gompertz model is described by 4 free parameters, each contributing to the characteristics of the curve. These parameters may be useful for describing biologically relevant metrics such as a developmental transitionary period (α), developmental lag, or (β) growth rates (γ, δ).

Non-linear regression via a Levenberg–Marquardt non-linear least squares (Levenberg, 1944) algorithm was used to determine the best fit for the free-parameters of each sigmoidal function (see Table 2). Fits of each model were compared using the Bayesian Information Criterion (BIC, Schwarz, 1978), a parsimony metric that compares the fit residuals while penalizing for the number of model parameters. The model that provided the lowest BIC measure consistently across the investigated brain regions was defined as the most representative.

Table 2

Fitting functions used to fit the myelination trajectories in children less than 2 years of age. The fitting routine first estimates the initial values for the free parameters and is then followed by non-linear regression that iteratively solves for the free parameter values by minimizing the sum of squares between the predicted and measured values. Free parameters are denoted by boldface Greek letters.

Functional form of modes
Model	Number of parameters	Functional form of modes
Gompertz	3	α ∗ exp[− exp(β − γ ∗ x)]
Modified Gompertz	4	α ∗ exp[− exp(β − γ ∗ x) + δ ∗ x]
Bleasdale and Nelder	4	(α + βxδ)− 1/γ
Richards	4	α ∗ (1 + β ∗ exp[γ ∗ (δ − x)])− 1/δ
Stannard	4	α ∗ (1 + exp[−(β + γ ∗ x)/δ])−δ
Hyperbolic tangent	4	α + β ∗ tanh(γ ∗ x − δ)
Logistic	3	α ∗ [1 + exp(β − γ ∗ x)]− 1
General logistic	6	α + (β − γ) ∗ (1 + ζ ∗ exp[−β ∗ (x − δ)])− 1/ε

Construction of a population-averaged MWF atlas

Using the BIC-selected ‘best’ model, voxel-wise fitting using wild-bootstrap with residual resampling (Efron, 1979) was performed to generate whole-brain, three-dimensional maps of the mean model parameters and their associated uncertainty. The resampling was performed 5000 times to provide accurate estimations of each parameter's distribution.

With this a priori knowledge of each parameter's mean, representative whole-brain MWF maps can be constructed for any age. Further, knowledge of the parameter uncertainty also allows calculation of the variance in this representative map aswhere represents the partial derivative of the growth model with respect to the i model parameter, is the uncertainty of the i parameter, and N is the number of free parameters in the model. Individuals can be directly compared to this population-derived model, e.g., using Z-statistics, without necessitating age-grouping or requiring substantive study sizes.

Comparison of MWF atlas to grouped and individually measured MWF data

To illustrate the similarity between model-derived and in vivo acquired MWF maps, two main analyses were performed on an independent sample of 63 male children. Subjects for these analyses were not included in the group of 122 participants used in generating the model and met the same inclusion/exclusion criteria as mentioned previously. Table 3 provides supplemental information on these additional subjects.

Table 3

Supplemental information on additional subjects used in comparing model-derived MWF maps to in vivo MWF maps. Comparison of model-derived MWF maps to group averaged MWF maps were broken into 9 separate age groups.

Model-derived MWF map comparison to group averaged MWF maps
	Group 1	Group 2	Group 3	Group 4	Group 5	Group 6	Group 7	Group 8	Group 9
Number of subjects	4	8	9	7	5	6	6	4	10
Age (range)	106–141	253–309	312–396	410–477	502–583	593–652	712–942	1263–1294	1351–1809
Age (mean)	121.75	280.625	350.89	445.86	536.6	624.83	815.5	1278	1645.3

First, we examined the model's ability to examine average group differences. In vivo comparison data was broken into 9 different age groups (Table 3) and compared with age-matched model-derived mean MWF and δMWF maps. In vivo mean and standard deviation MWF maps were generated by calculating the mean and standard deviation of the MWF for each age group (Table 3). Paired t-tests were performed using FSL's Randomise tool (www.fmrib.ox.ac.uk/fsl/) to compare these in vivo and model-derived data, with significant differences defined as p < 0.05, uncorrected for multiple comparisons. Permutation testing was restricted to a custom white matter mask for each group. This custom white matter mask was created using the mean in vivo MWF map (for each group) by thresholding voxels with a MWF value below 0.02.

Next, the utility of using model-derived MWF maps to investigate individual differences was examined. The developed MWF model was compared to four individuals from which three repeated datasets have been acquired. Repeated measurements were acquired from these individuals as part of the longitudinal study protocol (Deoni et al., 2012). For each repeated measurement, age-matched model-derived MWF and maps were calculated. Individually measured MWF maps were normalized into the common study template space as previously described. Z-statistic analysis was then performed to identify differences between the acquired and model MWF maps. The Z-statistic was calculated at each imaging voxel:where, x represents the i MWF voxel from the in vivo MWF map, u and s are corresponding voxels from the model-derived MWF and δMWF maps, respectively. Z-statistic analysis was also restricted to a white matter mask, created by thresholding voxels below with a MWF value below 0.02 from the individual's MWF map. Areas of significant deviation (p < 0.05, uncorrected for multiple comparisons) were defined as |Z| < 1.96.

Post-hoc construction of confidence interval for corpus callosum

To further demonstrate that the selected model best characterizes the underlying growth trajectory of the MWF, 95% confidence intervals were constructed for the regional trajectories of the genu, body, and splenium of the corpus callosum. Non-linear fitting using the wild-bootstrap with residual resampling was performed on these three regions to generate mean and uncertainty measures for the parameters of the growth model from 5000 resamples. From these mean and uncertainty measures, 95% confidence intervals of the growth trajectory were constructed. Regional trajectories of the independent sample of 63 additional subjects previously mentioned (Table 3) were also compared to these normative model fits and confidence intervals. Residuals between the predicted (mean model) and measured (in vivo) MWF values were then calculated for each region.

Results

Regional MWF trajectories for each investigated brain region are shown in Fig. 2. All trajectories follow the expected sigmoidal shape, with a period of lagged growth extending from birth through approx. 90–150 days, depending on the brain region, followed by rapid, exponential growth to approx. 400 days of age. Beyond 400 days, the trajectory begins to appear logarithmic in nature. We found transformed MWF values to be highly correlated with native-space values, with Pearson's r correlation coefficient ranging from 0.9225 to 0.9974 (Fig. 3). Thus, the normalization did not result in inhomogeneous changes to the MWF maps and consequently was found to be reliable for aligning individual MWF maps.

Fig. 2

Derived myelination trajectories for the left and right hemisphere frontal, temporal, occipital, parietal and cerebellar white matter, as well as body, genu, and splenium of the corpus callosum from 122 healthy male infants under 5 years of age. Right and left hemisphere measurements are denoted with black squares and gray circles, respectively. Error bars represent the standard deviation of the measurement. These developmental trajectories exhibit the characteristic ‘S’-shaped curve of a sigmoid, suggesting this class of models might be best to characterize the underlying pattern.

Fig. 3

Plot comparing representative native-space and template-space mean MWF values. High correlation between the native-space and template-space illustrate that MWF values were not significantly altered through the normalization process.

A representative sample of the best-fit results for each model to data from the body of the corpus callosum is shown in Fig. 4. A summary of the entire fitting results, including free parameter estimates and BIC values, for each region is shown in Table 4. Based on BIC measure comparisons, we found that the modified Gompertz model consistently provided the most faithful characterization of the MWF data. Thus, this function was subsequently used to construct the 4-D developmental atlas by fitting this function to each individual voxel.

Fig. 4

Representative fitting of the mean myelination trajectory for the body of corpus callosum. Blue points represent mean MWF values from the 122 male infants, while the red curve represents the best-fit curve for each investigated model. In total, 8 sigmoid models were examined. BIC metric values indicate that the modified Gompertz growth model provided the best fit.

Table 4

BIC values for investigated ROIs and fitting functions. Values highlighted in bold typeface indicate the model that best fit the data.

Region	Gompertz	Modified Gompertz	Bleasdale & Nelder	Stannard	Richards	Logistic	General logistic	Hyperbolic tangent
Body of CC	61.027	51.966	130.787	68.032	68.879	81.048	63.046	53.977
Genu of CC	70.769	55.949	252.652	82.122	76.479	88.235	62.128	74.276
Splenium of CC	67.410	54.001	340.922	74.421	75.885	86.653	68.075	58.786
Right cerebellar WM	49.022	53.036	68.268	54.375	52.338	52.119	61.976	52.321
Left cerebellar WM	65.687	70.389	87.534	70.371	70.774	65.401	80.390	70.787
Right frontal WM	35.208	28.558	102.033	42.239	43.058	54.718	41.890	33.543
Left frontal WM	40.532	29.578	105.180	47.943	44.163	61.924	44.343	35.627
Right occipital WM	43.941	43.018	101.119	50.556	48.005	60.463	45.040	35.595
Left occipital WM	43.446	42.028	107.705	49.585	50.898	56.949	52.541	43.130
Right parietal WM	39.359	36.017	109.783	46.222	48.493	57.173	45.331	36.154
Left parietal WM	41.062	34.827	114.331	48.321	45.502	62.224	46.280	37.341
Right temporal WM	35.452	36.739	108.376	42.192	54.712	55.017	44.458	36.346
Left temporal WM	41.543	40.311	117.204	48.892	57.261	62.864	47.770	39.517
Right cingulum	65.705	61.837	99.781	73.378	64.982	84.783	68.782	57.379
Left cingulum	64.300	48.183	107.511	71.280	65.806	83.456	67.624	58.438

The predictive validity of the model is shown in Figs. 5 and 6. Fig. 5 contains results from the comparison of in vivo group-averaged MWF maps and model-derived maps at the same ages. Qualitatively, the in vivo and model-derived MWF maps appear similar. Results from the paired t-test analysis yielded few statistical differences (p < 0.05 uncorrected for multiple comparisons) between the acquired and model derived MWF maps for each age group.

Fig. 5

Mean, acquired, model-derived MWF maps, and difference maps for 9 separate age groups (Table 3). Comparison of acquired and model derived MWF maps was done using a paired t-test for each age group. Few significant differences (p < 0.05, uncorrected) were detected between the in vivo and model-derived maps.

We also investigated how well the atlas predicted individual MWF estimates at variable ages by comparing derived MWF maps with individually acquired MWF data. Representative results of these analyses are shown in Fig. 6, which shows images from the acquired maps, model-derived images, and the Z-statistic maps. Z-statistic values are indicated at each image voxel of the representative axial slice. A small number of regions of significant difference (|Z| < 1.96, corresponding to p < 0.05) were found between the acquired and model-derived MWF maps in the examined subjects.

Fig. 6

Repeated individual MWF maps at three different time points. Atlas derived MWF maps were derived for each corresponding age and compared to the longitudinal data using Z-statistics. Few areas of significant differences (|Z| < 1.96) between the repeated data and corresponding model-derived MWF maps were found.

Bootstrap analysis of the regional MWF trajectories for the genu, body and splenium of the corpus callosum was performed to construct 95% confidence intervals using the modified Gompertz model. The mean modified Gompertz fit with the 95% confidence intervals are shown overlaid on the trajectories of the three regions in Fig. 7. Residual histograms (Fig. 7) illustrate the distribution of the residuals across the trajectories of these three regions. As expected, bootstrapping the parametric fit of these regional trajectories resulted in approximately 95% of the individual measurements to be contained within the bounds of the confidence intervals and the residuals normally distributed (Fig. 7).

Fig. 7

Bootstrapped mean Gompertz fit with 95% confidence interval for mean MWF values obtained from the genu, body and splenium of the corpus callosum.

Discussion

This work sought to model white matter myelination throughout the brain across the first five years of life in healthy male infants and toddlers. Investigating the overall shape of the myelination trajectories in the investigated regions revealed a non-linear sigmoidal developmental trajectory: with a lag period (0–150 days) preceding exponential growth (150–400 days) followed by a period of reduced growth thereafter. This growth trajectory is consistent with previous studies investigating this age-range (Deoni et al., 2012, 2013; Hermoye et al., 2006). Of the sigmoidal models investigated, we found a modified version of the Gompertz function to provide the best characterization of the regional MWF developmental data. This result is also consistent with a prior study that used the traditional, 3 parameter Gompertz function to model DT-MRI metrics across a younger age-range (Sadeghi et al., 2012). However, this study did not investigate alternative functions, used data from only 3 age points and focused on a more restricted age range (under 2 years of age). Furthermore, while the traditional Gompertz function may be best to describe the underlying growth trajectory at an early age, this function approaches an asymptote as age increases, and thus limiting its ability to account for the continued myelination that occurs throughout childhood and early adulthood (Kumar et al., 2011; Lebel and Beaulieu, 2011). The inclusion of an additional linear term inside the exponent of the Gompertz function corrects for the asymptotic behavior of the Gompertz function, allowing the model to account for continued growth.

The biological interpretation of the modified Gompertz function's free parameters may provide important insight into the developmental process. Fig. 1B illustrates the role that each parameter has on the overall shape of this function. The first parameter, α, sets the overall scale of the function's trajectory and corresponds to the amplitude of the model at the transition from rapid, exponential development to reduced, continued growth. In terms of the current model, this parameter corresponds to the MWF value at which the developmental trajectory transitions from rapid, exponential growth to slow steady growth. The second parameter, β, describes the initial lag period of the model that corresponds to the time before rapid myelination. The third parameter, γ, describes the growth rate of the rapid exponential developmental period; while the fourth parameter, δ, corresponds to the growth rate of the slower, continued myelination observed to take place after 24 months of age. These measures are descriptive of the underlying maturational process and may be useful for future investigations of intrinsic growth patterns of white matter microstructure. For example, investigating the growth rate among brain regions would simply entail comparing the parameter γ of the interested regions. Furthermore, the age at which the MWF trajectory transitions from rapid, exponential growth to slow, continued development (characterized by α) may represent a pivotal time of the neurodevelopmental process (i.e. transition from infancy to toddler).

The development of a parametric growth model of myelination is also important for understanding linkages between evolving neurodevelopment and malbehavior, as well as clinical-relevance in premature infants, infants diagnosed with multiple sclerosis, among other complications. For example, children with autism have been reported to exhibit widespread atypical patterns of early brain growth, including accelerated maturation of white matter (Ben Bashat et al., 2007). MWF growth trajectories of children with autism could be directly compared to trajectories generated by the developed model in order to see where, when and how the growth trajectory deviates from the normative population. While additional work is required to evaluate comparisons with such populations, the robustness of the method demonstrated here strongly intimates its feasibility.

Although model MWF maps may be reconstructed after performing a single whole brain fit of the in vivo MWF data, the resulting fitting parameters contain no information regarding their uncertainty. To obtain information about the fitting parameter distribution, we utilized wild bootstrap resampling of the residuals. Bootstrap resampling is not new to MR imaging, as it has been used to estimate the uncertainty of fiber orientations from DT-MRI (Haroon et al., 2009; Jones, 2008; Yuan et al., 2008) and uncertainty of neural brain activity from functional MRI (fMRI) (Darki and Oghabian, 2009; Kirson et al., 2008). This non-parametric, model-based resampling method is advantageous because it does not require redundant data (i.e. multiple image acquisitions) in order to estimate standard errors and is sufficient to use in cases of homoscedasticity (uniform variance among data) and heteroscedasticity (non-uniform variance among data) (Chung et al., 2006). One potential caveat to this approach is that the model of choice (in this context, the modified Gompertz growth model) needs to appropriately characterize the measured data. However, our initial analysis of comparing a wide variety of growth functions, using BIC parsimony measures, examined this issue. We consistently found the modified Gompertz growth model to provide a better representation of the regional MWF developmental trajectories and therefore concluded that this model best represents the underling developmental pattern.

MWF imaging has a long history in the field of known demyelinating disorders, such as multiple sclerosis (MacKay et al., 2009), however, its use in examining structural and functional development is new. Further, mcDESPOT differs from the conventional and established techniques that have been verified histologically (Laule et al., 2006; Webb et al., 2003). Similar verification of mcDESPOT has been limited to histological comparisons in the Shaking Pup model of dysmyelination (Hurley et al., 2010) and indirectly through comparison with the known histological time-course of myelination in human infants (Deoni et al., 2011, 2012), and demyelination studies in MS (Kitzler et al., 2012; Kolind et al., 2012). Thus, the specificity of mcDESPOT MWF measures as a reflection solely of myelin may be questioned. However, animal and in-vivo results garnered so far give confidence that if not specific to myelin, mcDESPOT provides novel information regarding white matter microstructure, and offers enhanced sensitivity to myelin changes relative to T1 and T2 relaxation times, or other measures.

Outliers, caused by intra-scan motion and other confounds, have the potential of influencing the calculation of the quantitative imaging maps (i.e. MWF maps) and therefore affecting the sum of squares fitting, which could lead to erroneous results (Chang et al., 2005). Further, these outliers may unreasonably inflate variability when performing wild bootstrapping. Motion artifacts were minimized as the scanning was performed during natural, non-sedated sleep. Images were further assessed for motion artifacts and none were found to be present. Other potential imaging-related confounds, such as B0 and B1 inhomogeneity, are accounted for within the mcDESPOT processing pipeline by mapping the B1 field and removing off-resonance effects by phase-cycling the bSSFP data (Deoni, 2010). A 3 mm full-width-at-half-maximum 3D Gaussian kernel was additionally used to spatially smooth the MWF maps of each participant. Although spatially smoothing the data reduces its effective spatial resolution, the kernel size chosen was reasonably conservative relative to comparable MRI studies, preserving image detail and minimally impacting MWF values. Qualitative inspection of the regional trajectories in Fig. 2 as well as the 95% confidence intervals and residual distributions constructed for the genu, body and splenium of the corpus callosum (Fig. 7) suggests that individual MWF measurements are well contained. Nonetheless, smoothing filter size has been shown to augment statistical neuroimaging findings (Jones et al., 2005) and therefore this parameter has an impact on the findings.

Quantifying the uncertainty of each parameter represents a critical step in the development of a model that describes MWF maturation. Mean and standard deviation estimates of each parameter can then be used to compute representative average and standard deviation MWF maps, for a given age between 70 and 1809 days, which can then be used to statistically compare against other typically developing and at-risk children in this age range. This aspect of the developed model was investigated using longitudinal measurements acquired in four subjects. The variation of the Z-statistic values observed in Fig. 6 suggests that there is individual deviation along the overall developmental trajectory. However, while this individual variation is known to exist (Giedd and Rapoport, 2010; Zilles and Amunts, 2013) and is therefore expected to be observed, a small number of these inter-individual deviations were found to be significant (|Z| < 1.96). In particular, examination of the first scan of subject 2 (Fig. 6) reveals areas of higher MWF in right hemispheric frontal and bilateral temporal white matter within the in vivo MWF map than the model-derived map. These results suggest that, at this point in time, the individual deviates from the normal myelination trajectory in these regions. The individual Z-statistic results suggest the developed MWF model to be representative of myelin maturation and sensitive to individual variation along the typical myelination trajectory. Moreover, although the paired t-test analysis of group averaged in vivo and model-derived MWF maps yielded few significant differences (p < 0.05, uncorrected), these differences did not remain significant after correcting for multiple comparisons using a cluster correction technique (FSL, cluster threshold of 1.96) in a post-hoc examination of these statistical results, suggesting the model-derived MWF maps to be representative of this normative population. This represents a significant new direction for investigating early neurodevelopment. Importantly it provides a means to identify where and when myelination deviates from the typical trajectory.

In this work, we have restricted ourselves to single gender data across a restricted age range. This was to ensure a homogeneous sample with which to explore the modeling functions as well as avoid potential gender-based variations in development. It should be noted, however, that while the goal of this work was not to identify potential gender differences in myelination, developmental trajectories have been reported to differ between males and females. For example, Lenroot et al. (2007) observed males and females to have distinct patterns of brain growth, specifically noting males to have an accelerated rate of white matter volume maturation. While volumetric measures, such as white matter volume, reflect broad microstructural composition, the observed gender differences are likely to be influenced by changes in myelin content. Hence, we would expect to find similar gender differences in the trajectories of MWF and a future study investigating the role that gender has on MWF development is of great interest.

Another aspect of the developed model not investigated in this work is generalizability of the model to data acquired from other imaging centers and scanners. Indeed, in order for model-derived MWF maps to have broad applicability and provide meaningful comparisons across data from different groups, values should be consistent across imaging centers and scanners. While the mcDESPOT imaging technique has been reported to have high intra- and inter-site reproducibility (Deoni et al., 2009), none of the data analyzed in this work has been separately acquired at a different imaging facility nor has any participant been imaged on a different MRI scanner. Scanning individuals at different imaging facilities as well as on multiple scanners would not only be of significant interest but necessary to investigate this aspect of the developed model.

In this work, we have outlined a framework that can be used to model myelination trajectories that are derived from MWF measurements in children under the age of 5 years. We have shown that a modified Gompertz growth model provides the most accurate representation of myelin maturation using the BIC parsimony metric and have developed the first 4D atlas of myelin maturation in healthy males under 5 years of age. This work demonstrates the ability to accurately model early myelin maturation trajectories and provides a normative template of typical myelination in healthy male infants across a large age span (between the age of 70 and 1809 days) from which atypical myelin development can be assessed. The resulting model thus provides an important step for understanding ‘typical’ myelin development as well as providing the ability to identify when and where white matter abnormalities occur in neurodevelopmental disorders.