SPHERIOUSLY? The challenges of estimating spherical pore size non-invasively in the human brain from diffusion MRI

The Soma and Neurite Density Imaging (SANDI) three-compartment model was recently proposed to disentangle cylindrical and spherical geometries, attributed to neurite and soma compartments, respectively, in brain tissue. The approach could also enable estimation of microstructure parameters such as the apparent size (radius) of the soma. There are some recent advances in diffusion-weighted MRI signal encoding and analysis (including the use of multiple so-called ‘b-tensor’ encodings and analysing the signal in the frequency-domain) that have not yet been applied in the context of SANDI. In this work, using: (i) ultra-strong gradients; (ii) a combination of linear, planar, and spherical b-tensor encodings; and (iii) analysing the signal in the frequency domain, three main challenges to robust estimation of soma size were identified: First, the Rician noise floor in magnitude-reconstructed data biases estimates of soma properties in a non-uniform fashion. It may cause overestimation or underestimation of the soma size and density. This can be partly ameliorated by accounting for the noise floor in the estimation routine. Second, even when using the strongest diffusion-encoding gradient strengths available for human MRI, there is an empirical lower bound on the spherical signal fraction and pore-size that can be detected and estimated robustly. For the experimental setup used here, the lower bound on the signal fraction was approximately 10%. We employed two different ways of establishing the lower bound for spherical radius estimates in white matter. The first, examining power-law relationships between the DW-signal and diffusion weighting in empirical data, yielded a lower bound of 7 μm, while the second, pure Monte Carlo simulations, yielded a lower limit of 3 μm and in this low radii domain, there is little differentiation in signal attenuation. Third, if there is sensitivity to the transverse intra-cellular diffusivity in cylindrical structures, e.g., axons and cellular projections, then trying to disentangle two diffusion-time-dependencies using one experimental parameter (i.e., change in frequency-content of the encoding waveform) makes spherical pore-size estimates particularly challenging. We conclude that due to the aforementioned challenges spherical pore size estimates may be biased when the corresponding signal fraction is low, which must be considered when using them as biomarkers in clinical/research studies.


