Spherical deconvolution with tissue-specific response functions and multi-shell diffusion MRI to estimate multiple fiber orientation distributions (mFODs)

In diffusion MRI, spherical deconvolution approaches can resolve local white matter (WM) fiber orientation distributions (FOD) to produce state-of-the-art fiber tractography reconstructions. In contrast, spherical deconvolution typically produces suboptimal results in grey matter (GM). The advent of multi-shell diffusion MRI data has seen a growing interest toward GM, but its modelling remains confined to that of an isotropic diffusing process. However, evidence from both histology and high-resolution diffusion MRI studies suggests a marked anisotropic character of the diffusion process in GM, which might be exploited by appropriate modelling to improve the description of the cortical organization. Indeed, improving the performance of spherical deconvolution in GM might allow to track WM bundles to their start/endpoint in the cortex, a crucial step to improve the accuracy of fiber tractography reconstructions. In this study, we investigated whether performing spherical deconvolution with tissue specific models of both WM and GM improves the FOD characterization in GM while also retaining state-of-the-art performances in WM. To this end, we devised a framework able to simultaneously accommodate multiple tissue response functions, to estimate multiple, tissue-specific, fiber orientation distributions (mFOD.) As proof of principle, we used the diffusion kurtosis imaging model to represent the WM signal, and the neurite orientation dispersion and density imaging (NODDI) model to represent the GM signal. The feasibility of mFOD is shown with numerical simulations, then the performances of the method were compared to gold-standard multi-shell constrained spherical deconvolution (MSCSD) with data of a subject from the Human Connectome Project (HCP). Results of simulations show mFOD can accurately estimate mixture of two FODs given sufficient SNR (≥ 50 at b = 0s/mm2) is achieved. With HCP data, mFOD provided improved FOD estimation in GM as compared to MSCSD, whereas the performances of the two methods were comparable in WM. When performing fiber tractography, tracts with mFOD entered the cortex with more spatial contiguity and for a longer distance as compared to MSCSD. In conclusion, mFOD allows to perform spherical deconvolution of multiple anisotropic response functions, improving the performances of spherical deconvolution in grey matter.


