Rutger H J Fick
email: rutger.fick@inria.fr

Demian Wassermann

Emmanuel Caruyer

Rachid Deriche

MAPL: Tissue Microstructure Estimation Using Laplacian-Regularized MAP-MRI and its Application to HCP Data

Keywords: Diffusion-Weighted MRI, Human Connectome Project, Microstructure Recovery, Laplacian-Regularized MAP-MRI, Axonal Diameter, Axonal Dispersion MAPL Extrapolation Isotropic Scaling Anisotropic Scaling Average peaks: 2.3 ∓ 0.6 Unregularized Positivity Constraint Laplacian Regularization Radial Order / Moment: 6 / 0 Average angle: 52.9 0 ∓ 51.6 0 Average peaks: 2.4 ∓ 2.3 Radial Order / Moment: 6 / 0 Average angle: 29.6 0 ∓ 53.7 0 Average peaks: 1.3 ∓ 2.4

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The recovery of microstructure-related features of the brain's white matter is currently a hot topic in diffusion MRI [START_REF] Le Bihan | Imagerie de diffusion in-vivo par résonance magnétique nucléaire[END_REF][START_REF] Taylor | The spatial mapping of translational diffusion coefficients by the nmr imaging technique[END_REF][START_REF] Merboldt | Self-diffusion nmr imaging using stimulated echoes[END_REF]. In many cases, the accurate estimation of these features requires the acquisition of data at large b-values. However, in practice the maximum b-value is limited and noise begins to dominate the signal at higher b-values. To robustly estimate these important features from noisy, sparsely sampled data we propose to analytically regularize the coefficient estimation of Mean Apparent Propagator (MAP)-MRI method ( Özarslan et al., 2013b) using the norm of the Laplacian of the reconstructed signal. We call this approach MAPL.

MAPL falls into the category of functional basis approaches in diffusion MRI, which in general assume the narrow pulse approximation. This allows them to conveniently interpret the signal as the EAP through a Fourier transform [START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF][START_REF] Callaghan | Principles of nuclear magnetic resonance microscopy[END_REF]. This Fourier relation lies at the heart of techniques such as q-space imaging [START_REF] Callaghan | Principles of nuclear magnetic resonance microscopy[END_REF], q-ball imaging [START_REF] Tuch | Q-ball imaging[END_REF] and diffusion spectrum MRI [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF]. However, the numerical implementations of these techniques require dense acquisition schemes to reconstruct the EAP. This makes them impractical for clinical applications, where scanning time is limited.

The introduction of functional bases to efficiently represent the dMRI signal partly overcame this restriction [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]Özarslan et al., 2013a;[START_REF] Assemlal | Efficient and robust computation of pdf features from diffusion mr signal[END_REF][START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF][START_REF] Caruyer | Diffusion mri signal reconstruction with continuity constraint and optimal regularization[END_REF]Özarslan et al., 2013b;[START_REF] Hosseinbor | Bessel fourier orientation reconstruction (bfor): An analytical diffusion propagator reconstruction for hybrid diffusion imaging and computation of q-space indices[END_REF][START_REF] Rathi | Multi-shell diffusion signal recovery from sparse measurements[END_REF]. That is, the fitting of these representations can be regularized using properties such as the smoothness and positivity of the EAP. This allows for more accurate EAP and tissue property recovery using fewer samples, resulting in faster acquisition schemes.

A number of microstructure-related scalar indices have been proposed based on functional basis approaches. They are typically expressed as integral operators of the diffusion signal and EAP, which are in turn expressed in terms of fitted functional basis coefficients [START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF]Özarslan et al., 2013a;[START_REF] Hosseinbor | Bessel fourier orientation reconstruction (bfor): An analytical diffusion propagator reconstruction for hybrid diffusion imaging and computation of q-space indices[END_REF][START_REF] Rathi | Multi-shell diffusion signal recovery from sparse measurements[END_REF]. Examples are the Return-To-Origin Probability (RTOP), Mean Squared Displacement (MSD) and Q-space Inverse Variance (QIV). However, the estimation of directional q-space indices, whose direction depends on the white matter orientation, have only been proposed for the MAP-MRI basis ( Özarslan et al., 2013b). Under the hypothesis that coherent white matter can be modeled as parallel cylindrical pores, the directional Return-To-Axis Probability (RTAP) can be used to recover moments of the axon diameter distribution ( Özarslan et al., 2013b). This was verified using parallel, water-filled microcapillaries in 1D-NMR experiments ( Özarslan et al., 2011, 2013a). Though, histology studies have shown that axonal dispersion is present even in coherent white matter, and one must be careful when assuming parallel axons [START_REF] Leergaard | Quantitative histological validation of diffusion mri fiber orientation distributions in the rat brain[END_REF][START_REF] Ronen | Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of n-acetylaspartate and post-mortem histology[END_REF] Still, For the EAP to be sensitive to particle movements as small as axon diameters -between 0.2 µm and 2 µm [START_REF] Aboitiz | Fiber composition of the human corpus callosum[END_REF][START_REF] Liewald | Distribution of axon diameters in cortical white matter: an electron-microscopic study on three human brains and a macaque[END_REF] -large diffusion gradients with q-values of over 200 mm -1 are needed. Using a typical diffusion time and pulse length of 30 ms and 10 ms, this corresponds to a b-value of over 40, 000 s/mm 2 . Yet, even the most recent data sets of the MGH Human Connectome Project have a significantly lower maximum b-value of 10, 000 s/mm 2 [START_REF] Setsompop | Pushing the limits of in vivo diffusion mri for the human connectome project[END_REF]. The accurate estimation of microstructure-related parameters in MAP-MRI therefore strongly depends on the smooth extrapolation of the diffusion signal beyond the largest measured b-value, where the SNR is lowest. To impose this smoothness, we propose to analytically regularize MAP-MRI's coefficient estimation using the norm of the Laplacian of the reconstructed signal. Throughout this work, we will refer to our Laplacian-regularized MAP-MRI as MAPL.

Figure 1: Visual abstract. On the left, we illustrate multi-shell dMRI, where diffusion-weighted images (DWIs) are measured along different gradient directions and b-values. DWIs at higher b-values (using strong gradients) provide information on smaller details in the tissue geometry, but also have lower SNR. We show that using Laplacian-Regularized MAP-MRI (MAPL) we can accurately estimate valuable microstructure features from noisy, subsampled data (bottom row).

We compare MAPL with the initially proposed positivity constraint on the EAP ( Özarslan et al., 2013b), which we refer to simply as "MAP", and the recently proposed modified Spherical Polar Fourier (mSPF) basis [START_REF] Caruyer | Diffusion mri signal reconstruction with continuity constraint and optimal regularization[END_REF]. We also include comparisons between MAP-MRI's Cartesian and spherical implementations, where the latter is equivalent to the previously proposed 3D-SHORE basis [START_REF] Merlet | Continuous diffusion signal, eap and odf estimation via compressive sensing in diffusion mri[END_REF]. In terms of signal fitting, our comparison will carefully disentangle the effects of using different basis orders (i.e. using different numbers of coefficients) and different regularization techniques and basis implementations. In addition, we also study the effect of estimating higher "radial moments" of the ODF in fiber crossing estimation. We perform these experiments on data from a physical phantom with known fiber configurations and gold standard data [START_REF] Ning | Sparse reconstruction challenge for diffusion mri: Validation on a physical phantom to determine which acquisition scheme and analysis method to use?[END_REF][START_REF] Moussavi-Biugui | Novel spherical phantoms for q-ball imaging under in vivo conditions[END_REF].

We then study the effect of data subsampling on the estimation of the scalar indices MSD, QIV, RTOP, RTAP, and RTPP using Cartesian MAPL. For this we use the WU-Minn Human Connectome Project data [START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF] as it has been acquired using an incremental sampling scheme [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion mri[END_REF]. We also illustrate an application of MAPL where we use RTAP to estimate the apparent axon diameter in the Splenium, Midbody, and Genu of the Corpus Callosum in six subjects of the MGH Human Connectome Project data [START_REF] Setsompop | Pushing the limits of in vivo diffusion mri for the human connectome project[END_REF].

Here "apparent" indicates that the axon diameter estimated through RTAP is only valid when the signal originates only from inside a set of parallel axons. In fact, contrary to RTAP's hypotheses, we know that the signal does not only originate from inside the axons and, in most parts of the white matter, axons are not parallel, but dispersed [START_REF] Leergaard | Quantitative histological validation of diffusion mri fiber orientation distributions in the rat brain[END_REF][START_REF] Ronen | Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of n-acetylaspartate and post-mortem histology[END_REF].

As an alternative, tissue models separate the signal contributions of different tissue types using biophysical models. This allows for the estimation of microstructural parameters such as tissue composition, axonal diameter [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF][START_REF] Alexander | A general framework for experiment design in diffusion mri and its application in measuring direct tissue-microstructure features[END_REF][START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion mri[END_REF] and axonal dispersion [START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF]. However, the variance of the estimated parameters in these models still depends on the noise and maximum gradient strength in the data. To reduce the variance of the parameter estimation, we propose to use MAPL's signal fitting and extrapolation as a preprocessing step before fitting a tissue model. This preprocessing consists of 1) fitting the signal using MAPL; 2) extrapolating signal from the fitted representation beyond the original maximum b-value; and 3) fitting the tissue model on the extrapolated signal. We study its effect on the estimation of axon diameter using a simplified Axcaliber model [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF] and axonal dispersion using the Neurite Orientation Dispersion and Density Imaging (NODDI) model [START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF]. We will show the effect of this preprocessing on the estimation variance by estimating these parameters with and without preprocessing in the Corpus Callosum of six different subjects of the MGH Human Connectome Project. Finally, we correlate the estimated axon diameter, dispersion and restricted volume fractions with Fractional Anisotropy (FA). We show that changes in FA correlate significantly with changes with all estimated parameters. All contributions of this work are illustrated in Figure 1. For reference, we write out all relevant acronyms used in this paper in the Appendix in Table A.4, including in which section we explain their meaning.

The paper is organized as follows: We first provide the theory of the MAP-MRI basis, our Laplacian regularization (MAPL), estimated scalar indices and competitive methods in Section 2. In Section 3 we explain the implementation of our method and the data we use in our experiments. Then in Section 4 we present our results on both the physical phantom and the HCP data. We then discuss our results in Section 5 and finally provide our conclusions in Section 6.

Theory

In this section, we first provide a brief explanation on the relationship between the measured dMRI signal and the ensemble average propagator (EAP) in Section 2.1. We then introduce the MAP-MRI basis, our proposed Laplacian regularization and regularization weight optimization in Section 2.2. In Section 2.3 we then describe how we estimate the radial moments of the EAP and recover microstructure related scalar indices. Finally, we give an overview of methods we compare against in Section 2.4.

The Relation between the diffusion signal and the EAP

In dMRI the EAP is reconstructed by first obtaining diffusion-weighted images (DWIs). In a pulsed gradient spin echo (PGSE) sequence, a DWI is obtained by applying two sensitizing diffusion gradients of length δ to the tissue, separated by pulse separation time ∆. The effective diffusion time is then given by τ = ∆ -δ/3. When the narrow pulse approximation is met, i.e., when the diffusion gradients are considered infinitely small (δ ≈ 0), the relation between the measured signal S(q, τ ) and the EAP P (R, τ ) is given by an inverse Fourier transform (IFT) [START_REF] Tanner | Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method[END_REF][START_REF] Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF]:

P (R, τ ) = R 3 E(q, τ )e i2πq•R dq with q = γδG 2π (1)
where E(q, τ ) = S(q, τ )/S 0 is the normalized signal attenuation measured at wave vector q and S 0 is the baseline image acquired without any diffusion sensitization (q = 0). The wave vector q on the right side of Eq. ( 1) is related to the nuclear gyromagnetic ratio γ and the applied diffusion gradient vector G. We denote q = qu and R = Rr, with u, r ∈ S 2 and q, R ∈ R + . The b-value is related to q and the effective diffusion time as b = 4π 2 q 2 τ . As a consequence of the Fourier relationship in Eq. (1), measuring E(q, τ ) at larger q makes one sensitive to smaller movements in P (R, τ ). Using a longer diffusion time τ gives the diffusing particles more time to interact with the tissue boundaries, making the EAP more specific to the tissue structure [START_REF] Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF][START_REF] Cory | Measurement of translational displacement probabilities by nmr: an indicator of compartmentation[END_REF].

