Introduction

The diffusion propagator $$$P(\mathbf{x}′,t|\mathbf{x})$$$ represents the probability that a particle located at position $$$\mathbf{x}$$$ travels to $$$\mathbf{x}′$$$ over a time interval of duration $$$t$$$. This quantity fully describes the diffusive motion in simple as well as complicated environments, which could feature reflective, relaxing, or semi-permeable walls, external forces, spatially varying diffusivity, etc. Such structural parameters imprint their trace on the propagator. The propagator, if available, could reveal a detailed static picture of the space within which diffusion is taking place.
For example, consider diffusion within a domain bounded by partially relaxing walls. If the diffusion propagator is available, the bulk diffusivity can be determined from its short-time behavior; the domain shape is given by its support; and the surface relaxivity can be obtained from its decay, most conveniently at long times.
The most commonly employed scheme for diffusion encoding is that introduced by Stejskal and Tanner¹ in 1965, see Figure 1a. The MR signal attenuation for this sequence has been related to an ensemble average propagator (EAP) since its early days². The EAP represents the distribution of net displacements, which is the diffusion propagator substantially compromised due to an integration over all possible initial positions.
More recently, another two-pulse sequence has been introduced³ (see Figure 1b) featuring one long pulse compensating another narrower pulse. This method provides a means to determine the pore geometry, or more adequately, the steady-state distribution of the particles, which is just the long-time limit of the diffusion propagator.
Here, a new experimental and data-analysis framework is introduced, which reveals the actual diffusion propagator.

Methods and Results

Pulse sequence and theory

One realization of the effective waveform proposed in this work is depicted in Figure 1c. The sequence comprises two independently-determined narrower gradient pulses, whose collective integral is compensated by one long pulse. Within a $$$d$$$-dimensional pore space, and in the ideal scenario of infinitely long first pulse followed by infinitesimally narrow pulses, the detected MR signal is simply the Fourier transform $$E_\Delta(\mathbf{q},\mathbf{q}')=\int \mathrm{d} \mathbf{x} \, \rho(\mathbf{x}) \int \mathrm{d}\mathbf{x}' \, P(\mathbf{x}', \Delta | \mathbf{x}) \, e^{-i (\mathbf{q}\cdot\mathbf{x}+\mathbf{q}'\cdot\mathbf{x}' )}\ ,$$
allowing the determination of the propagator via inverse Fourier transforms through the expression
$$P(\mathbf{x}',\Delta|\mathbf{x}) =\frac{\int\mathrm{d}\mathbf{q}\, e^{i\mathbf{q}\cdot\mathbf{x}}\int\mathrm{d}\mathbf{q}'\, e^{i\mathbf{q}'\cdot\mathbf{x}'}\, E_\Delta(\mathbf{q},\mathbf{q}')}{(2\pi)^d \int\mathrm{d}\mathbf{q}\,e^{i\mathbf{q}\cdot\mathbf{x}}\, E_\Delta(\mathbf{q},\mathbf{0})}\ .$$

Validation and beyond

The technique was validated for the problem involving one-dimensional diffusion near an attractive center² for which the analytical expression for the propagator is available⁴. The signal and its transform via the above expression were evaluated for arbitrary parameters. It was shown to yield the analytical form of the propagator for ideal pulses, i.e., in the limit $$$\delta_1\rightarrow \infty$$$, $$$\delta_2\rightarrow 0$$$, and $$$\delta_3 \rightarrow 0$$$.
Figure 2 shows a summary of our results. As the duration of the long pulse is shortened, the propagator obtained through the above procedure is blurred. On the other hand, the finite width of the narrow pulses changes the meaning of the propagator such that the displacements are defined between the centers-of-mass of the trajectories traversed during the application of the pulses⁵.
A signature of the presence of external forces is the violation of the reciprocity condition, i.e., $$$P(\mathbf{x}′,t|\mathbf{x})\neq P(\mathbf{x},t|\mathbf{x}’)$$$, which is revealed by the asymmetry of the estimated propagator about the $$$\mathbf{x}=\mathbf{x}’$$$ line.


Exchanging restricted compartments

By repeating the experiment with different values of $$$\Delta$$$, the method allows the study of diffusive dynamics. In Figure 3, we illustrate our results for diffusion within and between two restricted compartments featuring different diffusivities as well as sizes separated by a thin semi-permeable membrane. For this one-dimensional problem, the propagator can once again be displayed on a plane for each diffusion time (first three columns). At short times, the diffusivities can be estimated. In the intermediate time point considered, the molecules spread within their compartment and start leaking into the neighbouring one. At longer times, a significant portion of the molecules end up in the other compartment. In the last column of Figure 3, we show the results for parameters that can be employed in some preclinical scanners. The method provides interpretable results elucidating the structure within the pore space in great clarity. This is in sharp contrast to the case of EAP (bottom row) and the steady-state distribution (Figure 4), which are non-specific even for ideal experimental conditions.

Discussion and Conclusion

We introduced a new pulse sequence featuring two independent gradients and related the signal to the diffusion propagator through an expression involving Fourier transforms. Note that our method reduces to the Stejskal-Tanner sequence for $$$\mathbf{q}'=-\mathbf{q}$$$, and to Laun et al.'s sequence for $$$\mathbf{q}'=\mathbf{0}$$$. Our technique should not be confused with position exchange spectroscopy⁶, which also employs two independent pulses to encode positions making it suitable for flow measurements. For diffusion within a heterogeneous specimen composed of many pores, a long pulse is introduced to prevent large cancellations of signals for pores located at different positions.

The diffusion propagator is the key quantity that fully describes diffusion for Markov processes. Compared to EAP, the true diffusion propagator provides exquisite sensitivity to the fine characteristics of the pore space because of the extra dimension made available by our technique. As such, the method could be a powerful probe that could transform our capabilities of imaging tissue microstructure.

Acknowledgements

After conceptualizing the technique, the author had discussions with Cem Yolcu, who performed the simulations in Figure 3, Carl-Fredrik Westin, and Magnus Herberthson. Nicolas Moutal and Denis Grebenkov provided the code for their work (J Sci Comput, 81, 3, p. 1630-1654, 2019), which was helpful for the simulations in Figure 3.