Introduction

The cerebellum is a well-conserved structure across species and particularly across mammals. The cerebellar cortex is a highly folded structure organized in three layers: the molecular layer (containing the notable Purkinje cell dendritic trees), the Purkinje layer (containing the somas of Purkinje cells), and the granular layer (containing the granule cells). Projections to the deep cerebellar nuclei form the cerebellar white matter. Despite its small size (compared to the cerebrum), it contains more than half of the total brain neurons: this is mainly due to the granule cells, very small and densely packed cells, some of the smallest of the brain, presenting a very simple morphology¹. On the other hand, the Purkinje cells are one of the largest types of neurons, and most complex cells in the brain. Cerebellar development is protracted², starting from birth (or the last trimester of pregnancy) to a few weeks old in rodents (or 2 y.o. in humans)³. The cerebellum is involved in many basic and high-level functions, from motor control to social skills, and it is known to be affected in a number of neurodevelopmental disorders⁴.

During early development, neurons arborise. Their morphology becomes more complex and exhibits denser branching. In rodents, Purkinje cells’ dendritic trees start growing at birth. In 30 days, they evolve from a simple soma (∼10 µm diameter) to a large cell with a highly branched and planar dendritic tree (soma diameter ∼20 µm, area ∼10,000 µm² in the mouse⁵). Purkinje cell morphology can be altered in autism spectrum disorders (ASD), as shown, post-mortem, in humans⁶, and in multiple mouse models of ASD. Specifically, they show an increased spine density in Tsc1 mutant mice⁷, an overgrowth of the dendritic tree and abnormal branching in PTEN KO mice⁸, a reduction of the dendritic tree in AUTS2 mice⁹; and general alterations in developmental disorders^10–13. While cerebellar cell morphology is informative about brain development, there are very few translational techniques that allow for a non-invasive measure of their morphology.

All cerebellar outputs transit via the deep cerebellar nuclei (DCN), and a couple of them project to the thalamus before spreading to the cortex. The thalamus is considered as a relay station in the brain and contains about 60 nuclei on unique pathways. Some thalamic nuclei are involved in psychiatric disorders, and more specifically, the microstructure of the mediodorsal nucleus and the pulvinar nucleus are altered in schizophrenia and psychosis¹⁴. Some other nuclei are involved in neurodevelopmental disorders (e.g., the lateral reticular nuclei in ASD¹⁵). Generally, these disorders are associated with an altered thalamic connectivity¹⁴. Thalamocortical projections are crucial for the cytoarchitectural development of the cortex^16,17.

Rodents reach adulthood after a couple of months and are an ideal model to study mammal brain development. Manganese-enhanced MRI (MEMRI) is particularly suited to study whole-brain early development in rodents, providing high structural contrast and speed of acquisition in neonates^2,18, but the method is blind to the nature of the microstructural changes leading to the volumetric changes reported.

Magnetic resonance imaging (MRI) is the gold standard tool to access the properties of soft tissues in the body in a safe and non-invasive manner. Diffusion-weighted MRI (dMRI) measures the diffusion properties of water in the brain: constrained and hindered by cell membranes, water diffusion provides an indirect measure of the tissue microstructure¹⁹. However, water is ubiquitous, and can only partially disentangle the contribution from different tissue compartments (e.g. intra and extracellular space). Diffusion-weighted magnetic resonance spectroscopy (dMRS) measures the diffusion properties of metabolites in a large volume of interest. Brain metabolites are mostly intracellular, with some more specific to neurons (e.g., N-acetyl-aspartate, NAA, and glutamate, Glu) and others more specific to glial cells (e.g., choline compounds, tCho, and myo-inositol, Ins), therefore different cellular compartments can be resolved²⁰. However, metabolites are 10⁴-fold less concentrated than water, which necessitates large measurement volumes for dMRS (typically ∼8cm³ vs ∼8mm³ in dMRI in humans). Morphometric parameters can be extracted from dMRS experiments with multiple diffusion times and multiple diffusion-weightings, and they correlate with neuronal and glial morphologies extracted from histological measures in the healthy, adult mouse^21,22. More specifically, it has been shown that myo-inositol diffusion properties reflect astrocytic morphology in an induced model of astrocytic hypertrophy²³, and in a mouse model of demyelination where astrocytic activation was expected²⁴. Choline compound’s diffusion properties were associated with microglial activation in human disease²⁵, or in an LPS-induced brain inflammation model²⁶. However, to the best of our knowledge, very few studies have used dMRS to probe changes in neuronal morphology, and none have attempted to probe neuronal development in early life.

