Implicit neural representation of multi-shell constrained spherical deconvolution for continuous modeling of diffusion MRI Open Access

Diffusion-weighted magnetic resonance imaging (dMRI) uses directional magnetic gradients to measure the diffusion of water molecules in that direction (Novikov et al., 2019). These acquisitions are sensitive to diffusion at the micrometer scale. Therefore, one can gain insight into the (micro-)structure of the underlying tissue by combining many of these directional acquisitions. In dMRI, a b-value relates to the degree of diffusion weighting applied as a function of amplitude (G), time (⁠$δ$⁠), and duration between the gradients (⁠$Δ$⁠). Using multiple b-values, or shells, provides additional information, as different tissue types respond differently to increasing or decreasing b-value (Genc et al., 2020; Schilling et al., 2017).

The dMRI signal can be used to estimate the parameters of an ever-increasing selection of local (voxel-level) microstructural models (Alexander et al., 2019). Each model provides different insights and has its strengths and weaknesses. For white matter, it is particularly interesting to know how fiber bundles are oriented within a voxel. If there are multiple bundles, which is estimated to occur in 60 to 90 percent of voxels (Jeurissen et al., 2013), it is valuable to know how these bundles are angled and what their relative sizes are (Jeurissen et al., 2013). Once known, the fiber orientations can be used to evaluate the local properties directly or in downstream tasks such as tractography. Tractography aims to reconstruct macro-level fiber bundles by computing streamlines that follow the local fiber orientation (Jeurissen et al., 2019).

Constrained spherical deconvolution provides a method to calculate the fiber orientation distributions (FODs) for every individual voxel in a single-shell acquisition (Tournier et al., 2007). This method was later expanded to multi-shell multi-tissue CSD (MSMT-CSD), which uses multi-shell data to parameterize FODs for multiple tissue types (Jeurissen et al., 2014). Both methods require the estimation of response functions, which represent the expected signal for different tissues and b-values. CSD only requires a WM response function (obtained from voxels with a single fiber population), while MSMT-CSD requires a response function for every tissue and b-value. The response functions are deconvolved from the dMRI signal to obtain FODs. These FODs are most commonly used for tractography algorithms to reconstruct the structural architecture of the brain or to make inferences about the local microstructure.

Using (MSMT-)CSD often results in high-quality (i.e., sharp) FODs, an accurate estimation of partial volume effects, which will positively impact the results of tractography (Jeurissen et al., 2014). Therefore, this technique has become a staple in analyzing dMRI data. However, (MSMT-)CSD operates at the voxel level without using the structural continuity of the brain tissue. This means it is susceptible to noise, as a single noisy voxel can be completely different from and misaligned with its neighbors. This can negatively affect downstream tasks (e.g., fiber tracking, Maier-Hein et al., 2017). Fitting the FODs at locations where data quality is not optimal by using the information in the neighboring voxels (i.e., in noisy environments or low angular and spatial resolution) could help address this issue. Finally, for tractography, it is necessary to calculate FODs in locations that are not voxel center points. In the case of (MSMT-)CSD, these are not readily available, requiring interpolating the dMRI signal or the FODs. Interpolation is computationally costly, and interpolating FODs might produce artificial peaks in areas with strongly curved fibers (Nie & Shi, 2023).

Recently, many machine-learning models have been proposed to increase the quality of the fit of the local model (Karimi, 2024). For FOD estimation, most approaches are based on supervised learning, which requires a training dataset to learn the model parameters. While promising, these methods are not always ideal: they often require a large amount of data that might impose unwanted biases on the results upon inference, leading to erroneous predictions. To circumvent this problem, unsupervised methods, such as Elaldi et al. (2021) and Elaldi et al. (2024), are fit on a single dMRI dataset. Using a UNet-type architecture, the model has access to the spatial information of the dMRI signal to create accurate FOD estimations. However, this approach does not continuously represent the underlying volume as UNets operate in the discretized voxel space. These models, therefore, do not solve the interpolation problems.

Implicit neural representations (INRs) provide a way to use neural networks to represent a dataset continuously. Successful usage of INRs outside of the field of dMRI is seen in neural radiance fields, which capture 3D scenes using sparsely sampled 2D images as input (Mildenhall et al., 2020). The robustness to noise and accurate reconstruction of unseen parts of the scene makes INRs a promising avenue of research in dMRI, which deals with noisy, sparsely sampled data. INRs have already been used to compress dMRI data (Mancini et al., 2022), to model dMRI signal (Hendriks et al., 2023), and to model diffusion ODFs (Consagra et al., 2024). While showing promising results, Hendriks et al. (2023) have the limitation of modeling the diffusion signal directly, which limits the usability of the model output as further modeling is required. The diffusion ODFs output by the model from Consagra et al. (2024) are more informative but also require post-processing (i.e., sharpening) to be used in fiber tracking. Additionally, the approach is limited to using single-shell data.

