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 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.

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 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.

Neurodevelopmental disorders, or even psychiatric disorders with a neurodevelopmental origin such as schizophrenia¹⁸, generally need better neurobiological characterisation during development¹⁹.

Rodents reach adulthood after a couple of months and are an ideal model for mammal brain development. While brain cell developmental trajectories are key to assessing neurodevelopment, there is a lack of translational techniques allowing for a non-invasive measure of their properties. Typically, EEG and MEG signals contain brain activity cell-specific information. Still, it is hard to retrieve in practice, and the existing forward or inverse models notably use prior information on cell morphology. PET can be cell-specific and provide information about metabolism, but it is blind to cell morphology and relies on radioactive ligands, making it very expensive and potentially harmful. fNIRS and ultrasounds are too indirect to measure cell-specific information as they rely on the hemodynamic response. Last but not least, MRI is non-invasive and non-harmful. The MR signal can be sensitised to assess different tissue properties. In the following paragraphs, we will summarise the importance and limitations of some MR methods sensitive to brain structure and microstructure during development and introduce our primary method, namely diffusion-weighted MR spectroscopy (dMRS).

(Structural MRI: anatomy) Manganese-enhanced MRI (MEMRI) is particularly suited to studying whole-brain early development in rodents. It provides high structural contrast and a fast acquisition time in neonates^2,20. Still, the method is blind to the microstructural changes that lead to the volumetric changes reported.

(dMRI: microstructure) Diffusion-weighted MRI (dMRI) measures the diffusion properties of water in the brain: restricted and hindered by cell membranes, water diffusion provides an indirect measure of the tissue microstructure²¹. dMRI is versatile and applied in numerous applications, including the characterisation of the microstructural organisation during neurodevelopment in rodents and, to a lesser extent, in humans. Comparison with postmortem microscopy shows a correlation between a high cortical fractional anisotropy (FA) at birth and the presence of radial glia. The latter progressively disappears in the first week along with an FA decrease, illustrating the potential of dMRI to characterise neurodevelopment^22,23. In the cerebellum, the dense folding limits diffusion imaging to exploring white matter or ex vivo samples²⁴.

(dMRS: cell-specific microstructure) 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. 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, requiring large measurement volumes for dMRS (typically ~ 50-100mm³ in rodents), and long acquisition times. As for dMRI, some dMRS acquisition parameters, such as the diffusion weighting and diffusion time, enable the signal to be sensitised to different cellular morphological scales (e.g., cellular processes, cell bodies, cell extension). Using multiple diffusion times and a post-hoc model, cellular morphological parameters can be estimated from dMRS²⁶. They correlate with neuronal and glial morphologies extracted from histological measures in the healthy, adult mouse^27,28. 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 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. 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,33 (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). Normalising by the MM content is an unusual choice. Conventionally, metabolite absolute concentrations are normalised by the water content in the voxel or by the tCr absolute concentration, which is considered stable in the healthy brain. However, this study compares regions with different relaxation properties at different ages, with different metabolites profiles in adulthood. The water relaxation properties are different between the thalamus and cerebellum, and both water relaxation properties and content change with age^34,35. The tCr concentration changes substantially with age (about +50% in the cerebrum between P7 and P28) and is known to be higher in the cerebellum in adulthood (about +40-50% compared to the cerebrum)^35,36. The MM concentration increases by about 30% between P7 and P28. MM resonances mainly come from the cytosolic proteins and although their content might slightly vary between regions in the healthy brain, they make a reasonable intracellular reference if SNR allows. We report normalisation by water in Figure S2, but the results are noisier than those normalised by MM. 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 S3 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 for some metabolites, 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 S4 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 S5.

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 S6 for a region-to-region visual comparison and Figure S7 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) for tCr, 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. Early developmental cell-specific cerebellar and thalamic microstructural features 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 probing early cell structure development. f_sphere is notably higher in the cerebellum for metabolites compartmentalised fully (tNAA, Glu) or partially (tCr, Tau) in neurons, suggesting the influence of the cell-body dense granular layer of the cerebellum. f_sphere possibly reflects the relative volume of sphere-like structures (e.g., cell bodies) compared to the relative volume of fibrous structures (e.g., cell processes, fibres).

A study of human WM using diffusion MRI reported an increase in FA and a decrease in MD with age in early life.²² In mice, the same relationship was found in both WM and GM³⁷. We report similar dMRI trends, in both cerebellum and thalamus (Figure S8).

We also report dMRS changes with age in both regions. The general trend we find for lower signal attenuation and ADC decrease (for some metabolites) with increasing age aligns with these reported dMRI observations, and could reflect overall microstructural complexification. However, the different nature of the signal (metabolites: mostly intracellular and low compartmental exchange vs. water: ubiquitous and fast exchange) prevents close co-interpretation.

