The Diffusion Exchange Ratio (DEXR): A minimal sampling of diffusion exchange spectroscopy to probe exchange, restriction, and time-dependence

Abstract

Water exchange is increasingly recognized as an important biological process that can affect the study of biological tissue using diffusion MR. Methods to measure exchange, however, remain immature as opposed to those used to characterize restriction, with no consensus on the optimal pulse sequence(s) or signal model(s). In general, the trend has been towards data-intensive fitting of highly parameterized models. We take the opposite approach and show that a judicious sub-sample of diffusion exchange spectroscopy (DEXSY) data can be used to robustly quantify exchange, as well as restriction, in a data-efficient manner. This sampling produces a ratio of two points per mixing time: (i) one point with equal diffusion weighting in both encoding periods, which gives maximal exchange contrast, and (ii) one point with the same total diffusion weighting in just the first encoding period, for normalization. We call this quotient the Diffusion EXchange Ratio (DEXR). Furthermore, we show that it can be used to probe time-dependent diffusion by estimating the velocity autocorrelation function (VACF) over intermediate to long times (~ 2-500 ms). We provide a comprehensive theoretical framework for the design of DEXR experiments in the case of static or constant gradients. Data from Monte Carlo simulations and experiments acquired in fixed and viable ex vivo neonatal mouse spinal cord using a permanent magnet system are presented to test and validate this approach. In viable spinal cord, we report the following apparent parameters from just 6 data points: ${\tau }_k=17\pm 4ms$ , $f_{NG}=0.71\pm 0.01$ , $R_{eff}=1.10\pm 0.01\mu m$ , and ${\kappa }_{\text{eff}}=0.21\pm 0.06\mu m/ms$ , which correspond to the exchange time, restricted or non-Gaussian signal fraction, an effective spherical radius, and permeability, respectively. For the VACF, we report a long-time, power-law scaling with $\approx t^{-2.4}$ , which is approximately consistent with disordered domains in 3-D. Overall, the DEXR method is shown to be highly efficient, capable of providing valuable quantitative diffusion metrics using minimal MR data.

Keywords: diffusion exchange spectroscopy (DEXSY); exchange; low-field; static gradient; time-dependent diffusion; velocity autocorrelation.

Figure 1:

SG-DEXSY pulse sequence with timing parameters ${\tau }_1$, ${\tau }_2$, and $t_m$. Signal is acquired in a CPMG loop with ${\tau }_{\text{CPMG}}$. The effective gradient $G_{\text{eff}}(t)\in \{0,-\gamma g,+\gamma g\}$ and its modulation by RF pulses is shown below.

Figure 2:

Comparison of the SG-SE signal behavior in different regimes and in simulation and fixed spinal cord data w.r.t. ${\rho }^2={\left({\ell }_d/{\ell }_g\right)}^2\propto b^{1/3}$. (a) Signal curves $S/S_0$ plotted on a log-axis for the free (dash-dot line), motionally-averaged (red to black lines), and localized (blue) regimes compared to data from Monte Carlo simulations (magenta) and data acquired in fixed spinal cord (cyan). Curves are plotted up to ${\rho }^2=25$ using $D_0=2.15\mu {\mathrm{m}}^2/\mathrm{m}\mathrm{s}$, $g=15.3\mathrm{T}/\mathrm{m}$, which gives a dephasing length ${\ell }_g\approx 806\mathrm{n}\mathrm{m}$. Motionally-averaged signal is plotted for spherical radii from $R=0.4-1\mu \mathrm{m}$ or from ${\ell }_s=2R\approx 1-2.5{\ell }_g$, summing up to the first 5 roots in Eq. (2). Note the rapid approach of Eq. (2) towards the asymptotic, linear behavior predicted by Eq. (4) as ${\rho }^2$ exceeds ≈ 2, or $\rho \gtrsim 1.4$. Localized signal is plotted only for ${\ell }_d>1.5{\ell }_g$ using the prefactor $a_0=5.8441$, see Eq. (5). For the Monte Carlo simulation data, error bars indicate ±1 SD from 3 repetitions with different random seeds. For the spinal cord data, error bars indicate 95% confidence intervals estimated by bootstrapping 43 repetitions on the same sample. Fits to Eq. (7) yield $f_I\approx 0.44$, $⟨c_E⟩\approx 0.40$, $⟨c_I⟩\approx 0.21$ for the simulation data and $f_I\approx 0.16$, $⟨c_E⟩\approx 0.26$, $⟨c_I⟩\approx 0.18$ for the spinal cord data. (b) Zoomed plot up to ${\rho }^2=5$, highlighting the transition from Gaussian signal behavior to the characteristic non-Gaussian signal decay that is linear on this axis of ${\rho }^2\propto b^{1/3}$. Note the deviation from the fit in the spinal cord data, suggestive of potentially distributed non-Gaussian compartments.

