What is the exact diffusion MRI signal for water confined to a cylindrical surface?

An exact signal without propagator approximations

The question is what diffusion MRI measures when water molecules are restricted to the surface of a cylinder rather than filling its interior. The paper answers with a closed-form expression for the pulsed-gradient spin-echo (PGSE) signal that holds for arbitrary rectangular gradient durations and separations, and it does so without approximating the diffusion propagator or the distribution of spin phases.

The starting point is the Bloch-Torrey equation, which describes how magnetization evolves under diffusion and applied field gradients. Earlier surface-confined work by the same group treated the narrow-pulse limit exactly and extended it approximately to finite pulses; discrepancies appeared at large diffusivities and relatively high diffusion weightings, where finite-pulse effects matter. This paper removes that approximation by solving the equation directly in the eigenbasis of the cylindrical surface.

Why the exponentials cannot be merged

For the standard three-segment PGSE sequence, the signal becomes a product of matrix exponentials: one for the first gradient pulse, one for free evolution between pulses, and one for the reversed second pulse. These matrices do not commute, so they cannot be collapsed into a single exponential. The numerical evaluation must multiply them from right to left, propagating the initial state through the first pulse, then free evolution, then the second pulse.

Only the free-evolution matrix is diagonal, so its exponential is cheap to compute entry by entry. The gradient-pulse exponentials involve dense matrices, and that is where the computational cost concentrates. This structure explains why the paper's later speed-ups target the pulse operators rather than the whole expression.

A real cosine basis halves the matrix size

Truncating the angular eigenbasis at mode order M gives matrices of dimension (2M+1)×(2M+1). The paper shows that spectral coefficients for modes of opposite order stay identical throughout the evolution, c_n(t) = c_-n(t), for any gradient amplitude along a fixed transverse direction. That symmetry lets the +n and -n contributions be combined into a single normalized cosine basis function, leaving an (M+1)×(M+1) problem.

Because standard matrix-exponential algorithms scale cubically with dimension, this reduction implies an ideal speed-up of about eight in the large-M limit. The paper also proves a transformation identity that lets the compact-basis signal be written with the same three-exponential structure, so the reduction does not change the physics, only the size of the matrices. The symmetry would break under circumferential coherent motion or a time-varying transverse gradient direction, neither of which is considered here.

Where the Gaussian phase approximation drifts

The paper compares its exact spherical-mean signal with a recently proposed spherical-mean Gaussian phase approximation (GPA) for the same geometry, using Δ = 11 ms, δ = 3 ms, axial diffusivity 0.8 µm²/ms, and b-values from 0 to 40 ms/µm². At radius 1.0 µm the GPA follows the exact curve closely, with only minor deviations above roughly b = 20 ms/µm². At 2.0 µm and 3.0 µm the GPA departs from the exact solution from about b = 5 ms/µm² onward, and the gap widens with b.

At very high diffusion weightings the exact solution shows oscillatory diffraction-like features produced by angular confinement and interference among the discrete Laplacian eigenmodes. The GPA does not capture these features. The authors note that such large b-values are used to expose the diffraction patterns and are not necessarily reachable in standard clinical protocols, so the comparison is a stress test of the approximation rather than a clinical claim.

Validation across radii and a finite sensitivity window

Analytical signals computed with M = 50 modes were compared with Monte Carlo simulations for 50 discrete radii spanning 0 to 5 µm, at b-values of 1.0 to 6.0 ms/µm² and axial diffusivity 0.8 µm²/ms. The spectral model agreed with the simulations across the full range of radii and b-values, while the GPA reproduced the simulations at small radii and increasingly underestimated the signal at larger radii and higher weightings.

The same comparison defines a practical boundary: for this acquisition protocol the signal shows limited sensitivity to radii below about 0.5 µm and above about 4.0 µm. That window follows from the gradient strength and diffusion times of the PGSE sequence, so it is a property of the protocol, not a universal detection limit of the model.

Few modes and few nodes suffice here

Convergence with truncation order is rapid because the Laplace eigenvalues on the cylindrical surface grow quadratically with mode index, suppressing higher-order contributions. Using the M = 50 Voronoi-weighted spherical mean as reference, M = 10 was numerically indistinguishable from higher orders, and even M = 5 gave accurate approximations under the tested conditions. The GPA was far faster, with an acceleration factor of about 4.8×10³ relative to the reference, but its relative absolute errors exceeded those of even M = 3.

Two further accelerations were tested with Δ = 9.45 ms, δ = 4.61 ms, 92 diffusion-encoding directions, and b = 3, 6, and 20 ms/µm². Strang splitting approximates the pulse exponentials with p substeps and a second-order error that shrinks quadratically with p; diagonalizing the fixed matrix K removes the need for dense matrix exponentials. Gauss-Legendre quadrature replaces the angular average: with only 5 nodes the mean relative absolute error was 4.93×10⁻⁵, 1.37×10⁻⁴, and 6.23×10⁻² at b = 3, 6, and 20 ms/µm², and at 20 nodes the difference from a 100-node reference was on the order of 10⁻¹².

Combining M = 10, p = 10, and 5 quadrature nodes gave errors of 0.17%, 0.34%, and 1.11% at the three b-values, at 0.031 ms per radius and b-value. That is the same order of magnitude as the GPA's 0.023 ms, but with substantially smaller approximation errors. The authors stress that these settings are tied to the tested protocol and should be re-derived when acquisition parameters or radius ranges change.

An approximation-free reference for fitting and averaging

The paper's practical value is as a reference, not a replacement for fast models. Because the expression is exact up to numerical truncation, it can benchmark faster signal approximations and serve applications that need many model evaluations, such as fitting a radius from measured data or averaging over gradient orientations. The factorized structure separates motion along the periodic angular coordinate from unbounded motion along the cylinder axis, which is what makes directional and orientationally averaged signals tractable.

The distinction that matters most is between two kinds of angular averaging. The Voronoi-weighted average over 92 acquired directions estimates the spherical mean from finite sampling, with residual errors of 0.026%, 0.078%, and 0.653% at b = 3, 6, and 20 ms/µm² relative to a converged 100-node quadrature. Gauss-Legendre quadrature instead approximates the continuous, rotationally invariant spherical mean of the model. The gap between them reflects the acquisition scheme's finite angular sampling, not incomplete quadrature convergence.

For readers deciding what to trust, the paper's priority is clear: when restricted diffusion on cylindrical surfaces dominates and finite-pulse effects are significant, the exact spectral signal is the reference, and the GPA's convenience should be weighed against its growing deviation beyond roughly b = 5 ms/µm² at radii of 2 to 3 µm. The accelerated exact implementation narrows the cost gap enough that this trade-off no longer has to be settled by speed alone.

Sources used in this answer