Reconstruction of glymphatic transport fields from subject-specific imaging data, with particular emphasis on cerebrospinal fluid flow and tracer conservation

Abstract: The reconstruction of physically valid transport fields from subject-specific imaging data is a fundamental challenge in image-based computational modeling due to measurement noise, modeling uncertainties and discretization errors. Without a methodology to construct models that faithfully reflect the underlying physics, mechanistic understanding of complex biological systems is inherently limited. In this work, we address this challenge in the glymphatic system, the brain's waste-clearance network, where cerebrospinal fluid (CSF) is transported through perivascular spaces into the brain parenchyma to facilitate metabolic waste removal. We introduce a computational framework for the high-fidelity reconstruction of subject-specific glymphatic transport fields from spatiotemporal imaging data. The formulation utilizes an advection-diffusion model with a velocity decomposition that imposes mass conservation, enabling the recovery of solenoidal (divergence-free) velocity fields through the solution of a constrained inverse problem. The system is discretized using immersed isogeometric analysis with quadratic B-spline basis functions, providing smooth, high-continuity solutions and inherent regularization of imaging noise. We demonstrate the framework's utility by using contrast-enhanced magnetic resonance imaging of tracer transport in a mouse brain, obtaining spatially varying estimates of CSF velocity, diffusivity, and clearance parameters. Forward simulations using the recovered fields show close agreement with experimental observations, validating the framework's ability to characterize complex transport dynamics while preserving physical integrity. This approach provides a generalizable methodology for the robust inference of physically consistent transport fields from imperfect imaging data, with broad applicability to the image-guided modeling of biological and engineering systems.

• The paper presents a novel HW-modified advection–diffusion model that enforces mass conservation in reconstructing glymphatic transport from noisy MRI data.
• It employs an immersed finite element/isogeometric analysis framework with THB-splines to achieve robust, subject-specific reconstructions of CSF velocity, diffusion, and clearance fields.
• The approach accurately recovers divergence-free velocity maps and spatially varying transport parameters, offering new insights for studies on brain waste clearance and neurodegenerative diseases.

Physically-Constrained Inference of Glymphatic Transport from Imaging Data

Introduction and Motivation

Quantitative modeling of solute transport in the brain parenchyma is essential for elucidating the mechanisms underpinning brain waste clearance, particularly the glymphatic system—a subject of active controversy and deep clinical relevance for neurodegenerative diseases [mestre_brains_2020, bohr_glymphatic_2022]. Existing experimental and computational approaches have struggled with three technical bottlenecks: (1) obtaining anatomically faithful, subject-specific transport data; (2) inferring physically consistent (divergence-free) flow fields from inherently noisy and temporally/ spatially sparse MRI data; and (3) capturing the heterogeneous, multiscale nature of brain clearance. This work systematically addresses these challenges via a variational inverse modeling framework that enforces physical constraints on flow and transport fields extracted from mouse brain MRI (CE-MRI) data, with a specific emphasis on cerebrospinal fluid (CSF) velocity field reconstruction and conservation of mass.

The glymphatic system routes CSF from the choroid plexus through periarterial spaces, exchanges it with interstitial fluid (ISF) via aquaporin-4–rich astrocytic endfeet, and ultimately expels metabolic waste via perivenous and meningeal lymphatic pathways (Figure 1).

Figure 1: Schematic depicting CSF/ISF interchange, waste mobilization, and lymphatic drainage in the glymphatic system as modeled in this study.

Quantitative tools grounded in imaging data are required to advance both mechanistic understanding and subject-specific predictions for brain clearance, particularly in disease settings.

Model Formulation: Conservation with Non-Solenoidal Advection

The study's core contribution is an advection–diffusion model tailored for physical consistency, even when CSF velocity fields inferred from imaging are not inherently divergence-free. Standard formulations in which image-derived velocities directly parameterize advection operators lead to nonphysical mass loss/gain due to measurement noise and discretization error (Figure 2). To address this, the authors adapt the Hughes–Wells (HW) velocity decomposition [hughes_conservation_2005]. For any image-derived (potentially non-solenoidal) velocity $\bar{u}$, the true field is decomposed as $\bar{u} = \hat{u} + u'$, where $u'$ is constructed as the gradient of a scalar potential chosen to ensure weak solenoidality.

The numerical implementation leverages immersed finite element/isogeometric analysis, with quadratic B-spline basis functions on truncated hierarchical B-splines (THB-splines). This provides regularization, derivative continuity, and mesh flexibility over complex brain geometries (Figure 3).

Figure 3: Mouse brain geometry represented with immersed IGA mesh and THB-splines enabling flexible, noise-regularized field reconstructions.

The advection–diffusion–clearance equations are solved with Robin (flux-proportional) boundary conditions, capturing spatially variable clearance that reflects anatomical lymphatic heterogeneity. The measurable fields (velocity, diffusion, clearance) are themselves optimized in the inverse problem.

Experimental Image Data and Inverse Problem

The study utilizes dynamic CE-MRI data in mice following CSF tracer injection (Figure 4). Notably, signal changes at the voxel level are marginal relative to noise, severely constraining any velocity field estimates derived from finite differences (Figures 4, 5). The raw finite-difference approach produces checkerboard noise and spurious divergence.

Figure 4: Temporally resolved tracer intensity data, with window selection for calibration versus forward validation.

Figure 5: CE-MRI slices and voxel-level differencing, highlighting low SNR and the limitations for direct velocity estimation.