Here, we propose a non-invasive method to monitor cellular brain development. Neurodevelopmental disorders, or even psychiatric disorders with a neurodevelopmental origin such as schizophrenia²⁷, need better neurobiological characterisation during development²⁸. We focus on two regions with different developmental timelines: the cerebellum and the thalamus. As described above, the cerebellum expands greatly after birth, when the Purkinje cells dendritic trees start growing. This slow maturation exposes the cerebellum longer to developmental insults. In contrast to the cerebellum, the thalamus develops relatively quickly, and its neurons expand and complexify more rapidly^2,29 (Figure 1A).

Figure 1:

By contrasting these two regions, we want to assess the sensitivity of the proposed method. We scanned developing rats from postnatal day five (P5) to P30. We first used MEMRI to confirm the protracted cerebellar development compared to the thalamus. We then acquired diffusion properties of metabolites to probe short- and long-range restriction (i.e., soma restriction and branching). To interpret our data, we applied two biophysical models (one analytical and one computational) to characterise the specific changes in microstructure. We expected monotonic changes in microstructural parameters in both regions, reflecting cell growth, with thalamic properties plateauing earlier than cerebellar properties. This is mostly what we report in our observations. However, we report a U-shape trend for the cerebellar segment length (L_segment, the distance between two embranchments of a dendritic tree) that matches openly available 3D cell reconstructions. In contrast to common findings in the adult brain, we report that tNAA and Glu do not provide reliable morphological neuronal markers in neonates. We identify instead total creatine (tCr) and taurine (Tau) as more reliable markers for tracking early neuronal development. We believe that this is a promising method to track early neuronal development, particularly in the cerebellum, where microstructure changes drastically postnatally.

Results

Structural imaging confirms protracted cerebellar growth

Representative slices from intra-age non-linear registration are shown in Figure 2A, underlying the cerebellar unfolding. Both cerebellar and thalamic volumes are just below 50µl at P5 (Figure 2B). However, the segmentation at P5 is poor, as explained in Methods (Figure S1), hence data at this time point should be treated cautiously, although the cerebellum and the thalamus, which are medial line structures, are less affected by the atlas distortions.

Figure 2:

The cerebellum triples in volume, between P5 and P30, reaching 90% of its P30 volume at P20, while the thalamus reaches 90% of the P30 volume at P15.

Metabolic profile changes with age

The metabolic profile changes noticeably between P5 and P30 as shown in Figure 3A. The absolute concentrations of metabolites, estimated by LCModel using the b₀ averaged spectra, shown in Figure 3A, were normalised by the absolute concentration of macromolecules, MM (Figure 3B). Between P5 and P30, metabolite/MM ratios increase for tNAA, Glu, Gln and Ins, and decrease for Tau, Gly, PE and tCho. Note that the increase in tNAA signal with age is protracted in the cerebellum compared to thalamus and has not reached a plateau at P30. These results generally match the metabolic concentrations measured with gold-standard classical MRS in the striatum, cortex and hippocampus of developing rats reported in Reference ³⁰. We observe that tCr/MM increases strongly in the cerebellum between P5 and P30 but appears stable in the thalamus, while Reference ³⁰ reported a 30% increase between P7 and P28 in the 3 ROIs. However, they also reported a 30% increase of the macromolecular content during this developmental period, matching our tCr/MM observation in the thalamus.

Figure 3:

Metabolite diffusion properties differ between regions and change with age

Thalamic (Figure 4A) and cerebellar (Figure 4B) signal attenuations up to b=30 ms/µm² and ADCs at TM=100-1000 ms are shown for tNAA, Glu, tCr, Tau, tCho and Ins. These diffusion properties differ between both regions at each age (see Figure S2 for a region-to-region visual comparison). For most metabolites (except tCho), the signal attenuation is stronger in the thalamus, suggesting less restricted diffusion compared to the cerebellum. There is also a trend for lower signal attenuation and lower ADCs with increasing age, suggesting increasing restriction of metabolite diffusion. This is particularly striking for Tau in the thalamus. To give an estimation of the MRS fit performance, CRLBs are given for tNAA, Glu, tCr, Tau, tCho and Ins for the different diffusion conditions in Table S1.