Figure 3:

Signal difference $\mathrm{\Delta }S/S_0$ between non-Gaussian and Gaussian signal decay plotted as a function of $\rho ={\ell }_d/{\ell }_g$ for motional averaging in spheres of radii $R=0.4-1\mu \mathrm{m}$ using Eq. (2), compared to free diffusion, or Eq. (1). The difference between the components of the fits to Eq. (7) for simulated (magenta) and spinal cord (cyan) data shown in Fig. 2 are also plotted, i.e., $\mathrm{\Delta }S/S_0$ is calculated as the difference between the terms $\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\rho }^2⟨c_{NG}⟩\right)$ and $\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\rho }^6⟨c_G⟩\right)$. For these data, $\mathrm{\Delta }S/S_0$ is smaller and the maximum is farther to the right due to hindered diffusion. Overall, the optimal range to maximize $\mathrm{\Delta }S/S_0$ is about $\rho \approx 1.35-1.55$, with $\rho \gtrsim 1.4$ as a heuristic.

Figure 4:

Signal contour and difference maps between low and high exchange cases. (a) Plots for synthetic data generated using Eq. (13) and the fit parameters obtained by fitting Eq. (7) to the SG-SE spinal cord data in Fig. 2: $f_{NG}\approx 0.16$, $⟨c_G⟩\approx 0.25$, $⟨c_{NG}⟩\approx 0.18$. Exchanging signal fractions $f_{\text{exch}}=[0.02,0.13,0.27]$ are compared, where $f_{\text{exch},\text{ss}}\approx 0.27$. Exchange is seen to produce an inwards curvature in the signal contours around ${\rho }_1,{\rho }_2\gtrsim 1.25$ (see middle panel). The difference map $\mathrm{\Delta }S/S_0$ indicates that ${\rho }_1={\rho }_2\approx 1.55$ produces the most exchange contrast, which agrees with the optimum and range identified in Fig. 3. The parity axis ${\rho }_1={\rho }_2$ is marked with a dash-dot line. The heuristic optimum of ${\rho }_1={\rho }_2=1.4$ is marked by a cross. (b) The same plots as part (a) but using the maximal $f_{\text{exch}}=f_{\text{exch},\text{ss}}\approx 0.27$. Note that the exchange contrast $\mathrm{\Delta }S/S_0$ roughly doubles as $f_{\text{exch}}$ doubles between (a) and (b), proportional with the increase in $f_{\text{exch}}$. The peak value of $\mathrm{\Delta }S/S_0$ increases from ≈ 0.023 to ≈ 0.046 (see color bar values). (c) The same plots for data acquired in fixed spinal cord in a 6 × 6 grid at ${\rho }_1\approx {\rho }_2=[1.09,1.30,1.49,1.65,1.80,1.93]$ (recall that ${\tau }_2\ne {\tau }_1$, see Methods). In these data, the optimal point is shifted towards a smaller ${\rho }_1={\rho }_2\approx 1.37$ than in parts (a) or (b), and is also of a smaller peak amplitude than part (b), with a maximal $\mathrm{\Delta }S/S_0\approx 0.037$. This may be due to the model being an incomplete description of the distributed non-Gaussian microenvironments in tissue (see again the fit deviations in Fig. 2b) and/or larger compartments with a smaller expected optimum but larger volume dominating the exchange contrast (see Fig. 3). Nonetheless, the heuristic ${\rho }_1={\rho }_2=1.4$ remains a good choice. Despite the coarse sampling of this data and the shift in optimum, the qualitative similarity in shape and character to part (b) is evident.

Figure 5:

Comparison of $f_{\text{exch}}$ and $S_{\text{mid}}/S_0$ obtained from simulation data with ${\rho }_1={\rho }_2\approx 1.397$ and $t_m=[0.1,0.5,1,2,5,10,15,20,50]\mathrm{m}\mathrm{s}$. Here, the ground-truth $f_{\text{exch}}$ is quantified as the fraction of walkers that moved from inside/outside of a sphere between the start of the simulation and the start of the second diffusion encoding. The curve for $f_{\text{exch}}$ (left axis, black) is fit to the form $f\left(t_m\right)={\beta }_1\left[1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\beta }_2{t}_m\right)\right]+{\beta }_3$ whereas $S_{\text{mid}}/S_0$ (right axis, magenta) is fit to $f\left(t_m\right)={\beta }_1\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\beta }_2{t}_m\right)+{\beta }_3$. Note that a small intercept of ${\beta }_3\approx 0.02$ is estimated in $f_{\text{exch}}$ due to exchange during the encoding period. Error bars indicate mean ± SD from 3 repetitions. Solid lines are a fit to the mean. Fits to each repetition yield ${\tau }_k=1/{\beta }_2=21.7\pm 0.9$ and 20 ± 5 for $f_{\text{exch}}$ and $S_{\text{mid}}/S_0$, respectively. The values are in agreement, though noisier for $S_{\text{mid}}$.