Introduction
Diffusion Magnetic Resonance Imaging (dMRI) is a non-invasive technique sensitive to the microscopic motion (diffusion process) of water molecules 1 . In biologic tissues, the diffusion process is influenced by the presence of biologic membranes and macromolecules 2 , which can hinder or restrict the molecular random walk in both isotropic and anisotropic manners. In the '90s, Basser et al. 3 characterized the anisotropy of hindered diffusion processes, proposing the diffusion tensor imaging (DTI) framework. Later, it was shown that when the diffusion process exhibits sufficient anisotropy, its principal direction can be computed 3 and eventually tracked 4 . This approach can be applied, for instance, to reconstruct the configuration of brain white matter (WM) with fiber tractography 5,6 .
DTI-based fiber tractography is acquisition and computation efficient but suffers of a number of limitations. Firstly, the technique fails when multiple diffusion processes -i.e. multiple pathways -coexist in the same spatial location 7 , which is the case in over 90% of the human WM 8 and is not likely to be solved by technological improvements as higher imaging resolutions 9 . Secondly, DTI-based fiber tractography only applies to data acquired within the Gaussian diffusion regime, and thus cannot take advantage of the higher angular resolution provided by strong diffusion weightings 10 . For this reason, DTI is not likely not the method of choice in a heterogeneous structure as grey matter (GM), which is characterized by low fractional anisotropy and hence high uncertainty of the principal diffusion direction.
Spherical deconvolution approaches 11,12 effectively solve the first two limitations, and are the current gold standard to reconstruct fiber orientation distribution (FOD) and to perform fiber tractography in WM. However, their application to the study of GM remains limited, as FODs derived in the cortex are noisy and exhibit spurious peaks, a challenge when performing fiber tractography. Importantly, limitations in fiber tractography propagate to other analysis techniques making use of fiber reconstructions as, among others, structural connectivity approaches 13 . Indeed, fiber tracking is often seeded (or terminated) at the grey/white matter interface, but the inability to properly follow WM bundles in the cortex generates uncertainty onto the endpoints of the bundles. This has been partially addressed with geometric heuristics 14 , but these oversimplify the way axons traverse the cortex 15 .
The advent of specialized techniques for diffusion imaging, such as ultra-strong gradients 16 and simultaneous multi-slice 17 , are rapidly increasing the imaging quality and resolution achievable in dMRI, allowing to appreciate new details of the human cortex, as reviewed by Assaf 18 . In ex-vivo studies 15,19,20 , it has been consistently shown that the diffusion process observed in the cortex has a large anisotropic component perpendicular to the cortical folding, but also that tangential directions can be observed, for instance, in the central sulcus 21 . Such high-quality data acquisitions may soon be available for in-vivo applications, offering unprecedented research possibilities 22 in fiber tractography of GM, but a dedicated framework has not yet been investigated.
Recently, multi-tissue spherical deconvolution approaches 23 have been suggested to further increase the accuracy of FOD reconstructions in WM by accounting for partial volume effects, e.g. in spatial locations where WM overlaps with GM or cerebrospinal fluid (CSF). Although GM and CSF are merely considered nuisance factors 24 , such framework highlights the feasibility of explicitly accounting for multiple tissues in spherical deconvolution approaches. By separating multiple partial volume effects, multi-tissue approaches have shown improved FOD reconstructions in adults 23 , neonates 25 and in patients with Parkinson's disease 26 . In current state-of-the-art multi-shell spherical deconvolution, only WM is assumed to be isotropic, while GM is regarded as characterized by isotropic diffusion. To estimate the WM FOD, a WM fiber model -or response function -is estimated from voxels with single fiber population, which are typically located in large bundles 27 such as, for example, the corpus callosum or the corticospinal tract. The response function is then assumed to be valid throughout the brain. While the validity of this representation has been questioned even in WM itself 28 , it is highly unlikely that it can adequately represent the anisotropic signature of axons and dendrites in the cortex.
In this study, we investigate whether the spherical deconvolution problem can be revisited to take into account the different signal characteristics of WM and GM, by explicitly modelling both signal mixtures simultaneously. This novel formulation allows to reconstruct multiple FODs (mFODs) and is ultimately aimed at improving the FOD estimation and fiber tractography of GM while retaining state-of-the-art performances in WM. In mFOD, we disregard the existence of a single response function for the brain, and investigate whether it is feasible to i) simultaneously estimate multiple anisotropic FODs corresponding to tissue-specific response functions, and whether this improves ii) the fiber orientation characterization in GM, and iii) fiber tractography in GM as compared to existing spherical deconvolution methods.

Methods
The following paragraphs present the mFOD framework and explain the qualitative and quantitative analysis that we performed using both simulations and in-vivo data.

Theory: multi-shell Richardson Lucy (mRL)
In a previous work 29 , we introduced a multi-shell spherical deconvolution framework based on the Richardson-Lucy (mRL) framework to account for multiple tissue types, such as WM, GM and cerebrospinal-fluid (CSF), as shown in Eq. 1. [1] In the equation above, S is the signal collected at n diffusion weightings, S0 is the non-weighted signal, H represents the signal model associated to WM, GM and CSF, respectively, f is the associated signal fraction, and FOD is the fiber orientation distribution associated to WM. The columns of the deconvolution matrix shown in Eq. 1 represent the possible solution of the spherical deconvolution on the unit sphere, resulting in a total number of k+2 columns, as GM and CSF are both considered isotropic.

Theory: mFOD
In the mFOD framework, we redesigned the deconvolution matrix to allow for an arbitrary number of independent anisotropic deconvolution models. Considering to solve the spherical deconvolution problem once more for WM, GM and CSF, but now accounting for anisotropic diffusion in both WM and GM, we can write: [2] Eq. [2] can be solved using the iterative least-squares minimizer previously introduced for mDRL but replacing the Richardson-Lucy deconvolution scheme with a regularized non-negative least-squares solver 30 , producing 2 independent FODs.

Signal models
To underline the model independence of the mFOD framework, the H-matrices of WM, GM and CSF, shown in Eq. 2, were generated with three different models. HCSF was generated using the ADC model and the typical diffusion coefficient of free water at 37 degrees (3.0x10 -3 mm 2 /s). The columns of HWM and HGM were generated using instances of the DKI model 31 with isotropic kurtosis 32 , and of the NODDI model 33 , respectively, both in simulations and with in-vivo data. The parameters of the models were determined as explained in the following paragraphs. The DKI model was chosen as it can represent the anisotropy of the WM signal while also taking into account non-Gaussian diffusion effects. The NODDI model has been shown to characterize the GM microstructure better than tensor models 34 , and allows to take into account the larger presence of cell bodies in GM as compared to WM.