Figure 4:

Biophysical modelling underlines different developmental trajectories of cell microstructure between the cerebellum and the thalamus

Data were fitted with a spheres + “astrosticks” model and a morphometric model (more details in the Methods section). The spheres + “astrosticks” model captures diffusion within spheres (representing cell bodies) and randomly oriented sticks (representing cell processes/neurites): parameters extracted are sphere fraction, f_sphere (or equivalently neurite fraction as 1-f_sphere) and cell radius, R. The parameter space (sphere radius vs sphere fraction) is represented in the 3^rd column of Figure 4A (cerebellum) and Figure 4B (thalamus) for tNAA, Glu, tCr, Tau, tCho and Ins at all ages (see Figure S3 for a region-to-region visual comparison).

Mean values and standard deviations are reported in Table 1, along with the results from the linear mixed effect model ∼ Age + Region + Age:Region + (1|PupID). The fits on averaged data, for binned spectra by age and region, are shown in Figure S4.

Table 1:

The morphometric model, only applied to averaged ADCs at long t_d, interprets diffusion in branched structures: the segment length, L_segment is reported on the last column of Figure 4A/B (see Figure S5 for a region-to-region visual comparison and Figure S6 for all parameters estimated by the morphometric model: D_intra, N_Branch, L_segment).

Neuronal markers

tNAA and Glu are primarily concentrated in neurons in adulthood.

• Thalamus: The sphere fraction and the radius estimated from tNAA and Glu diffusion properties do not change much with age. For both metabolites, f_sphere is around 0.25 and R decreases to 5-6 µm. L_segment is stable with age for tNAA, but shows a reduction at P15-P20 for Glu.

• Cerebellum: The sphere fraction and the radius estimated from tNAA diffusion properties vary with age. f_sphere nears 0.5 in the cerebellum and R decreases below 5 µm for tNAA and Glu. L_segment decreases until P15 for both tNAA and Glu.

Between regions, sphere fraction differs significantly (Table 1).

Glial markers

tCho and Ins are primarily concentrated in glial cells in adulthood. The general trend is a high sphere fraction (∼0.7-0.9) for Ins and moderately high (∼0.5-0.6) for tCho. The radius estimated from tCho differs significantly between regions. Note the trend of decreasing sphere fraction in the cerebellum with age. Data at P15 for Ins in the thalamus were individually not reliable enough to apply the spheres + astrosticks model (c.f. Table S1 and the particularly higher CRLBs at P15 for Ins).

Non-specific markers

tCr and Tau are not known to be compartmentalised in a specific cell-type.

• Thalamus: Parameter estimates for tCr diffusion properties do not change substantially with age (f_sphere∼0.35, R∼6.5µm). Tau parameters estimation varies a lot with age: R shows a substantial decrease and f_sphere increases.

• Cerebellum: f_sphere decreases from 0.63 (P10) to 0.41 (P30), but R is stable. Tau follows a similar pattern, although with an overall lower f_sphere (0.43 at P10 to 0.29 at P30). L_segment increases in the cerebellum from P10 for tCr and increases P15 for Tau.

Between regions, the sphere fraction differs significantly for both metabolites (Table 1), but only Tau exhibits a significant change with age, and this change is different between regions (significant interaction term). For Tau, the radius is significantly different between regions, changes significantly with age and this change is different between regions (significant interaction term).

Discussion

As in Reference ² we used MEMRI to validate the protracted development of the cerebellum in rat neonates³, unfolding and tripling volume between P5 and P30. Cell-specific cerebellar and thalamic microstructural features during early development in rat neonates were assessed with dMRS. We specifically report diffusion properties for neuronal markers tNAA and Glu, glial markers tCho and Ins, and the non-specific markers tCr and Tau. Metabolites exhibit significantly different diffusion properties between regions, demonstrating the relevance of dMRS to probe early cell structure development. The sphere fraction is notably higher in the cerebellum for metabolites compartmentalised fully (tNAA, Glu) or partially (tCr, Tau) in neurons, suggesting the cell-body dense granular layer of the cerebellum. The sphere fraction possibly reflects the relative volume of sphere-like structure (e.g., cell bodies) compared to the relative volume of fibrous structures (e.g., cell processes), illustrated in Figure 5A.