Figure 6:

Estimation of the exchange time $t_k$ from two points per $t_m$ ($S_{\text{mid}}$ and $S_{\text{end}}$) in viable ex vivo spinal cord. The normalization and fitting approach shown here comprise the DEXR method. (a) Plots of the raw signal decay of $S_{\text{mid}}$ with ${\rho }_{\text{mid}}\approx 1.39(\tau \approx 0.59\mathrm{m}\mathrm{s})$ and $S_{\text{end}}$ with ${\rho }_{\text{end}}\approx 1.55(\tau \approx 0.74\mathrm{m}\mathrm{s})$ for 11 values of $t_m=[0.2,1,2,4,7,10,20,40,80,160,300]\mathrm{m}\mathrm{s}$. Data is plotted on a log y-axis to highlight the approximately linear decay at long $t_m$, indicative of diffusion-weighted $T_1$ relaxation. Error bars indicate mean ± SD from 3 repetitions on the same sample. Both points evolve by $T_1$, while $S_{\text{mid}}$ also evolves due to exchange. Fitting a monoexponential decay to $S_{\text{end}}$ yields an apparent diffusion-weighted $T_1\approx 600\pm 20\mathrm{m}\mathrm{s}$, which differs from the $T_1\approx 710\pm 10\mathrm{m}\mathrm{s}$ obtained by fitting an $S_0$ acquisition with ${\tau }_1={\tau }_2=0.05\mathrm{m}\mathrm{s}$ ms (fits and data not shown, see ref. [43]), highlighting the non-triviality of accounting for $T_1$. Note that because $S_{\text{end}}$ is not acquired precisely at ${\tau }_2=0$, but at ${\rho }_2\approx 0.81$, this point is also slightly exchange-weighted - see the non-linear behavior at short times. (b) Fit of Eq. (18) to the ratio $S_{\mathrm{m}\mathrm{i}\mathrm{d}}/S_{\text{end}}$, yielding ${\tau }_k=11\pm 3\mathrm{m}\mathrm{s}$. Fits to the mean using all $11t_m$ (solid line) or a minimal 3 values of $t_m=[0.2,20,160]\mathrm{m}\mathrm{s}$ (crosses, dashed line) are plotted. The minimal sampling yields a similar ${\tau }_k=17\pm 4\mathrm{m}\mathrm{s}$.

Figure 7:

Relationship between $f_{NG}$ and the fit-derived quantity $1-{\beta }_1/\left({\beta }_1+{\beta }_3\right)$ for various values of $\sigma =\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\rho }_{\text{mid}}^6⟨c_G⟩\right)/\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\rho }_{\text{mid}}^2⟨c_{NG}⟩\right)$ and $\varsigma =f_{\text{exch},0}/f_{\text{exch,ss}}$, which characterize the confounding effects of extant exchanged signal and exchange during the encoding, respectively. The solid black line indicates parity when $\sigma ,\varsigma =0$. Curves of Eq. (24) derived from an idealized $S_{\text{mid}}/S_0$ are plotted for $\sigma =[0.02,0.08,0.14,0.2]$ and $\varsigma =[0.02,0.06,0.1,0.14]$. The parameters are varied together (magenta circles) and independently (blue crosses, red diamonds), with deepening color representing increasing values. In all cases, the behavior manifests, roughly speaking, as a decrease in the linear relationship or slope between $f_{NG}$ and ${\beta }_1/\left({\beta }_1+{\beta }_3\right)$.

Figure 8:

Time-domain weighting $C(t)$ for an $S_{\text{end}}$ acquisition with ${\tau }_1=0.74\mathrm{m}\mathrm{s}$ compared to an $S_{\text{mid}}$ acquisition with ${\tau }_1={\tau }_2=0.59\mathrm{m}\mathrm{s}$; $t_m=2.36\mathrm{m}\mathrm{s}$ for both. These parameters are approximately consistent with the curvature method and $b_s=4.5\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^2$. The curves are plotted on a non-dimensionalized y-axis of $𝒞(t)/\left({\gamma }^2{g}^2\right)$, where $g=15.3\mathrm{T}/\mathrm{m}$. Gray vertical lines indicate the timing of RF pulses and echo formation for the $S_{\text{mid}}$ acquisition. The weightings over $0\le t<2\tau$ are similar between the two acquisitions, as expected for an equal total $b$-value or $b_s$. However, $C_{\text{mid}}(t)$ has an additional two lobes centered about $t=2\tau +t_m$ with peaks at $t=(4/3)\tau +t_m$and $(8/3)\tau +t_m$. These peaks arise from the autocorrelation of the separated diffusion encodings and integrate to $\mp b_s/2$, respectively, while the first lobe integrates to $+b_s$ (i.e., if the y-axis were multiplied by ${\gamma }^2{g}^2$).

