SANDI: A compartment-based model for non-invasive apparent soma and neurite imaging by diffusion MRI

This work introduces a compartment-based model for apparent cell body (namely soma) and neurite density imaging (SANDI) using non-invasive diffusion-weighted MRI (DW-MRI). The existing conjecture in brain microstructure imaging through DW-MRI presents water diffusion in white (WM) and gray (GM) matter as restricted diffusion in neurites, modelled by in ﬁ nite cylinders of null radius embedded in the hindered extra-neurite water. The extra-neurite pool in WM corresponds to water in the extra-axonal space, but in GM it combines water in the extra-cellular space with water in soma. While several studies showed that this microstructure model successfully describe DW-MRI data in WM and GM at b (cid:1) 3,000 s/mm 2 (or 3 ms/ μ m 2 ), it has been also shown to fail in GM at high b values (b ≫ 3,000 s/mm 2 or 3 ms/ μ m 2 ). Here we hypothesise that the unmodelled soma compartment (i.e. cell body of any brain cell type: from neuroglia to neurons) may be responsible for this failure and propose SANDI as a new model of brain microstructure where soma of any brain cell type is explicitly included. We assess the effects of size and density of soma on the direction-averaged DW-MRI signal at high b values and the regime of validity of the model using numerical simulations and comparison with experimental data from mouse (b max ¼ 40,000 s/mm 2 , or 40 ms/ μ m 2 ) and human (b max ¼ 10,000 s/mm 2 , or 10 ms/ μ m 2 ) brain. We show that SANDI de ﬁ nes new contrasts representing complementary information on the brain cyto-and myelo-architecture. Indeed, we show maps from 25 healthy human subjects of MR soma and neurite signal fractions, that remarkably mirror contrasts of histological images of brain cyto-and myelo-architecture. Although still under validation, SANDI might provide new insight into tissue architecture by introducing a new set of biomarkers of potential great value for biomedical applications and pure neuroscience.


Introduction
Mapping brain microstructure noninvasively using diffusionweighted MRI (DW-MRI) remains a formidable challenge due to the complexity of the underlying constituents and the relatively featureless diffusion-driven signal decay.Biophysical modelling can deliver more insight into the microstructure, thereby providing promising means for accessing MR-measurable parameters related to more specific features underpinning tissue or cellular structures.The existing conjecture or "standard model" of brain microstructure typically considers neural tissue as consisting of two compartments where endogenous water molecules diffuse (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2019;;Alexander et al., 2019;Panagiotaki et al., 2012): 1) a pool of water in neurites (axons, dendrites and neuroglial processes) thought to exhibit restricted diffusion and modelled by impermeable straight cylinders, or "sticks" if cylinder radius is assumed to be negligible (Panagiotaki et al., 2012); 2) another pool surrounding the neurites assumed to exhibit hindered diffusion and modelled as isotropic or anisotropic Gaussian diffusion (Fig. 1a).The extra-neurite pool in white matter (WM) corresponds to water in the extra-axonal space, but in gray matter (GM) it combines water in the extra-cellular space with water in cell bodies of any brain cell type: from neuroglia to neurons (collectively named soma).
This work introduces a biophysical model incorporating for the first time soma size and density in addition to neurite density, thereby enabling their joint estimation non-invasively using DW-MRI and a model-based approach.The model is motivated by recent studies that suggest the standard model of neural tissue microstructure (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2019;;Alexander et al., 2019) does not hold in GM at high b-values (McKinnon et al., 2017;Veraart et al., 2019Veraart et al., , 2020;;Palombo et al., 2018;Henriques et al., 2019a;Jespersen et al., 2019).We hypothesise that the observed departure of the standard model from the data at high b-values in GM can be largely explained by the breakdown of the assumption that water in soma exhibits similar diffusion properties as water in extra-cellular space (Palombo et al., 2018, Palombo et al., 2018).Furthermore, we propose that this unexplained signal can be accounted for by explicitly modelling the soma as one of the contributors to the intra-cellular signal (Fig. 1b).Most importantly, the resulting biophysical model enables us to estimate apparent soma size and density non-invasively using DW-MRI.
We tested our hypothesis using Monte-Carlo diffusion simulations in simplified digital models of neural cells (Palombo et al., 2019a).We show that soma size and density have indeed a specific signature on the direction-averaged (also known as powder-averaged (Callaghan and Soderman, 1983)) DW-MRI signal at high b-values that are consistent with the observed departure.Furthermore, using the same Monte-Carlo simulation framework, we show for the first time that, at reasonably short diffusion times (t d ) of few tens of milliseconds, the water exchange between neurites and soma can be ignored, supporting the design of a simple three-compartment model to separate and quantify the presence of soma (Fig. 1b).We note that it is still not clear whether using DW-MRI we can distinguish neuroglial from neuronal signal; therefore, we expect our model likely quantifies the presence of cell bodies of all cell types in the brain (e.g.neuroglia and neurons).We evaluate the resulting model with data from healthy ex-vivo mouse brain and in-vivo human brains, with results supporting the model as a promising tool for estimating apparent soma size and density.
The rest of the paper is organised as follows: in the Theory section we briefly recall the current standard model of brain microstructure and we formally introduce a new biophysical model of brain microstructure, explicitly accounting for the soma compartment and corresponding signal.In the Methods section we describe the numerical simulations and experiments used in this work to support such model, and in the Results section, we show promising results on how the proposed model enables us to characterize both cyto and myeloarchitectonic of the brain noninvasively using DW-MRI.We finally discuss the results, as well as the model limitations in the Discussion section.