Figure 5:

While soma morphologic maturation mostly happens in the first week³¹, dendritic expansion is rapid before P15 but then continues to mature up to the end of the 4^th week^31,32. Remarkably, the sphere fraction extracted from tCr and Tau decreases with age in the cerebellum only, potentially underlining the drastic Purkinje cell growth between birth and adulthood. The diffusion properties of glial metabolites, Ins and tCho, also display a trend to sphere fraction decrease with age. This could be explained by another cerebellar peculiarity, Bergmann glia develop long processes along Purkinje cells and this specific growth could play a role in the sphere fraction decrease of Tau and tCr.

We confirm in this study the observation from Reference ³⁰ that both tCr and Tau have high spectroscopic signals in neonates. They are the two metabolites most reliably quantified up to P15 in both regions (continuing for Tau in the cerebellum up to P30) and the estimates of their diffusion properties are the most confident. The notable decrease of Tau concentration in the thalamus, the extreme sphere radius decrease, and the sphere fraction increase with age extracted from Tau diffusion properties may indicate a change in Tau compartmentation with age in the thalamus. Note that Tau has a higher D_intra³³ compared to the other metabolites, potentially affecting the parameter estimation given for Tau. As described in Methods, the fit for taurine was substantially improved for D_intra > 0.6 µm²/ms. The results for D_intra=0.7 µm²/ms are similar for all the other metabolites (Figure S7), but it brings Tau sphere radius estimation to plausible values. However, in all cases, the increase in sphere fraction with age remains: at P5, f_sphere(Tau) nears f_sphere(tNAA), while at P30 f_sphere(Tau) nears f_sphere(tCho). We can reasonably hypothesize that Tau is relatively more concentrated in neurons at birth and loses its neuronal compartmentalisation with neuronal maturation, ending up concentrated in glial cells. Moreover, the volume taken by glial cells in the brain is lower than the volume taken by neurons: the decrease in Tau signal matches this proposed switch. Tau is crucial in early development, and it has been shown in cats (who are unable to synthetise Tau, and instead must source it from their diet), that Tau dietary deficiency delays cortical and cerebellar maturation³⁴. Older immunohistochemistry studies associated Tau with Purkinje cells specifically^35,36 and Tau levels in the cerebellum are still high at P30, matching the protracted development of Purkinje cells and the cerebellum. Moreover, Tau’s decreasing sphere fraction indicates that it stays relatively highly concentrated in neurons until P30, following the neuronal cerebellar growth and complexification (Figure 5A&C). Given that the low tNAA and Glu signal up to P10 (inclusive), particularly in the cerebellum, prevents confident inference of their diffusion properties, we suggest Tau as a cerebellar marker of early neuronal development, and possibly more specifically to Purkinje cells. Total creatine stands as a good marker too. However, there is no indication of a preferential neuronal compartmentation for this metabolite during early development, hence it is less specific to neurons. This is confirmed by the higher f_sphere(tCr) compared to f_sphere(Tau) up to P15 in the cerebellum.

We also observed a protracted increase in tNAA signal in the cerebellum compared to the thalamus. tNAA is usually described as a neuronal osmolyte, meaning it could reflect neuronal growth and maturation and be a marker of neuronal development. However, at the voxel level, its signal also reflects the neuronal volumetric contribution, i.e. differences in neuronal densities might mask the effect of development. In contrast, we want to reiterate that dMRS is a normalised measure and dMRS outcomes are not affected by the neuronal density: they represent the average morphological properties of the cell types the metabolite belongs to.

We chose to interpret the dMRS data with two biophysical models: an analytical spheres & “astrosticks” model and a computational morphometric model. The analytical model was applied simultaneously to the signal attenuation at high b-values and the ADCs at long diffusion times, and on individual animals. In addition to the valuable biophysical biomarkers, it provides a good signal representation to evaluate the significance of the differences between diffusion properties of individual datasets. The experimental design was not suitable for a randomly oriented infinite cylinders (“astrocylinders”) model as in previous studies^22,23. The mixing time used for the acquisition of signal decay at high b-values (TM=100 ms) is too long to enable sensitivity to fibre radius³⁷, hence the choice of “astrosticks” model. The sphere compartment was added to account for the high sphere fraction in the cerebellum, and for the nature of the neonatal brain (under-developed neurons). The RMSEs from the different models tested (Figure S2) tend to be more similar at P30, suggesting that it is less important to add a sphere compartment for modelling metabolites diffusion in the adult brain.