Figure 9:

Phase distributions $P\left(\varphi \right)\in [-\pi ,+\pi ]$ at (a) the first echo for $f_{NG}$ and $f_G$, and (b) at the final echo for $f_{NG,NG}$, $f_{\text{exch}}$, and $f_{G,G}$ in simulation data generated from a single repetition using ${\rho }_{\text{mid}}\approx 1.397$ and $t_m=50\mathrm{m}\mathrm{s}$. The signal $S=⟨\mathrm{c}\mathrm{o}\mathrm{s}\left(\varphi \right)⟩$ (i.e., real part) is shown as text. The non-exchanging, restricted signal $f_{NG,NG}$ is well-described by a Gaussian (green, dash-dot line) such that the GPA holds, as expected for this value of $\rho \gtrsim 1$. Quantitatively, an Anderson-Darling test yields a p-value ≈ 0.33. The Gaussian signal $f_{G,G}$ is fully dephased. The exchanged signal $f_{\text{exch}}$ is mostly, though not entirely dephased.

Figure 10:

MSD and other transport quantities from simulation. (a) The MSD along the gradient direction from one simulation up to $t\approx 52.36\mathrm{m}\mathrm{s}$ (blue dots), corresponding to $t_m=50$, $\tau =0.59\mathrm{m}\mathrm{s}$. Also shown is a linear relationship with $t$ (dotted) and an exponential fit (dashed) that describe the MSD at shorter and longer times, respectively. To analyze the first and second derivatives of the MSD, a piecewise, cubic Hermite polynomial was fit in the log-log domain (solid line), splitting the domain into 10 log-linearly spaced segements. For all panels, the bottom plot shows the short-time behavior $t\le 1\mathrm{m}\mathrm{s}$. Note the immediate deviation from $2D_{0}t$. (b) Same data and relationships on a log-log plot. The approach towards exponential behavior is clear. (c) Diffusivity quantities derived from the MSD. The time-dependent diffusivity $D(t)$ is the raw MSD divided by $2t$. The instantaneous diffusivity $D_{\text{inst}}(t)$ is estimated using a backward, first-order finite difference of the piecewise fit to the MSD with spacing $\mathrm{\Delta }t=5\times {10}^{-4}\mathrm{m}\mathrm{s}$, or twice the simulation time-step. Both quantities decay monotonically from $D_0$. $D(t)$ approaches a Gaussian limit described by an unknown $D_{\mathrm{\infty }}$ where as $D_{\text{inst}}(t)$ approaches 0, consistent with the bounding box in the simulation. (d) The VACF is estimated as half the curvature in the piecewise fit to the MSD obtained using a central, second-order finite difference with the same spacing $\mathrm{\Delta }t=5\times {10}^{-4}\mathrm{m}\mathrm{s}$. The VACF exhibits a sharp, initial decrease due to reflection before asymptotically approaching 0 as $t\to \mathrm{\infty }$ and the system loses its “memory” of the first interaction(s) with barriers via exchange.

Figure 11:

Apparent VACF from DEXR data viable spinal cord, calculated using Eqs. (45) and (18) to obtain ${\beta }_3$. Error bars = mean ± SD from repetitions with different seeds. For each repetition, ${\beta }_3$ was first estimated by fitting to Eq. (18) and the VACF was then calculated as $3/\left(2\tau b_s\right)\left[\mathrm{l}\mathrm{n}\left({\beta }_3\right)-\mathrm{l}\mathrm{n}\left(S_{\text{mid}}/S_{\text{end}}\right)\right]$ as in Eq. (45). The leading factor is calculated as $3/\left(2\tau b_s\right)\approx 0.55\mu {\mathrm{m}}^2/{\mathrm{m}\mathrm{s}}^2$ using $b_s\approx 4.59\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^2$ and $\tau =0.59\mathrm{m}\mathrm{s}$. Insets show a (natural) log-log plot obtained by first taking the absolute value of the VACF. The exchange time, ${\tau }_k=11.7\mathrm{m}\mathrm{s}$ is marked with a dotted line. A subset of the data is shown, omitting the longest mixing times at $t_m=[160,300]\mathrm{m}\mathrm{s}$ that are effectively 0. In the log-log inset, a linear fit to the points over $t_m=[20,40,80]$ (dashed line) indicates a power-law tail with $~t^{-2.4}$, or $\vartheta =1.4$, although this fit is highly sensitive to ${\beta }_3$.