In our work, we propose a method that benefits from the strengths of the CSD and the INR frameworks. We first introduce the methodology for creating a continuous parameterization of CSD and MSMT-CSD using INRs. This results in an unsupervised, response function-informed machine-learning model representing a single (multi-shell) dMRI dataset. This continuous model can be sampled at any coordinate to obtain an FOD directly usable for fiber tracking or microstructural analysis. We evaluate the correctness of the learned continuous representations and the robustness of the model to sparsity in both spatial and angular domains, noise, and a changing environment such as different acquisitions, scanners, and patients. Finally, we apply our model to a ‘real-world’ quality dataset to show its application in a realistic scenario.

In the upcoming sections, we first briefly overview the (MSMT-)CSD model (section 2.1). Then, we expand upon previous work using INRs to create a novel approach that fits the parameters of the (MSMT-CSD) model in a continuous volume and can leverage the spatial correlations for parameter estimation (section 2.2). The goal is to create an INR that outputs FODs ready for downstream analyses (i.e., to extract metrics (Raffelt et al., 2012) or perform tractography) at any desired location in the volume.

CSD estimates the distribution of fiber orientations present within each voxel. The distribution of these fiber orientations can be modeled with an FOD, a spherical function with an amplitude correlating to the density of the fibers in a given direction $(θ,ϕ)$⁠.

In diffusion MRI, these spherical functions are typically described using a truncated, real-valued, even ordered, spherical harmonics (SH) series. We use a basis identical to the one used in MRtrix3 (Tournier et al., 2019):

where $Ylm(θ,ϕ)$ is the SH basis function of order $l$ and phase¹$m$⁠. Any spherical function can then be represented as:

where $G(θ,ϕ)$ is a generic spherical function, $glm$ is the SH coefficient for the SH basis function with order $l$ and phase $m$⁠, and $lmax$ is the maximum order used in the SH series. The number of coefficients necessary to parameterize the SH-series depends on $lmax$ and can be found using $12(lmax+1)(lmax+2)$⁠. Given the dMRI signal of a single fiber population, called the ‘response function,’ the dMRI signal in a voxel can be reconstructed by doing a spherical convolution of the single fiber signal with the FOD:

where $S(θ,ϕ)$ is the measured signal, $F(θ,ϕ)$ the FOD, and $R(θ)$ the response function. The estimate for $R(θ)$ can be obtained from the dMRI data using various existing algorithms. Since $R(θ)$ is axis-aligned and symmetrical by construction, the spherical harmonic series for $R(θ)$ only depends on the harmonic phase $m=0$ terms of the series, the so-called zonal harmonics. The spherical convolution of a spherical harmonics series with a zonal harmonics series is reduced to a simple multiplication of coefficients:

where $slm$ and $flm$ are the SH-coefficients parameterizing $S(θ,ϕ)$ and $F(θ,ϕ)$⁠, respectively, $rl0$ are the zonal harmonics coefficients parameterizing $R(θ)$⁠. The scaling factor in front corresponds to the l-th order coefficients for the spherical harmonics series for a Dirac-delta function. The FOD can be retrieved by deconvolving the measured dMRI signal with the single fiber population signal. The retrieved FOD is constrained to be strictly positive by penalizing negative values during optimization, thereby filtering high-frequency noise that generates these negative values. Details about the CSD algorithm can be found in Tournier et al. (2007). MSMT-CSD expands this idea by using data with multiple b-values and response functions to get more accurate estimations of the FOD and separate volume fractions for other tissue types, such as grey matter (GM) and cerebrospinal fluid (CSF). The method then calculates an orientation distribution for every tissue type. These distributions are convolved with their respective response functions. The response functions in the case of MSMT-CSD are not only parmeterized by $θ$⁠, but also by the b-value. The convolution provides the expected signal for each shell. Since the GM and CSF are presumed to be tissues with (restricted) isotropic diffusion, an orientation distribution of a single coefficient suffices per tissue. These coefficients are often seen as a volume fraction of these tissue types. Details about the MSMT-CSD algorithm can be found in Jeurissen et al. (2014).

In our work, we represent a full multi-dimensional dMRI dataset using an INR that outputs FODs for any spatially encoded 3D input coordinate (Fig. 1). This borrows concepts from Neural Spherical Harmonics, previously used to model the dMRI signal directly (Hendriks et al., 2023). We extend the existing approach by adapting the model architecture to output FODs and introduce the (MSMT-)CSD framework to relate the FODs to the measured signal. This injects additional assumptions (e.g., the response functions) into the INR fitting process. In the upcoming sections, we expand upon the positional encoding (Fig. 1a), the adapted model architecture (Fig. 1b), and the fitting of the model (Fig. 1c).