Theory
The current paradigm in model-based microstructure imaging uses biophysical models, inspired by microscopy studies of tissue microarchitecture, to approximate the tissue microenvironment and estimate model parameters linked to specific tissue microstructure features from DW-MRI data (Novikov et al., 2018a(Novikov et al., , 2019;;Alexander et al., 2019).Among them, the most common class of models separate the contribution to the DW-MRI signal S measured in an image voxel into different parts that can be attributed to different "compartments" where water molecules diffuse (Panagiotaki et al., 2012).A common assumption of these so-called compartment models of diffusion is that there is no exchange between the compartments, i.e., water molecules do not move from one The total MRI signal is then given by the weighted sum of the signals from water molecules diffusing in each compartment, with relative signal fractions f in and 1-f in , respectively (a).We propose a new picture: the tissue component of an MRI voxel is subdivided into intra-cellular and extra-cellular non-exchanging compartments.The total signal is the weighted sum of the signal from water molecules diffusing in each compartment, with relative signal fractions 1-f ec and f ec , respectively.Furthermore, the intra-cellular compartment is itself divided into two non-exchanging sub-compartments: intra-neurite and intra-soma.The intra-cellular MRI signal is then given by the weighted sum of the MRI signal from water molecules diffusing within the two sub-compartments, with relative signal fractions f in and 1-f in , respectively (b).
compartment to the other.
In the simplest form, two compartments are typically used and are identified as intra-and extra-cellular.For neural tissue, it is common to assume the signal associated with the intra-cellular compartment (S intra ) is mostly due to water diffusing in elongated cellular fibres, generally called neurites (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2019;;Alexander et al., 2019), modelled as straight cylinders of zero diameter, namely sticks.Such models merge any signal contribution from soma or other large cellular domains with that of the extra-cellular compartment (S ex- tra ).Below, the detail of this standard model is first given, before we describe the proposed extension that models soma explicitly.

The standard model of neural tissue microstructure
Following (Zhang et al., 2012;Novikov et al., 2018b), we briefly recall here that in general the normalised DW-MRI signal associated with the intra-neurite compartment can be represented as a convolution between the fibre orientation distribution function (fODF) Рðb nÞ and the response kernel Κ from a perfectly aligned fiber (fascicle) pointing in the direction b n, such that (Zhang et al., 2012;Novikov et al., 2018b) where b is the diffusion weighting factor, measured along the direction b g, such that b ¼ b b g.The "stick" model assumes the functional form for the response kernel: modelled by axially symmetric Gaussian diffusion compartment, with the radial diffusion of the intra-neurite compartment D ?intra ¼ 0 (i.e.sticks).Starting from this common paradigm for the intra-neurite signal, many methods have been proposed to quantify neurite density and dispersion in both WM and GM (Jespersen et al., 2007(Jespersen et al., , 2010;;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2018bNovikov et al., , 2019;;Alexander et al., 2019;Lampinen et al., 2017Lampinen et al., , 2019;;Kroenke et al., 2004;Assaf and Basser, 2005;Sotiropoulos et al., 2012;Ferizi et al., 2015;Jelescu et al., 2016).They mostly differ in the way they model the fODF or in the way they model the contribution to the total signal from the extra-neurite compartment.Nevertheless, they all can be seen as multi-compartment models (Fig. 1a) where the total signal of the tissue component measured in an imaging voxel is given by: The signal from extra-neurite compartment A en ðbÞ is modelled as an isotropic or anisotropic diffusion tensor, with its principal direction of diffusion assumed to be parallel with the dominant direction of the fODF.Since intra-and extra-neurite compartments may generally have different T 2 values, the fraction f in is the relative signal fraction, not the absolute volume fraction (Dortch et al., 2013).Moreover, the myelin water contribution is assumed unobservable due to its short T 2 time compared to clinical DW-MRI echo time TE (MacKay et al., 1994).Also, further compartments, such as isotropic cerebrospinal fluid (CSF), can be added to Eq. (3) to accommodate partial volume contamination, such as in (Zhang et al., 2012).

SANDI: microstructure model for soma imaging
Here, we propose an intra-cellular model that consists of the intraneurite model accommodating an approximate description of the contribution from water spins diffusing within cellular soma (Fig. 1b).

Model assumptions
The proposed microstructural model is based on the same assumptions of the "standard" model and on the experimental evidence that at short t d ( 20 ms given a water bulk diffusivity of ~3 μm 2 /ms and estimated pre-exchange time !500 ms) the effect of cell (either neurons or glia) membrane permeability and corresponding water exchange between intra-and extra-cellular space is negligible (Yang et al., 2018).An additional assumption, whose validity is investigated in this work by numerical simulations, is that at short t d ( 20 ms), the two sub-compartments comprising the intra-cellular space: soma and neurites, can be approximated as two non-exchanging compartments.

General formulation
Under these assumptions, we propose the functional form for the new compartment model of brain tissue microstructure to be where f ic and f ec are the intra-cellular and extra-cellular relative signal fractions satisfying the condition f ic þ f ec ¼ 1; f in and f is are the neurite and soma relative signal fractions satisfying the condition f in þ f is ¼ 1; A in and A is are the normalised signals for restricted diffusion within neurites and soma, respectively and A ec the normalised signal of the extra-cellular space.Equation ( 4) represents a first attempt to provide a more general model to describe neural tissue microstructure that takes into account soma.Indeed, following a hierarchical decomposition of the tissue compartments similar to previous works (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Assaf and Basser, 2005;Assaf et al., 2008;Alexander et al., 2010;Stanisz et al., 1997;Yablonskiy et al., 2003), we identify an intra-cellular and an extra-cellular compartment, contributing to the total signal with relative signal fractions f ic and f ec that have to sum up to unity (Fig. 1b).Then, the intra-cellular compartment is comprised of intra-neurite and intra-soma compartment that contribute to the total intra-cellular signal with relative signal fractions f in and f is that have to sum up to unity (Fig. 1b).As such, f ic , f ec , f in and f is are not volume fractions of the corresponding constituent of the MRI voxel, but rather relative MRI signal fractions of the corresponding tissue compartment.In fact, the T2 of the intra-and extra-cellular compartments may be different (Dortch et al., 2013;Veraart et al., 2018), and here we are additionally neglecting myelin water assuming the echo time is sufficiently long to attenuate most of the myelin water contribution through relaxation.Moreover, in this first implementation, CSF contributions are not taken into account since; however, due to its quickly decaying signal with increasing b values we expect that its residual contribution would be simply captured by the extra-cellular compartment and would not significantly impact the estimates of the intra-cellular compartment model parameters.

Direction-average
Here, we focus on estimating orientation-independent features of microstructure by considering the direction-averaged signal SðbÞ (Callaghan et al., 1979).The direction-averaged signal, also known as the powder averaged, is defined as the average of the signals SðbÞ acquired along many uniformly distributed directions b g.The resulting signal takes the form: where Ãin Ãis , and Ãec are the direction-averaged normalised signals associated with their respective compartments and we used the relations f ic ¼ 1 -f ec and f is ¼ 1 -f in .The direction averaging eliminates the dependence on the fODF which is readily determined following the estimation of the orientation-independent microstructure features (Kaden et al., 2016;Callaghan et al., 1979;Celebre et al., 1992;Price, 1997).The forms of these direction-averaged signals are given below.
Extra-cellular compartment.The diffusion of water molecules associated with the extra-cellular compartment, Ãec , follows the assumptions made under the standard model.It is modelled as isotropic Gaussian diffusion with a scalar effective diffusion constant D ec : This approximation assumes that the extra-cellular space is fully connected for water molecules to sample during the course of the diffusion experiment.
Intra-neurite compartment.The signal contribution Ãin from neurites (dendrites and axons) also follows the standard model.On a voxel level, it is assumed that neuronal processes can be described as a collection of long thin cylinders, with a longitudinal apparent diffusion coefficient D k in D in and a transverse apparent diffusion coefficient D ? in e0.It is also assumed that on the timescale of our diffusion experiments (t d ~10 ms), the effects of branching and/or finite length of cellular processes can be neglected.These assumptions are appropriate under the considered experimental conditions.The root-mean-squared-displacement along the neurite would be ~5 μm, for typical longitudinal intra-neurite diffusivity which is half that of free water at body temperature.This distance is much smaller than the typical length of each cell fibre's branch, e.g.~55 μm for cerebral cortical pyramidal neurons (Jespersen et al., 2012;Palombo et al., 2016).Therefore, from the water diffusion standpoint, we can consider branching neurites as a collection of individual branches (or cylinders/sticks), randomly oriented in space (Jespersen et al., 2012;Hansen et al., 2013).Under these assumptions, the direction-averaged A in can be computed as powder average of randomly oriented sticks, such that (Jespersen et al., 2007;Kaden et al., 2016;Panagiotaki et al., 2012;Kroenke et al., 2004;Callaghan et al., 1979;Celebre et al., 1992;Price, 1997): Intra-soma compartment.The signal contribution Ãis from cell bodies is assumed to arise from a pool of diffusing water molecules restricted in spheres of radius r s .With this approximation, we are implicitly assuming negligible exchange between the pool of diffusing water molecules confined in the intra-neurite space and those in the soma.While this approximation is not valid in general, we will show, using Monte Carlo diffusion simulation in realistic models of neuronal cells, that under practical experimental conditions and typical soma size and volume fraction, the exchange between intra-neurite and intra-soma diffusing water is negligible.
Modelling soma as closed impermeable spheres, the normalised signal can be computed from the GPD approximation (Neuman, 1974;Balinov et al., 1993), such that where D is is the bulk diffusivity of water in somas, δ and Δ the diffusion gradient pulse width and separation, g the magnitude of diffusion gradient pulse, α m the mth root of the equation ðαr s Þ À1 J3 2 ðαr s Þ ¼ J5 2 ðαr s Þ, with J n (x) the Bessel function of the first kind.For simplicity, here we consider a single radius r s as representative for all the soma in a given MRI voxel.In reality, we would expect a distribution of radii P(r s ) in a given MRI voxel of few millimetres.In this case, the normalised signal could be computed from Eq. ( 8) following a volume average (Price, 2009): The r s 3 term is included to account for the spin volume, i.e., the increase in the number of spins as the radius increases.In principle, it is possible to use Eq. ( 9) to analyse experimental data.However, the inversion of Eq. ( 9), a Fredholm equation of the first kind, to provide P(r s ) is non-trivial, and to obtain an approximate estimate of apparent soma size, we prefer to use in this work Eq. ( 8).Total signal.Substituting Eqs. ( 6), ( 7) and (8) in Eq. ( 5), we get the approximated expression for the total direction-averaged signal: In general, the total free parameters to be determined from the direction-averaged data are thus six: f in , f ec , D in , D is , D ec and r s .However, from Eqn. (8), it is evident that it is challenging to disentangle D is from r s in many practical applications where data are usually acquired by varying only the magnitude of diffusion gradient; thus, it is possible to estimate only the ratio D is =r s , or D is from fixing r s or r s from fixing D is .Since the purpose of the proposed biophysical model is to characterize the microarchitecture of the brain tissue, an apparent MR estimate of r s would be more valuable than that of D is .Therefore, we will adopt the simplification of fixing D is to the value of the self-diffusion coefficient of free water, given the tissue temperature.Since all the experiments in this study were conducted in vivo, or ex vivo with the temperature kept constant at 37 C, for fitting purposes, we will fix D is e3 μm 2 /ms.However, in typical diffusion experiments, fixing D is to any value between 0.5 and 3 μm 2 /ms would not change the estimates of r s substantially, as suggested by previous studies using numerical simulations and PGSE experiments in murine erythroleukemia cancer cells (Li et al., 2017).

Numerical simulations
Monte-Carlo simulation of spin-diffusion in realistic digital models of dendritic structures were conducted to investigate the regime of validity of the assumption of non-exchanging intra-cellular compartments.Since the purpose of the simulations is to investigate when the assumption of non-exchanging neurite and soma compartments holds, only the intracellular component of the total MR signal is of interest.We first establish the regime of validity of our model in Eq. ( 10) and then further use the simulations to investigate the sensitivity to soma size and density (r s , f is ) within that regime.
Simulation setup.Detailed 3D geometries were constructed using our recently proposed generative model of complex cellular morphologies (Palombo et al., 2019a) that enables users to simulate molecular diffusion within realistic digital brain cells, such as neurons, in a completely controlled and flexible fashion (Fig. 2a and b).Here we use the generative model to mimic realistically connected neurites with different (r s , f is ) combinations.We assume cell fibres do not branch and simulate only the intra-cellular signal (Fig. 2c)hence, the generative model in Fig. 2a are not tested in simulation.Therefore, the experiments test the validity of only the assumptions of the intra-cellular compartment models: intra-neurites as randomly oriented sticks, intra-soma as sphere, and negligible exchange of diffusing spins between them.Specifically, we used 20 randomly oriented straight cylindrical segments of radius r n ¼ 0.50 μm and length L ¼ [50,200,500] μm, and spherical soma of radius r s ¼ [2, 4, 6, 8, 10] μm leading to a fraction of the total cell volume occupied by the soma, f is ¼ 4=3πr 3 s 20ÂLÂπr 2 n ranging from ~0.1 to 0.9, to simulate structures mimicking a range of possible brain cell types, from small microglia to large neurons.We note that radius of brain neurites (glial processes, dendrites, axons and neuropil in general) is typically 1.5 μm (Cardona et al., 2010).For such very thin fibres, we do not expect significant effects on the DW-MRI signal measured by sequences like Pulsed Gradient Spin Echo (PGSE) at the experimental conditions investigated here (Nilsson et al., 2017).The diffusion of 5 Â 10 5 non-interacting spins, initially uniformly distributed within the whole cellular volume, was simulated for each synthetic geometry with bulk diffusivity 2 μm 2 /ms and time step 20 μs, using CAMINO (Cook et al., 2006).The number of spins and the time step were chosen as the minimal values that guarantee stability of the simulated signal, according to previous work (Hall and Alexander, 2009).Using more spins or smaller time step would produce an identical simulated signal within <2% of error.To investigate the validity of the non-exchanging intra-cellular compartments approximation, a set of 3D digital models were created to explicitly prevent any exchange between soma and neurites.This is done by sealing off all the holes in the surface mesh of each sphere used to model soma (Fig. 2c), thus disconnecting each sphere from the neurites that extend from it, which are modelled by cylinders.Note that in these simulations we used simplified brain cell structures with non-branching neurites.This is adequate, because, as justified earlier, branching neurites can be approximated as a collection of individual non-exchanging branches under the considered experimental conditions.
Simulation 1 -Investigating the validity of the non-exchange approximation for different diffusion times.The purpose of this simulation experiment is to theoretically investigate using simulations when, in terms of chosen diffusion time, the non-exchanging neurite/soma compartments assumption used to build SANDI model holds.From the simulated spintrajectories, the normalised direction-averaged DW-MRI signal Ã was computed from a PGSE sequence with t d ranging from 1 to 240 ms, δ ¼ 1 ms and three b values: 500, 1,000 and 2,000 s/mm 2 (or, 0.5, 1 and 2 ms/ μm 2 ).This resulted in two sets of normalised direction-averaged signals: with exchange ( Ãw ) and without exchange ( Ãw=o ) that were used to compute the corresponding apparent diffusion coefficients ADC w and ADC w/o for different t d .The relative difference (ΔADC) between ADC w and ADC w/o was computed as a function of t d according to the following definition: A sensible regime where the non-exchanging-compartments approximation can be considered valid may be for those values of t d where ΔADC(t d ) < 10%.
Simulation 2 -Investigating the validity of the non-exchange approximation for different b values.The purpose of this simulation experiment is to theoretically investigate using simulations whether the non-exchanging neurite/soma compartments assumption holds for a wide range of (high) b values, when t d is fixed to a short or long value.From the simulated spins' trajectories, the normalised direction-averaged DW-MRI signal Ã was computed from a PGSE sequence with b-values ¼ [0 : 1,000 : 60,000] s/mm 2 (or [0:1:60] ms/μm 2 ) and 32 directions, uniformly distributed over the full sphere.Gradient pulse duration δ ¼ 3 ms and separation Δ ¼ 11 and 81 ms, were chosen according to the results of Simulation 1 (see Results section) and to match experimental data (see following section 3.2).Furthermore, in order to quantify the bias in model-parameter estimation due to the non-exchanging assumption, Eq. ( 10) without the extra-cellular compartment was fitted to the simulated signals with exchange between neurites and soma.The uncertainty in parameter estimation was evaluated with a Monte Carlo approach.Specifically, the residual sum of squares corresponding to the best initial fit for each (r s , f is ) configurations was used as standard deviation to randomly induce artificial Gaussian noise in our simulated signals before repeating the fitting operation.This process was performed 1,000 times.Then for each parameter, we computed its mean and standard deviation over the generated repetitions, and compared them with the groundtruth values.We performed this analysis for both conditions: t d ¼ 10 ms (when the non-exchange assumption should hold, according to Simulation 1) and t d ¼ 80 ms (when the non-exchange assumption should fail, according to Simulation 1).
Simulation 3 -Investigating the sensitivity of signal to soma size and density.The purpose of this simulation experiment is to investigate using simulation and comparison with real data the sensitivity of signal to soma size and density.Using the result in Simulation 2, we built a dictionary of simulated signals, corresponding to different microstructural scenarios, i.e. different (r s , f is ) configurations, and compared them to the experimental signals obtained from selected regions-of-interest (ROIs) in exvivo mouse brain (see following section 3.2).The direction-averaged DW-MRI signal from real data was computed for a GM ROI manually drawn in the cortex and a WM ROI in the corpus callosum.We chose these two brain regions because they are expected to have very different (r s , f is ) values.The average signal in each ROI was compared against the dictionary of simulated signals to determine whether different soma size and density could explain the non-exponential signal decay in experimental data.

Experimental data
We considered two datasets to evaluate our key modelling assumptiondiffusion of water molecules within soma has a non-negligible contribution to the normalised direction-averaged signal at high b values and can be modelled as restricted diffusion, separate from water diffusion in neuritesand to assess the proposed model's ability to estimate soma size and density.First, DW-MRI data of ex-vivo mouse brains were collected at ultra-high b values, with state-of-the-art preclinical hardware, to show that the specific signature of soma size and density on the DW-MRI signal, as predicted by Monte Carlo simulation (Simulation 3), is consistent with measured data.Second, DW-MRI data of in-vivo healthy human brains at high b values were analysed, producing maps of soma density that can be compared against published histological results, to show that the technique translates to in-vivo human studies and that it provides a novel contrast sensitive to neural tissue cytoarchitecture.

Ex-vivo mouse brain
All animal experiments were preapproved by the institutional and national authorities and were carried out according to European Directive 2010/63.A c57bl/6 mouse (N ¼ 1), male, 8 weeks old, was perfused intracardially with 4% paraformaldehyde.The brain was isolated and kept 48h in 4% paraformaldehyde and 5 days in PBS (changed daily), before being transferred to a 10 mm NMR tube filled with Fluorinert (Sigma Aldrich) for susceptibility matching.MRI experiments were performed using a 16.4 T MRI scanner (Bruker BioSpin, Karlsruhe, Germany) operating at 700 MHz for 1 H nuclei and equipped with a micro5 imaging probe (Bruker BioSpin, Rheinstetten, Germany) with maximum gradient strength 3000 mT/m isotropically.The brain was kept at constant temperature of 37 C using the probe's temperature controller.DW-MRI were acquired using a PGSE-EPI sequence with: TE/TR ¼ 20/2500 ms; δ/Δ ¼ 3/11 ms; 30 b values from 1 to 40 ms/μm 2 ; 40 gradient directions per b value, 30 b ¼ 0 images, slice thickness ¼ 0.250 mm; FOV ¼ 11.2 Â 11.2 mm; matrix dimension ¼ 224 Â 224; bandwidth ~250 kHz; resolution 50 Â 50 Â 250 μm 3 , 10 slices, 4 averages.The dataset was denoised using MRtrix3 (Tournier et al., 2012) (http://www.mrtrix.org)and corrected for Gibbs ringing (Kellner et al., 2016).No artifacts from eddy-current were observed.The direction-averaged DW-MRI signal was then computed for a GM ROI manually drawn in the cortex and a WM ROI in the corpus callosum.The average signal in each ROI was compared against the dictionary of simulated signals in Simulation 3 to determine whether different soma size and density could explain the non-exponential signal decay in experimental data.

In-vivo human brain
To provide proof-of-concept of translation to in-vivo human applications, we performed a retrospective analysis of 25 healthy young subjects (age between 25 and 35) from the MGH Adult Diffusion Dataset downloaded from the HCP data repository (https://www.humanconnectome.org).While this dataset was not acquired with the present application in mind, its sequence parameters turn out to be almost optimal for sensitivity of the direction-averaged signal to soma according to our model, i.e. t d < 20 ms and many b-values >3 ms/μm 2 .The dataset was acquired on a 3T Siemens Connectom scanner, customized with a 64 channel tight-fitting brain array coil (Keil et al., 2013) and consists of MPRAGE and diffusion scans with four levels of diffusion weighting.The b-values used were 1, 3, 5 and two acquisitions at 10 ms/μm 2 with respectively 64, 64, 128, 128 and 128 randomly distributed diffusion-encoding directions over a full sphere.The signal-to-noise ratio (SNR) of individual b ¼ 10 ms/μm 2 images was ~5 and the SNR of the direction-averaged images at b ¼ 10 ms/μm 2 was ~50.Every 14th volume was an image without diffusion weighting (b0) used for motion correction and normalisation.Other acquisition parameters were TE ¼ 57 ms, TR ¼ 8800 ms, δ ¼ 13 ms, Δ ¼ 22 ms, voxel size ¼ 1.5 mm 3 isotropic, field of view ¼ 210 Â 210 mm 2 , pixel bandwidth ¼ 1984 Hz/Px, echo spacing ¼ 0.63 ms and parallel imaging factor ¼ 3. Additional scan details can be found in (Setsompop et al., 2013).The DW-MRI data in the dataset were already pre-processed with software tools in FreeSurfer V5.3.0 (http://freesurfer.net/fswiki/FreeSurferWiki/)and FSL V5.0 (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/).Specifically, the distortion caused by the gradient nonlinearity was corrected based on the spherical harmonic coefficients (Glasser et al., 2013).For motion correction, the b ¼ 0 images interspersed throughout the diffusion scans were used to estimate the bulk head motions with respect to the initial time point (first b ¼ 0 image), where the rigid transformation were calculated with the boundary based registration tool in the FreeSurfer package V5.3.0 (Greve and Fischl, 2009).For each b ¼ 0 image, this transformation was then applied to itself and the following 13 diffusion weighted images to correct for head motions.The FSL's 'EDDY' tool was to correct for eddy current distortion.All 4 shells (b ¼ 1, 3, 5, 10 ms/μm 2 ) were concatenated (552 vol in total) and passed into the EDDY tool.After eddy current correction, the rigid rotational motion estimates obtained from both the motion correction step and the eddy current correction step were concatenated and applied to the original b-vectors for correction.

Model-parameter estimation
The five model parameters: f in , f ec , D in , D ec and r s are estimated by random forest regression (Nedjati-Gilani et al., 2017;Criminisi and Shotton, 2013), while f is is computed from the relation f is ¼ 1 À f in .To train the random forest regressor, 10 5 synthetic signals were generated using Eq. ( 10), with 10 5 random values of the five model parameters chosen uniformly distributed within the reasonable intervals: , 12] μm.For testing, we used 2 Â 10 4 previously unseen signals generated in the similar way.To match the SNR of the signals to be fitted, Rician-distributed noise was added to the synthetic data used for training and testing.We implemented a random forest regressor using the scikit-learn open source Python toolkit (Pedregosa et al., 2011).Following preliminary experiments, we built the final random forest regressor with 200 trees of maximum depth 20 and bagging as the setting that maximises the performance of the model.Further general implementation details can be found at http://scikit-learn.org/.
Accuracy and precision of intra-cellular model-parameter estimation.This section investigates the robustness of the machine-learning based fitting algorithm we use.We explored different (r s , f is ) parameters combinations, using the simulated intra-cellular signals in Simulation 2, with r s ¼ [2,4,6,8,10] μm and f is ¼ [0.01 0.02 0.05 0.15 0.30 0.45 0.60 0.65 0.85].For each combination we estimated the three free parameters f in , D in and r s (since we are not considering the extra-cellular contribution) by RF regression, and repeated the experiment 2500 times with different noise instances to estimate mean and variance of the parameter estimation and thus quantify accuracy (through bias) and precision (trough statistical dispersion or standard deviation).The f is was estimated by the relation f is ¼ 1 À f in .Different amount of Rician distributed noise was added to the simulated intra-cellular signals in Experiment 2 by adding complex Gaussian noise before computing the magnitude to simulate three SNR conditions: SNR ¼ 10 (worse scenario), SNR ¼ 50 (similar to our experimental SNR) and ∞ (i.e.no noise, ideal scenario).
Additionally, for the ideal scenario of SNR ¼ ∞, we also performed an ablation study to assess to what extent the accuracy and precision are compromised by using less and/or lower b values than those in Simulation 2. We tested four different combinations of b values, subsampled from Simulation 2 at t d ¼ 10 ms, that could be achieved by clinical scanners or more powerful human scanners such as the Connectom (Jones et al., 2018): b ¼ [0, 0.7, 1.5, 2, 3] ms/μm 2 ; b ¼ [0, 0.7, 1.5, 3, 10] ms/μm 2 ; b ¼ [0, 1, 3, 5, 10] ms/μm 2 ; b ¼ [0, 1, 2, 3, 5, 10, 25] ms/μm 2 .We explored the same set of (r s , f is ) parameter combinations as above.For each combination, we estimated the three free parameters f in , D in and r s in the same way as described above, for each b combination.We computed the mean squared error (MSE) to quantify the overall changes in accuracy and precision for each estimate compared to the ground-truth values, known by design.

Comparison with dot-compartment model
Because it has been reported (Panagiotaki et al., 2012;Alexander et al., 2010;Stanisz et al., 1997) that in fixed tissue a fraction of immobile water, known as the "dot-compartment", is not negligible for WM, we compared SANDI to a variant that replaces the sphere compartment with the simpler dot-compartment, to assess which one describes better the high b-value data, in both GM and WM.We fitted to the experimental data from ex-vivo mouse brain both SANDI (Eq.( 10)) and its dot-compartment variant where the sphere compartment in Eq. ( 10) has been substituted by a dot-compartment of relative signal fraction f dot ¼ 1f in ; specifically: The corrected Akaike's information criterion (AICc) (Burnham and Anderson, 2002) was used to compare the relative fit quality of the two models.Given a set of candidate models for the data, the preferred model is the one with the lowest AICc value.The models' degrees of freedom were 5 for SANDI (Eq.( 10)) and 4 for the dot-compartment variant (Eq.( 12)).

Comparison with histology from literature
To illustrate qualitatively that the contrasts in f in and f is maps mirror the myelo-and cyto-architecture of the brain, f in and f is maps for one representative subject were qualitatively compared against literaturederived histological images of myelin-and Nissl-stained sections of the human brain from https://msu.edu/~brains/brains/human. Furthermore, parametric maps of f is for each of the 25 analysed subjects were processed with FreeSurfer Software Suite (https://surfer.nmr.mgh.harvard.edu)and projected onto the inflated cortical surface extracted from each corresponding subject for visualisation.Projection of the average f is map across all the 25 subjects onto a common template (cortical surface-based atlas defined in FreeSurfer based on average folding patterns mapped to a sphere) was also computed and the parcellation of the cortical surface according to Brodmann areas (BA) was performed with FreeSurfer.Brodmann parcellation was chosen because it is based on differences in brain cytoarchitecture features, making it ideal to show correspondence between the proposed f is contrast and neural soma density in specific regions of the brain cortex.The particular Brodmann parcellation available on FreeSurfer and used in this study does not contain all the areas identified by Broadmann in his seminal atlas (Brodmann, 1909).However, the main BA, characterized by distinctive differences in neural soma densities, such as 1-3 (somatosensory areas), 4 (primary motor area), 6 (pre-motor area), 17 (primary visual area), 18 (secondary visual area), 44 (Broca's area, pars opercularis) and 45 (Broca's area, pars triangularis) are represented with high fidelity, following a rigorous probabilistic parcellation procedure performed by Amunts and Zilles (2015).For some of these areas, we provide also examples of histological images of cytoarchitecture from literature (Amunts et al., 1999(Amunts et al., , 2000;;Geyer et al., 2000), showing differences in neural soma arrangement and density.

Results
Regime of validity of the non-exchanging compartment model for different diffusion times and b values.Results from the first numerical simulation experiment (Simulation 1) are reported in Fig. 3. Specifically, Fig. 3a shows three different simulated ADC dependences on diffusion time t d : for the simulations where exchange between soma and neurites was considered ('exchange'); for the case of non-exchange ('no exchange') and the prediction of a simple compartment model (diffusion in randomly oriented sticks þ GPD approximation of restricted diffusion in spheres) with the relative diffusivities, f in , f is and r s known by construction ('compartments').The difference between exchange and nonexchange conditions ΔADC as a function of t d are reported in Fig. 3b, together with the threshold at 10% chosen as a reasonable level of approximation for modelling purposes.Considering different overall sizes of brain cell domains, ranging from 100 μm (approximating microglia) to 1000 μm (approximating big neurons), Fig. 3b suggests that t d 20 ms is the diffusion time regime where the exchange between soma and neurites can be neglected and a simple compartment model can be used to model the overall intra-cellular signal as a sum of two nonexchanging compartments, namely intra-neurite and intra-soma.This regime of validity for different b values is further demonstrated by the results in Fig. 4 (Simulation 2), where the direction-averaged signal as a function of b À1/2 is shown for two different t d : t d ¼ 10 ms < 20 ms (Fig. 4a) and t d ¼ 80 > 20 ms (Fig. 4b).While signals for exchange, no exchange and simple compartments perfectly overlap in almost all the simulated scenarios at t d ¼ 10 ms (Fig. 4a), they are clearly different at t d ¼ 80 ms (Fig. 4b).
Fig. 3 suggests that t d 20 ms is a suitable threshold that on one hand offers sufficiently long diffusion time to probe soma structures of diameter up to ~24 μm, and on the other introduces only a small error (≪10%) when using our multi-compartment analytical model in Eq. ( 10).This is true for cellular structures similar to large neurons and mediumsize neurons or glia (like astrocytes and oligodendrocytes) (first two rows in Fig. 3).However, for the typical size of microglia-like cellular structures, such as the panel of r s ¼ 4 μm and f is ¼ 0.31 in Fig. 3, the error expected from using the multi-compartment analytical model is too high (~20%) at t d ¼ 20 ms.This suggests that the suitable t d must be chosen according to the desired sensitivity to specific cell domains.For example, to support the study of microglia-like structures, t d < 10 ms must be used.However, this choice would reduce the sensitivity to soma of bigger celltypes, such as big neurons (e.g.r s ¼ 8 μm and r s ¼ 10 μm columns) because the characteristic length scale probed at t d < 10 ms is less than 10 μm.Fig. 4 also shows the same effect, but through the b dependence of the direction-averaged normalised signal.At t d ¼ 10 ms (Fig. 4a) the analytical model well describes the signal attenuation as a function of b (up to very high b ¼ 60,000 s/mm 2 , or equivalently 60 ms/μm 2 ), except for the case of very small cell domain and soma (panel r s ¼ 4 μm and f is ¼ 0.31 in Fig. 4a).However, when longer diffusion time is used (t d ¼ 80 ms, Fig. 4b), the analytical model fails to describe the signal attenuation as a function of b also for larger cell domains and soma (second and third rows in Fig. 4b).These results are confirmed by the direct comparison between ground-truth values and model-parameter estimates from fitting Eq. ( 10) without extra-cellular compartment to the simulated signals with exchange in Fig. 4, and reported in Fig. 5.The estimates for the case t d ¼ 10 ms are consistently close (within the error) to the ground truth for almost all the scenarios, except for the small cell domain (third row in Fig. 4a).In contrast, the estimates for the case t d ¼ 80 ms are similar to or worse (within the error) than those at t d ¼ 10 ms, especially for the small Fig. 3. Regime of validity of the compartment model.a) Comparison of the diffusion time dependence of the apparent diffusion coefficient (ADC) in cellular structures of different overall size and soma size/density, for three conditions: 1) fully connected cellular structures, simulating exchange between soma and neurites (exchange); 2) cellular structures where the connections between soma and neurites have been closed, simulating no exchange between soma and neurites (no exchange); 3) ADC computed from the compartment model Eq. ( 10) in the GPD approximation, without extra-cellular compartment (compartments).b) relative percentage difference between the ADC in the exchange and no exchange cases in a).The dashed lines show the 10% threshold used to define the diffusion time regime where the compartment model Eq. ( 10) is a reasonable approximation of cellular structures.Fig. 4. Direction-averaged normalised signal as a function of b À1/2 in (ms/μm 2 ) À1/2 for two diffusion times: 10 ms, where, according to the results in Fig. 3, we expect the compartment model to be a good approximation of the intra-cellular signal (a) and 80 ms where we expect the compartment model to fail (b).As in Fig. 3, cellular structures of different overall size and soma size/density were considered, for three conditions: 1) exchange allowed between soma and neurites (exchange); 2) exchange not allowed between soma and neurites (no exchange); 3) computed from the compartment model Eq. ( 10) in the GPD approximation, without extra-cellular compartment (compartments).
cell domains (second and third rows in Fig. 4b) scenarios.From Fig. 5, we also note that in some cases, such as f is ¼ 0.30 in cell domain 1000 μm and f is ¼ 0.32 in cell domain 400 μm, the estimates at t d ¼ 80 ms are closer to the ground truth than those at t d 10 ms.However, in these cases the standard deviation on the estimated parameters (error bars in Fig. 5) is higher, suggesting that higher fit instability and uncertainty in the parameter estimation may be the cause.
Finally, we note that these results do not change even if a much higher bulk diffusivity is used in our simulations, e.g. 3 μm 2 /ms (see supplementary Figure S1).
Sensitivity to soma size and density.Results from Simulation 3 are reported in Fig. 6 where a dictionary of pre-computed direction-averaged simulated signals as a function of b À1/2 is compared to real DW-MRI signal averaged across two ROIs representative of WM (SNR ¼ 50, at b ¼ 40 ms/μm 2 ) and GM (SNR ¼ 20, at b ¼ 40 ms/μm 2 ), collected in exvivo mouse brain.First, observe that in Fig. 6b, only the cases marked with # closely mirror the experimental data at b > 3 ms/μm 2 for WM, while only the cases marked with * for GM.The first set of cases corresponds to the cellular configurations f is ~1-5% and r s ¼ 2 μm in Fig. 6a, while the second set corresponds to f is ~60-65% and r s ¼ 6-10 μm.These configurations match our understanding of the neuroanatomy: in the WM ROI (corpus callosum) only a small volume fraction is occupied by oligodendrocytes, whose elongated soma has shorter axis of only a few μm; while in GM ROI (cortex) there is an abundance of large soma with typical size (r s ¼ 6-10 μm) compatible with cortical pyramidal neurons.
Second, the other cases demonstrate the specificity of our computational model: different cellular configurations can produce a range of distinct signal variations.For example, the panels (first row, second column) and (third row, third column) in Fig. 6b exhibit signal variations distinct from those of the two ROIs investigated.The corresponding panels in Fig. 6a show that they correspond to distinct cellular configurations: one containing large cellular domains with low volume fraction of soma of intermediate size, while the other containing small cellular domains with very high volume fraction of soma of large size.It is worth noting that these results also suggest that SANDI performs fairly well for characterizing WM, which is consistent with the standard "stick" model being a valid approximation for WM (Veraart et al., 2020).Furthermore, we see that for b-values as low as 3,000 s/mm 2 (or 3 ms/μm 2 ), the contribution from extra-cellular water is negligible, as it is not present in our simulation but present in the experimental data.
Accuracy and precision of model estimates.The study of accuracy and precision of model-parameter estimates (f is , D in and r s ) using the RF regression is reported in Fig. 7, showing that the proposed model can closely approximate (within 10% bias, or 90% accuracy) the connected cellular structure in the ideal case (SNR ¼ ∞), and maintains good accuracy and precision in more realistic case of SNR ¼ 50 and acceptable accuracy and precision in the worse-case scenario SNR ¼ 10.
The ablation study is reported in Fig. 8, showing that a minimum of five b values (or b shells), with two of them higher than 3,000 s/mm 2 (or equivalently 3 ms/μm 2 ) are required to produce parameter estimates of reasonable accuracy and precision.These results suggest that the in vivo human dataset used in this work (MGH Adult Diffusion Dataset) is adequate (third column in Fig. 8), but less and/or lower b values would be insufficient, which would lead to MSE values 2 to 30 times larger (first and second column in Fig. 7).
Model parameters maps in human brain.Parametric maps of all the five model parameters f in , f ec , D in , D ec , r s and f is for a representative human subject are reported in Fig. 9.
Comparison with histology.Qualitative comparison of f in and f is maps with histological images (from different subjects) of myelo-and cytoarchitecture (myelin-and Nissl-staining) from available human brain atlas (https://msu.edu/~brains/brains/human/)are reported in Fig. 10.
Comparison of the soma signal fraction map and brain cytoarchitecture in 25 healthy human subjects.Examples of f is maps projected onto the cortical surface (representing the whole cortical thickness) of 4 representative subjects together with the average map over the 25 subjects analysed in this study are reported in Fig. 11.Boundaries of Brodmann areas BA 1-6, 17, 18 and 44, 45 are also reported to show the remarkable match with boundaries where f is values change, in both individual subjects and average maps.Brodmann areas are known to identify regions of the brain cortex where cytoarchitecture differs in soma density and arrangement.Fig. 12 compares average MR estimates of soma density (f is map) among the 25 subjects with histological images from a few representative areas from literature showing a remarkable match of gradients in f is values with the gradient in soma density, one of the criterion used to parcellate the Comparison with the dot-compartment variant.In Fig. 13 we show the results of fitting SANDI model (Eq.( 10)) and the dot-compartment variant (Eq.( 12)) to the ex-vivo mouse data in Fig. 6, for both WM and GM ROIs.The values of the estimated fitting parameters are reported in the table, together with the AICc values for each model.We found that SANDI describes the data better than the dot-compartment variant, wich lower AICc (ΔAICc>2), for both WM and GM.

Discussion
In summary, this work proposes SANDI, a novel model to estimate apparent soma and neurite density non-invasively using DW-MRI.Our approach challenges the existing standard model (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2019;;Alexander et al., 2019) that considers water diffusion in WM and GM as restricted diffusion in neurites, modelled by "sticks" embedded in the hindered extra-cellular water (Fig. 1a).Motivated by recent studies that suggest this "stick" assumption fails in GM (McKinnon et al., 2017;Veraart et al., 2019;Palombo et al., 2018), we hypothesise that one plausible explanation for such failure is the abundance of cell bodies (namely soma) in GM relative to WM (Palombo et al., 2018).So far, the contribution from soma has not been directly modelled, but rather assumed to contribute to the overall extra-cellular compartment (for example see (Jespersen et al., 2007;Jespersen et al., 2012)) (Fig. 1a).Indeed, the underlying assumption has been that there is negligible restriction in soma because of the fast exchange rate with the extra-cellular space.However, a recent estimate of such exchange reports a water pre-exchange lifetime in neurons and astrocytes >500 ms (Yang et al., 2018).This suggests that for relatively short diffusion times t d ≪ 500 ms, restriction in soma may be not negligible.Here, we use advanced numerical simulations in realistically connected cellular structures to show that soma has indeed a specific signature on the DW-MRI signal at high b-values.The results from ex-vivo experiments in a mouse brain show that the signature predicted in simulation is both present and observable in measured signals (Fig. 6).The results from in-vivo experiments in a cohort of 25 healthy human subjects show that the proposed technique can provide maps of apparent soma density and size that meet expectations from histological imaging (Fig. 10) and current anatomical understanding (Fig. 12).These findings are also in agreement with other recent works that have challenged the validity of the standard model and its derived variants (like spherical mean technique (Kaden et al., 2016)), and showed that factors not considered by the underlying microstructural models, such as intercomponent and intracompartmental kurtosis, may cause misestimation of the model parameters (Henriques et al., 2019a;Jespersen et al., 2019).
SANDI as a first model for soma imaging.We propose a new microstructure model based on three non-exchanging compartments that explicitly includes the soma contribution to the intra-cellular signal as a pool of water diffusing in restricted geometries of non-zero size, i.e. not a dot-compartment (Panagiotaki et al., 2012;Veraart et al., 2019;Tax et al., 2019), but rather a restricted water pool, whose MR signal has a specific b and t d dependence (i.e.Eq. ( 8)) (Fig. 1b).Our numerical simulations (Figs.3-5) identify the time regime of validity for such a simple compartment model to be at relatively short diffusion times (t d 20 ms).In this time regime, the exchange between intra-and extra-cellular compartments is also negligible, for both neuronal and glial cell types, as suggested by a recent study from Yang et al. (2018).Furthermore, we show with numerical simulations and experiments in mouse brain that soma size and density have indeed a specific signature on the direction-averaged DW-MRI signal at high b values (i.e.b max > 3,000 s/mm 2 , corresponding to 3 ms/μm 2 ) and that the high b value regime can be used to increase sensitivity to geometrical restrictions of typical neural soma size ranging from a few microns (e.g. for microglia and glia) to a few tens of microns (e.g. for big neurons) (Zhao et al., 2006;Diaz-Cintra et al., 2004) (Figs.4-6 and 8).In the simulations we performed in this work, we did not include long axons because they have been shown to have small impact on the measured intra-cellular diffusion (see supplementary information of (Palombo et al., 2016)).Under these experimental conditions, the normalised direction-averaged (or powder averaged) DW-MRI signal can be expressed as the sum of three non-exchanging compartments: intra-cellular, comprised of intra-neurite and intra-soma compartments, and extra-cellular compartment (Fig. 1b).
We show that such a model, under such experimental conditions, approximates very well the expected intra-cellular signal (Figs.3-6) and provides reasonably accurate and precise estimates of neurite MR signal fraction and soma MR signal fraction and apparent size (Figs. 5,7 and 8).Note that the proposed model is very different from previously proposed models that include a sphere compartment to account for other extra-neurite compartments, such as in Stanisz et al. (1997).In fact, here we propose a model to disentangle the signal from cell-bodies of any cell type (modelled as a sphere compartment), from the signal from elongated cellular projections, such as neuroglial processes and neuronal dendrites and axons (modelled as sticks).In contrast, Stanisz et al. (1997) proposed Comparison with dot-compartment variant.The results reported in Fig. 13 show that SANDI describes the measured MRI signal better (lower AICc) than the model comprising a dot-compartment (ΔAICc>2), in both GM and WM.More importantly, we found that the signal fractions estimated using the dot-compartment model do not correspond to the neuroanatomy expected from histology and the apparent diffusivities are not in agreement with literature, especially for GM data.For WM, the dot-compartment model and SANDI provide similar estimates, suggesting that for WM, which has low soma density, both models are good approximations (Veraart et al., 2020), with SANDI still providing better fit (lower AICc).
Non-invasive cyto-and myelo-architecture maps by DW-MRI.Using the newly introduced microstructure model, we retrospectively analysed data from 25 healthy subjects from the MGH Adult Diffusion Dataset, freely available from the HCP data repository.This dataset was coincidently acquired with experimental conditions appropriate for our model assumptions: t d 20 ms and b max ¼ 10 ms/μm 2 .Parametric maps of the estimated model parameters in Fig. 9 show very reasonable and encouraging contrasts.Maps of f in have contrast highlighting the major WM tracts in the brain while f is values are higher in GM (Fig. 9).Specifically, f in values are higher in voxels mostly comprised of WM and they seem minimally affected by fibre orientation dispersion and crossing, as shown by the uniform contrast in all the WM regions, even in those characterized by high fibre orientation dispersion and crossing, e.g.regions where the radiation of the corpus callosum and the corona radiata cross.We also notice that f in values in cortical GM are between ~0.1 and 0.2, in good agreement with recent works focusing on estimating neurite density in GM using more advanced DW-MRI acquisition schemes such as spherical tensor encoding (Lampinen et al., 2017(Lampinen et al., , 2019)).Values of f is are consistently higher in all GM regions, from cortical to deep and cerebellar GM.The slightly lower values of f is in cerebellum may be due to higher partial volume in this region between WM and GM due to the large size of MRI voxels (1.5 Â 1.5 Â 1.5 mm 3 ).We note that more advanced DW-MRI acquisition schemes such as B-tensor encoding (Westin et al., 2016) seem to offer encouraging preliminary results on f is estimation using lower b values (Gyori et al., 2019), in good agreement with our estimates.Values of r s across the brain range from 2 to 12 μm, with a mean AE std values of 10 AE 3 μm.These values are in good agreement with the expected mean AE std radius of neural soma in human brain 11 AE 7 μm, as evaluated by a supplementary analysis we performed using about 3000 reconstructions of human brain cellular morphologies, available from the Neuromorpho Fig. 7. Correlation accuracy plot.Soma compartment signal fraction f is , soma apparent size r s and axial intra-neurite diffusivity D in estimated using relation Eq. ( 10) without the extra-cellular compartment and GPD approximation and labelled with the superscript "estimated" are plotted against the ground truth values labelled with the superscript "ground-truth".The perfect positive correlation line (solid line) and AE10% error (dashed lines) are plotted.In infinite SNR case, the correlation is very high (R 2 > 0.98) and bias within 10%.In the more realistic scenario of SNR ¼ 50 the correlation is still very high (R 2 > 0.85) and accuracy and precision are close to the ideal case of infinite SNR.In the worse-case scenario of SNR ¼ 10, the correlation is still high (R 2 > 0.75) and accuracy and precision acceptable.Error bars on data points indicate the statistical dispersion (standard deviation) in model parameter estimation as evaluated by Monte Carlo approach (2500 random drawn).database (neuromorpho.org).Finally, D in and D ec values are in good agreement with published values for intra-neurite diffusivity in WM being about 2.3 μm 2 /ms (Dhital et al., 2019) and extra-cellular diffusivity (Kunz et al., 2018).However, given the non-optimal experimental design of the human dataset, the maps of r s should be taken with care.In fact, r s estimation in this case is neither very accurate nor particularly precise because of the limited number of b values and only one t d .This is shown in supplementary Figure S2, where the random forest regressor predicts different and not always sensible values for r s when trained with an incorrect range, e.g.r s ¼ [1, 20] μm, instead of the correct one (r s ¼ [1, 12] μm).The model parameters f is and f in are MR signal fractions of the two compartments that in our model are linked to soma and neurites.As such, we expect them to correlate with soma and neurite densities or volume fractions within the MRI voxel.On the other hand, Nissl and myelin are two of the most used staining to characterize the brain cyto-and myelo-architecture.Although Nissl staining stains mostly the nucleus of neural cells rather than the whole cell body, we can reasonably expect that the contrast in Nissl staining correlates with density and arrangement of soma and thus the cyto-architecture of the brain.Similarly, myelin staining stains only the myelinated neurites, such as axons or myelinated dendrites, thus we can reasonably expect that the contrast in myelin staining highlights mostly WM tracts and the myelo-architecture of the brain.The qualitative comparison of f is and f in maps with Nissl-and myelin-stained histological images in Fig. 10 show a Fig. 8. Ablation study of accuracy and precision of model parameters estimation.Soma compartment signal fraction f is , soma apparent size r s and axial intra-neurite diffusivity D in estimated using relation Eq. ( 10) without the extra-cellular compartment and GPD approximation and labelled with the superscript "estimated" are plotted against the ground truth values labelled with the superscript "ground-truth" for four different b value combinations (subsampled from Experiment 2 with δ/Δ ¼ 3/11 ms).The perfect positive correlation line is plotted as solid line.Error bars on data points indicate the statistical dispersion (standard deviation) in model parameter estimation as evaluated by Monte Carlo approach (2500 random drawn) in case of infinite SNR.The mean squared errors (MSE) with respect to the groundtruth values are reported for each b value combination as metrics of accuracy and precision (lower the value, higher the accuracy and precision).remarkable similarity between the MRI maps and the contrast in the histological images, suggesting that f is and f in maps could be used to characterize in a non-invasive way the cyto-and myelo-architecture of neural tissue.However, we note that the concordance of SANDI maps and histology is not perfect.This could be due to several reasons: implicit contribution to f is from relaxation weighting; unavoidable differences between histology and imaging (e.g.thinner slice and higher in plane resolution in histology); histological and MRI images from different subjects.A particularly important difference is that, unlike its histological counterpart, f is is a fraction, thus reflecting the signal contribution of soma relative to the entire intra-cellular space.For instance, we note, unlike the Nissl staining image, f is values in the thalamus are lower than other neighbouring GM regions (e.g putamen and caudate).This is consistent with the fact that the thalamus, different from its neighbouring GM structures, consists of a large amount of WM and myelin (Whittall et al., 1997;Madler et al., 2008;Ganzetti et al., 2014) (e.g. the stratum zonale that covers the dorsal surface and the external and internal medullary laminae), leading to lower f is values (and higher f in correspondingly).Moreover, histological images and fitted SANDI parameter maps are from different subjects.On one hand, this represents a limitation of the present work that future validation work (Palombo et al., 2019) will aim to address.On the other, it makes the observed  Quantitative map of f is is expected to provide contrast related to soma density, while map of f in is expected to provide contrast related to neurite density.Since myelinated neurites (like axons) are the major constituent of white matter, f in is expected to provide contrast highly related to myelin in the white matter regions.
concordance between SANDI parameter maps and histology remarkable.
Non-invasive cyto-architectonics of the human brain by DW-MRI.Although a proper validation of the novel contrast in soma density and size introduced by the new model presented here is still in progress (Palombo et al., 2019), here we present individual and average results over 25 healthy subjects, showing remarkable similarity of f is values distribution on the cortical surface with Brodmann areas parcellation of the brain cyto-architecture (Figs. 11 and 12).As shown in Fig. 12, Brodmann areas are characterized by different soma density and arrangement.Changes in f is values on the cortical surface (representing the whole cortical thickness) follow very well the boundaries of Brodmann areas (within acceptable slight mis-alignment probably due to small errors in the co-registration procedure, performed using the automated toolbox within FreeSurfer), demonstrating that the contrast provided by this new MRI parameter could be used as non-invasive imaging marker of cyto-architectonics.Furthermore, the correspondence to Brodmann areas is also very good at the level of individual subject, as shown by Fig. 11.
Potential impact.A deeper understanding of cortical organization, including its complex fiber architecture and structural connectivity is still Fig.11.Projection onto cortical surface (representing the whole cortical thickness) of quantitative maps of MR soma signal fraction f is for 4 representative subjects and for the average over the whole cohort of 25 healthy subjects.The principal Brodmann's areas available on FreeSurfer are also reported.We notice a remarkable correspondence between the boundaries of Brodmann areas and the gradient in f is values.This is particularly evident in the average map, while in the maps of individual subjects we can also appreciate sensible inter-subject variability.
an open challenge in neuroscience.There has been considerable interest in the mapping of GM microstructure.Some examples include cortical laminar structure characterization using high-resolution DTI (Assaf, 2019), assessment of GM maturation in rodents with DKI (Cheung et al., 2009), and neurite density and dispersion quantification in human brain with NODDI (Fukutomi et al., 2018;Parker et al., 2018;Winston et al., 2014), all using PGSE experiments at intermediate diffusion times.Other works have used PGSE experiments at shorter diffusion times (t d 30 ms) (White et al., 2013;Taquet et al., 2019), or oscillating gradients at ultra-short diffusion times (Aggarwal et al., 2012), to characterize the GM microstructure.SANDI can help interpret and model in terms of soma size and density their results concerning, for example, DIAMOND   10) from the SANDI model, and the Eq. ( 12) from the dot-compartment variant.The results of the fit for the DW-MRI data from the WM and GM ROIs in Fig. 6 are reported in the table, together with the values of the AICc (lower AICc indicates better fit).(Scherrer et al., 2016) and restriction spectrum imaging (RSI) (White et al., 2013) metrics, and ADC frequency dependence in GM (Aggarwal et al., 2012).The additional contrasts provided by SANDI maps of soma MR signal fraction and apparent MR measured soma size could be useful to develop novel quantitative cytoarchitectonic and probabilistic mapping of cortical areas in a whole-brain fashion.For example, these could be used to define new parcellations of the brain based on cytoarchitectural features (e.g.Fig. 10 could be a proof-of-concept); and to improve the quality of currently available atlases of brain GM sub-regions that are notoriously difficult to delineate, such us the numerous nuclei comprising the brainstem (Bianciardi et al., 2015) (e.g.some encouraging preliminary results have been recently shown in (Bianciardi, 2020)).SANDI could also help provide metrics more specific to changes in the brain cyto-and myelo-architecture through development and due to the onset of disease.In longitudinal studies of developing brain (Neil et al., 1998;McKinstry et al., 2002;Deipolyi et al., 2005;Huang et al., 2006Huang et al., , 2009Huang et al., , 2013;;Yu et al., 2016;Ouyang et al., 2019aOuyang et al., , 2019b)), Ouyang et al. (2019b) has recently used DKI (Lu et al., 2006) to quantify the dynamic cortical microstructural signature of critical developmental stages.SANDI may help understand the exact neuroanatomical underpinning of observed cortical MK in terms of changes in neuronal soma density.In Multiple Sclerosis (MS), SANDI could provide more specific information about microglial and astrocytic activation during inflammation and astrocytic scarring, both processes involving increased accumulation of glial soma within the MS lesions (Holley et al., 2003) (e.g.some encouraging preliminary results have been recently shown in (Palombo, 2020)).In epilepsy, SANDI maps of apparent MR measured soma size and soma MR signal fraction could improve sensitivity and specificity to the remodelling of brain cytoarchitectonic occurring within the epileptic lesions (Aronica and Muhlebner, 2017).
Data acquisition requirement.In this work a comprehensive scheme was used a-priori so as to ensure high accuracy and precision of our results (Figs. 7 and 8).Further refinement and optimization will be required in the future to establish the limits of the methodology.In general, given the specific hardware characteristics of an MRI scanner (either preclinical or clinical) and acquisition time constraints, it is always possible to optimize the experimental protocol in order to achieve reasonable levels of precision in SANDI model parameters estimation.A possible general approach to perform such optimization has been previously proposed by Alexander (2008).
It is important to underline that conventional clinical DW-MRI data are in general not suitable for SANDI, for a number of reasons.First, b values are typically no higher than 3,000 s/mm 2 , or equivalently 3 ms/ μm 2 .SANDI requires several much higher b values.Second, diffusion times are typically much longer than 20 ms, the upper bound we have identified for SANDI modelling to be valid.The ablation study demonstrates that the MGH Adult Diffusion Dataset used in this work is an example of a minimal dataset that meets these requirements: short enough diffusion time, high enough b values, six b-shells to estimate five model parameters.
Similar considerations hold for ex-vivo experiments but as these experiments are commonly performed at room temperature, we must take into account the resulting lower diffusivity.Typically, apparent diffusion coefficient in ex-vivo brain at room temperature of 21 C can be ~4 times lower than in-vivo, suggesting that longer time scales (e.g.t d 80 ms) and higher b values (e.g.>12,000 s/mm 2 , or equivalently 12 ms/μm 2 ) have to be used.The ex-vivo mouse brain dataset used in this study is an example of a suitable ex-vivo experimental protocol.
Future directions and perspectives.The SANDI model presented here for apparent soma size and density estimation is used to analyse DW-MRI data acquired with simple classical MRI acquisition schemes like PGSE or Pulsed Gradient Stimulated Echo sequences.However, other more advanced acquisition schemes such as double diffusion encoding (DDE) or B-tensor encoding, may provide improved accuracy and precision of SANDI model parameters estimation.In fact, it has been shown that these acquisitions provide additional information that can help disentangling the relative contribution of the different compartments modelled by SANDI.For instance, B-tensor encoding has been successfully used to improve the estimation of neurite density based on the standard model Eq.(3) (Lampinen et al., 2017(Lampinen et al., , 2019)), while recent works using DDE showed that this acquisition scheme can help disentangling different sources of DW-MRI signal that can be linked to different features of the underpinning tissue microstructure (Coelho et al., 2019;Henriques et al., 2019b;Vincent et al., 2019;Shemesh et al., 2014).Future works will focus on harnessing the orthogonal information offered by these advanced acquisition schemes in order to maximize the sensitivity and specificity of the measured DW-MRI signal to the soma contribution (e.g.some encouraging preliminary results have been recently shown in (Afzali, 2020)).Together with improving the acquisition, another target of future work is the rigorous validation of the new parametric maps provided by SANDI.As already mentioned, SANDI provides apparent soma size and density maps in terms of MR measured spherical compartment size and relative signal fraction, respectively.By model design, these values are expected to correlate with the higher moments of the actual soma size distribution and the soma density in the brain tissue.However, histological validation is complex, expensive, and time consuming.Histological measurements of cellular content have their own inaccuracies and inconsistencies and appropriate metrics are difficult to fine tune.Preliminary investigation in ex-vivo mouse brain at ultra-high field and comparison with histological staining for cell bodies has already shown encouraging strong positive correlations between the soma size and density estimated from SANDI and those directly measured from histology (Palombo et al., 2019).Future work will extend such investigation to different mouse brain regions, such as cerebellum and olfactory bulb, and will provide more quantitative proof of the actual link between SANDI model parameter and the actual brain tissue microstructure.Nevertheless, the qualitative results (Figs. 10 and 12) showing trends consistent with known histological variation at macroscopic scale support the first demonstration of the value of SANDI.

