PRISM: Differentiable Analysis-by-Synthesis for Fixel Recovery in Diffusion MRI

Abstract: Diffusion MRI microstructure fitting is nonconvex and often performed voxelwise, which limits fiber peak recovery in narrow crossings. This work introduces PRISM, a differentiable analysis-by-synthesis framework that fits an explicit multi-compartment forward model end-to-end over spatial patches. The model combines cerebrospinal fluid (CSF), gray matter, up to K white-matter fiber compartments (stick-and-zeppelin), and a restricted compartment, with explicit fiber directions and soft model selection via repulsion and sparsity priors. PRISM supports a fast MSE objective and a Rician negative log-likelihood (NLL) that jointly learns sigma without oracle information. A lightweight nuisance calibration module (smooth bias field and per-measurement scale/offset) is included for robustness and regularized to identity in clean-data tests. On synthetic crossing-fiber data (SNR=30; five methods, 16 crossing angles), PRISM achieves 3.5 degrees best-match angular error with 95% recall, which is 1.9x lower than the best baseline (MSMT-CSD, 6.8 degrees, 83% recall); in NLL mode with learned sigma, error drops to 2.3 degrees with 99% recall, resolving crossings down to 20 degrees. On the DiSCo1 phantom (NLL mode), PRISM improves connectivity correlation over CSD baselines at all four tracking angles (best r=.934 at 25 degrees vs. .920 for MSMT-CSD). Whole-brain HCP fitting (~741k voxels, MSE mode) completes in ~12 min on a single GPU with near-identical results across random seeds.

• The paper introduces PRISM, an end-to-end differentiable framework that jointly fits multi-compartment tissue models, fiber orientations, spatial regularization, and intensity calibration across voxel patches.
• PRISM recovers crossing fibers with 3.5° angular error and 95% recall in MSE mode, improving to 2.3° error and 99% recall with learned-noise Rician likelihood fitting.
• The method improves phantom tractography connectivity to r=.934 and fits a 741,000-voxel human brain in about 12 minutes on one GPU, while remaining limited by fixed diffusivities, missing distortion modeling, and limited in-vivo validation.

Overview

PRISM is a differentiable analysis-by-synthesis framework for diffusion MRI (dMRI) microstructure fitting that replaces the conventional pipeline of voxelwise, per-stage estimation with a single end-to-end optimization over spatial patches. The method fits an explicit multi-compartment forward model — combining CSF, gray matter, up to $K$ white-matter fiber compartments (stick-and-zeppelin), and a restricted isotropic compartment — jointly with a nuisance calibration module comprising a smooth bias field, per-measurement affine intensity correction, and a learnable noise level $\sigma$. Two data-fidelity modes are supported: a fast MSE objective and a Rician negative log-likelihood (NLL) that recovers $\sigma$ self-supervisedly without oracle information. The work is positioned against scalar models (DTI, DKI), orientation estimators such as CSD and MSMT-CSD, deep-learning regressors, and forward-simulation/dictionary methods (ODF-FP, FORCE); its distinguishing regime is the multi-compartment inverse problem on clinical-style acquisitions with neighbor-aware regularization and joint calibration (2604.00250).

Forward model and parameterization

The synthesized measurement at voxel $\mathbf{x}$ and measurement $n$ is $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$, where the tissue signal is a softmax-weighted mixture of mono-exponential isotropic compartments ($D_{\mathrm{csf}}=3.0$, $D_{\mathrm{gm}}=0.9$, $$D_{\mathrm{res}}=0.2 \times 10^{-3}\,\mathrm{mm^2/s}$$) and $K$ stick-and-zeppelin fibers with fixed $\sigma$0, $\sigma$1, and learnable intra-axonal fraction. The restricted compartment captures slowly decaying high-$\sigma$2 signal from cellular structures in the spirit of SANDI-type models.

The calibration cascade applies an exponentiated trilinearly upsampled $\sigma$3 control-grid bias field, a per-measurement affine correction $\sigma$4, and a learned noise scale $\sigma$5. Because the product $\sigma$6 is scale-ambiguous, well-posedness rests on the frequency separation between the low-resolution bias grid and per-voxel $\sigma$7, plus $\sigma$8 penalties anchoring all calibration parameters to identity. This is an explicit assumption: on data where nuisance structure is not low-frequency or affine-per-volume, the decomposition may be ill-posed.

Regularization includes a Huber-Laplacian spatial smoothness prior over tissue fractions (6- or 26-connected; ablation shows 26-connectivity is not consistently better), direction repulsion $\sigma$9, $\sigma$0 sparsity on minor fiber fractions, orphan-WM suppression, directional continuity, and fiber ordering to break permutation symmetry. Optimization uses Rprop, selected empirically among eight optimizers for its sign-based updates that handle heterogeneous parameter scales; it attains the lowest final MSE ($\sigma$1) and fastest convergence (22 iterations to threshold) versus Adam-family baselines on Stanford HARDI slices. Convergence requires 50–300 iterations depending on dataset.