The model input is a 3D-coordinate $x=[x,y,z]$⁠. Using just the raw coordinates as input leads to a poor representation of the dataset as the model cannot capture the higher frequency details. Previous studies show that encoding the coordinates using a multitude of sine and cosine functions with different frequencies significantly improves the representation (Mildenhall et al., 2020). In this work, we leverage existing work based on frequencies randomly sampled from a Gaussian distribution to prevent an axis alignment bias during encoding, which has shown to be effective at representing different types of datasets (including conventional MRI) (Tancik et al., 2020). The encoding is shown in (5).

where $γ(x)$ returns the encoded coordinate vector $x$ (which is normalized to lie in range $[−1,1]$ prior to encoding), and $B$ is a size $np×d$ matrix with values sampled from $N(0,σ2)$⁠. The hyperparameter $np$ determines the number of positional encodings, $d$ is the number of dimensions in coordinate $x$⁠, and $σ$ is the standard deviation of the Gaussian. The hyperparameters $np$ and $σ$ can be tuned to generate the desired representation of the data, where $np$ is more related to dataset complexity, and $σ$ to the granularity of the representation (i.e., higher values for $σ$ capture finer details, but are more noisy) (Tancik et al., 2020).

We define a multi-layer perceptron (MLP) that maps positionally encoded input coordinate $γ(x)$ to an SH-coefficient vector $k$ (6)

where $ℳΨ$ is the MLP parameterized with weights $Ψ$⁠, and $k$ is given by:

where $flm$ is the vector of coefficient of the SH basis function of order $l$ and phase $m$⁠, and $fgm$ and $fcsf$ are the coefficients for the GM and CSF, respectively.

We use two similar but slightly different model architectures for the single-shell and multi-shell implementation. Both models take $np$ positionally-encoded 3D coordinates as an input, followed by an MLP consisting of four layers with ReLU activation functions. The size of each layer is a tunable hyperparameter. For the single-shell model, the MLP directly outputs SH-coefficients that parameterize the FOD at the input coordinate. For the multi-shell model, the output of the MLP is used as input for three output heads, one for each tissue type, consisting of one fully connected layer. These heads are optimized alongside the MLP during model fitting. The WM head outputs a full SH-series parameterizing the FOD, while the GM and CSF heads output a single coefficient. The complete model architecture can be seen in Figure 1b.

The model fits the dataset by calculating a mean-squared error (MSE) between the measured diffusion signal $S(θ,ϕ)$ and the diffusion signal reconstructed from the FODs $S^(θ,ϕ)$⁠. In single-shell CSD, $S^(θ,ϕ)$ is reconstructed by performing a convolution between the output FOD and the response function as shown in (4). This model will be referred to as INR-CSD.

For MSMT-CSD, the model outputs $12(lmax​+1)(lmax​+2)+2$ coefficients (a full SH-series for the WM and a zeroeth-order SH-series for both GM and CSF). The signal is obtained for each b-value individually as the sum over $T$ tissue types (three in the case of MSMT-CSD) of the convolution of the tissue FODs with the corresponding tissue response functions belonging to that b-value. The multi-shell output for a single measured coordinate, therefore, is:

where $B$ is the number of b-values and $S^i(θ,ϕ)$ is the signal for b-value $i$⁠. We obtain $S^i(θ,ϕ)$ from

where $Fj(θ,ϕ)$ is the FOD for tissue type j, parameterized by the corresponding coefficients in $k$⁠, and $Ri,j(θ)$ is the response function for b-value $i$ and tissue type $j$⁠. This process is shown in Figure 1(b, c). This model will be referred to as INR-MSMT.

Since $k$ and therefore $Fj$⁠, are dependent on 3D coordinate $x$⁠, we can define the signal in the entire volume as:

Finally, we constrain the FOD(s) to be strictly positive as described in Tournier et al. (2007). The (WM) FOD is sampled in 724 directions, and any values larger than 0.1 times the mean amplitude of the FOD are included in the loss. The GM and CSF are, by construction, non-negative.

Given a set of b-vectors with directions $(θ,ϕ)$ in $B$ b-values, where b-value $i$ contains $Qi$ b-vectors,

for a set of $N$ 3D coordinates $X={x0,…,xN}$ (one for every voxel in a dMRI dataset), we obtain the parameters for the INR $ℳΨ*$ by minimizing the MSE alongside the constraint:

where $Si(θi,q,ϕi,q,x)$ and $S^i(θi,q,ϕi,q,x)$ are, respectively, the measured signal and the reconstructed signal at coordinate $x$ for b-value $i$ in direction $(θi,q,ϕi,q)$⁠, and $Λx$ is the non-negativity constraint.