Simulations
The feasibility of simultaneously computing multiple anisotropic FOD corresponding to different models was evaluated numerically. HWM was generated with the DKI model by considering eigenvalues (1.7, 0.3, 0. 3) x 10 -3 mm 2 /s, and isotropic kurtosis equal to 0.6. HGM was generated with the NODDI model by simulating values of intra to extra-cellular ratio 0.4, parallel diffusivity 1.7x10 -3 mm 2 /s, the concentration parameter of the Watson distribution was set to 1.0, and no partial volume of free-water was considered as mFOD explicitly accounts for it. The parameters used for generating HWM and HGM matched the average values observed with in-vivo data.
Two set of simulations were generated. Simulation I aimed to evaluate the robustness to noise of mFOD, whereas simulation II characterized the ability of mFOD at separating partial volume effects for a given SNR level.

Simulation I
Three configurations were generated by considering partial volume effects between i) WM and CSF, ii) GM and CSF, iii) WM and GM. The simulated partial volume ratio was 0.8/0.2, with an angle of approximately 90 degrees between the main direction of WM and the main direction of GM. The effect of noise was simulated by adding 1000 Rician realizations in correspondence of SNR levels 10, 20, 30, 40, 50, 60, 70, 150. For each SNR level, we evaluated the error between the estimated and the simulated signal fractions, as well as the error between the simulated and the estimated main diffusion direction.

Simulation II
Three partial volume configurations were generated in analogy with simulation I but increasing the ratio between the two simulated tissues from 0 to 1 with step 0.1. Rician noise was added by considering 1000 realizations at SNR 50. The error between the simulated and the estimated signal fractions was determined.

In-vivo data processing
A subject from the Human Connectome Project (HCP) was included in this study. The dataset included a T1-weighted image acquired at a resolution of 0.7mm isotropic, and a dMRI dataset acquired at 1.25mm isotropic. The diffusion dataset consisted of 18 b = 0 s/mm 2 volumes in addition to 270 volumes acquired by sampling diffusion weightings b = 1000, 2000, 3000 s/mm 2 along 90 directions each. The T1-weighted data was rigidly registered to the pre-processed dMRI data, then the tissue type segmentations of WM, GM and CSF were derived with FSL FAST.

Determining the H-matrix
The dMRI data was fit with both the MSCSD and the mFOD approaches using custom implementations in MATLAB R2018b (The Mathworks Inc.) and functions from the ExploreDTI toolbox 35 . MSCSD was initialized as previously suggested by using the tissue maps derived from the T1-weighted image. For the generation of the H-matrices used for mFOD, a whole brain DKI and NODDI fit were performed. For the WM H-matrix, the eigenvalues of the diffusion tensor and the isotropic kurtosis were averaged within a white matter mask defined by fractional anisotropy values above 0.7. This resulted in eigenvalues [1.3, 0.5, 0.5]x10 -3 mm 2 /s and isotropic kurtosis value 0.64. For the GM H-matrix, the NODDI metrics were averaged within the GM segmentation obtained with the T1-weighted image, which resulted in intra to extra-cellular ratio 0.4, parallel diffusivity 1.7x10 -3 mm 2 /s and Watson's concentration parameter equal to 1.

FODs and signal fractions
The MSCSD fit produced signal fractional maps associated to WM, GM and CSF, and one FOD describing WM orientations. The fit of mFOD resulted in the same number of fractional maps but two FODs, mFOD1 describing WM geometry, and mFOD2 describing GM geometry, respectively.

Comparison of signal fractions
The signal fraction maps computed with MSCSD and mFOD were qualitatively compared to those obtained from the T1-weighted volume. Further, the FODs computed with MSCSD and mFOD were compared in selected cortical areas, to qualitatively determine whether mFOD2 provided additional information in GM as compared to mFOD1 and MSCSD.