f_sphere(tCr, Tau) decreases with age in the cerebellum, meaning that tCr and Tau are located in compartments where the relative contribution of sphere-like structures (e.g. cell-bodies) decreases. In other words, it reflects cerebellar cell growth and complexification (Figure 4A, Figure 5A). In the cerebellum, while soma morphologic maturation mostly happens in the first week³⁸, Purkinje cell dendritic expansion is rapid before P15 but then continues to mature up to the end of the 4^th week^38,39. Tau diffusion properties might reflect the development of Purkinje cells: (i) in cats, which cannot synthesise Tau and must source it from their diet, Tau dietary deficiency delays cortical and cerebellar maturation⁴³, meaning that Tau is also crucial in early development, and (ii) in the cerebellum, immunohistochemistry studies locate Tau in adult Purkinje cells^40–42. These observations align well too with the high Tau levels we still detect in the cerebellum at P30 (Figure 3B).

Figure 5:

The diffusion properties of glial metabolites, f_sphere(Ins) and f_sphere(tCho), also tend to decrease with age. This could be explained by another cerebellar peculiarity: Bergmann glia develop long processes along Purkinje cells and this specific growth could contribute to the f_sphere(tCr,Tau) decrease.

Validation would require an atlas of all types of cerebellar cellular processes and bodies with estimates of their relative contribution. Such a resource would be immensely valuable but does not exist. Using literature data, we use the relative thicknesses of the EGL, IGL, PL (“cell body-like” layers) and ML (“fibre-like” layer) to derive f_sphere,LT, a proxy of f_sphere, as illustrated in Figure 5A and reported in Table S2 and Figure S11. f_sphere,LT represents the contribution of sphere-rich layers compared to fibre-rich layers. The evolution of f_sphere,LT trend with age is comparable to f_sphere(tCr). Although f_sphere(tCr) and f_sphere(Tau) have similar trends, f_sphere(Tau) values are overall lower, and below f_sphere,LT values (Figure S11). This lends weight to the idea that Tau is preferentially located in cerebellar neurons during cerebellar development (and potentially less so in granule cells). In comparison, tCr, which appears more evenly spread between all cerebellar cell types, tracks the evolution of f_sphere,LT more closely. 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 Tau 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 S9), and generally improved the estimated sphere radius standard deviation. Additionally, the evolution with age of (R,f_sphere) is then very similar for tNAA and Tau, reinforcing the idea that Tau is neuronal in the cerebellum.

In the thalamus, the decrease in Tau content (Figure 3B) and the change in Tau diffusion properties (Figure 4B) point towards a change of compartmentation with age. Contrary to the trend in the cerebellum, f_sphere(Tau) consistently increases with age in the thalamus (Figure S9). This means that Tau diffuses in a space with fewer and fewer processes during development, which is incompatible with the hypothesis that Tau stays in the same compartment. Neurons, but also microglia⁴⁵ or astrocytes^46–48 develop with increasing ramification. Neurons are already branched at P10 in the thalamus³³ (shown by a low f_sphere(tNAA, Glu) at P10). f_sphere(Tau) at P10 nears f_sphere(tNAA, Glu), indicating that Tau might be located in neuronal compartments in the neonatal thalamus. Although this hypothesis is not verifiable in existing literature, some in vitro studies show that Tau contributes to neuronal proliferation⁴⁹ and neuronal growth⁵⁰. At P30, f_sphere(Tau) is closer to f_sphere(tCho,Ins), indicating that Tau might be located in glial compartments in adulthood. Indeed, immunohistochemistry studies locate Tau in adult thalamic glia⁴².

We suggest Tau as a cerebellar marker of early neuronal development, possibly more specifically of Purkinje cells. 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, contrary to the usual neuronal markers, tNAA and Glu. tCr 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, and probably even less to Purkinje cells. Because Tau seems to switch compartments in the thalamus between birth and adulthood, we do not think its diffusion properties reliably inform thalamic neuronal morphological development. In both human^51,52 and rhesus monkey⁵³ brains, Tau content increases at birth and then decreases with age, and is therefore a marker of the immature brain in both species. Although Tau concentration is low in the adult human brain, it is more concentrated in the cerebellum in adulthood^52,54. A similar finding was made in the adult rhesus monkey brain⁵³, indicating a conserved pattern across species of Tau properties in the cerebellum. These observations form a good basis for translation. Although dMRS in infants has substantial challenges, they are worth trying to overcome for a putative cerebellar marker of early neuronal development which could help understand neurodevelopmental disorders ³.

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 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^28,29. 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 S3) 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 S7). 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^27,29, 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. If we had robust DTI in all animals (Figure S8), we would have included the fibre dispersion in the modelling (notably for the cerebellum), as in Reference ²⁹. However, since we acquired enough directions for the dMRS data to be powder-averaged, we expect that the WM contribution in the cerebellum is mostly reflected by a slight f_sphere drop, as exemplified by the f_sphere,LT drop when including WM (Table S2).

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⁵⁹ 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. The macromolecular baseline was acquired in the adult rodent with the same STELASER parameters (TE and δ), at the shortest TM=100 ms, but the stimulated echo was preceded by a double-inversion recovery (TI₁=2000 ms, TI₂=710 ms). The b-value was set to 10 ms/µm² to reduce the residual contribution of bulk metabolites, and the tCr peak residual at 4ppm (due to the shorter T₁) was corrected by subtracting a Lorentzian.

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). The DIVE toolbox uses the Gaussian phase distribution approximation 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^44,64 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^27,28, 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.

Additional information

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).

Additional files

Supplementary Information