We use a synthetic high-resolution brain-like dataset to introduce our methodology and perform a quantitative analysis of its performance regarding sparsity and noise. The diffusion-weighted images are generated using Fiberfox (Neher et al., 2014), from curated tractograms that have been used previously in dataset generation of the 2015 ISMRM challenge (Maier-Hein et al., 2017; Renauld et al., 2023). Since we are interested in both single-shell and multi-shell CSD, the dataset we generate for this experiment has 67 volumes spread over three b-values: 7 for $b=0s​/​mm2$⁠, 30 for $b=1200s​/​mm2$⁠, and 30 for $b=3000s​/​mm2$⁠. For consistency, these are the same gradient directions and b-values used in real brain acquisition in the real-data experiments. We generate the dataset in a 1.25 mm isotropic resolution and add no artifacts or noise. The rest of the settings are identical to those used for generating the challenge data. Additionally, a dataset is generated at 2.5 mm isotropic resolution for the low spatial resolution experiment. Finally, three datasets are generated using Fiberfox at 1.25 mm isotropic resolution with simulated complex Gaussian noise, resulting in three datasets with signal-to-noise ratios (SNRs) of 15, 25, and 35 for the noise-related experiments.

Additionally, we use the computational diffusion MRI 2018 (cdMRI2018) challenge preprocessed dataset (Ning et al., 2020; Tax et al., 2019) to explore model performance and robustness in a realistic setting. Two randomly chosen subjects (B and K) were scanned on two scanners (3T Siemens Prisma and Connectom) with two acquisition settings per scanner: standard (ST) and state-of-the-art (SA). This provides a baseline for comparing the performance under varying acquisition parameters. Both scanners have ST acquisitions of 2.4 mm isotropic resolution with 67 volumes spread over three b-values: 7 for $b=0s​/​mm2$⁠, 30 for $b=1200s​/​mm2$⁠, and 30 for $b=3000s​/​mm2$⁠. Additionally, the Prisma scanner has a 1.5 mm isotropic resolution, 136 volume SA acquisition spread over three b-values: 16 for $b=0s​/​mm2$⁠, 60 for $b=1200s​/​mm2$⁠, and 60 for $b=3000s​/​mm2$⁠. Finally, the Connectom scanner has a 1.2 mm isotropic resolution SA acquisition, with volumes identical to the Prisma SA acquisition.

For the experiments using synthetic datasets described below, we perform a quantitative analysis of the resulting FODs by comparing them to the FODs generated on the clean, high-resolution multi-shell dataset using voxel-wise MSMT-CSD. We assume these FODs to be the ground truth (GT). We report the total apparent fiber density (AFD) as a measure of FOD scale, the angular correlation coefficient (ACC) as a measure of angular correctness, and the number of fiber orientations (NuFO) to compare the peaks used for tracking (Anderson, 2005; Dell’Acqua et al., 2013; Raffelt et al., 2012). Furthermore, we visualize the FODs in all experiments to make a qualitative assessment.

We use $lmax=8$ for every SH series in this paper. This choice is informed by previous work showing that higher orders have low contributions to the signal (Tournier et al., 2013). The spherical harmonic basis functions are obtained through scipy².

The response functions for single-shell CSD are retrieved from the datasets using the Tournier algorithm, and for the multi-shell algorithm using the Dhollander algorithm, both implemented in MRtrix3 (Dhollander et al., 2019; Tournier et al., 2013).

The INR-models in all experiments are implemented in PyTorch³. An Adam optimizer was used with a learning rate of 1 $×$10^-4 and default settings otherwise. Further hyperparameter settings are determined in the first experiment and used throughout the rest.

The AFD and NuFO are calculated using scilpy⁴, with an absolute threshold of 0.07 and default settings otherwise. AFD (total) is defined as being the first coefficient of the SH-series (Raffelt et al., 2012). The ACC is calculated between two SH-series with coefficients $ulm$ and $vlm$ as:

which corresponds with Anderson (2005) and Schilling et al. (2019). The metrics are computed in all voxels where the ground-truth FOD has a first coefficient greater than 0.05 for the WM.

The feasibility of our approach is first assessed by fitting the INR-CSD and INR-MSMT models for a full volume. Our approach is feasible if the model can reconstruct the original dMRI signal. The INR-CSD is fit to the $b=3000s​/​mm2$ shell of the 1.25 mm isotropic resolution synthetic dataset. The outputs are visualized for an area of interest and compared qualitatively and quantitatively to the GT. This experiment is repeated for INR-MSMT using all three b-values (⁠$b=0,1200,3000s​/​mm2$⁠) of the 1.25 mm synthetic dataset.

For the INR to handle spatial sparsity well, it must give valid outputs between the voxel positions for which it was fit. To assess this, we fit an INR-CSD model on the $b=3000s​/​mm2$ and an INR-MSMT on all three b-values of the 2.5 mm isotropic synthetic dataset. As a first indication of intervoxel correctness, we visualize the FODs resulting from supersampling the model at 1.25 mm isotropic resolution in two regions of interest. These FODs are compared to a set of FODs obtained by linear interpolation of the SH-coefficients of the FODs calculated at 2.5 mm resolution using voxel-wise CSD. The AFD and the ACC are computed between both models, and the GT and are reported alongside the NuFO. We then explore the intervoxel space by 1) selecting two voxels in each region, 2) sampling the model at multiple equidistant steps between these voxels, and 3) visualizing the FODs at those locations. Additionally, we show maps of the GM and CSF volume fractions for the ground truth, INR-MSMT, and linearly upsampled MSMT-CSD volume fractions.