Quantitative analysis of FODs
The FODs computed with mFOD can be either analyzed separately, or merged to combine accurate geometrical information in both WM and GM. In name of clarity and purely as proof of concept, we combined the two FODs by considering their normalized sum weighted by the respective signal fractions (mFOD-WS).
One of the difficulties when validating in-vivo dMRI techniques is the lack of ground truth. MSCSD can be regarded as the gold standard in WM, but not in GM. In a first analysis, we aimed to establish to which extent the principal diffusion directions of mFOD-WS agree with those of MSCSD in both WM and GM. Up to three peaks were voxel-wise identified for both FODs, then sorted per descending amplitude values. Subsequently, the angle between MSCSD and mFOD-WS was determined as the minimum dyadic angle between any pair of peaks.
If discrepancies are observed in GM, it remains impossible to ultimately judge which method is the most accurate, but comparisons with a third method could provide insights if not conclusions. Anatomical knowledge as well as previous MRI studies 22 suggest that the dominant diffusion directions in GM should be perpendicular to the cortical folding, with exceptions such as the primary visual cortex 15 . We therefore derived the normal vectors to the cortical folding by using an in-house method based on the structural tensor 36 and computed the minimum dyadic angle between the absolute value of the cortex normal vectors and the absolute value of the vectors dictated by the first three diffusion directions of MSCSD and mFOD-WS.

Tracking parameters
We investigated the results obtained by performing fiber tractography with MSCSD, mFOD1, mFOD2 and mFOD-WS both qualitatively and quantitatively, with intentional focus on GM areas. Deterministic fiber tractography was performed in ExploreDTI by using identical parameters for all methods: FOD threshold 0.1 (default), minimum/maximum fiber length equal to 20/500mm, angle threshold 45 degrees, step-size 0.6mm.

Quantitative analysis
Firstly, FODs computed with MSCSD and mFOD-WS were seeded in a frontal region and color encoded with a custom colormap assigning red, yellow and blue to tracts traversing WM, GM and CSF, respectively. Subsequently, whole brain tractography of all FODs was performed by seeding all voxels with an FOD value above the default threshold. The tracts were color encoded both with the conventional RGB color scheme 37 and the above-mentioned tissue-type color scheme, to convey both directional information and specificity of the FODs to different tissue types.
Tracts in a 3.75mm wide sagittal slab were color encoded by normalized distance between the GM/WM interface and the WM/CSF interface, obtained from the T1-weighted segmentation, as well as by tissue type, to appreciate differences in GM tracking distance / density between the compared methods.
Quantitatively, we evaluated how many tracts were generated with MSCSD and mFOD-WS, which is informative given the seed-points were kept constant and a deterministic tracking algorithm was used. Additionally, we counted the fraction of plausible tracts, i.e. tracts having at least one end-point in grey matter, as well as the spatial coverage of GM when considering only start/end points of tracts.

Results
Simulations Figure 1 shows the signal fractions estimated with mFOD in simulation I, where different mixtures of WM, GM and CSF were simulated for multiple SNR Levels. At SNR level equal or greater than 20, mFOD disentangles mixtures of WM and CSF with signal fraction errors below 2% and angular errors below 3 degrees. The separation of GM from both WM and CSF is more challenging and requires higher SNR to achieve an error below 10%. In both cases, the signal fraction errors at SNR 50 are about 9%, 11% and 2%, for WM, GM and CSF, respectively. At the same SNR level, the angular error on the GM FOD is 8 degrees for the GM and CSF mixture. With regards to the GM and WM mixture, the angular errors of the WM and GM FODs are 4 degrees and 10 degrees, respectively.   HCP data

Signal fractions
The WM, GM and CSF signal fractions estimated with MSCSD and mFOD, as well as those derived from the segmentation of the T1-weighted image are shown for an example axial slice in Fig. 3. The maps estimated with mFOD exhibit excellent agreement with both MSCSD and the T1-derived segmentation, despite the numerical ill-conditioning of the mFOD problem. The WM fractional maps estimated with both MSCSD and mFOD have non-zero values also in some GM areas, in contrast to the T1-derived segmentation.