The computational morphometric model was applied to attempt to extract cell extension-related parameters, but the study design and data quality are suboptimal, and we could only reliably extract L_segment and D_intra (Figure S6). D_intra is mostly stable with age in both regions. However, it sometimes drops or peaks, likely affecting the other parameters’ estimations. The longest mixing time (TM=1000 ms) is short compared to previous studies using a similar modelling approach^21,23, meaning that the diffusion path is shorter (i.e., about 30µm in this study compared to about 45µm in previous studies, assuming D_intra=0.5 µm²/ms). Note that the higher D_intra of Tau is therefore an asset for probing early neuronal cerebellar growth. In our study, the number of embranchments (N_branch) was very poorly evaluated because of the high level of noise in the data. However, even with a better dataset quality and a more optimal study design, it is unlikely that the N_branch estimation could reflect the Purkinje cell morphometry because of their size and complexity in adulthood.

Finally, we expected that L_segment estimated from the morphometric model in the cerebellum would increase regularly with age, matching the dendritic tree development. For neuronal, or partly neuronal metabolites, we report an L_segment decrease first. Except for tNAA, which consistently decreases in the cerebellum with age, L_segment decreases until P15 for tCr, Tau and Glu and then increases, forming a U-shape curve. To better understand these results, we extracted the morphometries of developing Purkinje cells and developing thalamic interneurons (both in the mouse) from NeuroMorph.org. Interestingly, from these 3D reconstructed cells, we observe a decrease in L_{segment, 3D-recon} up to P13 for Purkinje cells, whilst being quite stable in the thalamus (Figure 6), matching our non-invasive measurement trend. This could reflect the complexification of the dendritic tree following its expansion.

Figure 6:

Conclusion

This study proposes diffusion-weighted MRS as a method to longitudinally and non-invasively measure cell-specific cerebellar and thalamic microstructural features during early development in rat neonates. We demonstrate the sensitivity of Taurine to early neuronal cerebellar development. This is a stepping-stone enabling further studies in pathological conditions.

Methods

Animals

Experiments were approved by the United Kingdom Home Office. 18 Sprague-Dawley rats (9 females) were scanned at P5, P10, P15, P20 and P30. Pups came from 3 litters (6 pups scanned per litter). Animals were induced with 4% isoflurane in oxygen and anaesthesia was maintained with 1-2% isoflurane in oxygen. Respiration was monitored with a breathing pillow. Temperature was monitored with a rectal probe, which was used as a cutaneous probe from P5 to P20 to avoid any risk of harm. Isoflurane level was adapted to keep the animals at 50-70 rpm and a warm water-circulation heating pad was used to keep the animals at 36-37°C during scanning. Pups were gently held with foam and protected tape to minimise motion, because the neonatal skull is not suitable for the use of ear bars.

For logistic reasons we could not acquire MEMRI for all litters. Animals from the second litter had manganese-enhanced anatomical scans. To this end, they were intraperitoneally injected with a 50 mg/kg solution of MnCl₂ 24 hours prior to the scan. Under the age of P10, the dam received the intraperitoneal injection, and the manganese was transferred to the neonates via the maternal milk².

MRI acquisitions

Data were acquired on a 7T Biospec Bruker MRI scanner (Paravision 360.1, max gradient strength: 660 mT/m), equipped with a mouse quadrature cryoprobe, well-suited for rat neonates. An anatomical scan around the region of interest (thalamus or cerebellum) to place the voxel (T2-weighted, 2D RARE, R=8, 0.2 mm isotropic resolution, TE=33 ms, TR=3500 ms). For 6 animals, this scan was replaced by a whole brain manganese-enhanced anatomical scan (T1-weighted, 3D MGE, TE=2.4 s, TR=62 ms, 0.15 mm isotropic resolution).

dMRS acquisitions