Figure 2: Raw finite-difference velocity field (left) and HW-modified (right); original is highly noisy and non-divergence-free, corrected field is smooth and solenoidal.

By embedding the HW-modified advection–diffusion equation inside a Newton-based variational inverse problem, the framework simultaneously infers the spatial distributions of advection, diffusion, and boundary clearance coefficients to optimally reconcile model predictions with observed spatiotemporal MRI tracer data, while maintaining physical fidelity and mass conservation throughout.

The optimization variables are the spline coefficients for the field expansions. The parameter update at each Newton step uses linearizations of the model output with respect to all control variables (velocity field, diffusion, clearance).

The pseudocode for this inversion workflow is illustrated in Figure 6.

Figure 6: Inverse problem algorithm schematic for physics-constrained transport field reconstruction from MRI data.

Numerical Results and Comparative Analysis

Impact of Conservation-Preserving Reformulation

Forward simulations with various formulations (conservative, advective, HW-modified) using the noisy, non-solenoidal finite difference velocity recover the expected result: standard advective approaches exhibit nonconservation and spurious solute loss/gain, while the HW-modified scheme preserves global mass conservation and yields physically plausible (nonblow-up, nonnegative) solutions even on long time scales (Figures 7, 8).

Figure 7: Conservation error and min/max concentration evolution for three numerical schemes; only HW-modified maintains stability and conservation.

Figure 8: Simulated and experimental tracer distributions at 90 min; HW-modified achieves the lowest relative error.

Inverse Reconstructions and Sensitivity Analysis

Full inverse solutions initialized from literature-based parameter guesses produce subject-specific velocity, diffusion, and clearance maps (Figure 9). The reconstructed velocity fields exhibit peak CSF velocity in the inferior brain and divergence-free global structure. Spatial patterns of diffusivity and clearance closely reflect anatomical expectations (high clearance in the basal brainstem/neck, negligible superiorly). The residuals of the fit are insensitive to order-of-magnitude perturbations in initial guesses (Figure 10).

Figure 9: Convergence of residual during inversion (A), and reconstructed velocity, diffusion, and clearance fields (B); velocity divergence is minimized.

Figure 10: Inverse problem residuals and converged diffusivity with various initializations, demonstrating robustness to initial field guesses.

The reconstructed fields achieve strong qualitative and quantitative agreement with independent (non-calibration) MRI timepoints (Figures 11, 12).

Figure 11: Comparison of experimental and forward-simulated tracer distributions using converged inverse fields.

Figure 12: Time-resolved error between experiment and simulation; maximum error of ~0.1 across all timepoints.

Calibration with all available timepoints, as opposed to a five-point subset, reduces errors further at late times and recovers sharper, less diffusive solute gradients (Figures 13, 14, 15). However, increased calibration duration comes at added computational cost and may reflect perturbed physics during early post-injection transients.

Figure 13: Residual decrease and converged fields for "full-window" calibration; advective effects more localized to expected drainage anatomy.

Figure 14: Relative error across time for five-point versus all-timepoint calibrations, with image comparisons.

Figure 15: Detailed spatial comparison for final timepoint; all-window calibration recovers sharper features.

Localized Clearance Parameterization

A key theoretical contribution is the distinction between the intrinsic boundary clearance parameter $\gamma(\mathbf{x})$—a mechanistic, state-independent representation of periphery permeability—and emergent (state-dependent) global clearance rates $k_\mathrm{eff}(t)$, the latter of which varies with tracer spatial distribution and cannot directly capture anatomical drainage heterogeneity (Figure 16).

Figure 16: Time-dependent global clearance rate $k_\mathrm{eff}(t)$ (left) versus inferred, anatomy-specific $\gamma(\mathbf{x})$ (right, superior and inferior projection), highlighting inferior "hotspots."

Theoretical and Practical Implications

This study provides the first advection–diffusion inference framework that rigorously enforces conservation with noisy, non-solenoidal velocity fields extracted from imaging, while simultaneously reconstructing spatially varying diffusion and boundary clearance in realistic brain anatomy. The framework advances beyond previous optimal mass transport approaches by (i) endowing velocity fields with correct physical units and conservative properties, (ii) enabling subject-by-subject interpretation, and (iii) producing interpretable field maps for each transport parameter.

Numerical results indicate that diffusion and advection both contribute to glymphatic transport, but that advection—localized near basal efflux routes—plays a critical, nonuniform role. The field regularization provided by B-splines and the immersed isogeometric framework ensures robust convergence and mitigates overfitting to noise-limited data.

Future Directions

Extension of this approach to human imaging data is an immediate practical target, given the surge in high-resolution MRI studies quantifying glymphatic transport in health and disease [vinje_human_2023-1, bohr_glymphatic_2022]. The methodology is compatible with future inclusion of time-varying parameter fields, explicit coupling to amyloid deposition/damage models, and integration with alternative imaging modalities (e.g., PET). Generalization to other biological transport networks—where incompressible flows must be reconstructed from experimental data—remains an open avenue.

Methodological synergies with neural Green’s operator learning [melchers_neural_2025] or divergence-free kernel networks [ni_representing_2025] could further enhance scalability and generalizability in large-scale inference and surrogate modeling.

Conclusion

The additive, conservation-enforcing velocity decomposition coupled with a variational inverse modeling framework offers a powerful, interpretable, and robust foundation for subject-specific inference of glymphatic and related transport networks in the brain (2605.00730). This approach bridges the gap between noisy, indirect imaging observations and mechanistically grounded, field-resolved transport models—opening new pathways for understanding and potentially modulating waste clearance in normal and diseased states.