FODs
The FODs estimated with MSCSD and mFOD in two cortical areas are shown in Fig. 4 and Fig. 5, for an example axial slice and an example coronal slice, respectively. The WM fiber orientation estimated with mFOD (mFOD1) show remarkable similarity to the FOD estimated with MSCSD in WM. Both FODs have very sharp and comparable orientations in the vicinity of WM, but they tend to zero in deeper cortical GM. In contrast, the GM fiber orientation with mFOD (mFOD2) show dominant and physiologically plausible peaks in most GM voxel, with orientation mostly perpendicular to that of the underlying WM. When mFOD1 and mFOD2 are linearly combined (mFOD-WS), a continuous vector field from WM to deep cortical GM is revealed.  and mFOD2 refer to the FODs computed when using response functions meant to represent WM and GM, respectively. mFOD2 is non-zero only in proximity of GM, and shows dominant orientations also in the deeper cortex, in contrast to both MSCSD and mFOD1. mFOD1 provides directional information in strong agreement with MSCSD but occasionally shows more peaks. Fig. 6 shows an analysis of the peak orientations obtained with MSCSD and mFOD-WS as compared to the spatial direction normal to the cortical folding derived from the T1-weighted image with a structural tensor approach. When comparing the peak direction of MSCSD and mFOD-WS in WM, we observe that more than 70% of the WM voxels have an angular deviation smaller or equal than 10 degrees. In GM, the angular deviation between the two methods is larger, and mFOD-WS result in smaller angular deviation with the cortex normal than MSCSD.

Fiber tractography
The results we obtained by performing deterministic fiber tractography on the FODs computed with MSCSD and mFOD-WS in a frontal lobe region are shown in Fig. 7. The tracts, which are color encoded according to the tissue type they are spatially traversing, are more densely reconstructed with mFOD-WS. Further, the tracts with mFOD-WS cover the GM more extensively and uniformly, and traverse the tissue with a more plausible shape to a deeper extent. Subsequently, we performed whole brain fiber tractography of the FODs derived with MSCSD and mFOD, both individually (mFOD1, mFOD2) and after a weighted sum procedure (mFOD-WS). The tracking of mFOD1 produces a whole brain reconstruction that is largely comparable to that of MSCSD, suggesting preservation of the WM tracking performances. The tracking of mFOD2 results mostly in cortical reconstructions and short U-fibers, in-line with the tissue-specific nature of this FOD. The deterministic fiber tractography of mFOD-WS results in 694752 reconstructed tracts, 183821 more than those reconstructed with the FOD derived with MSCSD. After filtering out implausible tracts -i.e. tracts not starting or ending in GM, mFOD-WS reconstructed 594199 valid tracts, 130414 more than MSCSD. Among these, 42.1% have the first or the last tracked point in GM, as compared to 34.1% when performing tractography with MSCSD. This is qualitatively visible in the tissue-type encoded reconstruction of Fig. 8, which shows that mFOD-WS covers more extensively the GM as compared to MSCSD, although introducing some extra spurious fibers. An example sagittal slab of the tracts reconstructed with MSCSD and mFOD-WS was obtained by cropping the whole brain fiber tractography, as shown in Fig. 9. In the first column, the tracts are color encoded according to their distance from the WM/GM interface (blue) to the GM/CSF interface (red). The results confirm that the tracts reconstructed with mFOD-WS consistently travel a longer distance within the cortical GM than with MSCSD, in-line with observations in Fig. 7. Further, the frontal, parietal and occipital side of the cortex appear to be better reconstructed with mFOD-WS. The first column shows the tracts color encoded according to their normalized distance from the GM/WM and the GM/CSF interfaces. The distance was determined in 3D as flow between the two surfaces using the same method we previously used to determine the cortex normal. The central and right column show the same sagittal slab color encoded by tissue type, and its zoom in a cortical region. The WM-GM flow color encoding shows that a deeper cortical tracking was achieved with mFOD-WS as compared to MSCSD, especially in the frontal and occipital areas. When looking at the tissue type encoded images, mFOD-WS resulted in a more complete and contiguous reconstruction of the cortical folding.

Discussion
In this work, we have introduced a framework to perform spherical deconvolution of multiple anisotropic response functions which allows to reconstruct multiple fiber orientation distributions (mFOD) for each voxel. With simulations and data from the Human Connectome Project, we have shown that i) it is possible to reconstruct two anisotropic FODs specific of WM and GM while also accounting for CSF partial volume, given sufficient SNR (greater or equal than 50 with respect to the b = 0 s/mm 2 data) is achieved. Further, we have shown that ii) mFOD improves the reconstruction of FODs in GM as compared to existing stateof-the-art methods, and that iii) this results in improved fiber tractography reconstruction in GM.