A STELASER (Ligneul et al., MRM 2017) sequence (TE=34 ms, diffusion gradient duration, δ=5 ms) was used for dMRS acquisitions. The voxel was placed in the cerebellum for 9 pups, and in the thalamus for the other 9 pups. Voxel size was increased with age to follow brain growth. Diffusion-weighted spectra were acquired at 4 mixing times (TM) and 2 b-values (c.f. Table 2 for detailed parameters): a small “b₀” used as a crusher for the stimulated echo, and b₀+3 ms/µm². Higher b-values were also acquired at the shortest diffusion time to probe restriction within cell bodies and enable inference of cell fraction and radius.

Table 2:

Processing

dMRS

Processing was achieved with a custom Matlab (The Mathworks) routine (https://github.com/clemoune/dmrs-neonates/tree/main/A_dMRS_processing):

• Signals from different coils were combined, using the water signal at the lowest b-value and shortest diffusion time.

• Outliers were identified and discarded, as described in⁴².

  • For each transient i, the mean absolute upfield signal (range 1.9-4.5 ppm), M_up(i) was computed.

  • Transients affected by motion are characterised by a clear drop of M_up(i). Transients falling below the median(M_up) - 3xMAD threshold were discarded.

  • Experiments at high b-values are more susceptible to changing gradient direction (particularly if the voxel contains white matter tracts). To predict its effect across b-values and overcome the small number of acquisitions per unique condition (i.e. a (b-value, direction) pair), a diffusion kurtosis tensor was computed based on M_up. The MAD calculated at b=0 was used as a “ground truth” MAD and was reported to higher b-values (per direction). High b-values acquisitions are more easily corrupted by motion and therefore present an artificially high MAD.

• Individual transients were corrected for phase and frequency drifts and averaged across all diffusion directions (powder-averaged) via an arithmetic mean. In general, metabolites signal was visually detectable on single transients (with a line-broadening of 5 Hz), except at b=20,30 ms/µm² at P5, rendering these acquisitions more susceptible to motion artefacts.

• Any water residual was removed using a singular value decomposition method centred around the water region. If the baseline was heavily distorted despite this step due to very poor water suppression, the acquisition was discarded from the analysis.

Spectra were subsequently fitted with LCModel. The simulated basis set contained the following metabolites: Ace, Ala, Asp, Cr, GABA, Gln, Glu, Gly, GPC, GSH, Ins, sIns, Lac, NAA, NAAG, PCho, PCr, Tau, PE, and included a macromolecules (MM) baseline acquired experimentally (double-inversion recovery and b=10 ms/µm²).

Diffusion properties of NAA+NAAG (tNAA), Glu, Cr+PCr (tCr), Tau, GPC+PCho (tCho) and Ins are reported. In order to model both signal compartmentalisation and restricted diffusion by cell membranes, we investigated metabolites diffusion through two sets of measurements: i) the direction-averaged signal decay as a function of the diffusion weighting b, keeping the diffusion time constant at 104.25 ms (time between diffusion gradients Δ fixed at 105.91 ms and δ fixed at 5 ms) and varying the strength of the diffusion-sensitizing gradients; and ii) the apparent diffusion coefficient (ADC) as a function of increasing diffusion time (keeping δ fixed at 5 ms). For the measurements i), data are reported as the ratio S(b)/S(b=0). For the measurements ii), ADCs diffusion time dependencies are reported. ADCs are calculated following Negative ADCs were systematically removed as they point to a quantification error.

In the cerebellum, metabolites’ signal attenuation at P5 is always much stronger and quite far off the other time points. The macromolecules’ signal attenuation at high b-values is supposed to remain low⁴³ and can be used as an internal probe for motion artefact³³. However, the macromolecules’ signal attenuation is quite strong at P5 in the cerebellum (Figure S8), pointing out to a residual motion artefact despite careful processing. The spectroscopic signal in individual transients was too low at high b-values, notably in the cerebellum: it is possible that the spectral realignment failed.

MEMRI

The three first echoes of each MGE acquisition were averaged and denoised using the DenoiseImage function from the MINC toolbox, based on⁴⁴. For each time point, an image with no obvious artefact and a good visual CNR was taken as the initial model and transformed into the Fischer atlas⁴⁵ space. Images were then processed using registration_tamarack.py with the MAGeT segmentation in the available Pydpider pipelines⁴⁶. Due to the large deformations between P30 and P5, the P30 timepoint could not be used as a reference. An adequate atlas from P15 was generated in a couple of iterations (using a P20 atlas as intermediary). It was then used as a reference atlas and reference timepoint for the pipeline. The olfactory bulb and the more frontal parts were poorly registered because the image often had an intensity drop there, but it did not affect the registration of caudal regions, such as the cerebellum and the thalamus. Image with strong intensity differences, or corrupted by motion were removed from the analysis. The registration did not work well at P5, even when attempting to use an atlas generated at P10 (that is already slightly distorted).