To evaluate the performance of the INR in settings with low angular resolution, we iteratively remove b-vectors from the dataset. At every step, we remove the b-vectors with the smallest effect on the hemispherical coverage. Since we expect model failure only to occur at considerably sparse data, the starting point for our evaluation is a single shell of 30 b-vectors at $b=3000s​/​mm2$ in the 1.25 mm isotropic synthetic dataset. For every subset of b-vectors, we fit an INR-CSD and a voxel-wise CSD model. The AFD and ACC are calculated between both models and the GT for different numbers of b-vectors. The ACC maps, FODs, and NuFO are also reported for a single set of b-vectors of interest.

We evaluate the robustness of our approach to noise by fitting the INR-CSD to synthetic datasets with an SNR of 15, 25, and 35. The resulting FODs are compared visually to the FODs of voxel-wise CSD. For the dataset with SNR = 25, we report the AFD and the ACC calculated between both models and the GT, and NuFO.

We aim to show the performance of the model across different scanners and acquisitions when using identical hyperparameters. We fit our INR-MSMT model to the cdMRI 2018 dataset (Ning et al., 2020). Two subjects were randomly selected, each with two acquisitions on two scanners, resulting in 8 different INR-MSMT fits. To assess the robustness of our model, we qualitatively analyze the consistency of the FODs within the subject and the difference between subjects by visualizing them. The FODs are compared to FODs obtained through voxel-wise MSMT-CSD.

Having a multi-shell high-resolution dataset is not common in day-to-day practice. To test the performance of the INR on a ‘real-world quality’ dataset, we fit an INR-CSD on a single $b=3000s​/​m2$ shell of the ST of the Connectom scanner and sample it at SA resolution. The resulting FODs are compared to a set of FODs obtained by linear interpolation to SA resolution of the SH-coefficients of the FODs calculated at ST resolution using voxel-wise CSD. We show the MSMT-CSD fit on the SA data as an additional comparison, which corresponds to FODs obtained in a ‘best-case scenario’ in current practice.

The following hyperparameter settings were found using a grid search and used for fitting both models: 1024 neurons per layer, four layers, $σ=4$⁠, $np=5000$⁠. The influence of these parameters is detailed in Appendix A. The single-shell and multi-shell INRs converge after approximately 650.000 iterations over 30 epochs with batch size 500 on a 16 GB NVIDIA RTX 4080. Training times are 2 m 58 s and 5 m 24 s respectively. Figure 2 shows the ACC map and FODs in the centrum semi-ovale (CSO) for the single-shell model alongside the GT FODs. Quantitative results for both models are shown in Table 1. We observe minimal differences between the angles of the fiber populations and the relative sizes for this feasibility experiment, compared to the ground truth. Quantitative analysis confirms that the errors of the INRs in representing the target data are low.

We used hyperparameter settings identical to those used in the previous experiment to fit the INRs. Figure 3a, b, c shows the error maps of the AFD for both models alongside the difference in error. Our method has a smaller error in areas where fibers cross but larger in areas with a single-fiber population. Figure 3 also shows two regions of interest: an area of crossing fibers and an area of strongly curved fibers. We observe that INR-CSD can correctly interpolate crossing and curving fibers at subvoxel levels. Linear interpolation does not correctly interpolate the FODs in curved fibers, creating spurious peaks as shown in the yellow box in Figures 3 and 4 that shows maps of the GM and CSF volume fractions of the MSMT-CSD and INR-MSMT fits. Both show a loss of detail compared to the ground truth. Notably, the INR-MSMT has a higher estimate for the GM volume fraction and a lower estimate for the CSF volume fraction in areas where brain tissue is present. The quantitative results in Table 1 show identical AFD values and a slightly larger ACC for the INR-CSD for this low spatial resolution experiment, compared to the ground truth. The INR-MSMT shows a minor underestimation of the AFD, with similar ACC scores. The distribution of the NuFO looks similar across methods, showing an underestimation of single-fiber populations and an overestimation of dual-fiber populations. The overestimation of dual-fiber populations is more pronounced in the interpolated voxel-wise models.

The mean ACC and the mean absolute error (MAE) in AFD for both the CSD and INR-CSD model for an increasing number of directions are shown in Figure 5. One can see that the INR-CSD model shows higher ACC and lower MAE in AFD compared to CSD, although differences are marginal. In Figure 6, where the outputs of the models fit on eight directions are shown in an area of low crossing angle (top row), we can see that in regions of low crossing angle, INR-CSD does not find two distinct peaks where CSD (mostly) does (NuFO = 2). Meanwhile, in the area of high curvature (bottom row), we see spurious peaks in the CSD fit, which do not appear in INR-CSD. The ACC maps show an overall higher ACC for the INR-CSD fit, although both models perform sub-optimally in areas of crossing fibers. Table 1 shows an almost identical AFD, while the ACC is marginally higher for INR-CSD for this angular sparsity experiment, compared to the ground truth. NuFO shows a slight shift toward higher numbers for the CSD fit. The first row in Figure 7 shows how spurious peaks show up throughout the WM for CSD.