In the rest of this work we assume the diffusion time is constant throughout dMRI acquisitions and will omit the τ dependence in our notation. In the next section we explain our method for reconstructing the EAP from the signal.

MAPL: Laplacian-Regularized MAP-MRI

We improve EAP and q-space index estimation by regularizing the coefficient estimation of the MAP-MRI functional basis ( Özarslan et al., 2013b) with the Laplacian of the reconstructed signal. We call our Laplacian-regularized MAP-MRI approach MAPL for brevity. In dMRI, this type of regularization has successfully been applied to several other techniques [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Caruyer | Diffusion mri signal reconstruction with continuity constraint and optimal regularization[END_REF]. We provide a brief summary of the MAP-MRI basis in Section 2.2.1, present our proposed Laplacian regularization in Section 2.2.2, and explain the generalized cross-validation (GCV) method to find optimal regularization weights in Section 2.2.3.

MAP-MRI Functional Basis Description

MAP-MRI is a functional basis that reconstructs the EAP from the dMRI signal. It has the convenient property that its basis functions are eigenvectors of the Fourier transform [START_REF] Walter | Properties of hermite series estimation of probability density[END_REF]. For this reason, it can describe both the signal and the EAP with the same coefficients as

E(q) = i c i Φ Ni (q) and P (R) = i c i Ψ Ni (R) (2)
where Ψ = IFT(Φ) and {c i } are the basis coefficients. The basis functions for the signal attenuation and the EAP (left and right side of Eq. ( 2)) are given as products of three orthogonal onedimensional basis functions, which are known as Simple Harmonic Oscillator based Reconstruction and Estimation (SHORE) functions ( Özarslan et al., 2013a,b). For the signal attenuation they are given as

Φ n1n2n3 (A, q) = φ n1 (u x , q x )φ n2 (u y , q y )φ n3 (u z , q z ) with φ n (u, q) = i -n √ 2 n n! e -2π 2 q 2 u 2 H n (2πuq). ( 3
)
and for the EAP as

Ψ n1n2n3 (A, R) = ψ n1 (u x , R x )ψ n2 (u y , R y )ψ n3 (u z , R z ) with ψ n (u, R) = 1 √ 2 n+1 πn!u e -R 2 /(2u 2 ) H n (R/u) (4)
where H is a Hermite polynomial of order n and ψ = IFT(φ). The fitting of the MAP-MRI signal basis is then not performed in the regular "image space" -that does not change from voxel to voxel -but in voxel-dependent "anatomical space" ( Özarslan et al., 2013b). In this anatomical space the q-space vectors of the data (and therefore the data itself) has been rotated such that the main eigenvectors of an estimated DTI tensor [START_REF] Basser | Mr diffusion tensor spectroscopy and imaging[END_REF] coincide with the axes of the Cartesian coordinate system {x, y, z}. The eigenvalues of the DTI tensor are then used to set the stiffness of the oscillator in each direction (i.e. {u x , u y , u z }). In this way the anisotropy of MAP-MRI's basis functions depends on the anisotropy of the data as estimated by a DTI tensor.

From this point on we define image space variables with an added accent ( ) and their corresponding anatomical space variables without it. To fit the MAP-MRI basis we first obtain the anatomical space DTI tensor D from the image space DTI tensor D as D = VD V T where V is the orthonormal rotation matrix. We recover the image space q-space vectors as q = Vq . We then estimate the covariance matrix of displacements as A = 2Dτ with τ the diffusion time [START_REF] Basser | Relationships between diffusion tensor and q-space mri[END_REF], and obtain the scaling factors as A = Diag(u 2

x , u 2 y , u 2 z ). The basis series is then fitted up to a maximum radial order N max . For a real, symmetric propagator (as is the case in dMRI), N max is even-valued and the total number of basis coefficients is N coef = 1 6 (F + 1)(F + 2)(4F + 3) with F = N max /2. Using N max = 0 reduces the MAP-MRI basis to a DTI fit of the data, while using higher N max 'corrects' the initial DTI approximation to the true shape of the data.

The data-dependent anisotropic scaling from A lies at the heart of MAP-MRI. However, if the scaling is chosen isotropic (u x = u y = u z ) the MAP-MRI basis corresponds to the 3D-SHORE basis [START_REF] Merlet | Continuous diffusion signal, eap and odf estimation via compressive sensing in diffusion mri[END_REF] ( Özarslan et al., 2013b, Appendix). The 3D-SHORE basis is thus a particular case of MAP-MRI where the bases in Eqs. ( 3) and ( 4) have been rewritten as a product of a radial oscillator and angular spherical harmonics. In the experiments of this work we will compare the Cartesian (anisotropic) and spherical (isotropic) implementations of MAP-MRI, but in the theory here we will only consider the 'regular' Cartesian implementation. We provide the derivations of scalar measures and Laplacian regularization for the spherical implementation in Appendix C.

Closed-Form Laplacian Regularization

In this section, we provide the formulation of our proposed Laplacian-regularized MAP-MRI. We call this approach MAPL for brevity. In MAPL, we fit MAP-MRI's basis functions to noisy data by first casting the coefficients into an N coef -dimensional vector c and the signal values in an N datadimensional vector y. Design matrix Q ∈ R N data ×N coef then has elements Q ij = Φ Ni (A, q j ). The coefficients c (recall Eq. ( 2)) are found by minimizing the quantity c = argmin c y-Qc 2 +λ ∆ U (c) where λ ∆ weights our Laplacian regularization functional

U (c) = R 3 ∆E c (q) 2 dq (5)
with ∆E c (q) = i c i ∆Φ Ni (q) the Laplacian of the reconstructed signal. We then express U (c) as a summation of MAP-MRI basis functions

U (c) = R 3 i c i ∆Φ Ni (q) 2 dq = i k c i c k R 3 ∆Φ Ni (q) • ∆Φ N k (q) dq (6)
where the subscripts i and k indicate the basis order of the i-th or k-th basis function Φ Ni (q) = Φ n1n2n3(i) (q). We write the summations in quadratic form such that U (c) = c T Uc where regularization matrix U has elements

U ik = R 3 ∆Φ Ni (q) • ∆Φ N k (q)dq. ( 7
)
It follows that the Laplacian of Φ Ni (q) with respect to q corresponds to the Laplacian of its orthogonal components φ n (q, u) along q x , q y and q z . The equation for the elements of U can therefore be solved using the general differential equation whose solutions form the functional basis functions φ n of the MAP-MRI basis ( Özarslan et al., 2012, Eq. (17))

- 1 (2πu) 2 ∂ 2 ∂q 2 + (2πuq) 2 φ n (q, u) = λ n φ n (q, u) (8)
with λ n = 2n + 1. Inverting this equation we show that ∆φ n (q, u) = (2πu) 2 (2πuq 2 ) -λ n φ n (q, u).

We solve Eq. ( 7) by inserting Eq. ( 9) and using the orthogonality of H n with respect to the weighting function e -x 2 on [-∞, ∞]. This allows us to describe every entry in U as a function of the anisotropic scaling factors (u x , u y , u z ) and basis orders (x i , y i , z i ) = (n x(i) , n y(i) , n z(i) ) as

U ik = u 3 x u y u z S x k xi U y k yi U z k zi +2 u x u y u z T x k xi T y k yi U z k zi + u 3 y u z u x S y k yi U z k zi U x k xi +2 u y u z u x T y k yi T z k zi U x k xi (10) + u 3 z u x u y S z k zi U x k xi U y k yi +2 u x u z u y T x k xi T z k zi U y k yi
with functions S m n , T m n and U m n given as

S m n = 2(-1) n π 7/2 δ m n 3(2n 2 + 2n + 1) + δ m n+2 (6 + 4n) m!/n! + δ m n+4 m!/n! + δ m+2 n (6 + 4m) n!/m! + δ m+4 n n!/m! (11)
T m n = (-1) n+1 π 3/2 δ m n (1 + 2n)

+ δ m+2 n n(n -1) + δ m n+2 m(m -1) (12) U m n = δ m n (-1) n /(2π 1/2 ) (13)
with δ m n the Kronecker delta. Note that regularization matrix U is symmetric, mostly sparse and its elements depend only on the data-dependent scale factors {u x , u y , u z } and basis orders. Also note that S m n , T m n and U m n do not depend on the scale factors and can be precomputed for a given N max . Using this formulation we obtain the regularized MAP-MRI coefficients using penalized least squares with unique minimum

c = (Q T Q + λ ∆ U) -1 Q T y. (14)

Optimal Weighting Parameter Choice

To find optimal regularization weights λ ∆ in Eq (14) we use the Generalized Cross Validation (GCV) algorithm [START_REF] Craven | Smoothing noisy data with spline functions[END_REF][START_REF] Koay | A signal transformational framework for breaking the noise floor and its applications in mri[END_REF]. GCV is based on an N data -fold cross validation. Fortunately, the estimation of λ can be calculated as the minimum argument of the GCV function

GCV (λ, y) = y -ŷλ N data -T r(S λ ) ( 15
)
where

S λ = Q(Q T Q + λU) -1
Q T is the smoother matrix and ŷλ = S λ y.

Figure 2: The effect of increasing the radial moment s on a 90 degree and 35 degree fiber crossing. Here s = -2 represents Tuch's ODF and s = 0 the marginal ODF. As s increases further both profiles become sharper, where it can be seen that the 35 degree crossing can only be resolved at higher radial moments.

Estimation of EAP-based Microstructure Parameters

We can estimate EAP-based features of the tissue structure directly from the fitted MAP-MRI coefficients ( Özarslan et al., 2013b). To estimate the orientation of the white matter tissue we resume the formulation of the radial moments of the EAP in Section 2.3.1. To relate the EAP to microstructure properties such as the apparent axon diameter, we also provide the formulation for the estimation of several boundary cases of the EAP in Section 2.3.2. Finally, we describe two multi-compartment tissue models that we use to estimate the actual axon diameter and axonal dispersion in Section 2.3.3.

Radial Moments of the EAP

The orientation of the underlying tissue can be estimated by marginalizing the EAP P (Rr) with respect to its radius R. Any maxima in the resulting angular profile are likely to coincide with the orientation of the underlying axon bundles. The general equation to marginalize the EAP is given as ODF s (r) = R 2+s P (R r)dR.

(16)

This marginalization can be expressed in MAP-MRI coefficients in closed form by inserting the right side of Eq. (2) into Eq. ( 16) and evaluating the integral. Here s denotes the s th radial moment of the propagator. Setting s = -2 gives Tuch's Orientation Distribution Function (ODF) [START_REF] Tuch | Q-ball imaging[END_REF] while setting s = 0 gives the marginal ODF [START_REF] Aganj | Reconstruction of the orientation distribution function in single-and multiple-shell q-ball imaging within constant solid angle[END_REF] where the integral over the sphere is unity. Increasing s induces a "sharpening" effect on the ODF along the fiber directions, as can be seen in Figure 2. As we will illustrate in Section 4.1.2, using higher radial moments can be used to more precisely recover the directionality of the white matter tissue.

Boundary Cases of the EAP

The EAP can be related to the mean tissue sizes, e.g. mean volume, mean cross-sectional area and mean length, by modeling the tissue as pores [START_REF] Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF]. Examples of these are parallel cylinders for aligned axon bundles and spherical pores for cell bodies and astrocytes. Given the diffusion time is long enough for spins to travel the longest end-to-end space of the tissue, the relation between the EAP and these mean tissue sizes is given through so-called boundary cases of the EAP ( Özarslan et al., 2011, 2013b).