Analysis

Signal decays at high-b values (at TM=100 ms) and ADCs at varying diffusion times were fitted together with a spheres and “astrosticks” (corresponding to randomly oriented sticks) model using the functions implemented in DIVE (https://git.fmrib.ox.ac.uk/fsl/DIVE). We fitted the signal to:

Free parameters were R (sphere radius) and f_sphere (sphere fraction). D_intra (i.e. the metabolite intracellular free diffusivity) was fixed for both compartments at 0.5 µm²/ms. Previous studies measuring brain metabolites diffusion properties at ultra-short diffusion-times^33,47 estimated the intracellular free diffusivity at 0.5-0.6 µm²/ms for most metabolites (except Tau). However other studies using different (less direct) approaches reported lower D_intra^21,22, hence we tested a range of reasonable D_intra. D_intra was fixed 0.3, 0.4, 0.5, 0.6, 0.7 and 0.8 µm²/ms. RMSEs between model fit and data were equivalent for all D_intra, except for Tau in the thalamus where increasing D_intra improved the fit, particularly until D_intra=0.5µm²/ms. Moreover, parameters estimation was unrealistic for low D_intra. Increasing D_intra between 0.5 µm²/ms and 0.7 µm²/ms did not change the trend of the results, except for Tau in the thalamus (see Figure S0 and Discussion). The fit was made robust to missing data (when excluded during processing).

To test for the model accounting best for our data, we tested different analytical models. RMSE between model fit and data were calculated for “astrosticks” only (1 free parameter, D_intra), “astrocylinders” (corresponding to randomly oriented infinite cylinders, 2 free parameters: D_intra and cylinder radius), spheres and astrosticks without constraint (3 free parameters) and spheres and astrosticks with D_intra fixed at 0.5 µm²/ms (Figure S0). RMSE were the lowest for spheres and astrosticks, and were very similar with and without the D_intra constraint. Hence the model with the least free parameters was chosen. Moreover, the results from the constrained spheres and astrosticks model were more stable and physically meaningful with the constraint on D_intra. The code used for the analysis can be found on: https://github.com/clemoune/dmrs-neonates/tree/main/B_Astrosticks_Spheres_model.

Long diffusion times

Average ADCs at varying diffusion times were fitted with the morphometric model introduced in²¹. Modeling was performed using the dictionary generated for²³ with TM = 100, 750 ms predicted by interpolation, and a random forest regressor (‘treeBagger’ function in MatlabR2020b, with tree depth = 20; number of trees = 200 and other model hyper-parameters as per default) trained in a fully supervised fashion, using the ADCs time dependencies as input and the ground truth morphological parameters as targets; as loss function the mean squared error between predicted and ground truth parameters values was used. Parameters fitted were D_intra, N_branch, L_segment, SDN_branch and SDL_segment. Boundaries for fitting the model to data were respectively [0.1, 1], [2, 25], [5, 100], [2, 3] and [5, 10]. Given the data noise, only L_segment and D_intra are estimated robustly. Fits for all parameters are shown in Supplementary Information.

Statistical analyses

Statistical analyses were run with R, using the “lmerTest⁴⁸ package. To assess whether diffusion properties were significantly different in the 2 regions during development, we applied a linear mixed effect model to the parameters extracted from the spheres + “astrosticks” model (sphere fraction and sphere radius), using age and region as fixed effects, and a random intercept for pups ID. An analysis of variance enabled to identify whether the inclusion of an interaction term age^*region was necessary or not. Results were corrected for multiple hypotheses testing with a Bonferroni correction (6 metabolites,2 parameters = 12 hypotheses). The threshold corrected p-value is 4.2e^-3 for α=0.05 and 8.3e^-4 for α=0.01.

Funding

This work was funded by Wellcome. The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z and 203139/A/16/Z). For the purpose of Open Access, the corresponding author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission. MP is supported by UKRI Future Leaders Fellowship (MR/T020296/2).