Introduction
Diffusion magnetic resonance imaging (dMRI) is a noninvasive technique widely used to study brain microstructure in vivo. Most dMRI methods are based on the conventional Stejskal-Tanner experiment (Stejskal and Tanner, 1965) that applies a pair of pulsed field gradients along a single axis for each signal preparation, which we refer to here as 'linear' encoding. Using linear encoding, disentangling different microstructural properties such as their size, shape, and orientation is far from trivial (Lampinen et al., 2017;Novikov et al., 2019). Such features may be entangled in the encoding process resulting in low specificity in their estimation. This is particularly problematic in dMRI where the image voxel is on the scale of a millimeter, and can therefore contain multiple microenvironments. Biophysical modeling is often used to tackle the inverse problem of inferring relevant tissue features (such as cell size, shape, and orientation) from the measured dMRI signal (Mitra et al., 1992;Wiegell et al., 2000;Stanisz et al., 1997;Zhang et al., 2012;Assaf et al., 2008). Most contemporary dMRI models for neural tissue share some common assumptions and features. First, they separate the tissue into intra-and extra-neurite compartments. Second, the exchange between the compartments is considered to be negligible, such that each compartment has a fixed and time-invariant signal fraction f i , where i f i = 1. Third, most models treat the intra-neurite compartment as a 'stick' -that is a compartment in which the diffusivity perpendicular to the long axis of the compartment is assumed to be effectively zero. This assumption is based on the lack of sensitivity to the neurite diameter . Different models have been used for the orientational dispersion of the compartments, represented by an orientation distribution function (ODF). Some models consider only purely parallel orientations (i.e., a delta function on the sphere ODF) (Stanisz et al., 1997;Assaf et al., 2008) while others use spherical harmonics (Jespersen et al., 2007) or a function such as the Watson distribution (Zhang et al., 2012) to characterise orientation dispersion. A Gaussian anisotropic representation is most often used for the extra-neurite compartment. Its orientation is determined by the mean of the fiber ODF and is characterized by axial and radial diffusivities. In white matter, the extra-neurite compartment largely corresponds to the extra-axonal space (and possibly extraastrocytic-process space). However, in gray matter, the extra-neurite compartment is comprised of the extracellular space (where the Gaussian representation is likely valid) and an additional contribution from the water restricted inside isotropic cellular structures, such as the cell bodies (soma) . Recent studies have shown that the two-compartment model is not a good representation of the signal in gray matter (Afzali et al., 2020a,c;McKinnon et al., 2017;Veraart et al., 2019;Palombo et al., 2018a;Henriques et al., 2019;Jespersen et al., 2019). This can be due to non-negligible water exchange processes occurring between intra-and extra-cellular compartments and between different intracellular compartments Jelescu and Novikov, 2020), or the assumption that the water inside the soma behaves the same as water in extracellular space. (Palombo et al., 2018a,b). Addressing this model insufficiency, Palombo et al. (Palombo et al., 2020) first demonstrated with non-trivial numerical simulations that, under specific experimental conditions, the contribution of soma to the total intracellular dMRI signal can be disentangled from that of neurite, and then introduced a three-compartment model called Soma And Neurite Density Imaging (SANDI), which decomposes the measured dMRI signal into three main sources: extra-cellular space, neurite and soma where, if there is any sensitivity to the size of that compartment, the diffusion MRI signal in each of these compartments will have a time-dependence. This enabled both MR apparent soma size and soma signal fraction to be estimated non-invasively using dMRI for the first time.
However, the inverse problem of using complex multicompartment models to infer microstructural information from the diffusion-weighted signal can be highly ill-posed and give rise to degeneracies in the model-parameter estimation, i.e., in the forward sense, completely different sets of model parameters predict the same dMRI signals (Jelescu et al., 2016;Novikov et al., 2018b;Jones et al., 2013;Lampinen et al., 2017Lampinen et al., , 2020Lampinen et al., , 2019. SANDI, like any other multi-compartment model (Jelescu et al., 2016), may suffer from the same degeneracy problems. In SDE, the MR signal is sensitized to diffusion using a pair of gradient pulses that encode the position of the spins along the axis defined by the diffusion gradients. Double Diffusion Encoding (DDE) contains two pairs of pulsedfield gradients that are separated from each other with a mixing time τ (Cory et al., 1990;Callaghan, 2011). This approach has been utilized by several groups for extracting microstructural information (Özarslan et al., 2009;Jespersen et al., 2013;Benjamini et al., 2014;Ianuş et al., 2016;Yang et al., 2018b;Coelho et al., 2019). A framework called q-space trajectory imaging (QTI) was recently introduced by (Westin et al., 2016) to probe tissue using different gradient waveforms. The traditional, pulsed field gradient sequences attempt to probe a point in q-space but in q-space trajectory encoding, time-varying gradients are used to probe a trajectory in q-space, and the b-matrix is defined as an axisymmetric second order tensor (Topgaard, 2017). In this framework, SDE is just a special realization of linear tensor encoding (LTE) where the b-tensor has only one non-zero eigenvalue as all gradients are along the same axis. DDE is a special case of planar tensor encoding (PTE) as all gradients lie on a plane and the b-tensor has two non-zero eigenvalues. In spherical tensor encoding (STE) the gradients may point in all directions giving rise to a rank-3 b-matrix. Recently, b-tensor encoding has been used to resolve degeneracy problem Coelho et al., 2019;Fieremans et al., 2018;Gyori et al., 2019;Lampinen et al., 2020). While these studies have shown improvement in the accuracy of parameter estimates, a consideration of timedependence of the diffusion-weighted signal can be used as another feature to add information. In particular, Gyori et al. (Gyori et al., 2019) have recently proposed a method based on a three-compartment model to estimate neurite and soma features (e.g. signal fractions and intracompartment apparent diffusivities) from combined LTE and STE data. Notably, Gyori et al. treated the signal coming from the spherical compartment as a simple monoexponential with a fixed, time-invariant small diffusivity. This implicitly assumes that the spherical compartment would show the same signal behavior for linear and spherical tensor encoding for a given b-value. However, Lundell et al. (Lundell et al., 2019) demonstrated that this is only true in restricted geometries when the LTE and STE waveforms have the same frequency power spectra. In this study, we used b-tensor encoding with variable power spectra (including LTE and STE waveforms that were not spectrally matched to each other) to investigate whether and how it improves fitting of the SANDI model. Here, we exploit all three forms of b-tensor, i.e., LTE, STE, and PTE. We adopted van Gelderen's model of the spherical compartment (Vangelderen et al., 1994), which explicitly includes both ∆ and δ. The challenge in using the free gradient waveforms, however, is that ∆ and δ are poorly defined, and so the time-dependency of the obtained signal is not well-defined in the time domain. Therefore, to find a closed-form for the diffusion-weighted signal decay in the spherical compartment, we adopted this model to the frequency domain Stepišnik, 1993;Lundell et al., 2019). The main findings of this paper are as follows: • Noise Sensitivity: Even when complementing LTE with STE-and PTE-data, fitting the spherical pore properties remains challenging (when the signal fraction is small, i.e. < 10%), with simulations showing biases in parameter estimates. Here we demonstrate that it is predominantly the Rician distribution of the noise (and associated noise floor) that impacts the estimation of spherical compartment properties. Such biases disappear when simulating purely Gaussian noise). However, if the Rician noise floor is accounted for in the model-fitting (albeit naively)much of the noise-floor induced bias is ameliorated.
• Lower Bound on Sphere Signal Fraction: Here, by using the F-statistic to compare nested models (i.e., those that do or do not include a sphere fraction) in simulated data where the spherical pore properties are varied systematically, it was possible to identify a lower bound on the detectable MRI signal fraction limit. This was around 10% for data with SNR = 50.
• Lower Bound on Spherical Pore Size: The empirical lower limit on spherical pore size in brain tissue was estimated by comparing exponents in powerlaw relationships between the dMRI signal and bvalue fitted to simulated data, with exponents observed empirically in vivo. We employed two different ways of establishing the lower bound for spherical radius estimates in white matter. The first, examining power-law relationships between the DW-signal and diffusion weighting in empirical data, yielded a lower bound of 7µm, while the second, pure Monte Carlo simulations, yielded a lower limit of 3µm. In addition, there is little differentiation in signal attenuation for low radii spherical pores (e.g. R sphere < 4µm, Fig. .3 (a)).
• Sticks + Ball + Sphere vs Cylinders + Ball + Sphere: In addition to the challenge of estimating sphere size, the fitting becomes even more challenging if we have cylinders instead of sticks. As shown by Veraart et al. (Veraart et al., 2020), at 300 mT/m, and with appropriate diffusion times, we have sensitivity to the internal perpendicular diffusivity in cylindrical pores (demonstrated by a break from a power-law relationship between signal intensity and b-value). A challenge then arises when trying to disentangle two time-dependencies by varying the same experimental parameter (i.e, changing the frequency content of the gradient-encoding waveform). With currently-available pipelines, this prevents reliable estimates of spherical pore properties in white matter when there is sensitivity to intraaxonal radial diffusivity, and indeed may plague grey matter modelling if there is sensitivity to water in the astrocytic processes. We should note that most of the cellular projections are smaller than 3 microns in radius, while the majority of soma are above 3 microns (Savtchenko et al., 2018;Di Benedetto et al., 2016;Zhang et al., 2016;Papageorgiou et al., 2011;Fannon et al., 2015;Mohamed et al., 2020). Therefore, the ambiguity here is more relevant to WM voxels given a low soma density there.