The FODs in the centrum semi-ovale for different levels of SNR for INR-CSD and CSD are shown in Figure 8. A marked difference can be observed between both methods: the INR-CSD appears to be more robust to noise, showing similar results from SNR = 35 to SNR = 15. In contrast, the FODs obtained through CSD are affected by data quality, especially with lower SNRs. Table 1 shows the quantitative results for SNR = 25, which support the qualitative analysis. The mean ACC for INR-CSD is 0.87, while voxel-wise CSD scores 0.61. NuFO counts show the distribution of fiber orientations of INR-CSD, which more closely resembles the GT but is shifted toward higher population counts. This shift is much more pronounced for the CSD fit, where nearly all voxels have three or more fiber populations. The second row in Figure 7 exemplifies this difference.

The lack of ground truth for a real-world dataset like cdMRI makes quantitative analysis difficult. The qualitative comparison between subjects, scanners, and acquisitions is shown in Figure 9. We see comparable size and sharpness of FODs for both INR-CSD fits between subjects. Lower spatial resolution (in the standard acquisitions) also leads to decreased sharpness of the FODs. The voxel-wise CSD fit on subject K appears to have sharper FODs but displays spurious, less aligned peaks (green arrow in Fig. 9) compared to the INR-CSD fit.

The FODs for INR-CSD, linearly interpolated CSD, and MSMT-CSD are shown in Figure 10. In the CSO (A), we see fewer spurious peaks in the INR fit than in the voxel-wise CSD, especially in the ascending/descending corticospinal tract. Compared to MSMT-CSD, the INR and voxel-wise CSD show less variability between voxels. In area B, the INR fit has a sharper border between the corpus callosum and the CSF in the ventricle than the voxel-wise fit, as indicated by the yellow line in Figure 10B. Area C again shows a reduction in spurious peaks from the INR-fit to the voxel-wise fit, best seen below the apex of the sulcus (yellow box in Fig. 10C). Compared to MSMT-CSD, there is increased coherence between voxels in the INR-CSD but decreased FOD sharpness.

In this work, we propose INR-CSD and INR-MSMT, an approach to create a continuous representation of FODs from a dMRI dataset by combining (MSMT-)CSD with INRs. Our approach differs from recently proposed INR approaches (Consagra et al., 2024; Hendriks et al., 2023). Most importantly, we apply a model-based approach, which allows us to apply dataset and outcome-specific priors (e.g., the response functions) to the training process using multiple shells. The INR outputs the SH-coefficients of the FOD at the desired location directly, bypassing the need for interpolation. The INR can produce FODs representing the underlying structures, even in noisy environments, and when sampled at resolutions significantly higher than the input. Fitting the model to the data can be done in seconds to minutes on consumer-grade hardware based on the input resolution and the number of gradient directions. This, alongside the displayed robustness of the hyperparameter settings, shows that our method can be feasibly used in daily practice. Our method produces similar or better quantitative results in all experiments on synthetic datasets for which the ground truth is available. The visualizations show improvement compared to the baseline in the spatial interpolation of highly curved fibers, where fewer spurious peaks are produced. Noisy (real-world) acquisitions are modeled more consistently in our approach, showing more coherence between voxels and sharper interfaces between bundles and tissues. These results indicate that the proposed method increases the abilities of the model to estimate the coefficients of the FOD based on the underlying structure. The differences between our suggested method and the baseline are the smallest when dealing with angular sparsity: areas with low crossing angles produce sharper FODs using the baseline method, but areas of high curvature produce more spurious peaks. In general, the FODs produced by the INR appear slightly less sharp but with more intervoxel coherence than the FODs from CSD.

The methods in this paper are based on early INR literature (Mildenhall et al., 2020; Tancik et al., 2020), after which many alternatives have been suggested (Molaei et al., 2023). The positional encoding, model architecture, and training process can all be subject to improvements. An example of a possible improvement would be the sinusoidal networks used by Consagra et al. (2024) and Mancini et al. (2022). We should not only consider applying the existing state-of-the-art but also look into dMRI-specific improvements (e.g., positional encodings based on bundle shapes or dMRI-specific loss functions). Additionally, the model outputs FODs directly, which limits the usability of this approach to scenarios where downstream FOD metrics or fiber tracking is the end goal. In other scenarios, we can revert to the models introduced in previous work that output dMRI signal or use an INR to parameterize different dMRI models (Alexander et al., 2019).