Here we consider three particular cases: (1) The Return-To-Origin probability (RTOP), which is related to the mean pore volume as RT OP = V -1 , (2) the Return-To-Axis probability (RTAP), which is related to the mean pore cross-sectional area as RT AP = A -1 and (3) the Return-To-Plane probability (RTPP), which is related to the mean pore length as RT P P = L -1 ( Özarslan et al., 2013b). Though, as the modeled cylinders are not capped, the long diffusion time condition is hard to fulfill in the case of RTPP, meaning its value should be close to that of unrestricted diffusion. These three boundary cases are defined on the EAP as

RTOP P (0) (17) RTAP R P (R r |r ⊥ = 0)dR (18) RTP P R 2 P (R r ⊥ |r = 0)dR. ( 19
)
Unlike RTOP, the RTAP and RTPP are directional scalar indices that assume the white matter tissue is modeled by parallel cylinders, with r parallel and r ⊥ perpendicular to the cylinder axis.

Because of this assumption, their values are only related to the mean tissue sizes when the assumed cylinder axis coincides with the orientation of the underlying axon bundle. In this work, this orientation is estimated using the principal eigenvector of the diffusion tensor [START_REF] Basser | Mr diffusion tensor spectroscopy and imaging[END_REF], but can also be estimated using any other method that estimates this feature. The closed-form expressions for RTOP, RTAP, RTPP can be found in ( Özarslan et al., 2013b), but we newly provide the Mean Squared Displacement (MSD) and q-space Inverse Variance (QIV) [START_REF] Wu | Computation of diffusion function measures in-space using magnetic resonance hybrid diffusion imaging[END_REF][START_REF] Hosseinbor | Bessel fourier orientation reconstruction (bfor): An analytical diffusion propagator reconstruction for hybrid diffusion imaging and computation of q-space indices[END_REF] in Appendix D. Of these scalar indices, the RTAP is especially interesting as its related tissue size, the mean cross-sectional area A , can be related to the mean diameter as D = 2 1/(RT AP π). However, The validity of the relation between RTAP and the mean diameter D depends on several important assumptions ( Özarslan et al., 2013b;[START_REF] Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF]:

1. The tissue can be modeled as an ensemble of parallel cylinders. 2. The diffusion signal originates only from the intra-axonal compartment. 3. The applied gradient pulse length is infinitesimally small (δ ≈ 0). 4. The diffusion time τ is long enough for diffusing particles to traverse the longest end-to-end space (τ R 2 /D) with D the diffusivity.

Though, we know that the first two tissue-based assumptions are not met: The in-vivo diffusion signal originates from both the intra-and extra-axonal space [START_REF] Novikov | Revealing mesoscopic structural universality with diffusion[END_REF][START_REF] Burcaw | Mesoscopic structure of neuronal tracts from time-dependent diffusion[END_REF] and even in the Corpus Callosum we must account for axonal dispersion [START_REF] Leergaard | Quantitative histological validation of diffusion mri fiber orientation distributions in the rat brain[END_REF][START_REF] Ronen | Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of n-acetylaspartate and post-mortem histology[END_REF]. We are therefore careful to call the estimated value from RT AP the "apparent axon diameter" (AAD). Furthermore, when (3) and ( 4) are not exactly met an underestimation of the AAD can be expected [START_REF] Bar-Shir | The effect of the diffusion time and pulse gradient duration ratio on the diffraction pattern and the structural information estimated from q-space diffusion mr: experiments and simulations[END_REF]. Despite these practical limitations, we demonstrate an application of this approach by estimating the AAD in the Corpus Callosum in Section 4.3, where we do find AAD values characteristic of trends also found in histology [START_REF] Aboitiz | Fiber composition of the human corpus callosum[END_REF][START_REF] Liewald | Distribution of axon diameters in cortical white matter: an electron-microscopic study on three human brains and a macaque[END_REF]. We interpret our findings using AAD in the discussion.

Signal Extrapolation as Preprocessing for Multi-Compartment Models

To estimate properties of the white matter tissue -such as the axon diameter or dispersion -it is necessary to separate the signal contributions that originate from the intra-or extra-axonal parts of the tissue. To this end, multi-compartment models have had a profound impact on microstructure imaging in diffusion MRI, starting from a simple ball-and-stick model [START_REF] Behrens | Characterization and propagation of uncertainty in diffusion-weighted mr imaging[END_REF] to more complex models such as Neurite Orientation Dispersion and Density Imaging (NODDI) [START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF], AxCaliber [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF] and many others [START_REF] Panagiotaki | Compartment models of the diffusion mr signal in brain white matter: a taxonomy and comparison[END_REF].

In this work we propose to use MAPL as preprocessing to improve microstructure estimation using multi-compartment models. The preprocessing consists of (1) fitting the signal using MAPL, (2) extrapolating signal from fitted representation beyond the original maximum b-value and ( 3) fitting the tissue model on the extrapolated signal. Similarly as the smoother matrix in Section 2.2.3, we generate an 'extrapolation' matrix as S extra = Q extra (Q T Q + λU) -1 Q T and recover the extrapolated signal as ŷextra = S extra y.

We study the effect of our preprocessing on two different compartmental models that estimate either the mean axon diameter or axonal dispersion. To estimate the mean axon diameter we use a simplified Axcaliber model [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF], adapted for data that was acquired with one diffusion time (See Appendix E). To estimate axonal dispersion in terms of the Orientation Dispersion Index (ODI) we use the Neurite Orientation Dispersion and Density Imaging (NODDI) model [START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF]. In both models we fit an intra-axonal signal E ic (q), extra-axonal signal E ec (q) and isotropic signal E iso (q) with volume fractions ν iso + ν ic + ν ec = 1. Also in both cases E iso (q) is an isotropic Gaussian with free water diffusivity and E ec (q) is an axially symmetric, anisotropic Gaussian that in the case of NODDI is also dispersed. Following the same formulation as in [START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF], we describe the tissue using three components as

E(q) = (1 -ν iso )(ν ic E ic (q) + (1 -ν ic )E ec (q)) + ν iso E iso (q) (20)

Comparison with State-of-the-Art

We compare MAPL with competing techniques. In Section 2.4.1 we describe the previously proposed positivity constraint for MAP-MRI. in Section 2.4.2 we describe a competing functional basis approach to reconstruct the EAP.

Positivity Constrained MAP-MRI

We compare our regularization method with an approach that enforces positive-definiteness of the EAP ( Özarslan et al., 2013b). Here the signal fitting is cast in a quadratic programming framework which imposes a positivity constraint on the EAP on a discrete set of points while constraining the discrete integral of the EAP to unity. This approach is elegant as positivity and unity of the EAP are intrinsic properties of the EAP, but as the constraints require repetitive dense sampling of the EAP it is computationally expensive.

Laplacian-Regularized mSPF

As a functional basis approach that has a similar regularization technique, we compare our approach with the recently proposed Laplacian-regularized modified Spherical Polar Fourier (mSPF) basis [START_REF] Caruyer | Diffusion mri signal reconstruction with continuity constraint and optimal regularization[END_REF]. The mSPF basis, unlike MAP-MRI, uses decoupled radial and angular basis functions. This means the maximum radial order N max and angular order A max can be chosen independently, and the total number of estimated coefficients is given by N coef = (N max -1)(A max + 1)(A max + 2)/2. It also regularizes the fitting of its coefficients using the Laplacian of the reconstructed signal. We use GCV to find optimal weighting parameters for the mSPF Laplacian functional.

Materials and Methods

Implementation

Our numerical implementations of MAPL, MAP, and multi-compartment tissue models are written to be entirely compatible with the open source DiPy framework [START_REF] Garyfallidis | Dipy, a library for the analysis of diffusion mri data[END_REF]. All results were generated on an Intel(r) Core(TM) i7-3840QM CPU 2.80GHz computer with 32GB RAM. For brevity here we describe the algorithmic fitting procedure for both MAPL and the compartmental models in Appendix B.

Data Set Descriptions 3.2.1. SPARC Phantom Data

In our first experiment, we compare our MAPL approach to the other methods with respect to signal reconstruction and crossing angle recovery. We use the data of the sparse reconstruction challenge hosted at the 2014 CDMRI workshop on computational diffusion MRI [START_REF] Ning | Sparse reconstruction challenge for diffusion mri: Validation on a physical phantom to determine which acquisition scheme and analysis method to use?[END_REF]. This data was acquired from a physical phantom with known fiber configurations providing a gold standard reference data. The phantom is made of polyfil fibers of 15µm diameter and was developed along the same lines as in [START_REF] Moussavi-Biugui | Novel spherical phantoms for q-ball imaging under in vivo conditions[END_REF]. It consists of two fiber bundles crossing at a 45 degree angle with isotropic diffusion outside of the fiber bundles. The gold standard of the diffusion signal was obtained by acquiring 81 different gradient directions at b-values {1000,2000,3000,4000,5000} s/mm 2 averaged over 10 repetitions. The challenge data we use has 30 gradients directions per b-shell at b-values {1000,2000,3000} s/mm 2 , resulting in only 90 samples plus one b 0 , with an average SNR of 9.5 over all directions and voxels. The in-plane resolution was 2 × 2 mm 2 with slice thickness of 7 mm. Other acquisition parameters were T E/T R = 141/3400 ms and δ ≈ ∆ = 62 ms.

WU-Minn Human Connectome Project Data

