In dMRI, functional basis approaches have been used to efficiently represent the diffusion signal with little assumptions on its shape. Following this methodology, we represent the measured multi-spherical signal E(q, r) in terms of a continuous functional basis E(q, r; c), where the signal is now represented in terms of coefficients c 2 R^{Nc} .An effective representation E (q, r; c) should be able to

1. closely approximate the measured multi-spherical dMRI signal,

2. smoothly interpolate between and outside the measured fq, r} points,

3. have a sparse representation in c,

4. be able to reconstruct the EAP from the fitted signal.

Requirements 1-3 are described in Eq. (2), while the fourth will follow by choosing a functional basis that is also a Fourier basis.

Note that the integrals over three-dimensional q have limits [—1,1 ] and those over x have limits [0,1]. As stated in Sect. 2.1, the boundary constraints are important to respect the Fourier relationship between the fitted signal attenuation and the EAP.

Functional Basis Signal Representation

We represent the multi-spherical signal using an orthogonal basis that allows for the implementation of all our previously stated requirements. As we assume an infinitely short gradient pulse (1 ! 0), we follow Callaghan et al.’s description of time-dependent diffusion in pores and assume separability in the dependence of the dMRI signal to q and x [12]. Following this hypothesis, we can independently choose any representation for these two spaces. We represent the combined space E(q, x; c) using the cross-product between the spatial basis Ф(я) and temporal basis T(x) as

where N_{q} and N_{x} are the maximum expansion order of each basis and c_{ik} weights the contribution of the ikth basis function to E(q, x; c).

A plethora of functional bases to represent q have been proposed, e.g. [8, 9, 13, 15]. Of these bases, we use the Mean Apparent Propagator (MAP) basis [13] as it neatly fulfills all four previously stated requirements; (1) being an orthogonal basis, it can accurately represent any signal over q using few coefficients; (2) it allows to impose smoothness using analytic Laplacian regularization [14]; (3) the isotropic MAP implementation was successfully used to obtain sparse signal representation [8] and (4) MAP is a Fourier basis. It is worth noting that this basis is different than Fick et al.’s basis, who used the isotropic implementation (3D- SHORE) to represent q [4].

MAP’s signal basis is a product of three orthogonal Simple Harmonic Oscillator- based Reconstruction and Estimation (SHORE) functions ф„(ы) [16]:

with its Fourier transform, the EAP basis as

where H is a physicist’s Hermite polynomial of order n and u is a data-dependent scale factor. As in MAP [13], before fitting, the data is rotated such that the DTI eigenvectors are aligned with the coordinate axis and we can use the data-dependent scaling matrix A = Diag(uZ, Uy, u^{2}z) to scale the MAP basis functions according to the anisotropy of the data. The zeroth order is a purely Gaussian function while higher orders use the Hermite to correct this approximation to the true shape of the data. For a given radial order N_{rad} the number of coefficients is N_{q} = (N_{rad} + Z)(Nr_{a}d C 4)(ZNrad C 3)/24.

Our functional basis to describe r was introduced in Fick et al. [4]. As a limiting case the diffusion signal dependence on r is exponential for pure Gaussian diffusion and constant for diffusion in restricted geometries. To represent r we, therefore, choose a product of the negative exponential and a Laguerre polynomial L, which together form an orthogonal basis over r

with basis order p and temporal scaling factor u_{t}. The zeroth order is a pure exponential function and higher orders use the Laguerre polynomials to correct this approximation to the true shape of the signal.

For the rest of this work we will linearize the ordering of our multi-spherical basis such that we use one basis index i with notation

where the total number of fitted coefficients is N_{c} = (N_{r} C 1)(N_{q} CZ)(N_{q} C4)(2N_{q} C 3)/24. Using this notation, the fitted signal E(q, r; c) in the Data Fidelity term in Eq. (2), with measured signal y e R^{Ny} and N_{y} the number of samples, can be represented as y = Фc with Ф e R^{NyxNc} with values Фу = Ej(q_{;}, r, A, u_{t}).

The multi-spherical EAP can be reconstructed using MAP’s Fourier properties [13]. The Fourier transform only concerns the q-space, so the EAP is found simply by switching Ф^, A) in Eq. (7) by its Fourier transform in Eq. (5).

Analytic Laplacian Regularization

We impose smoothness in the multi-spherical signal reconstruction by using the squared norm of the Laplacian of the reconstructed signal as a regularizer. We define the Smoothness term in Eq. (2) as Laplacian functional U(c) as

where, due to our choice of basis, the Laplacian of the reconstructed signal can be estimated as V^{2}?_{c}(q, r) = _{;} c,V^{2}S;(q, r, A, u_{t}). Equation (8) can be further

rewritten in quadratic form as U(c) = c^{T} Uc with

where the subscript ik indicates the ikth position in the regularization matrix. We use the orthogonality of the basis functions (standard inner product on [0,1]) to compute the values of the regularization matrix to a closed form depending only on the basis orders and scale factors. We provide U in Appendix.

Coefficient Estimation from Multi-Spherical Data

We represent the multi-spherical signal ?(q, r) in terms of a sparse coefficient vector c as y = Фc + e with Ф the observation matrix, y the signal values and e the acquisition noise. We frame the numerical implementation of our approach in the same way as we did continuously in Eq. (2):

where we described the Data Fidelity and Smoothness term in Sect. 2.2, and the Sparsity term and constraints are imposed by framing our problem as a convex optimization using the open-source package CVXPY [17]. We find optimal values for regularization weights a and using cross-validation and implemented the surrounding code infrastructure inside the DiPy framework [18].

Estimation of т -Dependent q-Space Indices

Once coefficients c are known, our basis allows us to freely explore, for any diffusion time, all previously proposed scalar metrics for the three-dimensional EAP [13, 14], also known as q-space indices. We can do this because our basis reduces to the MAP basis when the temporal basis is evaluated for a particular diffusion time. In this work we illustrate this using the r-dependent Return-To- Origin Probability (RTOP) and Mean Squared Displacement (MSD):