Feasibility of mFOD
The reconstruction of multiple anisotropic FODs is a hard numerical problem that requires the inversion of a matrix whose rank (number of measurements) is much lower than the number of possible solutions (FOD amplitudes on the unit sphere), a recurring challenge in deconvolution problems. Previous studies, as Tournier et al. 27 and Jeurissen et al. 23 used an L2 regularization term enforcing contiguity (smoothness) in the FOD amplitudes. In this work, we used a non-negative least-squares solver with L2 regularization encapsulated in an iterative solver. Taking into consideration that the mFOD formulation (eq. [2]) hosts multiple FODs, this might be suboptimal and introduce discontinuities in the solution. Nevertheless, our results suggest that mFOD does not propagate smoothness between multiple FODs, likely due to the attenuation of small perturbations operated by the median operator 29 . Indeed, the results of simulations, shown in Fig. 1 and Fig. 2, suggest that this approach allows to effectively disentangle the partial volume effects between mixtures of isotropic (n=1) and anisotropic (n=2) signal mixtures, if sufficient SNR is provided. In particular, the partial volume effects between WM and CSF are the easiest to disentangle, and an SNR equal or greater than 20 is sufficient to obtain accurate signal fractions separation and minimal angular error on the WM FODs. In contrast, the separation of GM from both WM and CSF is a harder task, probably as a result of the low anisotropy of this tissue type. Nevertheless, if an SNR equal or above 50 is achieved, mFOD can excellently preserve the ratios between different tissue types and estimate their mixtures with minimal biases (Fig. 2). While an SNR of about 50 might be difficult to achieve in clinical settings, accurate pre-processing and denoising 38 , or even deep-learning based reconstructions 39 might mitigate this need. Furthermore, in both simulation I and simulation II the diffusion direction of GM was set perpendicular to that of WM, but in-vivo the two are likely to coincide as white matter bundles travel to their cortical endpoint perpendicularly 40 , with exceptions such as the primary visual cortex 15 or the Ufibers adjacent to the cortex 41 . If the direction of the axons leaving the cortex (to the WM) wouldn't match that of intra-cortical axons, the results of simulations suggest that the GM FOD would be estimated with angular error between 10 (SNR 30) and 5 degrees (SNR 50), which are likely to perturbate the tracking trajectory. Furthermore, Fig. 3 shows that with in-vivo data the signal fraction map associated to GM is zero-valued in all non-GM voxels, and that it assumes a value of 1 in almost all GM voxels, which suggests excellent separation of the two tissues.

Why modelling GM?
The unique characteristic of mFOD is indeed the ability to reconstruct multiple FODs specific of the different employed signal models, as shown in Fig. 4 and Fig. 5. The figure shows that WM fiber orientations estimated with mFOD (mFOD1) and MSCSD are very similar in GM, and that both rapidly assume a value of zero in GM with the default FOD thresholds. In contrast, the GM fiber orientations estimated with mFOD (mFOD2) show consistent and well-defined diffusion directions for a large extent of the cortical GM. The performance improvement over MSCSD of mFOD in GM is explained by considering the improved tissue modelling in mFOD. Indeed, approaches as MSCSD essentially describes the WM and GM mixture as the sum of an anisotropic and an isotropic signal, similarly to the ball-and-stick approach 42 , where the anisotropic signal is modelled with properties extracted from single fiber population WM 27,43,44 . In practice thus, by applying any constrained spherical deconvolution method in GM, one assumes the axonal component of GM to produce a similar signal to that of a WM bundle without crossing fibers. However, assuming cortical axons and dendrites to be identical to axons in the corpus callosum or part of the corticospinal tracts might be a stretch given their different diffusion signature 45 . Differently, in this work we represent the GM signal with a special case of the NODDI model 33 , which represents a first attempt to take into account additional effects as tortuosity 46 and intra to extra neurite signal ration into the deconvolution problem.
Multiple FODs: how do we use them?
One of the open challenges in the mFOD framework regards the choice of how to combine the multiple FODs computed by the method. In this work, we adopted a pragmatic approach and considered their weighted sum (mFOD-WS). This choice produced convincing and contiguous FODs, as shown in Fig. 4 and Fig. 5. However, this might be sub-optimal for fiber tracking approaches, as it might blur the peaks detected in the individual FODs. Future work is needed to establish whether dynamically tracking the locally larger FOD, or dynamically choosing the FODs most aligned to the tracking direction might represent better alternatives. Additionally, sparse solvers as those including L1 or L0 norm 47 could be taken into consideration to improve the FOD separation and minimize overfit. Nevertheless, when we compared the main peaks of mFOD-WS to those obtained with MSCSD (Fig. 6), we observed remarkable similarity between the two in WM, where an angle error smaller than 10 degrees was observed in over 70% of WM voxels. This suggests that mFOD-WS does not worsen the performances of spherical deconvolution in WM. In GM, mFOD-WS resulted in smaller angular deviations with respect to the cortex normal as compared to MSCSD. While the cortex normal might not be always representative of the direction with which axons are organized in the cortex, they represent an accepted assumption 14 and might serve as a comparison term between the two methods. Further, we observe that the average of the angular deviation between the cortex normal and mFOD-WS is 20 degrees, which is in-line with the angular error predicted in case of WM and GM mixtures at SNR 50 ( Fig. 1).