Crossing-fiber recovery

On the synthetic 16-angle benchmark (SNR = 30, 3 shells, 3400 voxels), PRISM achieves $\sigma$2 overall best-match angular error with 95% recall in MSE mode — a $\sigma$3 reduction relative to the strongest baseline, MSMT-CSD ($\sigma$4, 83% recall). In NLL mode with learned $\sigma$5, error falls to $\sigma$6 with 99% recall and F1, resolving crossings down to $\sigma$7. Notably, CSD and MSMT-CSD largely detect only one fiber at crossings $\sigma$8, while ODF-FP's dictionary discretization degrades at wider angles ($\sigma$9). Two aspects of this comparison deserve emphasis: PRISM uses a deliberately mismatched forward model ($\mathbf{x}$0 vs. true $\mathbf{x}$1), whereas CSD derives its response from ground-truth single-fiber voxels and MSMT-CSD receives oracle response functions — so the reported advantage is obtained under a handicap in model specification. However, the benchmark is synthetic with known ground truth, and results at SNR 10 and 50 are mentioned but not detailed in the provided content.

Under per-measurement gain perturbations ($\mathbf{x}$2–$\mathbf{x}$3), the calibration module halves angular error at $\mathbf{x}$4 ($\mathbf{x}$5) and removes 85% of reconstruction MSE, while self-regularizing to identity on clean data — confirming the module activates only under corruption rather than absorbing biological signal.

Tractography and in-vivo validation

On the DiSCo1 digital phantom (16 ROIs, 120 ROI pairs, SNR = 50), identical deterministic tracking shows PRISM (NLL mode) outperforming both CSD baselines at every tracking angle, reaching best connectivity correlation $\mathbf{x}$6 at $\mathbf{x}$7 versus $\mathbf{x}$8 for MSMT-CSD at $\mathbf{x}$9, with the largest gains ($n$0 pp) at tighter angles ($n$1–$n$2). The shift of PRISM's optimum to a narrower tracking angle indicates its recovered orientations support more stringent angular thresholds.

Ablation on DiSCo1 identifies the restricted compartment as the dominant contributor: removing it drops $n$3 from $n$4 to $n$5 ($n$6 pp), compared with $n$7 pp for the topology prior and $n$8 pp for the spatial prior. This is a strong claim about model composition — high-$n$9 restricted diffusion, not fiber-count priors, drives most of the tractography gain on this phantom.

In-vivo whole-brain fitting on HCP subject 100307 (~741k voxels, 288 measurements) completes in ~12 minutes on a single RTX A6000 in MSE mode, reaching MSE $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$0. Tissue fractions agree with FreeSurfer $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$1 segmentations at Dice 0.79 (WM) and 0.80 (GM) with Pearson $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$2, without any $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$3 input to the fit, and all summary metrics vary by less than 0.003 across three random seeds. These agreement figures are moderate rather than near-ceiling, which is expected given that dMRI-derived fractions and histological-style segmentations measure related but distinct quantities.

Limitations and open questions

The paper concedes several boundaries explicitly. PRISM models only intensity-domain nuisance factors; geometric artifacts (motion, susceptibility distortion) remain the responsibility of upstream preprocessing. Convergence relies on random initialization within a nonconvex landscape, and the authors identify warm-starting from dictionary methods (ODF-FP, FORCE) or a self-supervised encoder as untested accelerations. Fiber dispersion is not modeled (no Bingham distribution), and fixed diffusivities across all compartments impose biophysical rigidity whose sensitivity is quantified only through the single deliberate mismatch experiment. Whether the DiSCo1 ablation ranking — restricted compartment dominating — transfers to in-vivo data at clinical SNR remains open, as does scaling behavior beyond $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$4–$$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$5 fibers on real acquisitions.

Conclusion

PRISM demonstrates that treating microstructure estimation as differentiable analysis-by-synthesis — with explicit fiber parameters, a restricted compartment, learned-$$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$6 Rician fitting, and joint nuisance calibration — yields measurable improvements in crossing-fiber resolution ($$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$7 vs. $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$8 angular error; $$\hat{y}(\mathbf{x},n) = \mathcal{A}(\mathcal{S}(\mathbf{x},n;\boldsymbol{\theta}_{\mathrm{micro}});\boldsymbol{\theta}_{\mathrm{cal}})$$9 in NLL mode), phantom tractography agreement ($D_{\mathrm{csf}}=3.0$0), and anatomically consistent whole-brain maps at clinically practical compute cost (~12 min per brain on one GPU). The evidence base is strongest on synthetic and phantom data with ground truth; the in-vivo validation is limited to segmentation agreement and reproducibility rather than downstream task performance, leaving quantitative in-vivo fiber-recovery accuracy as the principal open question.