While fiber tracking is the most obvious use case for our proposed models, we have chosen not to include actual fiber tracking experiments in this work. We believe that tracking on a continuous space requires a novel set of tracking methodologies—that goes beyond traditional step-wise streamline methods. The evaluation of the results would also be challenging: many comparisons to existing methods can be made (i.e., different fiber tracking algorithms and different interpolation methods), for which the computational speed and accuracy depend on the implementation. Nonetheless, we have performed a preliminary analysis of fiber tracking using our methods, which is discussed in Appendix B.

We designed the synthetic experiments to separate the different difficulties a model can encounter while fitting a dataset. This is both a strength and a weakness. It allows for quantitative analysis of the model and its comparison in a controlled environment. However, in a real dataset, low spatial and angular resolution and noise artifacts do not present themselves separately, nor are the issues limited to these three. The real data experiments overcome this limitation somewhat but do not allow us to analyze the results quantitatively. Whether the FODs obtained through MSMT-CSD in these experiments are the ground truth is debatable since imperfections in response function and partial volume estimations can lead to inaccuracies. We can see a higher estimation of the GM fraction by the INR compared to MSMT-CSD. This could be an area of improvement for the INR fit if we consider MSMT-CSD to be the ground truth for the volume fractions. Possibly, tuning the hyperparameters of the fitting process or changing the model architecture can generate a similar output to MSMT-CSD. If, however, the true volume fraction, indeed, lies higher, fitting using an INR (in this case) would be an improvement over the existing method. Furthermore, during one of the synthetic experiments, we applied CSD to a noisy dataset, while it would be common to first denoise the data before doing so (Veraart et al., 2016). This is done to showcase the denoising abilities of the INR and make a comparison between methods when fit on ‘raw’ data. Finally, although all three quantitative metrics are commonly used, they do not provide a perfect comparison between FODs. For example, NuFO can vary considerably based on the threshold values for a peak. This is seen in the angular sparsity experiment, where minor differences in FODs (small peaks or slightly different lobe orientations) change the NuFO from 1 to 2 or vice versa.

Everything considered, the usage of CSD model-based INRs for dMRI data is feasible, and it shows great promise for handling spatial sparsity and noise in particular. There are many highly effective upsampling and denoising techniques available (Alexander et al., 2017; Dyrby et al., 2014; Tax et al., 2022), but using INRs has specific advantages compared to these methods. For one, the output of the INR can be tailored to the desired diffusion model output, that is, after fitting once, the inference immediately produces the requested outcome. Future work would be necessary, comparing the quality and speed of computation for existing interpolation and denoising methods and the inference quality and speed of an INR model, which are highly dependent on the encoding and model sizes. Secondly, the ability to sample in a continuous space is a significant advantage. One can imagine probabilistic fiber tracking using our models, benefiting from the computational efficiency and high-quality representation of the INRs. Thirdly, the INR and spatial encoding provide the user with hyperparameters that can be fine-tuned to fit their specific needs. We have analyzed the effects of different hyperparameters in Appendix A, but further research might be necessary to pinpoint exactly how they can influence the results. Some examples could be using lower frequency encodings (by lowering sigma) if we are only interested in the more extensive tracts or coarse structures. Or decreasing the model or encoding sizes to create lower complexity representations in scenarios with poor-quality datasets. Finally, and most interesting, is the possibility, as suggested by Consagra et al. (2024), that the spatial encoding allows the model to learn spatial correlations in the dataset. The results in the spatial resolution, noise, and cdMRI experiment especially point toward the INR having a notion of structural coherence it can exploit to create a better representation. In our model specifically, we see that the INR can capture curved fibers accurately at a resolution much higher than the input. Curvature, fanning, crossing, and kissing fibers are difficult to deal with at a voxel level. This has resulted in the suggestion of using asymmetric FODs, which make use of neighborhood information to more accurately represent the fiber populations in the voxel (Barmpoutis et al., 2008; Bastiani et al., 2017; Poirier & Descoteaux, 2024; Reisert et al., 2012; Wu et al., 2020). Theoretically, a well-fit INR that correctly captures the underlying structure could circumvent the usage of asymmetrical FODs by sampling the regions at a sufficiently high resolution. More interestingly, we could apply INRs to the work done on asymmetrical FODs to further improve the quality of the INR. The improved representation through spatial correlations is not limited to INR-CSD and INR-MSMT, however, but could be used to fit the parameters of any dMRI model. This could increase fit quality, especially when fitting more data-demanding models (e.g., the standard model, Novikov et al., 2019).

This work introduces the INR-CSD and MSMT-CSD models, which integrate the models of CSD and MSMT-CSD with recent research on INRs in dMRI, highlighting both their promising capabilities and limitations. Notably, the denoising abilities and precise on the fly interpolation of FODs offer valuable starting points for further investigation. Ultimately, this research paves the way for applying INRs to other dMRI models, which may benefit from the spatially informed parameter estimation provided by an INR approach compared to a voxel-wise method.