Improved fiber tractography in GM
Fiber tractography is one of the main applications of spherical deconvolution techniques. When seeding deterministic tractography in a frontal region (Fig. 7), we obtained more tracts with mFOD-WS than with MSCSD with identical FOD thresholds. Further, the tracts reconstructed with mFOD-WS covered the cortical structure more extensively, and their geometrical configuration appeared more similar to previous reports of high resolution cortical dMRI 22 . Eventually, a similar tractography result might be achievable also with MSCSD by tuning the FOD threshold 48 . However, this is likely to introduce not only false positives in WM, but also to worsen the tracking in GM itself. Indeed, being the model used for GM not optimal, classic deconvolution methods are likely to produce a large number of spurious peaks 49 . When performing whole brain tractography, the use of mFOD-WS resulted in more tracts having a start or endpoint in the GM (+8%), and generally in more extensive GM coverage, as shown in Fig. 8. This was even more apparent when looking at a sagittal slab (Fig. 9), which showed how tracts could be followed in GM up to the GM/CSF interface, at the price of more spurious tracts as compared to MSCSD. Practically, this could have a major impact when performing fiber tractography from/to the cortex is critical as, for instance, pre-surgical planning 50 or structural connectivity studies 51 . Indeed, existing fiber tractography methods do not seed the tracking procedure in the cortex, but rather at the GM/WM interface or in the adjacent WM, potentially introducing a number of false positives 52 . For similar reasons, connectivity studies are often performed by first seeding WM or the WM/GM interface, then filtering the tracts with anatomical constraints [cit. ACT], and finally assigning the tracts to the closest GM node after spatial dilation, potentially introducing uncertainty in the assignment [cit.]. mFOD may contribute addressing these short comes, allowing to follow WM bundles with higher accuracy than existing methods from/to their cortical start/endpoint.

Limitations
To follow, we discuss some limitations of this work. The models used to describe WM and GM were arbitrarily chosen among popular methods, but their combination might not be the optimal representations to jointly estimate multiple FODs. The DKI model, here used in its simplified version with isotropic kurtosis 32 , efficiently represents non-Gaussian signals with a limited number of parameters but is sub-optimal at accounting for direction dependent signal restrictions. Alternatively, geometric representations of the WM signal 53 , as well as data driven methods 44 could also be explored. The choice of models here suggested for both WM (DKI) and GM (NODDI) is not meant to be final but is rather a proof of concept of the feasibility of mFOD with tissue specific models.
In Fig. 6, we used the cortex normal direction as a comparison for FODs computed with both MSCSD and mFOD. While it is reasonable to assume that the most axonal projects will traverse the cortex with a perpendicular trajectory, this is not always the case. Further, the method here employed to compute the normal trajectory is not accurate when the sulci are separated by less than 2 voxels and should therefore interpreted in sight of his limitations. The performances of mFOD have here been shown on a single subject but should be investigated on a larger cohort and with more clinically achievable protocols to generalize the observed results. Finally, while we have shown that mFOD improves the performances of spherical deconvolution and fiber tractography in GM, it remains to be proved to which extent such improvements will benefit clinical studies and structural connectivity investigations.
In conclusion, we have presented a framework (mFOD) to perform spherical deconvolution of multiple user chosen models, to simultaneously reconstruct multiple fiber orientation distributions. By accounting for WM and a GM specific models, we have shown that mFOD can improve the performances of spherical deconvolution in GM while retaining state-of-the-art methods in WM.