Conclusion
The current conjecture in brain microstructure imaging (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2019;;Alexander et al., 2019;Panagiotaki et al., 2012) envisions the brain tissue component in an MRI voxel as subdivided into two non-exchanging compartments: intra-neurite and extra-neurite space.The total MRI signal is then given by the weighted sum of the signals from water molecules diffusing in each compartment.Although very successful in describing the DW-MRI signal in WM and GM at relatively low b values (b 3,000 s/mm 2 , or 3 ms/μm 2 ) in both healthy and diseased conditions (Jespersen et al., 2007(Jespersen et al., , 2010;;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2018bNovikov et al., , 2019;;Alexander et al., 2019;Lampinen et al., 2017Lampinen et al., , 2019;;Kroenke et al., 2004;Assaf and Basser, 2005;Sotiropoulos et al., 2012;Ferizi et al., 2015;Jelescu et al., 2016), this microstructure model fails in describing DW-MRI signal at high b values (b≫3,000 s/mm 2 , or 3 ms/μm 2 ) (McKinnon et al., 2017;Veraart et al., 2019Veraart et al., , 2020;;Palombo et al., 2018).Here we introduce a new picture: the tissue component of an MRI voxel is subdivided into intra-cellular and extra-cellular non-exchanging compartments.The total signal is the weighted sum of the signal from water molecules diffusing in each compartment.Furthermore, the intra-cellular compartment is itself divided into two non-exchanging sub-compartments: intra-neurite and intra-soma.The intra-cellular MRI signal is then given by the weighted sum of the MRI signal from water molecules diffusing within the two sub-compartments.This new microstructure model that we call SANDI (Soma And Neurite Density Imaging) directly accounts for one of the major differences between WM and GM: soma abundance in GM compared to WM, enabling the non-invasive estimation of apparent soma density and size trough MRI.
Using advanced numerical simulations, we identified the regime of validity of the assumption of non-exchanging intra-cellular subcompartments (neurites and soma) and propose SANDI as a new method for non-invasive soma imaging.We demonstrate it in ex-vivo DW-MRI mouse data and in-vivo cutting-edge human acquisitions.We showed that the new microstructure model for soma imaging defines new contrasts, dissimilar to the simple tensor analyses, representing new complementary information on the brain cyto and myeloarchitecture.
Although still under validation (Palombo et al., 2019), the maps here reported already show some interesting contrast that might provide new insight into tissue architecture and provide markers of pathology, as well as a new set of biomarkers of potential great value for biomedical applications and pure neuroscience.With the availability of powerful human scanners like the Connectom (Jones et al., 2018), this technique has the potential for translation into the clinic, opening a promising avenue for more in-depth assessment of cellular microstructure in-vivo in human brain.
CRediT authorship contribution statement The authors regret that due to an error in the code used to process the data, the estimated parametric maps reported in Figure 9 are incorrect.The corrected Figure 9 is included here.
The primary differences to the original Figure 9 are in the estimated parametric maps D in , D ec and f ec .However, the observed contrast in f in , f is and r s is minimally affected and remains compatible with the histological counterparts shown in Figure 10-12, hence the main conclusions still hold.
Furthermore, it has been brought to the authors' attention that our statement about the adequacy of the MGH Adult Diffusion Dataset is incorrect.The statement drew inference from the results of the ablation study shown in Figure 8, which considers only the simulated intra-cellular signal (i.e.we assume f ec = 0; thus, f ec and D ec are not estimated, leaving only three free parameters).The in vivo human data have additional signal contributions from the extra-cellular space (f ec > 0 so f ec and D ec become relevant and the number of free parameters increases five), which our simulations do not reflect.Thus, indeed, the conclusions may not extend to the in vivo human data and we would like to correct the following statement: • The original paragraph on page 10 (second column, third paragraph) reads: ' The ablation study is reported in Fig. 8, showing that a minimum of five b values (or b shells), with two of them higher than 3,000 s/mm 2 (or equivalently 3 ms/ m 2 ) are required to produce parameter estimates of reasonable accuracy and precision.These results suggest that the in vivo human dataset used in this work (MGH Adult Diffusion Dataset) is adequate (third column in Fig. 8), but less and/or lower b values would be insufficient, which would lead to MSE values 2 to 30 times larger (first and second column in Fig. 7).' • The corrected paragraph should read: ' The ablation study is reported in Fig. 8, showing that having four or more nonzero b values (or b shells), with two of them higher than 3,000 s/mm 2 (or equivalently 3 ms/ m 2 ), produces parameter estimates of reasonable accuracy and precision; but having less than two b values (or b shells) higher than 3,000 s/mm 2 (or equivalently 3 ms/ m 2 ) may be insufficient and can lead to MSE values 2 to 30 times larger (first and second column in Fig. 8).' This highlights further that the results in Figure 9 come ostensibly from fitting five free parameters to a data set with only four independent data points (four nonzero b-shells).Despite the apparent underdetermination, the fitting procedure still provides robust estimates of four parameters (D ec , f ec , f in , and r s ) because the data has little or no sensitivity to D in and the random-forest regressor implicitly fixes it to a value close to the average of the settings in the training data.An additional sensitivity analysis reported here in Figure X1

Fig. 1 .
Fig. 1.Schematics of current standard model of brain microstructure (a) and the novel model proposed in this work (b).Current conjecture envisions the tissue component in an MRI voxel as subdivided into two non-exchanging compartments: intra-neurite and extra-neurite space.The total MRI signal is then given by the weighted sum of the signals from water molecules diffusing in each compartment, with relative signal fractions f in and 1-f in , respectively (a).We propose a new picture: the tissue component of an MRI voxel is subdivided into intra-cellular and extra-cellular non-exchanging compartments.The total signal is the weighted sum of the signal from water molecules diffusing in each compartment, with relative signal fractions 1-f ec and f ec , respectively.Furthermore, the intra-cellular compartment is itself divided into two non-exchanging sub-compartments: intra-neurite and intra-soma.The intra-cellular MRI signal is then given by the weighted sum of the MRI signal from water molecules diffusing within the two sub-compartments, with relative signal fractions f in and 1-f in , respectively (b).

Fig. 2 .
Fig. 2. Summary and a few examples (a) of the 12 morphological features used in the generative model of brain cells generation introduced in (Palombo et al., 2019a) to simulate realistic cellular structures like Purkinje cells, motor neurons and pyramidal spiny neurons (b).Here, the generative model is used to investigate simplified cellular structures (c) comprised of straight long cylindrical fibres connected to a spherical soma structure, with and without the possibility for diffusing spins to exchange between neurites and soma.

Fig. 5 .
Fig. 5. Comparison between ground-truth model parameters and their estimates from fitting Eq. (10)without extra-cellular compartment to the simulated signals with exchange between neurites and soma in Fig. 4. Uncertainty (error bars) in parameter estimates was quantified by a Monte Carlo approach (see Methods section 3.1 for details).

Fig. 9 .
Fig. 9.An example of the parametric maps of the proposed compartment model for brain microstructure, obtained by fitting Eq. (10) to the normalised directionaveraged DW-MRI data from a representative subject.

Fig. 10 .
Fig. 10.Novel contrasts in apparent neurite and soma density of human brain.Qualitative comparison of MR soma signal fraction f is and MR cell fibers signal fraction f in maps with histological images (from different subjects) stained for brain cytoarchitectonic (Nissl staining for cell nuclei, left side) and myeloarchitectonic (myelin staining, right side), respectively.Brain histological images from the human brain atlas at https://msu.edu/~brains/brains/human.The contrast in the MRI maps show remarkable similarity to the contrast from histological staining.Quantitative map of f is is expected to provide contrast related to soma density, while map of f in is expected to provide contrast related to neurite density.Since myelinated neurites (like axons) are the major constituent of white matter, f in is expected to provide contrast highly related to myelin in the white matter regions.

Fig. 12 .
Fig.12.Projection onto cortical surface (representing the whole cortical thickness) of the quantitative map of MR soma signal fraction f is for the average over the whole cohort of 25 healthy subjects.In correspondence of the principal Brodmann areas we report histological images from literature, showing the typical soma arrangement and density used as criteria to delineate Brodmann areas.We find a very good correspondence between f is values and the expected pattern of soma density from histology (high, intermediate and low soma density as indicated).

Fig. 13 .
Fig. 13.Comparison of the SANDI model and the dot-compartment variant to describe measured DW-MRI signal decay at high b values.Two microstructure models were fitted to the data in Fig. 6: the Eq.(10) from the SANDI model, and the Eq.(12) from the dot-compartment variant.The results of the fit for the DW-MRI data from the WM and GM ROIs in Fig.6are reported in the table, together with the values of the AICc (lower AICc indicates better fit).