This work was conducted per approval from the Eindhoven University of Technology Ethical Review Board (ERB2023MCS5).

The ISMRM dataset was obtained from: https://tractometer.org/ismrm2015/home/. The CDMRI dataset was obtained from: https://www.cardiff.ac.uk/cardiff-university-brain-research-imaging-centre/research/projects/cross-scanner-and-cross-protocol-diffusion-MRI-data-harmonisation. The code will be made available to the community via GitHub upon publication at: https://github.com/tomhend/MSMT-CSD_INR.

T.H. conceived and implemented the method, performed the experiments, analyzed the results, created the figures, and wrote the first draft of the manuscript. A.V. and M.C. provided guidance in data use and assisted in experimental design and results interpretation. All authors contributed to revisions of the manuscript prior to submission.

This research was funded in part by the Dutch Research Council (NWO) grant number OCENW.M.22.352 awarded to M.C.

The authors declare no conflict of interest.

The CDMRI data were acquired at the UK National Facility for In Vivo MR Imaging of Human Tissue Microstructure located in CUBRIC funded by the EPSRC (grant EP/M029778/1), and The Wolfson Foundation. Acquisition and processing of the data was supported by a Rubicon grant from the NWO (680-50-1527), a Wellcome Trust Investigator Award (096646/Z/11/Z), and a Wellcome Trust Strategic Award (104943/Z/14/Z). This database was initiated by the 2017 and 2018 MICCAI Computational Diffusion MRI committees (Chantal Tax, Francesco Grussu, Enrico Kaden, Lipeng Ning, Jelle Veraart, Elisenda Bonet-Carne, and Farshid Sepehrband) and CUBRIC, Cardiff University (Chantal Tax, Derek Jones, Umesh Rudrapatna, John Evans, Greg Parker, Slawomir Kusmia, Cyril Charron, and David Linden).

We also thank Ruben Vink for his input and interesting discussions.

This appendix provides some preliminary results and quantitative measures concerning fiber tracking for our proposed models. The initial fit of the INR on a considerably high resolution (1.25 mm isotropic) dataset with an average number of gradient directions (30 or 60 non $b=0s​/​mm2$ images) was acquired in 3 to 5.5 minutes for INR-CSD and INR-MSMT respectively. This is considerably faster than previous work using INRs for diffusion MRI (Consagra et al., 2024). After the initial fit, a batch of 100,000 FODs can be produced at arbitrary locations in the volume in 0.8 ms once the data are on the GPU using the Pytorch library. Moving 100,000 coordinates from the CPU to the GPU, running the model for the coordinates, and moving back the FODs from the GPU to the CPU increases this time to 150 ms.

While a complete evaluation of tracking in the continuous space is too elaborate and out of scope for this work, we want to provide insight into the potential qualitative improvements with a simple experiment. We fit the INR-CSD and voxel-wise CSD on the 2.5 mm resolution noiseless synthetic dataset using 30 non $b=0s​/​mm2$ directions. The INR-CSD is then sampled at 1.25 mm resolution, and the CSD FODs are upsampled to 1.25 mm resolution using linear interpolation. Both FODs are then input into the MRtrix implementation of the iFOD2 algorithm to perform full brain fiber tracking, using GM-WM interface seeding regions, stopping when 2 M valid streamlines are reached (Tournier et al., 2010). This implementation uses linear interpolation between the FODs. This implies that any difference in output (excluding the stochasticity inherent to the iFOD2 algorithm) is because of differences in the upsampling step from 2.5 mm to 1.25 mm isotropic resolution. We compare the output to the ground-truth streamlines used to generate the synthetic dataset. The results can be seen in Figure B1. We can see that both INR-CSD and CSD show partial voluming of the fiber bundles into regions that have no bundles in the ground truth. In both regions, this appears to be less profound in the INR-CSD than in the CSD fiber tracking. Possibly, these differences can be more significant when tracking in the continuous space, as we are still relying on linear interpolation for the sub 1.25 mm FODs.

Finally, we use the Tractometer created for the ISMRM 2015 challenge to score both tractograms (Maier-Hein et al., 2017; Renauld et al., 2023). Metrics include the number of valid bundles found (VB, out of the 21 ground truth bundles), the number of valid and invalid streamlines (VS, IS), the overlap (OL, the percentage of voxels covered by the ground truth which are also covered by the tractogram), the normalized overreach (ORn, the percentage of voxels covered outside of the ground truth, normalized by the ground truth size), and the F1-score. An in-depth explanation of these metrics can be found in Renauld et al. (2023). Table B1 shows the INR-CSD and the voxel-wise CSD approach results. CSD finds more valid streamlines but produces the same OL and a higher ORn. This confirms the findings of the qualitative comparison, as these scores indicate that more streamlines are going outside the ground truth bundles. The resulting F1-score is, therefore, higher in INR-CSD compared to CSD.