In our second experiment we investigate the effects of subsampling the data on the estimation of the microstructure-related scalar indices described in Section 2.3.2. We use the WU-Minn Human Connectome Project data [START_REF] Moeller | Multiband multislice ge-epi at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fmri[END_REF][START_REF] Feinberg | Multiplexed echo planar imaging for sub-second whole brain fmri and fast diffusion imaging[END_REF][START_REF] Setsompop | Blippedcontrolled aliasing in parallel imaging for simultaneous multislice echo planar imaging with reduced g-factor penalty[END_REF][START_REF] Xu | Highly accelerated whole brain imaging using aligned-blipped-controlled-aliasing multiband epi[END_REF][START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF][START_REF] Glasser | The minimal preprocessing pipelines for the human connectome project[END_REF]. In this dataset the diffusion directions were obtained such that every subset of the first M directions is still isotropic (HCP-Manual, 2014, p.46) [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion mri[END_REF]. The data was sampled on 3 shells with b-values {0,1000,2000,3000} s/mm 2 , with {14, 90, 90, 90} directions, respectively. Following the sampling scheme, we subsample this data by truncating last N diffusion weighted images from the data, producing data sets with either 180, 90 or 60 gradient directions. The diffusion time and pulse separation time in this data are δ/∆ = 10.6/43.1 ms with 2 × 2 × 2 mm resolution and T E/T R = 89.5/5520 ms.

MGH Adult Diffusion Human Connectome Project Data

In our final experiments we test the robustness of our methods when estimating the AAD, axon diameter and axonal dispersion in different parts of the Corpus Callosum for different maximum b-values. We use the MGH Adult Diffusion Data of the Human Connectome Project [START_REF] Greve | Accurate and robust brain image alignment using boundary-based registration[END_REF][START_REF] Andersson | A comprehensive gaussian process framework for correcting distortions and movements in diffusion images[END_REF][START_REF] Keil | A 64-channel 3t array coil for accelerated brain mri[END_REF][START_REF] Setsompop | Pushing the limits of in vivo diffusion mri for the human connectome project[END_REF]. This data set was acquired at particularly high b-values {0, 1000, 3000, 5000, 10000} s/mm 2 with {40, 64, 64, 128, 256} directions, respectively. We use 6 different subjects of this data, having the HCP database numbering {1007, 1010, 1016, 1018, 1019, 1030}, which we from now on call subjects 1 through 6.

(c) Optimal λ ∆ [-]
Figure 3: The estimated Mean Diffusivity (MD), Fractional Anisotropy (FA) and optimal regularization weight for the physical phantom. It can be seen that white MD differs for different fiber configurations, the FA is consistent in fibrous areas. The optimal regularization weight using GCV also shows similar values in fibrous areas, but much higher in the empty areas.

The diffusion time and pulse separation time in this data are δ/∆ = 12.9/21.8 ms with 1.5 × 1.5 × 1.5 mm resolution and T E/T R = 57/8800 ms.

Experiments and Results

Signal Fitting and ODF Reconstruction on Physical Phantom Data

We compare our MAPL approach with MAP on the SPARC phantom data in two ways: (1) to reconstruct the normalized signal values at the gold standard q-space positions and (2) to accurately recover the known 45 degree crossing angles using the ODF. An illustration of the ODFs of the gold standard data is given in Figure 7. For completeness, we also compare the MAPL and MAP implementations with both the 'regular' anisotropic scaling factors and with isotropic scaling factors.

Figure 3 shows scalar maps of the Fractional Anisotropy (FA), Mean Diffusivity (MD) and optimal Laplacian regularization weight λ ∆ . The MD shows clearly different diffusivities in crossing fiber areas (blue-green), single fiber areas (red) and empty areas (dark blue). In fibrous areas (single and crossing) the FA shows very similar FA (0.86 ± 0.03) and lower in empty areas (0.6 ± 0.19). Finally, for a radial order of 8, we find fairly consistent Laplacian weights λ ∆ in fibrous areas (0.04 ± 0.02) and much higher in the empty areas (0.45 ± 0.11).

Effect of Regularization and Radial Order on Signal Reconstruction

We compare signal reconstruction quality between the different methods over different radial orders of the MAP-MRI basis. In Figure 4 we show the mean squared error (MSE) of the signal recovery for single, crossing and isotropic voxels, along with the computation time for every technique. The dashed lines indicate the use of isotropic scaling, the solid lines indicate anisotropic scaling and the gray dotted line is the unregularized reference. We use the acronyms MAPL and MAP to indicate the use of either Laplacian Regularization or Positivity Constraint, respectively.

Figures 4a and4b show that the MSE values in single and crossing voxels are similar. For radial orders lower than 8 the anisotropic scaling (purple and magenta) results in lower MSE values than isotropic scaling (green and blue). For higher radial orders both approaches have similar MSE. In terms of regularization, MAPL has significantly lower MSE than MAP for isotropic scaling, especially in crossing voxels, while this difference is more subtle for anisotropic scaling. However, we will show in the next section that ODF reconstruction differs significantly between MAPL and MAP, despite the similarity in MSE. We also study the fitting error of the interpolated points on b ={0, 1000, 2000, 3000} s/mm 2 compared to the extrapolated points on b ={4000, 5000} s/mm 2 in the fibrous areas of the phantom. We do this by separately estimating the MSE for either the interor extrapolated points. Without showing the graphs, we find for any method and any radial order that the MSE on the extrapolated points is generally between 2-3 times higher than the MSE on the interpolated points, which corresponds to similar ratios found in the original challenge [START_REF] Ning | Sparse reconstruction challenge for diffusion mri: Validation on a physical phantom to determine which acquisition scheme and analysis method to use?[END_REF].

Figure 4c shows that in empty areas the differences in MSE depend only on the choice of using either MAPL or MAP, where the latter has lower MSE values for higher radial orders. Although, none of the combinations give a lower MSE than a simple DTI fit (radial order 0) as diffusion is expected to be isotropic and Gaussian here.

Figure 4d shows the computation time of each method. We show that our MAPL is 3 to 4 times faster than MAP because of the iterative nature of the positivity constraint. We also show that the isotropic implementations are faster than the anisotropic ones as no rotation of the q-space vectors is required.

As a functional basis approach that has a similar regularization approach as ours, we compare MAPL with the recent Laplacian-regularized modified Spherical Polar Fourier (mSPF) basis [START_REF] Caruyer | Diffusion mri signal reconstruction with continuity constraint and optimal regularization[END_REF], which we outlined in Section 2.4.2. We now consider both the single and crossing fiber voxels together. Figure 5 shows the MSE and number of estimated coefficient for different maximum radial and angular orders of the mSPF basis (in shades of red) and our best MAPL result (in green), which was using MAPL at radial order 8. It can be seen that for no matter how many basis functions mSPF uses, it cannot approximate the signal as well as MAPL.

Effect of Radial Order, Radial Moment and Regularization on ODF Reconstruction

In this section, we compare the recovery of the phantom's 45 degree crossing over both radial order of the basis N max and radial moment s of the EAP. Recalling Eq. ( 16) and Figure 2, increasing the radial moments of the EAP yields sharper angular profiles from which to recover fiber bundle orientations. As in the original challenge [START_REF] Ning | Sparse reconstruction challenge for diffusion mri: Validation on a physical phantom to determine which acquisition scheme and analysis method to use?[END_REF], we compute the crossing angle by first finding the three largest peaks in the angular profile of the radial moment and compute the angles between them. The crossing angle closest to 45 degrees is registered as the detected crossing angle. When no crossing is found the voxel is ignored. We only include voxels where we know there is a crossing in the comparison. With respect to these metrics, an optimal regularization technique should have a good estimation of both the angle and the number of peaks.

Figure 6 shows intensity plots indicating the average estimated crossing angle (top row) and average found number of peaks (bottom row). Starting with the anisotropic methods (left two columns), we see that the crossings start to be recovered from a N max of 4, but the crossing angle is significantly underestimated. Increasing N max reduces this bias for our MAPL, but barely improves or even worsens using MAP for N max higher than 6. On the other hand, MAP does regulate the number of peaks better than MAPL for higher N max . In both cases increasing the radial moment improves the number of peaks, but typically worsens the underestimation of the angle.

For the isotropic methods (right two columns) the crossings start to be recovered after an N max = 6 for MAPL, and after N max = 8 for MAP. Here we do not find a consistent angle underestimation compared to the anisotropic methods. Increasing the radial moment for MAPL improves the peak estimation, but going over 2 increases the number of spurious peaks. For MAP, we find a similar trend, but a higher radial order must be chosen to accurately estimate the crossing angle with the correct number of peaks.

We visualize the radial moments of the EAP in Figure 7. We choose the optimal radial order and radial moment for every method based on Figure 6 and report the averages and standard deviation of the estimated angle and number of peaks in the top-right corners. On the top row, we show the reconstructed golden standard data using anisotropic scaling (left) or isotropic scaling (right). In these reconstructions we still needed to slightly regularize the ODF estimation (Laplacian weight 0.02) to eliminate some spurious peaks. Comparing the fibrous areas between these two high-quality reconstructions (green arrow A) it can be seen that anisotropic scaling yields sharper ODFs than isotropic scaling. But, for anisotropic scaling, we find an underestimation of over 13 degrees in the crossing angle while isotropic scaling nearly exactly finds the 45 degree crossing. Both methods find the correct number of peaks.

In the middle two rows of Figure 7 we compare MAPL and MAP in the challenge data, whose details we describe in Section 3.2.1. We see for anisotropic scaling (left column, blue arrow B) that using MAPL results in a lower underestimation of the crossing angle compared to MAP (37 • versus 29.4 • ), but has slightly higher standard deviations. For isotropic scaling (right column, red arrow C) we find a similar trend, where the Laplacian finds the peak angle almost exactly (45.9 • ) and the positivity constraint finds a slightly lower angle (42.3 • ), but with a slightly lower standard deviation. The shapes of the ODFs using MAP also seem to be 'pushed together' causing some of the crossing profiles to appear like a single population. In the bottom row, we show an unregularized least squares estimation for reference. In both implementations we used a radial order of 6. Using higher radial order than that resulted in extremely noisy profiles. This is also accompanied by the steep rise in MSE as shown in Figure 4. Still, in both implementations, the standard deviation of the estimated angle and number of peaks is much higher than in MAP and MAPL. Furthermore, in the anisotropic implementation the crossings are mostly pushed together to appear as one population and in the isotropic implementation many spurious peaks are estimated.

In Figure 7, we show the radial moments on human data of the WU-MINN HCP data, described previously in Section 3.2.2. This data has a maximum b-value of 3, 000 s/mm 2 with 3 shells and a total of 288 samples. We selected a small region of interest where projection, association and commissural tracts intersect each other. Comparing reconstructions of this crossing area between our laplacian regularization and the positivity constraint using anisotropic scaling (left column, red arrow D) shows that the positivity constraint produces sharper profiles, but slightly smaller secondary peaks. Comparing the isotropic scaling implementations (right column, green arrow E) similarly shows that the positivity constraint produces slightly sharper profiles. Finally, comparing anisotropic and isotropic implementations (bottom row, blue arrow F) shows an area where the isotropic implementation consistently finds a small crossing, but the anisotropic one only finds a single fiber population. This is possibly a result of the 'pushing together' effect also found in Figure 7.

The phantom data results in Figures 6 and7 show that the laplacian-regularized isotropic implementation performs better for recovering crossing angles. However, on human data, it is debatable which of the compared techniques performs better. Nonetheless, for the signal fitting itself the results in Figure 4 show that the anisotropic implementation with laplacian regularization is more efficient. For this reason, we will omit comparisons of between regularization techniques on our next experiments on human data and always use anisotropic MAPL with the radial order set to 8.

Effect of Subsampling on Scalar Index Estimation

In this section, we investigate how robustly MAPL estimates the scalar indices RTOP, RTAP, RTPP, MSD, and QIV (recall Section 2.3.2) when we reduce the acquisition protocol with respect to the number of samples. The estimation of these indices is of interest as they, under certain assumptions, are related to microstructural features of the brain tissue. For this experiment, we use again the WU-Minn HCP data (description in Section 3.2.2), whose samples are distributed in such a way that a subset of the first M samples of the acquisition scheme is still uniformly distributed on the sphere [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion mri[END_REF]. To give an idea of how important regularization is on real data we show in Figure 9 on the top row the regularized and on the bottom row the unregularized reconstruction of RTOP, RTAP, RTPP, MSD and QIV on a subsampled 60 sample data set. It can be seen that the regularized versions produce very smooth and detailed scalar maps for all indices. On the other hand, without regularization the directional RTAP and RTPP and QIV indices only produce noise, while in the RTOP and MSD images the brain structure can be seen, though the images are very noisy. In Figure 10 we show the MAPL estimations of RTOP, MSD, and QIV, together with a visualization of the optimal Laplacian weighting parameter obtained using GCV. On the left, we show a coronal slice of these indices using the full data. In the middle, we show reconstructions of the same indices in a section near the Corpus Callosum and centrum semiovale (indicated by the red box) using 270, 180, 90 or only 60 samples. On the right, we present histograms of the reconstructed values in the chosen subsection for the different numbers of samples. Both the images and histograms show that the reconstructed values remain very similar even though we only use 60 samples instead of 270. The only thing that does change is the optimal regularization weights -we find that lower regularization weights are required when less data is used. The optimal regularization weight also appears to be tissue specific, finding lower values in anisotropic white matter and CSF and higher in gray matter.

In Figure 11 we show a similar experiment on the scalar indices RTAP and RTPP. However, as we explained in Section 2.3.2, these values are directional and only valid in highly coherent white matter. For this reason, we show the reconstructions of RTAP and RTPP on a coronal and sagittal slice for the whole brain but construct the histograms only using values recovered from a manually drawn region of interest (ROI) in the Corpus Callosum. The segmented area can be seen in green

Effect of Maximum b-value on Apparent Axon Diameter Estimation

The reconstruction of scalar indices depends not only on the number of samples, but also on the number of shells and the b-or q-values (relation given in Section 2.1) at which they were acquired. Given that the diffusion time is not changed, using higher b-values makes the acquisition more sensitive to smaller particle motions. This is important when estimating microstructural features such as the axon diameter. In this section, we study how using different numbers of shells and bvalues influences the estimation of the apparent axon diameter (AAD) estimated using RTAP (recall Section (2.3.2)). We estimate RTAP using MAPL with preset regularization weight to 0.2, which was found to be appropriate in white matter in Section 4.2. We use 6 subjects of the MGH-HCP data (description in Section 3.2.3), which has 4 shells with a maximum b-value of 10, 000 s/mm 2 . As Figure 12 shows, we manually segmented and estimated the AAD in the Genu (Red), Midbody (Green) and Splenium (Blue), as well as the whole Corpus Callosum for each subject. We performed the estimation while we incrementally removed the outer b-shells such that the maximum b-value changes from b max = 10, 000 s/mm 2 (4 shells, 552 samples) to b max = 5, 000 s/mm 2 (3 shells, 296 samples) to b max = 3, 000 s/mm 2 (2 shells, 168 samples). The bar charts in Figure 12 show that, in accordance with literature [START_REF] Aboitiz | Fiber composition of the human corpus callosum[END_REF][START_REF] Liewald | Distribution of axon diameters in cortical white matter: an electron-microscopic study on three human brains and a macaque[END_REF], we find higher AAD in the midbody than in the genu and splenium for all subjects. However, we do estimate significantly larger axon diameters (> 4µm) compared to what is found in histology (< 1.5µm). As we reduce b max we estimate larger diameters, but the relative differences between the different parts of the Corpus Callosum remain similar as shown in the heat maps next to every bar chart.

The decrease in AAD at higher b max is a result of an increase in RTAP, which is the integral of the (interpolated and extrapolated) signal attenuation perpendicular to the bundle orientation. As b max increases, assuming only the gradient strength changes, data points measuring higher qvalues along the perpendicular direction are included. Along this direction the signal decays more slowly, i.e., the diffusion becomes more restricted as the separation between the walls restricting the diffusion becomes smaller. To estimate the amount of diffusion restriction, i.e., the deviation of the measured signal from a DTI approximation (assuming mono-exponential decay), we can use the normalized non-Gaussianity (NG) metric ( Özarslan et al., 2013b, p.22), defined as

N G = 1 -R 3 E DTI (q) 2 dq R 3 E MAPL (q) 2 dq where E MAPL (q) = E DTI (q) + E MAPL\DTI (q). ( 21
)
As is shown on the left side of Eq. ( 21), the NG is estimated through the ratio of the DTI approximation of the data and the entire MAPL approximation. The right side of Eq. ( 21) emphasizes that Figure 11: Effect of subsampling on RTAP and RTPP estimation. As these indices are defined directionally, we construct the histograms only using manually segmented voxels in the Corpus Callosum (in green), where the white matter is highly coherent. Again we find that there is little difference between the reconstructed values between using 270 samples and 60 samples.

the MAPL approximation already conveniently includes the DTI approximation as its first basis function, and the non-Gaussian elements of the signal in all the others, as indicated by E MAPL\DTI . It can be seen that NG ranges from zero when the signal is completely Gaussian -which means E MAPL (q) = E DTI (q) -to one when the signal is completely non-Gaussian.

By assuming that the underlying tissue is axially symmetric, we can treat the signal as separable along the perpendicular and parallel axes [START_REF] Assaf | New modeling and experimental framework to characterize hindered and restricted water diffusion in brain white matter[END_REF]. As a consequence, we can study the perpendicular and parallel NG separately (NG ⊥ and NG ) by taking advantage of the orthogonality of the MAP-MRI basis ( Özarslan et al., 2013b, p.22). In Figure 13 we show the entire (3D) NG, NG ⊥ , and NG estimated in voxels coming from the same corpus callosum ROI we used in subject 2 shown in Figure 12, where we can reasonably assume axial symmetry of the tissue. We study the NG at different b max . We see that as b max increases, the mean of each metric for that b max also increases, but much more in NG and NG ⊥ than in NG . This shows that indeed most of the non-Gaussian (restricted) signal comes from the perpendicular signal component. This result supports the hypothesis as to why the higher-order basis functions in MAPL are useful to estimate these microstructure-related indices that depend on an accurate approximation of the restricted diffusion signal at high b-values. That being said, the AAD values still do not take into account the different tissue compartments, which is required to actually estimate microstructural tissue parameters. In the next section, we therefore propose a novel approach where we combine MAPL and multi-compartment tissue models to estimate the axon diameter and axonal dispersion.

Using MAPL as a Preprocessing For Multi-Compartment Tissue Models

In this section, we propose a novel approach where we combine MAPL with the multi-compartment tissue models described in Section 2.3.3. We consider two models that either estimate axon diameter or axonal dispersion. We estimate the axon diameter using a simplified version of Axcaliber [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF] that uses the Callaghan model [START_REF] Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF] to simulate the intra-axonal diffusion. We will refer to this model simply as "Callaghan". To estimate axonal dispersion we use the In every subject we segmented the whole Corpus Callosum and separately the splenium (red), midbody (green) and genu (blue). Using MAPL we then estimated the AAD in these sections while removing the outer shells of the data set from the full 4-shell data (bmax = 10, 000 s/mm 2 ) to 2-shell data (bmax = 3, 000 s/mm 2 ). The bar chart shows the mean and standard deviation of the recovered AAD in each section, while the heat maps show the recovered AAD in the whole Corpus Callosum. As our results in Section 4.1.1 show, without regularization the AAD cannot be accurately approximated. We estimate larger diameters in the midbody than in the genu and splenium in all subjects. Then, as we remove the outer shells we estimate larger diameters, but as the heat maps and bars show, the relative differences in diameters between the different sections remain quite similar.

NODDI model [START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF]. The idea of our combined approach is to use MAPL as a preprocessing technique to improve the compartment estimation. The idea is to reduce the noise and to include higher b-value data to Callaghan fitting than was originally available by extrapolation the MAPL representation, improving the reliability of the estimation at lower b-values. As MAPL is a functional basis approach, the extrapolation should, to some extent, be guided by the shape of the data itself. We illustrate our approach using the same 6 subjects from the MGH data sets as in the previous section with the same segmentations. In our hybrid approach, we fit the truncated data with either b max = 5, 000 s/mm 2 or b max = 3, 000 s/mm 2 using MAPL and always resample the data on the same signal points as the original b max = 10, 000 s/mm 2 data set. Figure 13: The difference between a DTI and MAPL signal approximation at different bmax, here given as the Non-Gaussianity (NG). The data comes from the same Corpus Callosum ROIs of subject 2 shown in Figure 12, where it is reasonable to assume the tissue is axially symmetric. The left graph shows the NG of the entire 3D signal and the middle and right graphs show the NG along the perpendicular (restricted) and parallel (free) direction (NG ⊥ and NG ). As bmax increases NG and especially NG ⊥ also increase, but NG remains fairly similar, indicating that Non-Gaussian behavior is more dominant primarily in the restricted perpendicular direction.

Estimating Axon Diameter using MAPL + Callaghan

In Figure 14 we show the mean axon diameter distributions of subject 2 for different parts of the Corpus Callosum using either Callaghan or our MAPL + Callaghan hybrid. The black dashed curve represents the 'gold standard' where we fitted Callaghan directly to the b max = 10, 000 s/mm 2 data. It can be seen that directly fitting Callaghan to the lower b max data results in distributions with higher means and standard deviations. Our combined approach, however, results in a mean and standard deviation much closer to that of the gold standard, though the mean shifts slightly towards higher diameters as b max decreases.

Table 1 shows the mean and standard deviation of all 6 subjects for all combinations of the two methods, maximum b-value (b max ) and Corpus Callosum sections. For every subject, the first column shows which method is used, the second column shows which b max was in the original data, and the third through the fifth column show the means and standard deviations of the estimated axon diameter. In this case the b max = 10, 000 s/mm 2 result for the MAPL + Callaghan means that the data was fitted and resampled at exactly the same points. As expected, when b max decreases the standard deviation increases when using Callaghan. The dagger symbols indicate instances where our MAPL + Callaghan approach results in a lower standard deviation of the results compared to the regular Callaghan approach. As can be seen, this is nearly always the case. Though, when just refitting the b max = 10, 000 s/mm 2 data we recover slightly higher axon diameters than just using Callaghan. For the rest, in almost all cases MAPL + Callaghan results in axon diameters that are closer to the is b max = 10, 000 s/mm 2 result than Callaghan.

Estimating Axon Dispersion using MAPL + NODDI

In Figure 15 we show the estimated orientation dispersion index (ODI) using the same setup as Figure 14. Here ODI is dimensionless between 0 (completely parallel) and 1 (completely dispersed). When just using NODDI we find that the ODI increases as b max increases. Similarly as before, MAPL + NODDI has smaller standard deviations than NODDI, but for all b max finds a mean ODI that is close to the b max = 5, 000 s/mm 2 NODDI result. Table 2 shows the mean and standard deviation of the estimated orientation dispersion index (ODI). The layout of the results is the same as in Table 1. When only using NODDI, for all subjects, increasing b max results in an increase in the estimated ODI. On the other hand, The Figure 14: The estimated axon diameters using either only Callaghan (magenta) or our MAPL + Callaghan method (orange) on subject 6 of the MGH-HCP data. The estimation is done for different maximum b-values and sections of the Corpus Callosum. The black dotted line is the reference reconstruction using the full data with bmax = 10, 000 s/mm 2 . Our combined approach typically has lower standard deviations than the regular Callaghan approach.

Corpus Callosum Section

estimated ODI using MAPL + NODDI remains stable over b max , and is in all cases similar to the b max = 5, 000 s/mm 2 ODI of using only NODDI. Furthermore, in contrast to the results on axon diameter, the ODI standard deviations for NODDI remain approximately the same over all b max , but MAPL + NODDI still has lower standard deviations in all cases.

Genu Splenium Midbody

Figure 15: The estimated axonal dispersion using either only NODDI (red) or our MAPL + NODDI method (green) on subject 6 of the MGH-HCP data. The estimation of done for different maximum b-values and sections of the Corpus Callosum. The black dotted line is the reference reconstruction using the full data with bmax = 10, 000 s/mm 2 . Our combined approach always has a lower standard deviation than the regular NODDI approach.

Table 1: Tables of estimated axon diameter values for 6 subjects of the MGH-HCP data set. Estimates are made on the same segmentations of the Corpus Callosum as Figure 11, using either the regular Callaghan method or our MAPL + Callaghan hybrid. The daggers ( † ) symbols after the MAPL + Callaghan standard deviations indicate instances where our hybrid approach reduces the variance of the estimated axon diameters compared to regular Callaghan in the same section and maximum b-value (bmax). It can be seen that this is the case in nearly all instances. Current diffusion MRI research applications still rely heavily on the interpretation of DTI measures such as Fractional Anisotropy (FA) [START_REF] Yeatman | Development of white matter and reading skills[END_REF][START_REF] Jolles | Plasticity of left perisylvian white-matter tracts is associated with individual differences in math learning[END_REF]. To investigate whether microstructure-related indices can be used to better interpret white matter changes, we correlate FA with the estimated axon diameter, dispersion and restricted volume fraction v r estimated using Callaghan. Figure 16 shows the correlation of all three values with FA for subject 1. It can be seen that FA correlates significantly with all values -negatively with axon diameter and dispersion and positively with restricted volume fraction. This means that an increase in FA could mean either (or at the same time) a decrease in axon diameter or dispersion or an increase in restricted volume fraction. There is no way to know which one, if not all, from just from looking at FA. We correlate these values in the Corpus Callosum for all 6 subjects of the b max = 10, 000 s/mm 2 MGH-HCP data, where we used the same segmentations as in Sections 4.3 and 4.4. The results for all 6 subjects are given in Table 3, which all confirm this finding. This last result underlines the importance of looking for more specific biomarkers that describe specific changes in the white matter microstructure.

Discussion

In this section, we discuss the results of our experiments using laplacian-regularized MAP-MRI (MAPL) in Section 4.

MAPL is an Efficient Way to Reconstruct the Ensemble Average Propagator We pose MAPL's fitting process as a regularized least squares problem using generalized crossvalidation (GCV) [START_REF] Craven | Smoothing noisy data with spline functions[END_REF][START_REF] Koay | A signal transformational framework for breaking the noise floor and its applications in mri[END_REF] to find its optimal regularization weight. The GCV procedure was found to be well-suited for our approach, finding consistent values in different white matter tissues (Figure 10). In practice, we can preset the regularization weight to an appropriate value for white matter and omit the GCV without significantly impacting the results, making the estimation as fast as ordinary least squares. On human data of the human connectome project (HCP), we find that MAPL estimates the boundary cases of the EAP contextually smooth and with positive values in almost all cases (Figures 10 and 11). Though, the Laplacian does not explicitly guarantee positive-definiteness of the EAP.

MAPL Improves Signal Reconstruction and Reduces ODF Estimation Bias A key feature of MAP-MRI is its use of anisotropic basis functions. We show in Figure 4 that this feature indeed reduces the fitting error for lower radial orders compared to the isotropic version of MAP-MRI (also known as 3D-SHORE) and outperforms the mSPF basis [START_REF] Caruyer | Diffusion mri signal reconstruction with continuity constraint and optimal regularization[END_REF] (Figure 5). MAPL also has a great speed advantage over positivity constrained MAP-MRI (MAP), especially when presetting the regularization weight (Figure 4(d)). Isotropic MAPL is still over 4 times faster than the anisotropic implementation due to the omission of basis rotations.

Interestingly, Figures 6 and7 show that the anisotropic scaling also causes an underestimation of fiber crossing angles, which we do not see for the isotropic implementation. We previously showed that that this underestimation is directly related to the anisotropy of the basis functions (Fick et al., 2015a), but MAPL in Figure 6 (top-left) shows that this bias lessens as the radial order increases. This indicates that the higher order basis functions correct the shape of the first basis function, representing a DTI tensor, which is aligned with the 'average' direction of the crossing.

When comparing the influence of regularization technique on ODF reconstruction, we find that MAPL has a lower underestimation of the 45 degree crossing than MAP (Figure 7, blue arrow B). The values of the estimated coefficients reveal that the positivity constraint ( Özarslan et al., 2013b) suppresses the higher order coefficients more than MAPL does, reducing the ability of the corresponding basis functions to correct the initial bias. In Figure 7 (right-bottom) we also see this effect in the isotropic implementation, where within the crossing area some voxels are reconstructed to look like a single fiber. The same experiment on ODF reconstruction in human data of the WU-MINN HCP dataset (Figure 8) produces results in which it is harder to choose the optimal estimation approach. We find that both the anisotropic and isotropic implementations of MAP produce slightly sharper angular profiles than those of MAPL. We also find that both anisotropic implementations seem to reduce the crossing angle between crossing fibers, or even merge secondary peaks completely, but this effect appears slightly less for MAPL. Together, these results indicate that when the accurate recovery of ODF peaks is important (e.g. for tractography applications) we must use either isotropic MAP or MAPL, where we find that a radial order of 8 with a radial moment of 2 works better. When only signal fitting is important the 'regular' anisotropic MAPL is more efficient.

Microstructure-Related Scalar Indices Estimated With Less Samples Using MAPL We show on human data of the WU-MINN Human Connectome Project (HCP) that MAPL can estimate the boundary cases of the EAP with the same quality despite subsampling the data from 270 samples to only 60 samples (Figures 10 and11). Note that we keep the same maximum b-value in this experiment, and subsample data points using an incremental scheme [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion mri[END_REF] that guarantees an 'as isotropic as possible' sampling on the sphere. We find that presetting the regularization weight to 0.2 regularizes the reconstruction over the entire brain well, meaning it is possible to omit the GCV optimization without significantly changing the results.

We also find that the optimal regularization weight that GCV finds has tissue-specific contrast, finding lower values in white matter and CSF and higher values in gray matter (Figure 10). The lower values in white matter show that the Laplacian should not impose too much 'smoothness' into the reconstruction, or the reconstruction will become too smooth to represent the anisotropic nature of the data. In this sense, the diffusion profiles in both gray matter and CSF can be seen as isotropic, but the signal decays more quickly in CSF (free water) and more slowly in gray matter (hindered diffusion). This explains the higher regularization weights in gray matter and lower ones in CSF areas.

MAPL Can Distinguish Axon Populations at Lower b-values Using The Apparent Axon Diameter A current focus in diffusion MRI is the estimation of the mean axon diameter [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF][START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion mri[END_REF][START_REF] Huang | The impact of gradient strength on in vivo diffusion mri estimates of axon diameter[END_REF]. The Apparent Axon Diameter (AAD) (see Section 2.3.2) ( Özarslan et al., 2013b) estimates a value that is related to the axon diameter, but is called "apparent" because it does not separate the effects of tissue composition or axonal dispersion on the signal.

In the experiment described in Section 4.3, we used 6 subjects of the MGH-HCP data set to study the effect of using a different number of shells and maximum b-value on the estimation of the AAD. We investigated if we could distinguish AAD populations in the Genu, Midbody, and Splenium of the Corpus Callosum of each subject. For most subjects in Figure 12 we indeed observe the characteristic 'small-big-small' trend that is found in literature [START_REF] Aboitiz | Fiber composition of the human corpus callosum[END_REF][START_REF] Liewald | Distribution of axon diameters in cortical white matter: an electron-microscopic study on three human brains and a macaque[END_REF], which is mostly preserved even though we truncate the data from a maximum b-value of b max = 10, 000 s/mm 2 to b max = 3, 000 s/mm 2 . This shows that the signal extrapolation, on which the estimation of the AAD is based, is well preserved with MAPL at lower b-values. Although, we do see that the mean AAD and its standard deviation increase as b-value decreases, which corresponds to what was found in [START_REF] Huang | The impact of gradient strength on in vivo diffusion mri estimates of axon diameter[END_REF]. The increase in AAD (i.e. decrease in RTAP) at lower b-values is accompanied by a decrease in the non-Gaussianity (NG) of the signal (Figure 13), which shows the difference between a Gaussian approximation of the data and the estimation using MAPL ( Özarslan et al., 2013b). Note that for NG to be physically meaningful, MAP-MRI's scaling parameters must reflect the mean-squared displacements of tissue water spins undergoing Gaussian diffusion. For this reason, only in this part of our study, we estimated the scaling matrix by fitting a DTI model only to the lowest b = 1000 s/mm 2 shell of the data [START_REF] Avram | Clinical feasibility of using mean apparent propagator (map) mri to characterize brain tissue microstructure[END_REF]. As expected, we find that the lower b-value data is better explained by mono-exponential, Gaussian signal decay, but as higher b-value data is included the signal becomes increasingly non-Gaussian (restricted). Indeed, this result indicates that the higher-order basis functions in MAPL improve the estimation of all q-space indices whose estimation depends on the extrapolation of the signal at high q-values. This does not necessarily benefit the MSD as it is estimated as the Laplacian of the signal at the origin q = 0, see Appendix C.3, but certainly proves the added value of using MAPL when estimating RTOP, RTAP, RTPP, and QIV. It should be noted that both DTI measures and q-space indices characterize properties of the entire signal and EAP. Therefore, we do not expect one to have higher specificity to tissue changes than the other. However, it is likely that q-space indices will provide an increased sensitivity to microstructural changes that affect tissue boundaries, i.e., tissue properties that affect the amount of diffusion restriction. In recent work, we found evidence that supports this hypothesis by comparing the sensitivity of both DTI and q-space indices to changes in axonal diameter, dispersion and tissue composition in simulated data [START_REF] Fick | A Sensitivity Analysis of Q-space Indices With Respect to Changes in Axonal Diameter, Dispersion and Tissue Composition[END_REF].

Nonetheless, as histological values of axon diameters are in the range of 0.5-2µm [START_REF] Aboitiz | Fiber composition of the human corpus callosum[END_REF][START_REF] Liewald | Distribution of axon diameters in cortical white matter: an electron-microscopic study on three human brains and a macaque[END_REF], we significantly overestimate the true axon diameter for any b-value. This is the opposite of what we would expect when the short gradient pulse and long diffusion time assumptions are violated [START_REF] Bar-Shir | The effect of the diffusion time and pulse gradient duration ratio on the diffraction pattern and the structural information estimated from q-space diffusion mr: experiments and simulations[END_REF]. Instead, we find an overestimation that is mainly caused by two things: (1) the lack of separation between intra-and extra-axonal compartments [START_REF] Novikov | Revealing mesoscopic structural universality with diffusion[END_REF][START_REF] Burcaw | Mesoscopic structure of neuronal tracts from time-dependent diffusion[END_REF] and (2) the use of low gradient strengths, which results in a resolution limit for small axon sizes [START_REF] Lätt | Accuracy of-space related parameters in mri: Simulations and phantom measurements[END_REF]. Moreover, we find that AAD values typically correlate with FA values. This makes sense because both values globally describe the diffusion signal -the AAD perpendicularly and the FA as a ratio between parallel and perpendicular DTI eigenvalues. Nonetheless, the AAD is a straightforward marker that characterizes the non-Gaussian behavior in the diffusion signal, and its clinical value should be explored further.

Reducing The Variance in Multi-Compartment Tissue Model Estimation Using MAPL To estimate microstructural tissue properties such as the axon diameter or axonal dispersion, the signal contributions from different tissue compartments must be considered [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF][START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion mri[END_REF][START_REF] Zhang | Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF][START_REF] Huang | The impact of gradient strength on in vivo diffusion mri estimates of axon diameter[END_REF]. In this work, we used MAPL's signal extrapolation as a preprocessing to estimate either axonal diameter using a simplified AxCaliber model [START_REF] Assaf | Axcaliber: a method for measuring axon diameter distribution from diffusion mri[END_REF] or axonal dispersion and the NODDI model. It is noteworthy that these parameters cannot be estimated at the same time as we would need at least two diffusion times to disentangle the effects of dispersion and restricted diffusion [START_REF] Nilsson | The role of tissue microstructure and water exchange in biophysical modelling of diffusion in white matter[END_REF]. Furthermore, for axonal diameter estimation, we could also have used the Van Gelderen model [START_REF] Vangelderen | Evaluation of restricted diffusion in cylinders. phosphocreatine in rabbit leg muscle[END_REF] which was used by [START_REF] Huang | The impact of gradient strength on in vivo diffusion mri estimates of axon diameter[END_REF] and [START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion mri[END_REF], but chose Callaghan because it also reconstructs the diffusion propagator just as MAPL.

With respect to axonal diameter, it can be seen from Figure 14 that our combined approach significantly reduces the variance of the mean axon diameter estimation of the Callaghan model. The preprocessing also reduces the increase in estimated mean axon diameter as b-value decreases, with notable exception the b max = 3, 000 s/mm 2 graph of the splenium (right bottom). In this case, the average axon diameter is closer to the gold standard reference for the Callaghan model than for MAPL + Callaghan. Though, this is because Callaghan's axon diameter distribution flattens out between axon diameters of 2 µm and 7 µm, while MAPL + Callaghan maintains a smaller standard deviation.

The fact that the MAPL preprocessing actually results in the smaller estimated axon diameters shows that the restricted nature of the data is still represented to some extent in the signal extrapolation. If this was not the case, and the extrapolated signal quickly decays after the maximum b-value, we would expect the estimated axon diameter to increase as a fast signal decay indicates larger axon diameters.

Table 1 further underlines the results in Figure 14, where the result for all 6 subjects of the MGH-HCP data is given. The dagger sign ' † ' indicates instances where our combined approach has lower variance than the regular Callaghan approach. This is the case for 48 out of 54 instances, where for the remaining 6 the standard deviations are very similar. Interestingly, just fitting and resampling the b max = 10, 000 s/mm 2 data using MAPL also reduces the variance in the result, but consistently slightly increases the estimated axon diameter.

Our results on the estimation of axonal dispersion (ODI) using NODDI and MAPL + NODDI are given in Figure 15 and Table 2. We observe that increasing the maximum b-value increases the average ODI estimation (by about 1-1.5 standard deviations per jump in b-value) while leaving its standard deviation relatively stable. At the same time in restricted volume fraction also increases (result not shown). When we repeat the same experiment using synthetic data with known ground truth ODI and restricted volume fraction we find the same effect with increasing b-value, but that the results at b max = 3, 000 s/mm 2 are closest to the ground truth (result not shown). It is possible that the combination of lower SNR at higher b-values or the sampling scheme cause this effect, but from the results here we cannot point to the exact cause.

When we enhance the estimation using MAPL + NODDI the results change in two ways: (1) The ODI estimation stabilizes over b-value around b max = 5, 000 s/mm 2 result for NODDI and (2) the standard deviation decreases between 25% to 50%. The mean ODI between different sections often comes closer together, but this is in proportion to the decrease in standard deviation. We also find a slightly decreased estimation of the restricted volume fraction compared to NODDI (result not shown).

Overall, these results are meant as an initial exploration into the use of functional basis approaches such as MAPL as a complementary technique to multi-compartment tissue models, rather than a competing technique. Reducing the variability of tissue parameter estimation in this way could possibly increase the statistical power of population tests and the clinical applicability of this approach should be investigated in further studies. It should be noted here that MAPL cannot be used as a preprocessing of the original version of Axcaliber as the MAP-MRI basis is not compatible with data with multiple diffusion times ( Özarslan et al., 2013b). Furthermore, [START_REF] Burcaw | Mesoscopic structure of neuronal tracts from time-dependent diffusion[END_REF] showed that when estimating the axon diameter using multiple diffusion times it is necessary to take into account the time-dependence of the extra-axonal diffusion signal, which biases the estimation even at the long diffusion time limit. Mesoscopic effects, i.e. effects of spatial arrangement, as well as correlations and structural disorder in axon packing should also be considered. To overcome the limitation of MAP-MRI with respect to fitting multiple diffusion times, we have recently proposed an extension of MAP-MRI that is able to simultaneously represent the diffusion signal over both q-space and diffusion time (Fick et al., 2015b).

Compartment Models Better Characterize Tissue Changes Than Fractional Anisotropy In Section 4.5 we investigated whether the microstructure-related indices we estimate add information in terms of describing white matter changes compared to Fractional Anisotropy (FA) [START_REF] Basser | Mr diffusion tensor spectroscopy and imaging[END_REF]. We correlated FA the estimated mean axon diameter, dispersion and restricted volume fraction. Figure 16 proves that FA correlates significantly for all metrics -negatively with axon diameter and dispersion and positively with restricted volume fraction. This means that a change in FA can mean a change in all parameters, but without further information there is no way of knowing which one. The Pearson correlations and p-values in Table 3 underline that this result is consistent for all 6 subjects. Interestingly, axon diameter and axonal dispersion only correlate significantly with each other when b max = 10, 000 s/mm 2 data is used. For lower b max we do not find significant correlations between the estimated parameters of the different models.

This last result underlines the importance of looking for more specific biomarkers that can add more insight to the information that DTI measures already provide.

Conclusion

We proposed a fast, analytic and robust Laplacian regularization of the MAP-MRI basis, which we call MAPL, that allows us to robustly estimate microstructural contrast with fewer samples and at lower b-values. We showed on phantom data with reference golden standard data that MAPL outperforms previously proposed methods with respect to signal fitting and EAP reconstruction in noisy, sparsely sampled data. We also demonstrated the robustness of our method with respect to subsampling on the WU-Minn HCP data set. We showed that using MAPL we find an almost identical reconstruction of microstructure-related scalar indices between the full 270 sample data and subsampled 60 sample data. We also showed that we can use the apparent axon diameter to differentiate between the axon diameter populations of different parts Corpus Callosum. We retain this ability even when we remove the outer shells of the MGH-HCP data from a maximum b-value of 10, 000 s/mm 2 to only 3, 000 s/mm 2 . However, we do find that the value of the apparent axon diameter correlates with FA, meaning that its clinical value should be more carefully studied in the future. We also proposed a new combined approach to estimate the axon diameter and axonal dispersion, where we use MAPL as a preprocessing for multi-compartment models. We showed that this hybrid approach consistently reduces the variability of the estimated metrics. Overall, we showed that using MAPL it is possible to accurately recover important parameters related to the tissue microstructure with much fewer acquisitions. We plan to release our implementation of Laplacian-Regularized MAP-MRI in the open-source dipy framework [START_REF] Garyfallidis | Dipy, a library for the analysis of diffusion mri data[END_REF] and make it publicly available.

The positivity approach for MAP ( Özarslan et al., 2013b) was implemented by casting the coefficient estimation in the cvxopt convex optimization software (http://cvxopt.org/). This approach imposes positivity of the EAP by sampling the EAP P (Rr) in 21 × 21 × 11 grid, resulting in 4851 points, with a maximum sampling distance R max of 20 µm. Here the last dimension is only sampled on its positive axis as the EAP is antipodally symmetric. Lastly, when we use MAP-MRI with isotropic scaling we use the implementation using spherical harmonics (Appendix C), which does not require rotation of the q positions.

Figure 4 :

 4 Figure4: The mean squared error of the reconstructed signal with respect to the ground truth signal. A separation is made for voxels containing one fiber (a), two crossing fibers (b) or isotropic diffusion (c). In (d) the computation time in seconds for every regularization method is given.

  Figure 5: (a) MSE Comparison between our best MAPL result (in green) and mSPF (in red). The mSPF results are given for different radial orders Nmax and angular orders Amax. (b) The number of coefficients in mSPF and our best MAPL (radial order 8). It can be seen that MAPL always has lower errors.

Figure 6 :

 6 Figure 6: Intensity plots of the average estimated crossing angle (top row) and number of peaks (bottom row) of the recovered crossings both anisotropic and isotropic scaling, using either Laplacian Regularization (MAPL) or the Positivity Constraint (MAP). We use different color tables for each row: For the crossing angle, green indicates a correct average recovery of the 45 degree angle, while white indicates the correct average recovery of 2 peaks. The color in each intensity plot is a function of both radial order of the MAP-MRI basis and radial moment of the ODF reconstruction, where higher moments yield sharper ODFs.

Figure 7 :Figure 8 :

 78 Figure 7: The radial moments of the EAP on the phantom data for anisotropic scaling (left column) and isotropic scaling (right column). The top row is the reconstruction on the close-to-noiseless golden standard data, and the bottom three rows show the noisy challenge data. In the second row we use Laplacian regularization (MAPL), in the third the positivity constraint (MAP) and in the last unregularized least squares. In every image we report in the top-right the radial order, radial moment, average and standard deviation of the estimated angle and number of peaks. The rectangles indicate a crossing area in the phantom. The arrows are meant to facilitate the explanation of the results.

Figure 9 :

 9 Figure 9: Regularized (top row) and unregularized (bottom row) reconstruction of RTOP, RTAP, RTPP, MSD and QIV with only 60 samples on a coronal slice. It can be seen that the regularized reconstructions produce smooth scalar maps, while the unregularized ones are very noisy and contain negative values.

Figure 10 :

 10 Figure 10: Effect of subsampling on RTOP, MSD, QIV and optimal regularization weight as found by GCV. Left: a coronal visualization of these indices. Middle: Visualizations of a subset of the coronal slice (in red) as we progressively remove samples from the whole 270 samples to only 60 samples. Right: Histogram representations of these quantities for the different amounts of samples. It can be seen that the values of RTOP, MSD and QIV change very little even when only 60 samples are used. Only the GCV regularization weight changes: as fewer samples are used we find lower optimal regularization weights.

Figure 12 :

 12 Figure 12: Estimation of apparent axon diameter (AAD) in the Corpus Callosum for 6 subjects of the MGH-HCP database.In every subject we segmented the whole Corpus Callosum and separately the splenium (red), midbody (green) and genu (blue). Using MAPL we then estimated the AAD in these sections while removing the outer shells of the data set from the full 4-shell data (bmax = 10, 000 s/mm 2 ) to 2-shell data (bmax = 3, 000 s/mm 2 ). The bar chart shows the mean and standard deviation of the recovered AAD in each section, while the heat maps show the recovered AAD in the whole Corpus Callosum. As our results in Section 4.1.1 show, without regularization the AAD cannot be accurately approximated. We estimate larger diameters in the midbody than in the genu and splenium in all subjects. Then, as we remove the outer shells we estimate larger diameters, but as the heat maps and bars show, the relative differences in diameters between the different sections remain quite similar.

Figure 16 :

 16 Figure16: Correlations between FA and axon diameter (left), axon dispersion (middle) and restricted volume fraction (right). It can be seen that they are all significantly correlated, meaning that a change in FA can mean either a change in diameter, dispersion or volume fraction, but it is impossible to know one from the other.

Table 2 :

 2 Tables of estimated axonal dispersion values for 6 subjects of the MGH-HCP data set. Estimates are made on the same segmentations of the Corpus Callosum as Figure11, using either NODDI or our MAPL + NODDI hybrid.

	Subject 1			Section		Subject 4			Section
			Genu	Midbody	Splenium			Genu	Midbody	Splenium
	Method	bmax(s/mm 2 )	diameter (µm)	diameter (µm)	diameter (µm)	Method	bmax(s/mm 2 )	diameter (µm)	diameter (µm)	diameter (µm)
		3,000	5.95±0.92	6.83±0.62	5.43±0.69		3,000	5.83±0.99	6.51±0.77	5.07±0.92
	Callaghan	5,000	5.21±0.60	6.08±0.40	4.93±0.46	Callaghan	5,000	4.94±0.57	5.55±0.45	4.79±0.50
		10,000	4.22±0.27	4.76±0.30	3.99±0.36		10,000	3.91±0.35	4.39±0.33	4.06±0.30
		3,000	5.55±0.47 †	6.04±0.42 †	5.22±0.54 †		3,000	5.61±0.53 †	5.81±0.46 †	4.97±0.45 †
	MAPL + Callaghan	5,000	4.48±0.45 †	5.14±0.41	4.12±0.57	MAPL + Callaghan	5,000	4.43±0.42 †	4.78±0.49	4.16±0.52
		10,000	4.23±0.20 †	4.61±0.29 †	3.90±0.14 †		10,000	4.31±0.28 †	4.71±0.31 †	4.24±0.28 †
	Subject 2			Section		Subject 5			Section
			Genu	Midbody	Splenium			Genu	Midbody	Splenium
	Method	bmax(s/mm 2 )	diameter (µm)	diameter (µm)	diameter (µm)	Method	bmax(s/mm 2 )	diameter (µm)	diameter (µm)	diameter (µm)
		3,000	5.13±0.99	6.58±0.81	4.28±0.93		3,000	5.58±1.09	5.51±1.35	5.47±0.87
	Callaghan	5,000	4.96±0.57	5.81±0.55	4.22±0.60	Callaghan	5,000	5.06±0.54	5.16±0.95	4.75±0.59
		10,000	4.09±0.31	4.57±0.40	3.76±0.34		10,000	4.51±0.24	4.46±0.62	3.90±0.39
		3,000	4.95±0.42 †	5.64±0.53 †	4.71±0.34 †		3,000	5.63±0.40 †	5.67±0.47 †	5.06±0.41 †
	MAPL + Callaghan	5,000	4.48±0.45 †	5.15±0.56	3.93±0.44 †	MAPL + Callaghan	5,000	4.58±0.37 †	4.76±0.65 †	4.01±0.50 †
		10,000	4.35±0.30 †	4.90±0.35 †	3.97±0.32 †		10,000	4.70±0.25	4.62±0.60 †	4.22±0.32 †
	Subject 3			Section		Subject 6			Section
			Genu	Midbody	Splenium			Genu	Midbody	Splenium
	Method	bmax(s/mm 2 )	diameter (µm)	diameter (µm)	diameter (µm)	Method	bmax(s/mm 2 )	diameter (µm)	diameter (µm)	diameter (µm)
		3,000	4.35±1.47	5.24±1.16	4.68±1.12		3,000	4.37±1.53	5.71±1.05	3.59±1.46
	Callaghan	5,000	4.50±0.92	5.07±0.74	4.53±0.74	Callaghan	5,000	4.24±1.00	5.22±0.65	3.73±0.86
		10,000	4.10±0.45	4.35±0.53	3.96±0.44		10,000	3.99±0.43	4.25±0.49	3.39±0.47
		3,000	5.11±0.40 †	5.32±0.42 †	4.98±0.51 †		3,000	5.25±0.42 †	5.51±0.45 †	4.84±0.44 †
	MAPL + Callaghan	5,000	4.50±0.40 †	4.63±0.55 †	4.11±0.58 †	MAPL + Callaghan	5,000	4.21±0.48 †	4.59±0.51 †	3.84±0.44 †
		10,000	4.23±0.39 †	4.52±0.48 †	4.09±0.43 †		10,000	4.15±0.40 †	4.46±0.48 †	3.54±0.45 †

Table 3 :

 3 Tables showing the correlations between Fractional Anisotropy (FA) and the estimated axon diameter, dispersion and restricted volume fraction on the full data of all 6 subjects of the MGH-HCP data. In this case the voxels for the segmentations of the Corpus Callosum were joined together. It can be seen that FA always correlates negatively with axon diameter and dispersion and positivity with restricted volume fraction.

	Subject 1	Pearson r p-value	N	Subject 4	Pearson r p-value	N
	Diameter	-0.60	3.1e-42 417	Diameter	-0.43	4.4e-19 402
	Dispersion Index	-0.59	4.9e-40 417	Dispersion Index	-0.36	5.2e-14 402
	Restricted Volume Fraction	0.76	8.2e-79 417	Restricted Volume Fraction	0.79	1.4e-88 402
	Subject 2	Pearson r p-value	N	Subject 5	Pearson r p-value	N
	Diameter	-0.69	1.7e-61 422	Diameter	-0.38	9.5e-14 354
	Dispersion Index	-0.86	3.5e-59 422	Dispersion Index	-0.59	1.4e-34 354
	Restricted Volume Fraction	0.89	2.7e-94 422	Restricted Volume Fraction	0.81	5.1e-85 354
	Subject 3	Pearson r p-value	N	Subject 6	Pearson r p-value	N
	Diameter	-0.52	1.9e-30 427	Diameter	-0.69	9.7e-69 476
	Dispersion Index	-0.47	5.3e-25 427	Dispersion Index	-0.64	5.2e-57 476
	Restricted Volume Fraction	0.76	1.6e-81 427	Restricted Volume Fraction	0.85	1e-132 476

Acknowledgments

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. This work was partly supported by the French ANR "MOSIFAH" under ANR-13-MONU-0009-01. We would also like to thank Marco Pizzolato and Mauro Zucchelli for their constructive input.

Appendix A. Acronym Glossary

In this section, we provide in detail the algorithmic steps involved in fitting the basis coefficients in MAPL and the multi-compartment models.

Appendix B.1. MAPL

Our implementation of Laplacian-regularized MAP-MRI consists of the following steps:

1. For every voxel, we first estimate the scale factors {u x , u y , u z }, which we find by fitting a diffusion tensor to the data [START_REF] Basser | Mr diffusion tensor spectroscopy and imaging[END_REF]. In our experiments, we did not impose any positivity constraints on the DTI estimation. 2. Using {u x , u y , u z } and the rotated q positions of the data points, we compute the design matrix Q and Laplacian regularization matrix U as outlined in Section 2.2.2. 3. We find the optimal regularization weight λ ∆ (recall Eq. ( 14)) using the generalized crossvalidation (GCV) approach outlined in section 2.2.3. The minimum of the GCV cost function is found from a preset number of regularization weights, which we adjusted for different datasets. 4. Using Q, U and λ ∆ we estimate the unknown coefficients using the regularized least-squares equation in Eq. ( 14). 5. The estimated coefficients were fed to Eq. ( 16) to estimate the radial moments of the EAP and to Eqs. ( 17), (D.1) and (D.2) to estimate the boundary cases of the EAP, mean squared displacement and q-space inverse variance. The maxima of the radial moments were found by projecting the angular profiles on a tessellation of 2172 spherical points and finding the maximum amplitudes with a minimum crossing angle of 10 degrees and a minimum peak amplitude of 30% of the largest peak amplitude.

Appendix B.2. Multi-Compartment Tissue Models

For the NODDI model we used the open-source NODDI toolbox for MATLAB (http://mig. cs.ucl.ac.uk/). In our simplified Axcaliber implementation we use free water diffusivity D = 3µm 2 /ms for the isotropic compartment and a hindered and restricted diffusivity of D = 1.7µm 2 /ms [START_REF] Huang | The impact of gradient strength on in vivo diffusion mri estimates of axon diameter[END_REF]. Our stepwise implementation is as follows:

1. For every voxel in the region of interest we fit a tensor and use its principal eigenvector as the axon direction, and rotate the q-space vectors into the anatomical frame of reference.

we use a brute force algorithm to find an initial set of parameters between

We then finetune the found parameters using a quasi-Newton algorithm to solve the nonlinear optimization problem [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] with bounds

The finetuned parameters are then reported.

Appendix C. Isotropic MAPL

The isotropic implementation of MAP-MRI ( Özarslan et al., 2013b, Appendix A), which is equivalent to 3D-SHORE [START_REF] Merlet | Continuous diffusion signal, eap and odf estimation via compressive sensing in diffusion mri[END_REF], describes the signal and EAP as

where Υ = IFT(Ξ) and c i are the basis coefficients. These basis functions are given as

where j = (n + 2 -l)/2 is related to the radial order n and angular order l of the basis where j ≥ 1, l ≥ 0. The scale factor u 0 is related to the diffusivity of the measured data as u 0 = √ 2Dτ with D the mean diffusivity and τ the effective diffusion time. The real spherical harmonic basis Y m l [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] has angular order l and phase factor m such that -l ≤ m ≤ l, and L l+1/2 j-1 is the generalized Laguerre polynomial.

Appendix C.1. Laplacian Regularization for Isotropic MAPL

Exactly the same as for the anisotropic MAP-MRI basis in Eq. ( 7), the values of the Laplacian regularization matrix are given as

(C.4)

The equation for the elements of U can again be solved using the general differential equation whose solutions form the functional basis functions Ξ jlm of the isotropic MAP-MRI basis ( Özarslan et al., 2013b, Eq. ( 56))

with Λ jlm = 2l + 4j -1. Inverting this equation we show that

Inserting Eq. (C.6) into Eq. (C.4), using the fact that Y m l is an orthonormal basis with respect to the dot product on S 2 and L α n (x) is orthonormal with respect to the weighting function x α e -x on [0, ∞), we find the general equation for U as

with δ the Dirac delta function. Note that regularization matrix U is symmetric, mostly sparse and its elements depend only on the ordering of the basis functions and their radial and angular indices j, l and m. Also note that U(u 0 ) = u 0 U(1), meaning that the regularization matrix can be entirely precomputed and later multiplied by the data-dependent scale factor u 0 .

Appendix C.2. Radial Moment Computation

The fitted isotropic MAP-MRI coefficients analytically describe the EAP in Eq. (C.1). They therefore also describe the ODF at an arbitrary radial moment in Eq. ( 16). Inserting the right side of Eq. (C.1) into Eq. ( 16) and evaluating the integral gives

with v an orientation on the unit sphere and

the ODF basis function with

The scaling factor u 0 is just a multiplication or division of Ω jlm s , depending on the radial moment. Given a set of orientations v the ODF can be computed as ODF s (u 0 , v) = u s 0 cI with design matrix

(1, v k ), allowing for the precomputation of I for the whole data set after the radial moment has been chosen.

Appendix C.3. Scalar Indices for q-space imaging

In this work we computed the boundary cases of the EAP -the RTOP, RTAP and RTPP ( Özarslan et al., 2013b) -and also the Mean Squared Displacement (MSD) and q-space Inverse Variance (QIV) [START_REF] Wu | Computation of diffusion function measures in-space using magnetic resonance hybrid diffusion imaging[END_REF]. Here we provide these quantities in terms of isotropic MAP-MRI coefficients, which were not previously given:

with δ (l,0) is nonzero when l = 0, c {j,l,m} are the isotropic MAP-MRI coefficients and

with u f iber the estimated direction of the underlying fiber. The MSD has a Fourier relation to the diffusion signal and is given as [START_REF] Cheng | Estimation and Processing of Ensemble Average Propagator and Its Features in Diffusion MRI[END_REF]. Evaluating this equation finally gives

The QIV is similarly given as (Hosseinbor et al., 2013, Eqs. (21, 22)), which gives

Appendix D. MSD and QIV for Anisotropic MAP-MRI

We also newly provide the equations for MSD and QIV in terms of anisotropic MAP-MRI coefficients (recall Section 2.2.1):

For both these quantities only the contributions of terms where n 1 , n 2 and n 3 are even are nonzero.

Appendix E. Callaghan Model

In Section 4.4 we use the Callaghan model [START_REF] Callaghan | Pulsed-gradient spin-echo nmr for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF] to estimate the axon diameter in the Corpus Callosum. The Callaghan model describes the intra-axonal diffusion signal as E r (q, τ ) = where J n are the derivatives of the n th -order Bessel function, β nk are the arguments that result in zero-crossings and the cylinders are of radius a. This model describes the signal in parallel cylinders for any diffusion time, but keeps the short gradient pulse assumption. As in [START_REF] Assaf | New modeling and experimental framework to characterize hindered and restricted water diffusion in brain white matter[END_REF], we model the three-dimensional intra-axonal signal as a product of the restricted perpendicular and free parallel diffusion, i.e. E ic (q, D ) = E restricted (q ⊥ , D )E free (q ).