Theory
Multi-compartment models express the diffusionweighted signal as the sum of several compartments.
where f k is the signal fraction ( k f k = 1) and S k is the signal from the kth compartment. For a general B-matrix, the diffusion-weighted MR signal is modeled as: in (t))nn T + D ⊥ in (t)I, D ball and D sphere are the cylinder, ball and sphere signal fractions and diffusivities, respectively (Murday and Cotts, 1968). W (n) is the Watson orientation distribution function (ODF) and κ is the dispersion parameter. The diffusion weighting matrix B is given by (Eriksson et al., 2015;Westin et al., 2016), and γ is the gyromagnetic ratio. Axial and radial elements in the diagonal axisymmetric b-tensor are b || and b ⊥ respectively, b-value, b is the trace of B-matrix and b ∆ = (b || − b ⊥ )/b. For linear (LTE), planar (PTE) and spherical (STE) tensor encoding, b ∆ = 1, −1/2, and 0 respectively (Eriksson et al., 2015). To remove the effect of fiber orientation dispersion (Jespersen et al., 2013;Lasič et al., 2014), the acquired signal is averaged over all diffusion directions for each shell. This so-called 'powder-averaged' signal (Callaghan et al., 1979;Edén, 2003) has less complexity than the orientation-dependent signal, and yields a signal whose orientationally-invariant aspects of diffusion are preserved but with an orientation distribution that mimics complete dispersion of anisotropic structures. Compartmental diffusion is represented with axisymmetric diffusion tensors which are described by isotropic diffusivity, D I = 1/3D || + 2/3D ⊥ , and anisotropy, where D || and D ⊥ are the axial and radial diffusivities, respectively. D ∆ changes between −1/2 for a planar tensor to 1 for a stick. The signal attenuation from the kth compartment is given by (Lampinen et al., 2019;Eriksson et al., 2015): and erf(.) is the error function (Callaghan et al., 1979). Diffusion inside the sphere and ball is isotropic (D ∆;sphere = 0, D ∆;ball = 0) while for the cylinder and stick it is anisotropic (D ∆;cylinder > 0, D ∆;stick = 1). Therefore, the full signal equation is given by: 2.1. Two and three-compartment models 2.1.1. Cylinder + Ball + Sphere (Extended SANDI Model) In the original SANDI framework, the dMRI signal in brain tissue is assumed to arise from three main nonexchanging compartments: (i) intra-neurite (modeled as diffusion in sticks); (ii) intra-soma (modeled as diffusion constrained to a sphere); and (iii) extra-cellular (modeled as isotropic Gaussian diffusion). Here we extend this model to consider the perpendicular diffusivity in the intra-neurite compartment, D ⊥ cylinder (t), thereby modeling it with cylinders instead of sticks. We additionally explore the feasibility of modeling the intra-axonal perpendicular diffusivity and the additive spherical, D sphere (t), sensitivity simultaneously.
For complex gradient waveforms, the diffusion time is illdefined. We therefore consider the diffusion spectrum D ⊥ cylinder (ω), D sphere (ω) (Stepišnik, 1993;Lundell et al., 2019) in our analyses of compartment size. The restricted DW-signal inside the sphere and cylinder is S = exp(−ρ) where ρ is (Stepišnik, 1993;Lundell et al., 2019): is the gradient waveform. D(ω) can be expressed with a rotation matrix R as D(ω) = RΛ(ω)R −1 where Λ(ω) is the diagonal matrix containing diffusion spectra λ j (ω) along the restriction principal axes. The analytical expression for λ j (ω) in the case of restricted diffusion in planar, cylindrical and spherical geometries, is the weighted sum of negative Lorentzians: where for a cylinder, λ 1 (ω) = λ 2 (ω) = λ(ω) and λ 3 (ω) = 0 and where µ i are the roots of J 1 (µ i ) = 0 and J 1 (.) is the Bessel function of the first kind and order (Stepišnik, 1993;Nilsson et al., 2017;Lundell et al., 2019) and R c is the cylinder radius. For a sphere, λ 1 (ω) = λ 2 (ω) = λ 3 (ω) = λ(ω) and where µ i are the roots of the derivatives of the first order spherical Bessel function j 1 (µ i ) = 0 and R s is the sphere radius. D 0 is fixed at 3µm 2 /ms for the sphere, as proposed in  and D 0 = D || in for cylindrical geometry.
2.1.2. Stick + Ball + Sphere (R cylinder = 0 (Original SANDI Model) We define a three-compartment model, Stick + Ball + Sphere to investigate the sensitivity of the diffusion signal to the sphere radius and signal fraction. This model is the same as the original SANDI model with the difference that here we use b-tensor encoding and frequency-domain analysis. (Behrens et al., 2003) (f sphere = 0 and R cylinder = 0) In this section, we compare the Stick + Ball + Sphere model with a Stick + Ball model and provide the range of sphere signal fractions and radii that make these two models significantly different. This latter model is the simplest model and does not have any time dependency.

Method
Using both numerical simulations and in vivo experiments in healthy volunteers, we explore estimation of the soma size and density R sphere and f sphere using the combination of efficient gradient waveforms for LTE, PTE, and STE. Note: as we analyze the powderaveraged/orientationally-averaged signal, we do not need to estimate orientational dispersion. We study the challenges of the fitting landscape, the effect of noise, the lower limit on detectable signal fraction, the empirical lower limit on detectable pore size, and the challenge of disentangling two time-dependent properties (cylinder and sphere radius) of the model.

Fitting landscape
There is little differentiation in signal attention for low radii (R s < 4µm) ( Fig. .3(a)). When there is low sensitivity to some parameters, the numerical optimization algorithm terminates prematurely and therefore the estimates are not accurate.

Noise sensitivity
To explore the sensitivity of parameter estimation to noise perturbations, we simulated three different scenarios: (i) addition of Gaussian noise to the magnitude of the signal; (ii) addition of Gaussian noise to the real and imaginary channels which results in Rician-distributed magnitude signal; and (iii) addressing the noise-floor problem in case (ii) with a (simple) correction. In general, when there is Gaussian noise in the signal, averaging improves the signal to noise ratio (SNR) and because of the orientationalaveraging used in this work, we expect some improvement in the SNR in the first scenario and therefore better estimates of the model parameters. In the second scenario, the signal is corrupted by Rician-distributed noise, and therefore the orientational-averaging that improved the SNR in case (i), does not remove the non-zero positive 'noise-floor' bias in Rician-distributed noise and therefore we expect some bias in the parameter estimates. In the third scenario, we use a simple correction for the Rician bias in case (ii). To estimate the standard deviation of the noise, we include the noise floor in the model so that the predicted signal is S n = √ S 2 + σ 2 where S is our original model prediction and S n is the prediction after accounting for the noise floor (Jones and Basser, 2004;Koay et al., 2009;Pieciak et al., 2016bPieciak et al., , 2018Pieciak et al., , 2016a. We expect some improvement in the parameter estimates in case (iii) compared to case (ii) but the results of case (i) are expected to best out of all three cases.

Lower Bound on Resolvable Spherical Pore Signal
Fraction Here we considered that the framework had sensitivity to the sphere signal fraction if the inclusion of a sphere component to a stick + ball model was statistically supported by an F-test. To determine the lower bound on the spherical signal fraction and the pore size that can be detected using a diffusion-weighted signal, we systematically varied both parameters, while comparing the fit from two models: (i) the stick + sphere + ball model (threecompartments including a spherical component); and (ii) a stick + ball model (two-compartment without a spherical component). To test whether inclusion of the spherical compartment was needed to describe the signal (thereby showing sensitivity to this component) we considered the stick + sphere + ball model justified if the p-value from the F-test was less than 0.05 (Nilsson and Alexander, 2012;Panagiotaki et al., 2012). Here, the F-statistic is calculated where SSR is the sum of squared residuals, M is the number of fitted parameters of the simplified (1) and full SANDI model (2), and N is the number of measurements. The pvalue is estimated using p where fcdf is the cumulative distribution function of Fdistribution.

Stick + Ball + Sphere vs Cylinder + Ball + Sphere
If ultra-strong gradients and diffusion-time settings are such that we do, indeed, have sensitivity to the intraneurite perpendicular diffusivity, then the intra-neurite compartment should be more correctly modeled using cylinders instead of sticks (Veraart et al., 2020). This introduces an additional challenge, as we now have two compartments (sphere and cylinder) with a diffusion-time dependence. To explore this, we conducted further simulations to investigate the impact of including a non-zero perpendicular intra-neurite diffusivity (or, rather, a cylinder with a finite radius of R cylinder = 4 µm) on estimation of spherical pore-size.

Empirical Lower Bound on Spherical Pore Size
To identify the empirical lower bound on spherical pore size, we simulated signals for experiments with fixed diffusion time, ∆ = 37.05ms, and δ = 29.65ms (matching our in vivo experimental set-up), diffusivities D || in = 2 µm 2 /ms and D ball = 2 µm 2 /ms, but with variable sphere signal fractions f sphere = (0.01 : 0.01 : 0.1, 0.15, 0.2 : 0.1 : 1), f ball = f cylinder = (1 − f sphere )/2, and sizes, R sphere = 1 : 0.5 : 10 µm. For each set of signals, we fitted a power-law (Veraart et al., 2019;McKinnon et al., 2017;Lampinen et al., 2017) to the direction-averaged signal from the LTE measurements for b = 6, 7.5, 9, 10.5ms/µm 2 according to (S/S 0 = βb −α ) and then compared the values of the exponent, α, with values observed empirically in vivo to establish a lower bound on the size of spherical pores. The rationale behind the choice of using the power-law to drive an empirical conclusion is that it is free of any model assumptions, and simply considers the rate of signal decay versus b-value. We know that the α value for a pure stick-like geometry is 0.5, and thus any deviation from this value is indicative of sensitivity to an additional compartment (a deviation from the stick-like geometry could be due to any shape that is not stick-like). The compartment that we choose to change is, in fact, a spherical compartment. By systematically increasing the size of the spherical compartment until such a deviation is detected, we can obtain an empirical lower bound on the spherical compartment. Any sensitivity to the intra-axonal perpendicular diffusivity would make the signal decay faster (Veraart et al., 2020) but with the timing parameters used here, we do not expect any such sensitivity.

In vivo Data
Two healthy participants who showed no evidence of a clinical neurologic condition were scanned with the approval of the Cardiff University School of Psychology Ethics Committee. Magnetization-prepared rapid gradient echo (MPRAGE) images were also acquired for anatomical reference. 192 sagittal slices with TE = 2.3 ms, TR = 1900 ms, TI = 900 ms, and a voxel size of 1 mm isotropic and 256×256 matrix size were acquired in 5 minutes. Diffusion-weighted images were acquired with the protocol detailed in the simulation section 3.6 on a 3T Connectom MR imaging system with 300 mT/m gradients (Siemens Healthineers, Erlangen, Germany). Forty-two axial slices with 3mm isotropic voxel size and a 78×78 matrix size, TE = 88 ms, TR = 3000 ms, partial Fourier factor = 6/8, and heat dissipation limit = 1, were obtained for each individual. The total acquisition time was around one hour. To take full advantage of q-space trajectory imaging, it is imperative to respect the constraints imposed by the hardware, while at the same time maximizing the diffusion encoding strength. Sjolund et al. (Sjölund et al., 2015) provided a tool for achieving this by solving a constrained optimization problem that accommodates constraints on maximum gradient amplitude, slew rate, coil heating, and positioning of radiofrequency pulses. The gradient waveform is optimized and Maxwell-compensated (Szczepankiewicz et al., 2019) based on a framework that maximizes the b-value for a given measurement b-tensor and echo time. Substantial gains in terms of reduced echo times and increased signal-to-noise ratio can be achieved, in particular as compared with naive planar and spherical tensor encoding. Duration of the first, pause, and the second waveform in Fig. .

Preprocessing
The diffusion weighted images were corrected for Gibbs ringing (Kellner et al., 2016). We acquire some interleaved b0 images between the diffusion-weighted images (DWIs) to use for motion correction. In PTE and STE data, we registered the interleaved b0 images to the first b0 image and used the corresponding transformation to correct the motion in the DWIs. In LTE data the eddy current and subject motion were corrected by FSL EDDY (Andersson and Sotiropoulos, 2016) and finally the gradient nonlinearity was corrected by the method proposed by Rudrapatna et al. . We applied a 3D Gaussian filter with a standard deviation of 0.5 to the preprocessed data to make the images smooth. We normalized the direction-averaged signal based on the b = 0 s/mm 2 signal in each voxel.

Regions of interest
We defined five regions of interest (ROIs) in the splenium and internal capsule as white matter regions and putamen, ventrolateral thalamus, and mediodorsal thalamus, as gray matter regions. To define anatomical regions of interest, publically-available atlases were coregistered to each participant's T 1 -weighted MPRAGE image. The white matter ROIs were obtained from the JHU atlas (Mori et al., 2005), while the putamen ROI was obtained from the work of Tziorti et al. (Tziortzi et al., 2011) and ventrolateral thalamus and mediodorsal thalamus from (Danos et al., 2003). The MPRAGE image was co-registered to the diffusion-weighted data, and the resulting transform applied to the ROIs to translate them to the native diffusion-weighted space. The five ROIs are illustrated in Fig. .1.

Goodness of fit for in vivo data
To check the stability of the model fit and that the global minima of the cost function had been found, we first fixed the signal fraction of the spherical compartment and estimated the remaining parameters (The same procedure was used by Lampinen et al. (Lampinen et al., 2019) for finding the stick fraction that can be detected reliably). To assess the precision of parameter estimation, a metric of the goodness-of-fit (see below) was plotted for different values of sphere signal fraction, which was varied systematically between zero and one in 40 equal steps. If the model determined all parameters unequivocally, a clear optimum in the goodness-of-fit would be seen for some signal fractions. Conversely, a flat plot of the goodness-of-fit over a wide range of signal fractions would indicate degeneracy in the fitting, i.e. two or more sets of solutions yield a similarly good fit. For each ROI, all the voxels contained therein were concatenated to provide a sufficient number of data points for our estimation. Goodness-of-fit was determined using reduced chi-square, χ 2 red (Andrae et al., 2010). The reduced chi-square or normalized residual variance (NRV) (Lampinen et al., 2019) was obtained by dividing the residual variance (σ 2 R ) by the variance of noise (σ 2 noise ): where n i is the number of directions for ith b and b ∆ . S i and S i are the direction-averaged measured and predicted signals, n is the number of samples and k is the number of free parameters in the model. The standard deviation of the noise (σ 2 noise ) was estimated for each ROI.

Simulations
Fig. .2 (a) shows the gradient waveforms used for the linear, planar and spherical tensor encoding and their corresponding frequency spectra. Fig. .2 (b) shows the signal decay inside the spherical and cylindrical compartments using different encoding schemes. It is important to note that, for a given b-value, STE results in the most signal loss, followed by PTE and then LTE. LTE appears relatively insensitive for small radii spherical compartments. Fig. .3 (a) shows the changes in apparent diffusivity of the spherical compartment (D sphere ) as a function of radius of the spherical compartment for the three b-tensor shapes. Clearly, there is a large difference in the sensitivity to pore size, with LTE being the least sensitive and PTE and STE tracking each other closely in the plot of D sphere vs R sphere . Notably, for all wave-forms, there is little differentiation in sphere signal attenuation for low radii, (e.g. R sphere < 4µm). For small radii, the noise dominates over measurable effects. The error bars in Fig. .3 (b) show the confidence interval. In the fitting, we use lsqcurvefit in MATLAB which returns the jacobian at the solution, and then nlparci command is used to find the 95% confidence interval. Fig. .4 shows the results of fitting the sphere radius for different sphere signal fractions under different noise simulations. The figure also shows the p-value of the Ftest between the two and three-compartment models in the presence of Gaussian, Rician, and corrected Rician noise. Here, we take p < 0.05 as an indication that the full model (three-compartment) is preferred over the simplified model (two-compartment). When the sphere radius or the signal fraction of the sphere is small (R sphere < 2µm and f sphere < 0.05) the simplified model is preferred. Fig. .4 shows that the estimates are largely positively biased in the Rician-noise data, whereas their errors are symmetrically-distributed about the line of identity in Gaussian only. After Rician noise correction, it looks more like the Gaussian -in terms of symmetrical distribution around the line of identity.

Signal fraction resolution limit
The examination of the F-test measures indicates a lower bound on the signal fraction that can be detected or reasonably modeled from the diffusion-weighted signal. This is around 10% for SNR = 50 but this limit changes at different noise levels (Fig. .4). The changes of sphere radius estimates versus ground truth for a wider range of SNR values are shown in this link: https: //bit.ly/StickBallSph. The lower limit scales inversely with SNR, i.e. as SNR increases, smaller sizes are detectable.

Stick + Ball + Sphere vs Cylinder + Ball + Sphere
In addition to the challenge of estimating sphere size (given the above 4 points), the fitting becomes more challenging if we have cylinders instead of sticks because in this model the time-dependence of the signal can now arise from two independent sources (the cylinder and the sphere). Fig. .5 shows the estimated sphere size when there is a non-zero diameter cylinder, (SNR = 200). For signal fractions smaller than 0.2 and sphere radius larger than 5 microns, the estimated sizes deviate from the ground truth. This behavior is observed where the noise floor is very low (SNR = 200) in the simulated signal. The deviation from the ground truth in the presence of noise, for different noise levels, is shown in the following link: https://bit.ly/CylBallSph.

Empirical lower limit on pore size
In addition to the lower bound imposed by the noise floor on the sphere signal fraction and size, there is an empirical lower limit on the pore size which can be derived from power-law measurements. Fig. .6 shows the effect of sphere signal fraction and size on the estimated exponent, α, in the power-law fit (S/S 0 = βb −α ) experiments. As will be discussed below, in our in vivo data, we observe power-law exponents greater than 0.5 in the white matter. Fig. .6 shows that to observe such values of α the sphere radius needs to be larger than 7 microns. We will discuss the implications of this key result in the context of the in vivo data below. A good 'sanity check' is that after reconstructing the signal from the estimated parameters the exponent α in the power-law fit (S/S 0 = βb −α ) of the original and the modeled signal should not change considerably (see the last row of Fig. .7), otherwise, the parameters are representing the signal incorrectly. This approach has several limitations that should be acknowledged; The first limitation is that we used a pre-determined and fixed set of parameters in the simulation. If we had used a smaller diffusivity for the extra-cellular compartment, for the range of b-values used here (6 < b < 10.5 ms/µm 2 ), there might have been some residual contribution from the extracellular compartment and therefore the exponent α and the behaviour of the signal would change. The second limitation is that we fixed the intra-neurite and extracellular signal fractions to be equal which may not reflect reality. Despite these limitations, we consider the result a useful empirical benchmark because the fixed diffusivities used in this experiment are close to the range of diffusivities estimated in previous works . As such, the contribution from the extracellular compartment at high b-values will be negligible meaning that the only remaining signal contributions come from spheres and sticks.

In vivo results
Here we provide in vivo results from fitting the Stick + Ball + Sphere model. Fig. .7 shows the results of fitting the model to the signal by fixing the sphere signal fraction to different values. Results are shown for five different ROIs: Splenium; internal capsule; mediodorsal thalamus; ventrolateral thalamus; and putamen. Note that, for the fitting, the data points from all three voxels in the ROI are concatenated to provide enough data points for a stable fit, but in the figure, the average of these three voxels is shown. The flat valleys in χ 2 red correspond to plausible sphere signal fractions in different ROIs of the brain, for example, 0-0.125 for Splenium and Internal Capsule and 0.2-0.3 for the Putamen (Fig. .7). When the valley is flat for a large range of sphere signal fractions, it means the data do not provide a unique solution and a range of parameters can represent the signal equally well. This may be related to a large range of acceptable parameters as shown by the second to the fourth row of Fig. .7. For example, acceptable D || in values for the stick compartment ranged between 2.2 and 2.5 µm 2 /ms in the splenium and between 1.2 and 1.5 µm 2 /ms in the putamen. If there was no degeneracy in the estimation of parameters and the data could provide useful information about the underlying microstructure, then the plot would have a sharp valley at the local optimum. Among the five ROIs we selected in this experiment, we see a quite sharp minimum for the putamen and mediodorsal thalamus (with f sphere around 0.2-0.3), and the minima for white matter (splenium and internal capsule) clearly puts the sphere fraction in the sub 10% regime, which is where it is expected. The first column of Fig. .7 shows the parameter estimates for a synthetic signal generated with the following parameters; f sphere = 0.5, f ball = f stick = 0.25, D || in = 2 µm 2 /ms, D ball = 0.6 µm 2 /ms, R sphere = 5 µm, and SNR = 100, Rician distributed signal. There is a sharp minimum in χ 2 red which shows there is only one set of parameters that fits the signal accurately. Fig. .8 shows estimated parameter maps in vivo. When fitting on a voxel by voxel-level, the fitting is unstable and the resulting maps are not smooth. Besides, the large voxel size (3 mm) used here and the resulting problems with partial volume may affect the estimated parameters as well. Maps of the estimated parameters in Fig. .8 show a reasonable contrast that matches the results in ). The f stick map has higher values in white matter tracts in the brain while f sphere values are higher in the gray matter. The f stick values in cortical GM range from 0.1 to 0.2 which is in agreement with recent works on estimating neurite density in GM using b-tensor encoding (Lampinen et al., 2017(Lampinen et al., , 2019. .10 shows the results of fitting a power-law (S/S 0 = βb −α ) to the diffusion-weighted signal. CSF has the fastest decay and therefore no signal remains from CSF in the high b-value data and the alpha is close to zero. The decay in the GM is faster than the WM and therefore the estimated α values are correspondingly larger in GM than WM. Within WM, α is usually larger than 0.5. As noted by Veraart et al. (Veraart et al., 2019) and McKinnon et al. (McKinnon et al., 2017), if there is only a stick compartment, then α = 0.5. The decay faster than the stick model may come from the exchange between compartments (Stanisz et al., 1997), sensitivity to the axon diameter (Veraart et al., 2020) or the presence of a non-negligibly-sized third compartment  that makes the signal decay faster than 0.5. If the radius of the spherical compartment is less than approximately 7 microns then the exponent is smaller than 0.5. As will be seen in the map of α in Fig. .6, we do not observe values of α less than 0.5 in the white matter. This places a lower bound on the size of spherical pores in white matter of around 7 µm.

Discussion
The SANDI model  extends existing multi-compartment models that only consider two pools of water in brain tissue (Jespersen et al., 2007;Zhang et al., 2012;Fieremans et al., 2011;Kaden et al., 2016;Novikov et al., 2018aNovikov et al., , 2019Alexander et al., 2019). Palombo et al. (Palombo et al., 2018a) suggested that the failure of the stick model in gray matter (McKinnon et al., 2017;Veraart et al., 2019) can be due to the abundance of cell bodies (namely soma). In previous works, the contribution from soma was considered as part of the extracellular space (Jespersen et al., 2007;Jespersen, 2012) because the exchange between the restricted water in soma and the hindered water in the extracellular space was assumed to be fast. However, recent studies (Yang et al., 2018a) suggest that the pre-exchange time of intracellular water in neurons and astrocytes is 500 ms. These findings prompt the conclusion that for diffusion times much smaller than 500 ms (e.g.
10-20 ms) the exchange between intra and extra-cellular water may be negligible, supporting the idea that the signal from spins restricted in soma may be non-negligible. The typical size of soma ranges from a few microns (for microglia and glia) to a few tens of microns (for big neurons) (Zhao et al., 2006;Díaz-Cintra et al., 2004). In this study we found that based on the power-law fit to the signal in voxels containing white matter, the sphere radius is larger than 7µm (Fig. .6). While at first glance, one might not expect spherical compartments of this size in white matter in vivo, there are studies that report the size of astrocytes (Savtchenko et al., 2018;Di Benedetto et al., 2016;Zhang et al., 2016;Papageorgiou et al., 2011) and oligodendrocytes (Fannon et al., 2015;Mohamed et al., 2020) and confirm that having effective soma radius > 7µm is possible (data also publicly available on Neuromorpho.org). This suggests the presence of glial cells in the central nervous system whose soma can have such a size. Although in healthy conditions glial volume fraction in white matter has been estimated to be low (5% − 20% (Veraart et al., 2020;Coelho et al., 2018;Duval et al., 2016)), the presence of large glial cells may become more relevant in pathological conditions, where glia may migrate and react undergoing hypertrophy of the whole cellular structure (Ligneul et al., 2019;Genovese et al., 2020). In this work, we studied the minimal sphere signal fraction and radius that can be detected from diffusion MRI of water inside spheres. Our finding is in agreement with the results reported by Dhital et al. (Dhital et al., 2018) (2% of the unweighted signal for moderate diffusion times using Prisma scanner) and Tax et al. (Tax et al., 2020) (isotropic signal fraction of 9.7% for the apparent diffusivity of 0.12µm 2 /ms). It should be noted that the protocol in this study is not specifically optimized for size estimation. For example, intentionally varying the frequency spectra (e.g. to include high frequency components), might result in better sensitivity to smaller pore-sizes (Molendowska et al., 2020). The estimation of stick signal fraction in (Lampinen et al., 2019) was challenging while in (Lampinen et al., 2020) Lampinen et al. provided reliable estimates of the two-compartment model parameters because the sequence was optimized in the latter. Comparing the estimation of sphere signal fraction in this work with the stick signal fraction in (Lampinen et al., 2019(Lampinen et al., , 2020 we conclude that to estimate the parameters accurately, the acquisition should be optimized toward the estimation of the model parameters. Here we used a combination of linear, planar, and spherical tensor encoding to ameliorate the degeneracy problem that exists in the fitting of multi-compartment models. Nilsson et al. (Nilsson et al., 2017) reported that for the estimation of diameter in complete orientation dispersion (which we effect by powder-averaging the signal), from an SNR perspective it is advantageous to use oscillating gradient spin echo (OGSE) compared to standard SDE. The benefits of double diffusion encoding (DDE) for size estimation have been presented in other studies (Benjamini et al., 2014;Katz and Nevo, 2014;Vincent et al., 2020). Our results suggest using multiple waveforms provides the best estimates. In the case of cylinder diameter estimation, the resolution limit is determined by the amount of attenuation due to radial diffusion. This attenuation is estimated by the integral of the gradient squared and can be maximized by either a fat-pulse SDE or a rectangular oscillating pulse. However, when the long axis of the cylinder is not perpendicular to the direction of the applied gradient, the high b-values should be avoided because of the signal attenuation and decrease in SNR. To improve the SNR and the resolution limit for cylinder diameter estimation, waveforms should have more oscillations and hence lower b-values . Here, we are targeting the sphere diameter, and therefore OGSE or SDE can be both useful. In the estimation of spherical pore size from the diffusionweighted signal, different confounding factors have to be considered. One of the challenges is that for small sphere radii (< 3 micron), the fitting landscape is flat and there is a negligible change in the signal for small sphere sizes (Fig. .3). Noise is another confounding factor that affects the estimation of parameters in both model-based and signal representation based techniques. Parameters obtained from multi-compartment models, (the stick + ball + sphere model in this paper) applied to noisy data are biased because of the effect of noise. Three different noise scenarios were simulated here: Gaussian noise, Rician noise, and corrected Rician noise. If the data were corrupted with purely Gaussian noise, then this could be removed to some extent through the act of computing the orientationally-averaged signal. However, as we invariably use the magnitude-reconstructed data, the noise has a Rician distribution, which presents a more challenging scenario because averaging does not remove the bias. This effect is more pronounced when there is a small contribution from the spherical compartment. We wish to stress that the challenges we identified are mainly relevant to WM. The SANDI model was developed mainly for soma imaging in GM. Nothing in our results suggests that SANDI is unreliable in GM, and will indeed, benefit from the multi-waveform and frequency-domain approach presented here. The challenge in the gray matter is to determine whether the deviation from the 'stick' model comes from the soma compartment, exchange between compartments (Jelescu and Novikov, 2020) or both. We note that the extracellular signal fraction obtained from the fit in the gray matter areas is around 0.45 which is higher than expected. This discrepancy can be explained by three factors: first, T 2 relaxation is not explicitly accounted for in our model, and second, we consider the same proton density in both intra-axonal and extra-axonal spaces. And also non-negligible partial volume with CSF, particularly problematic for cortical GM at 3 mm 3 resolution. The model presented here is only sensitive to relative signal amplitudes while differences in T 2 relaxation can impose different weights to the amplitudes of the pools (Szafer et al., 1995;Callaghan, 1995). The specific assignment of nerve water T 2 components by simultaneously considering compartmental diffusion and transverse relaxation properties was already studied by Peled et al. (Peled et al., 1999) in myelinated nerves of the frog sciatic nerve tissue. More and more evidence is given for the T 2 relaxation time constants of intra-and extra-axonal water to be different from each other in case of slow exchange between the intra-axonal and extra-axonal pools (Peled et al., 1999;Dortch et al., 2010). In the case of fast exchange (Szafer et al., 1995), the signal loss in all the pools would be weighted in the same manner. Although it is relatively easy to incorporate relaxation into the model and fit the experimental data with additional parameters (relaxation rates), such a model results in an unstable fit with the current protocol. To incorporate additional parameters (i.e., relaxation rates) in the model we need to obtain additional information in our experiment. For example, repeating the acquisition at different echo times as suggested by Lampinen et al. (Lampinen et al., 2020). We highlight six limitations of the current study in the estimation of spherical compartment signal fraction and size. First, we assume that the intra-axonal water comes from straight axons, which is not the case in most of the white matter voxels Nilsson and Alexander, 2012;Jeurissen et al., 2013;Lee et al., 2020;Özarslan et al., 2018). Second, the extra-cellular component is modelled as a ball with isotropic diffusion. This assumption is not valid when there are coherently-oriented fibers (such as in the midline of the corpus callosum), where diffu-sion in the extra-axonal space can have a high anisotropy (Özarslan et al., 2011). Third, at low frequencies, the timedependency of diffusivity in the extracellular space can dominate over the time-dependency in the intra-axonal space (Burcaw et al., 2015;Nilsson et al., 2017). Fourth, it is assumed that the exchange between water environments is negligible. This assumption might be valid since previous studies have shown the exchange times in the white matter are of the order of seconds or longer (Nilsson et al., 2013) which is much larger than the time-scales that the effects of restricted diffusion can be observed. Fifth, in our tissue model, we have neglected the distribution of restricted dimensions (e.g. range of soma sizes, axon diameters). However, adding extra parameters to the model to account for this will make the fitting unstable. Six, the effect of T 2 relaxation is not considered in our model which may result in bias in the estimation of the other parameters of the model. In addition, the lack of cerebrospinal fluid (CSF) component in the model is another limitation of this work. Future directions. The fitting method used in this work is a nonlinear least-square fit that can be replaced with new deep learning approaches to improve the quality of fit (Gibbons et al., 2019). From the acquisition perspective, the protocol that is used in this study is not optimized for the estimation of small sizes. Using a range of frequency spectra will help (Molendowska et al., 2020). The protocol, used in this work, imposes a long acquisition time which can be minimized by optimizing the directions as well as the number of shells. In this paper, simple arithmetic averaging is used for powder averaging which can be replaced with some other techniques such as weighted averaging (Knutsson et al., 1999;Szczepankiewicz et al., 2017;Afzali et al., 2020b) to obtain a better orientationinvariant signal that improves the parameter estimates.

Conclusion
In this work, we have demonstrated key challenges and limitations in estimating spherical pore size non-invasively in the human brain from diffusion MRI. Our simulations show the effect of Rician bias on the estimation of pore size and identified the lower bound limit of the sphere signal fraction and size that can be detected from the diffusion-weighted signal from both an SNR and empirical perspective. We showed that for small signal fraction of soma, i.e. < 10%, this is a problem. However, we know from detailed microscopy of brain cortex (Beul and Hilgetag, 2019;Collins et al., 2010), that in GM the soma signal fraction is on average > 20%. Therefore, reliable estimation of soma properties in GM is possible, while in WM it presents several challenges. The flat landscape of the fitting was also investigated. Using the ultra-strong gradients of the Connectom scanner, the diffusion signal in the white matter can be made sensitive to the axon diameter, and therefore the three-compartment model of stick + ball + sphere changes to cylinder + ball + sphere which has two time-dependencies, one for the diffusivity in the sphere and the other one for the diffusivity in the cylinder. Disentangling these two time-dependencies using only one sequence parameter (i.e., changing the frequency content of the encoding waveform) in the acquisition is challenging. Studying all these challenges prevents misinterpretation of the biased estimated parameters as the potential biomarkers in clinical studies.

Acknowledgements
The data were acquired at the UK National Facility for In Vivo MR Imaging of Human Tissue Microstructure funded by the EPSRC (grant EP/M029778/1), and The Wolfson Foundation. MA and DKJ are supported by a Wellcome Trust Investigator Award (096646/Z/11/Z) and a Wellcome Trust Strategic Award (104943/Z/14/Z). MN was supported by grants from Swedish Research Council (2016-03443) and Random Walk Imaging AB (grant no. MN15). The authors would like to thank Filip Szczepankiewicz for providing the pulse sequences for btensor encoding. We thank Lars Mueller for setting up the protocol for b-tensor encoding. MP is supported by UK EPSRC EP/N018702/1 and UKRI Future Leaders Fellowship MR/T020296/1.   ). The third row shows the reduced chi-square, χ 2 red , values for two scenarios where the sphere radius is fixed to 5 and 8 µm, blue and red curves respectively. Note that there is a sharp minimum for R = 5µm while the χ 2 red curve has several local minima or even flat in the R = 8µm case (GT = ground truth and E = estimated).   Figure .7: The results of fitting the stick + ball + sphere model to the diffusion-weighted signal by fixing the sphere signal fraction to different values. Five different ROIs of the brain are used here; putamen, internal capsule, mediodorsal thalamus, ventrolateral thalamus, and splenium. The mean value of the direction-averaged signal for each ROI is represented in the first row (in different columns). The second row shows the estimated signal fraction of stick (f stick ) and ball (f ball ) for different predefined sphere signal fractions (We used f as y-label here to show the signal fraction of both ball and stick in one plot.). The third row illustrates the parallel diffusivity of the stick (D || in ) and the diffusivity of the ball (D ball ) for different ROIs. The estimated radius of the sphere is illustrated in the fourth row. And finally, the last two rows show how well this model can explain the signal for different predefined sphere signal fractions (f sphere ) in terms of reduced chi-square and power-law. The first column in the figure shows the results of fitting a synthetic signal. Note that we do not estimate diffusivity of the compartment when the signal fraction is estimated as zero, this is the reason for discontinuity in the plots of estimated diffusivities. , ball (f ball ), and sphere (f sphere ) signal fractions, intra-axonal parallel diffusivity (D || in (µm 2 /ms)), extracellular diffusivity (D ball (µm 2 /ms)), and sphere radius (R sphere (µm)) on axial, sagittal, and coronal views of the smoothed brain image (A 3D Gaussian kernel with standard deviation of 0.5 is used for smoothing). , and sphere (f sphere ) signal fractions, intra-axonal parallel diffusivity (D || in (µm 2 /ms)), extracellular diffusivity (D ball (µm 2 /ms)), sphere radius (R sphere (µm)), and standard deviation of the noise (σ) on axial, sagittal, and coronal views of the smoothed brain image (A 3D Gaussian kernel with standard deviation of 0.5 is used for smoothing).   Figure S1: Estimated stick (f stick ), ball (f ball ), and sphere (f sphere ) signal fractions, intra-axonal parallel diffusivity (D || in (µm 2 /ms)), extracellular diffusivity (D ball (µm 2 /ms)), and sphere radius (R sphere (µm)) on axial, coronal, and sagittal views of the in vivo brain image . : Estimated stick (f stick ), ball (f ball ), and sphere (f sphere ) signal fractions, intra-axonal parallel diffusivity (D || in (µm 2 /ms)), extracellular diffusivity (D ball (µm 2 /ms)), sphere radius (R sphere (µm)), and standard deviation of the noise (σ) on axial, coronal, and sagittal views of the in vivo brain image .