Complete fourier direct magnetic resonance imaging (CFD-MRI) for diffusion MRI

1. Introduction

Since the conception of mathematical models for the effect of the magnetic moment diffusion in nuclear magnetic resonance (NMR) experiments by Hahn (1950), Carrand Purcell (1954), and Torrey (1956), several methods have been proposed for analysis of diffusion-weighted (DW) magnetic resonance imaging (MRI) signal. These advancements endowed with the non-invasive, in vivo nature of the technique, have propelled the initial utilization of DW imaging measures, e.g., apparent diffusion coefficient in early detection of ischemia (Moseley et al., 1990; Baird and Warach, 1998), to many highly crucial areas in research and clinical imaging: for example in cancer diagnosis (Song et al., 2002; Turkbey et al., 2009; Xu et al., 2009), follow-up on treatment, pre- and post-operative assessment for different organs [e.g., fiber tracking (Conturo et al., 1999; Mori and van Zijl, 2002)] white matter integrity assessment (Budde et al., 2007; Correia et al., 2008) as in monitoring of neurological diseases such as multiple sclerosis (Song et al., 2002) and disorders (Ciccarelli et al., 2008) like schizophrenia (Seal et al., 2008; Voineskos et al., 2010) and Alzheimer's disease (Mielke et al., 2009), as well as neonatal development (McKinstry et al., 2002) and traumatic brain injury (Mac Donald et al., 2011).

In brief, diffusion weighted magnetic resonance imaging (DW–MRI) has become an indispensable and versatile technique playing an important role in several applications by its ability to estimate diffusion. The abundance of DW–MRI models is an indicator of room for improvement as well as the necessity for unification [see Özcan et al. (2012) for a detailed account of the partial differential equation (PDE) based adaptation's implications as well as a thorough mathematical analysis and a description of the background of existing methods].

DW–MRI's aim is to obtain measures and characterization of microstructure by investigating the diffusion process. Several methods and models have been proposed, all originating from the seminal work of Stejskal and Tanner (1965). Therein, under the influence of the additional motion sensitizing magnetic field gradients, the self-diffusion PDE of the magnetic moments is included in the Bloch PDE to model the attenuation in the DW–NMR spectroscopy signal. The result is the estimation of the scalar diffusion coefficient of the entire sample. In a sense, DW–NMR added another dimension, i.e., the magnetic moment motion, to the spectroscopic information even before the introduction of magnetic moment position later by the invention of MR imaging.

Accordingly, herein, DW–MRI model is naturally progressed to a higher dimensional construct that jointly presents magnetic moment position and motion. This is achieved by carefully re-examining the first principles of DW–MRI signal formation using particle methods in the spirit of the work of McCall et al. (1963). The mathematical model constructed in section 2.1 is specifically obtained for DW–MRI signal (rather than DW–NMR) by including all of its elements (e.g., imaging gradients) using complex values without taking the signal's magnitude.

The approach reveals that for an l_mr-dimensional MRI slice, the DW–MRI complex valued signal that comes out of the scanner is the (l_mr + 3)-dimensional Fourier transform of the joint distribution function of the number of magnetic moments (that are at a given position at the initial time) and their displacement integrals. In other words, the first l_mr dimensions correspond to the usual MRI k-space with position information and the remaining three dimensions constitute the frequency space of displacement integrals. The values of imaging and motion sensitizing magnetic field gradient vectors together define in the (l_mr + 3)-dimensional Fourier space the sampling points of the joint distribution function's Fourier transform. The distribution function is recovered by taking the Fourier transform of the complete data directly (i.e., without any scaling or use of magnitudes), giving the method its name: Complete Fourier Direct (CFD) MRI (Özcan 2010b).

2. Materials and Methods

2.1. Complete Fourier Direct MRI Signal Formation

The MR signal is generated by the vectorial sum of transverse magnetization of magnetic moments (Haacke et al., 1999):

$\begin{array}{cc}M(t)={\displaystyle \sum _i{m}_i(t)}.& (1)\end{array}$

By neglecting the effect of spin–spin relaxation, the evolution of the transverse magnetization of the ith magnetic moment is described in a standard manner by a rotating magnetization vector according to Bloch equations (see Appendix B):

Here, γ is the gyromagnetic ratio, the transverse magnetization vector, m_i, is written in complex number form with m_i(t₀) denoting the initial magnetization tipped to the transverse plane,

$\begin{array}{cc}{\Omega }_i={\displaystyle \underset{t_0}{\overset{t}{\int }}G}(x_i,\tau )·x_i(\tau )d\tau & (3)\end{array}$

describes the phase (when multiplied by γ) as a function of the magnetic field gradients $G(x,t)\in {ℝ}^3$, and the position of the magnetic moment is $x_i\in {ℝ}^3$. The time-dependent position, x_i, in Equation (3) affects the phase, Ω_i, thereby also affecting the total signal in Equation (1).

Any kind of displacement (such as Brownian motion, molecular movement in biological tissue with different medium and obstacles, coherent motion or any combination thereof) is incorporated straight into the model, by modeling the position in a general and direct form herein without any stochastic assumptions [such as Markovian property used in Wedeen et al. (2005)] on the motion

$\begin{array}{cc}{x}_i(t)=x_i(t_0)+w_i(t),& (4)\end{array}$

where $w_i(t)\in {ℝ}^3$ represents the displacement of the magnetic moment from its initial position, x_i(t₀), [i.e., w_i(t₀) = 0]. The only physical requirement is the continuity of w_i(t) since a magnetic moment cannot disappear at a given point and reappear at another.

The signal is calculated using Equations (1– 4) in reverse order. Following the motion described by Equation (4), the phase of the ith magnetic moment in Equation (3) during the digital acquisition period of the two dimensional imaging (l_mr = 2) pulsed-gradient spin-echo (PGSE) sequence of Figure 1, is obtained after tedious but routine derivations [see Appendix B for a brief exposition of the derivations and Özcan (2012)] using the definitions of the variables in Figure 1 with $G_{*}\in {ℝ}^3$ denoting the magnetic field gradient vectors labeled as read-out, ro; phase encode, pe; slice select, ss; and diffusion, D. Omitting routine calculations for trapezoidal shapes for clarity, the derivation is carried out assuming ideal gradient amplifiers providing rectangular shaped gradient pulses. The initial time, t₀, is the end time of π/2 radio frequency (RF) pulse when the magnetization is fully tipped to the transversal plane. The resulting phase of the transverse magnetization is a function of time, t, and, the imaging and diffusion gradients (see Appendix B):

$\begin{array}{l}\begin{array}{cc}\begin{array}{l}{\Omega }_i(t,G_{\text{pe}},G_{\text{ss}},G_{\text{D}})\\ =\left((t-t_{\text{acq}}-\Delta t_{\text{rw}})G_{\text{ro}}-\Delta t_{\text{rw}}{G}_{\text{pe}}\right)·x_i(t_0)-\Delta t_{\text{rw}}{G}_{\text{ss}}·x_i(t_0)\end{array}& (5)\end{array}\\ \begin{array}{cc}+G_{\text{D}}·W^{\text{d}}+((t_{d4}-t_{d3})-(t_{d2}-t_{d1}))G_{\text{D}}·x_i(t_0)+{\Phi }_{\pi i}& (6)\end{array}\\ \begin{array}{cc}+G_{\text{ro}}·W_i^{\text{acq}}(t)-(G_{\text{ro}}+G_{\text{pe}}+G_{\text{ss}})·W_i^{\text{rw}}.& (7)\end{array}\end{array}$

FIGURE 1

Figure 1. The pulsed-gradient spin-echo (PGSE) pulse sequence and the definition of the variables used in the calculations. SS–EX is for the slice select gradient during the excitation (π/2) pulse, RO for read out, PE for phase encode, SS for is the slice select gradient, Diff marks the motion sensitizing pulses, SS–PI is the slice select gradient during the π-pulse and ACQ stands for digital acquisition period. In practice, the MR pulse sequences implement the rewind (rw) gradients such that the amplifiers are turned on and off at the same times.

The second term in Equation (6) removes the injection of the initial position into the DW signal because of equal pulse duration times, δ = t_d4 − t_d3 = t_d2 − t_d1. The term Φ_πi describes the systematic phase (see Appendix B1) created by the π-pulse slice select gradient (SS–PI) and in this work it will be automatically taken out by the phase correction algorithm in section 2.2. Equations (6) and (7) incorporate three integrals of the displacement w_i(t): (W^d_i, W^acq_i, W^rw_i) corresponding to the displacement integrals for diffusion (d), analog to digital conversion acquisition (acq), and initial rewind (rw) gradient time periods, respectively

$\begin{array}{l}\begin{array}{cc}{W}_i^{\text{d}}={\displaystyle \underset{t_{d3}}{\overset{t_{d4}}{\int }}{w}_i}(\tau )d\tau -{\displaystyle \underset{t_{d1}}{\overset{t_{d2}}{\int }}{w}_i(\tau )d\tau },& (8)\end{array}\\ \begin{array}{cc}{W}_i^{\text{acq}}(t)={\displaystyle \underset{t_{\text{acq}}}{\overset{t}{\int }}{w}_i}(\tau )d\tau ,W_i^{\text{rw}}={\displaystyle \underset{t_0}{\overset{t_{\text{rw}}}{\int }}{w}_i}(\tau )d\tau .& (9)\end{array}\end{array}$

First term in Equation (5) is the definition of the k-space in regular MRI, $k_{\text{mr}}\in {ℝ}^2$:

$\begin{array}{cc}{k}_{\text{mr}}(t,G_{\text{pe}})=(t-t_{\text{acq}}-\Delta t_{\text{rw}})G_{\text{ro}}-\Delta t_{\text{rw}}{G}_{\text{pe}}& (10)\end{array}$

with the additional term,

$\begin{array}{cc}{\phi }_{\text{slice}}=-\Delta t_{\text{rw}}{G}_{\text{ss}}·x_i(t_0),& (11)\end{array}$

which is constant because the slice select axis component of x_i(t₀) is the slice position¹.

Without loss of generality, by adopting the imaging coordinate frame defined by the directions of the read-out, phase encode and slice select gradients, G_ro = [g_ro1, 0, 0], G_pe = [0, g_pe2, 0], G_ss = [0, 0, g_ss3], time and g_pe2 become functions of k_mr = [k_mr1, k_mr2, k_mr3] using Equation (10):

$t=k_{\text{mr}1}/g_{\text{ro}1}+t_{\text{acq}}+\Delta t_{\text{rw}}\text{ and }{g}_{\text{pe2}}=-k_{\text{mr2}}/\Delta t_{\text{rw}}.$

Accordingly, W^acq_i(t) becomes a function of k_mr1 and the coefficients of W^rw in Equation (7) are written as a vector which is an affine function of k_mr:

$\begin{array}{cc}{k}_{\text{rw}}=\left[-g_{\text{ro}1},k_{\text{mr}2}/\Delta t_{\text{rw}},-g_{\text{ss}3}\right].& (12)\end{array}$

Consequently, the phase, Ω_i, of Equation (3) is expressed concisely by defining k_D = G_D:

$\begin{array}{cc}\begin{array}{l}{\Omega }_i=k_{\text{mr}}·x_i(t_0)+k_{\text{D}}·W_i^{\text{d}}+k_{\text{rw}}·W_i^{\text{rw}}\\ +W_{i,1}^{\text{acq}}(k_{\text{mr}1})g_{\text{ro}1}+{\phi }_{\text{slice}},\end{array}& (13)\end{array}$

reflecting the effect of initial position and displacement integrals on the phase² of each magnetic moment. Since ϕ_slice is constant for all i, it is taken out of Equation (13) with the appropriate rotation of the magnetization coordinate frame on a slice by slice basis.

Finally, assuming that all of the magnetic moments have the real valued initial magnetization m_i(t₀) = m₀, Equation (1) can be re-written using Equation (13) to reveal a Fourier relationship,

A more efficient way to evaluate the sum in Equation (14) is first to group the magnetic moments with the same position-displacement properties and then to count the numbers elements in the groups.

Definition 1. The joint position-displacement integral distribution function, P^total_cfd(x, W), is defined as the number of magnetic moments with the initial position $x\in {ℝ}^3$ at time t₀, possessing the displacement integral values of $W=(W^{\text{d}},W^{\text{rw}},W_1^{\text{acq}})\in {ℝ}^7$.

The signal in Equation (14) is calculated by integrating over the whole position-displacement space (absorbing m₀ into P^total_cfd for ease of notation):

Equation (15) is by definition the Fourier transform of P^total_cfd with non-linearities added by W^acq₁.

Among the elements of W, the focus is on the most descriptive MRI observable, W^d. Its distribution is obtained by marginalizing W^rw from Equation (15)

$\begin{array}{cc}\begin{array}{l}{P}_{\text{cfd}}(x,W^{\text{d}})={\displaystyle \int P_{\text{cfd}}^{\text{total}}(x,W^{\text{d}},W^{\text{rw}})d{W}^{\text{rw}}}\\ ={ℱ}_{(k_{\text{mr}},k_{\text{D}})}^{-1}\{S_{\text{cfd}}(k_{\text{mr}},k_{\text{D}},0)\}.\end{array}& (17)\end{array}$

However, the affine dependence of k_rw in Equation (12) makes it impossible to fix k_rw = 0 and to sample in (k_mr, k_D, 0) subspace. The following example demonstrates how the affine dependence affects the measurements by using a two dimensional Gaussian function, exp(−(k²₁ + k²₂)), with the “undesirable” variable k₂ sampled on a line k₂ = a k₁ + b:

The result is complex valued in comparison to the real valued Fourier transform of exp(−k²₁) which can be obtained by setting a = b = 0.

2.2. Phase Corrections for the Estimation of P_cfd

In addition to inherent affine dependence and non-linearities, different experimental factors, noise, hardware imperfections etc., affect the DW–MR signal adversely. CFD–MRI addresses these issues by adopting a pivotal, physically meaningful standpoint originating from the following Fourier transform property (Bracewell, 2000):

(real valued function) ⇆ ℱ ⇆ (Hermitian symmetric function).

Accordingly CFD–MRI reconstruction is based on the following:

Furthermore, an immediate implication of the transform property and the Hermitian symmetry of S_cfd is that theoretically, taking its magnitude before Fourier transforming will result in a symmetric real valued distribution under ideal conditions. As noted above, in practice the experimental S_cfd is never Hermitian symmetric resulting in an asymmetric magnitude. Consequently, the magnitude's Fourier transform used in existing methods (Callaghan, 1991; Wedeen et al., 2005), results in a complex valued (Hermitian symmetric) distribution function. The difficulty of a physical interpretation forced those methods to take the magnitude of the transform as well to obtain a real valued function.

Herein, in order to obtain a real valued P_cfd, the re-establishment of Hermitian symmetry in the signal during the computation of the inverse transform for Equation (17), is realized by phase corrections. The strategy is similar in principle to the correction of the k_mr-space center's (echo time) shift during the read-out period of acquisition in MR imaging. The resulting linear phase shift in the physical read-out axis uniformly and systematically appears in all of the data. The shifts in both phase-encode (e.g., due to sample shaking) and read-out directions are corrected by first determining from the Fourier transform in k_mr,

$\begin{array}{cc}{I}_{k_{\text{mr}}}^{\text{complex}}((x_{\text{ro}},x_{\text{pe}}),k_{\text{D}})={ℱ}_{k_{\text{mr}}}^{-1}\{S_{\text{cfd}}(k_{\text{mr}},k_{\text{D}},k_{\text{rw}})\},& (18)\end{array}$

the angle (∠ I_k^complexmr) from the signal regions along the center lines of each physical direction at k_D = 0,

$\begin{array}{cc}{r}_{\text{ro}}(x_{\text{ro}})\approx \angle I_{k_{\text{mr}}}^{\text{complex}}((x_{\text{ro}},0),0),r_{\text{pe}}(x_{\text{pe}})\approx \angle I_{k_{\text{mr}}}^{\text{complex}}((0,x_{\text{pe}}),0).& (19)\end{array}$

The phase corrections are then applied systematically at each value of k_D (see Figure 3):

The Fourier transform in the remaining variables,

$\begin{array}{cc}{P}_{\text{cfd}}(x,W^{\text{d}})={ℱ}_{k_{\text{D}}}^{-1}\{I_{k_{\text{mr}}}(x,k_{\text{D}})\},& (21)\end{array}$

is evaluated sequentially in each k_D-dimension with the aim of re-establishing the Hermitian property, I_kmr(x, −k_D) = I_k*mr(x, k_D), using the following steps.

CFD Phase correction algorithm:

1. Pick a pixel at location x_c, preferably near the center of the image where tissue or a good signal area is located.
2. Starting from the first direction, l = 1 of the k_D space calculate the phase on the line passing through the origin (i.e., [k_D1, 0, 0], [0, k_D2, 0], [0, 0, k_D3], respectively for l = 1, 2, 3), e.g., ∠I_kmr(x_c,(k_D1, 0, 0)).
3. Investigate the plot of the phase versus k_Dl. Pick as many as possible consecutive values of k_Dl near 0 without sudden changes to assure high signal to noise value.
4. Construct a polynomial of degree m (with m less than the number of points) that approximates the phase at the selected points. The polynomial's constant term must be set to be 0 to guarantee that I_kmr(x_c, 0) remains unchanged. For example, for the first direction, at selected values of k_D1 the polynomial looks like
   $\begin{array}{cc}\begin{array}{l}\angle I_{k_{\text{mr}}}(x_c,(k_{\text{D}1},0,0))\approx r_{\text{D}1}(k_{\text{D}1})\doteq a_m{\left(k_{\text{D}1}\right)}^m\\ +a_{m-1}{\left(k_{\text{D}1}\right)}^{m-1}+\dots +a_1{k}_{\text{D}1}\end{array}& (22)\end{array}$
   as demonstrated in Figure 2.
5. Apply the same phase correction systematically to the entirety of the data along the lth direction at each of the other dimensions, at all of the pixel locations. For example in the first direction, k_D1, update I_kmr to be equal to for all k_D2, k_D3 and x.
6. Repeat steps 2–5 for the remaining directions.

FIGURE 2

Figure 2. Özcan (2012): Top row: The Nyquist plots of uncorrected (left) and corrected (right) S_cfd obtained from the experimental data described in section 2.4. The plots show data acquired at each diffusion gradient value k_D1 on the complex plane. Bottom row: The magnitude and phase plots of the data. Uncorrected data (bottom row left, second column) exhibit a linear phase shift around 0 frequency, indicative of coherent motion. After the phase corrections obtained using the polynomial 0.266k_D1 estimated from the points k_D1 = ‒6, 0, 6, 12 Gauss/cm, the magnitude is unchanged but the signal's imaginary part is smaller for the corrected values visible by the difference between the vertical axis spans of Nyquist plots and the phase plots.

The algorithm transfers the signal to the real channel by preserving its energy as the phase corrected spin-echo image without diffusion gradients, I_kmr(x, 0). The distribution of the magnetic moments with low mobility in all three directions [i.e., P_cfd(x, 0)] shows the result of the transfer in Figure 3 for the sample described in section 2.4.

FIGURE 3

Figure 3. On the left, the fixed baboon brain images acquired in the anatomical-coronal plane. On the top row, the real and imaginary parts of I_kmr(x, 0) and on the bottom row, P_cfd(x, 0) are displayed. The imaginary part is approximately 10% of the real part in both cases. On the right, full representation of P_cfd with isosurfaces (${\overline{P}}_{\text{cfd}}=0.17$, see section 3.1) around CC and EC junction. Starting from left bottom going clockwise, the sample pixels are from cerebro-spinal fluid (CSF), CC, white matter (WM) and CC, and EC junction, respectively (see also Figure 5).

In P_cfd(x, 0), areas with high level of organization inducing low mobility, such as the corpus callosum (CC), the external capsule (EC), the mid-brain and the pons, appear brighter. The image is not an anisotropy map, e.g., mineral oil would appear brighter than water due to a smaller diffusion coefficient despite both liquids being isotropic. Spin-echo image is more blurred because it is a low pass filtered version of P_cfd:

$\begin{array}{cc}\begin{array}{l}{I}_{k_{\text{mr}}}^{\text{complex}}(x,k_{\text{D}})={ℱ}_{k_{\text{mr}}}^{-1}\{S_{\text{cfd}}(k_{\text{mr}},k_{\text{D}},k_{\text{rw}})\}\Rightarrow I_{k_{\text{mr}}}(x,0)\\ =\frac{2\pi }{\gamma }{\displaystyle \int P_{\text{cfd}}^{\text{total}}}(x,W)dW\end{array}& (23)\end{array}$

(see Appendix C).

CFD phase correction algorithm outperformed the fitting of the phase values up to the fourth degree multinomials in ${ℝ}^3$. The reasons behind this outcome, which will provide information about DW–MR signal artifacts, as well as inclusion of different functions for corrections will be investigated in the future.

2.3. CFD-MRI Sampling and Windowing

Whereas the standard MRI field of view (FOV) calculations (Haacke et al., 1999) are used for k_mr-space, the infinite bandwidth in k_D-space due to P_cfd's finite support in W^d-space (originating from finite length displacements) falls beyond the reach of the gradient hardware's limits for small diffusion gradient duration and separation times (δ and Δ, respectively in Figure 1). Even with a powerful gradient system, a large magnitude of k_D causes substantial signal uncertainties due to an increasing performance deterioration as the power requirements push the hardware to its limitations.

With such a hardware constraint, in order to reduce ripple effects caused by truncation, P_cfd's bandwidth (i.e., S_cfd's support) is shrunk by increasing δ and Δ causing the dispersion (covariance) of the displacement integral W^d (and therefore P_cfd's support) to increase. This is directly visible in the special case of Brownian Motion characterized by the diffusion tensor D in Figure 4. P_cfd and S_cfd are zero mean Gaussians with covariances, respectively equal to (see Appendix A)

FIGURE 4

Figure 4. Effects of fitting the bandwidth within the field of view (FOV) in the Fourier space on the left, with the covariance of the Gaussians given as $\widehat{D}={({\delta }^2(\Delta -\delta /3)D)}^{-1}$. The sampling is done in the fixed interval, [‒5, 5], of the Fourier Space followed by reconstruction on the right using the discrete Fourier transform. Theoretically this is equivalent to convolution with the sinc function (see Brigham, 1988) in physical space. The ripple effects created by sinc lobes diminish on the left as the Gaussian falls into the FOV with increasing δ, Δ but constant D.

$\begin{array}{cc}E[W^{\text{d}}{(W^{\text{d}})}^T]=b_{t}D\text{ and }{(b_{t}D)}^{-1}\text{ where }{b}_t\doteq {\delta }^2(\Delta -\delta /3)& (24)\end{array}$

(see Özcan (2009, 2010a) because the Fourier transform of a Gaussian with a covariance matrix $\widehat{D}$ is also a Gaussian with covariance ${\widehat{D}}^{-1}$:

$\begin{array}{cc}ℱ\left\{\exp \left({\left(W^{\text{d}}\right)}^T{\widehat{D}}^{-1}{W}^{\text{d}}\right)\right\}~\exp \left(k_{\text{D}}^T\widehat{D}{k}_{\text{D}}\right).& (25)\end{array}$

The procedure is graphically displayed in Figure 4 also emphasizing the effect of ripples on the small values of P_cfd which are especially important in revealing microstructure as explained in section 3.1.

The second sampling criterion is an appropriate sampling rate i.e., sufficient number of points in k_D-space to prevent aliasing artifacts on P_cfd. This is constrained by the time available for acquisitions as each point in k_D-space requires the scan time of the entire k_mr-space.

2.4. Experimental Setup and Analysis Methods

A fixed baboon brain immersed in 4% paraformaldehyde was used for the experiments. The primate was prematurely delivered on the 125th day and sacrificed on the 59th day after delivery. All animal husbandry, handling, and procedures were performed at the Southwest Foundation for Biomedical Research, San Antonio, Texas. Animal handling and ethics were approved to conform to American Association for Accreditation of Laboratory Animal Care (AAALAC) guidelines. Further details of the preparation are described in Kroenke et al. (2005).

The experiments were carried out on a 4.7 Tesla MR scanner (Varian NMR Systems, Palo Alto, CA, USA) with a 15 cm inner diameter gradient system, 45 Gauss/cm maximum gradient strength and 0.2 ms rise time using a cylindrical quadrature birdcage coil (Varian NMR Systems, Palo Alto, CA, USA) with 63 mm inner diameter. CFD–MRI data were obtained using the standard pulsed-gradient spin-echo multi-slice sequence. The k_mr-space was sampled to result in images of 128 × 128 pixels with a FOV 64 × 64 mm² and 0.5 mm slice thickness. The k_D-space was sampled in a uniformly spaced Cartesian grid in a cube [−30, 30 Gauss/cm]³ with 6 Gauss/cm sampling intervals at each dimension resulting in 11× 11 × 11 voxels. The repetition time T_R = 1 s, echo time T_E = 56.5 ms, diffusion pulse separation time Δ = 30 ms and diffusion pulse duration δ = 15 ms were used.

The data were transferred to a two quad core 2.3 GHz Intel Xeon^® cpu and 8 GB memory Dell Precision Workstation 490 running Windows XP^® 64-bit operating system. The DWI data were placed in a 5-dimensional array in the computer memory and the discrete Fourier transform (DFT) was computed along with the phase corrections. In-house Matlab^® (Mathworks, Natick, MA, USA) programs were used for all of the computations and to display the graphics and maps.

3. Results

3.1. Visualization of the CFD Distribution

The joint distribution's high-dimensionality [e.g., two dimensions for position in regular MR images (l_mr = 2), plus three more for displacement integrals] creates a visualization challenge which is addressed herein by using P_cfd(x, 0) as the background image. Furthermore, the isosurface³ of normalized P_cfd,

${\overline{P}}_{\text{cfd}}(x,W^{\text{d}})=\frac{P_{\text{cfd}}(x,W^{\text{d}})}{P_{\text{cfd}}(x,0)}$

is overlayed on the pixel at location x, as in Figure 3. For the sake of an objective assessment, the isosurfaces are defined using a common level value c (0 < c ≤ 1),

$\left\{W^{\text{d}}\in {ℝ}^3:{\overline{P}}_{\text{cfd}}\left(x,W^{\text{d}}\right)=c\right\},$

over all locations. The key point is the choice of an appropriate c-value that will reveal the outskirts of P_cfd corresponding to the small number of “scout” magnetic moments that travel further away thereby portraying the microstructure of the environment. In summary,

1. Too high values do not provide enough structural information (see first rows in Figure 5).
2. The appropriately informative value depends on the properties of the motion (thus of the microstructure) at a given location (compare columns of Figure 5, right side).
3. Too low values force the isosurfaces to become extremely noisy (see last row of Figure 5).

As the motion in highly organized tissue is less dispersed (i.e., a smaller support for P_cfd which implies a larger support for S_cfd thereby causing bigger truncation effects), increasing the diffusion gradient times δ and Δ in section 2.3 will create almost “flat” displacement integral distributions at an isotropic medium like CSF. In this case, the small valued distribution [caused by constant integral value ∫ P_cfd(x, W^d) dW^d = number of particles] is susceptible to noise, creating noisy isosurfaces of Figure 5 for all level values. In contrast, for the experiments conducted with an isotropic (water) phantom at much lower δ and Δ values, the isosurfaces were spheres for a wide range of c-values (not shown).

FIGURE 5

Figure 5. On the left, one dimensional graphical representation of the choice of c-value. The drawing below the horizontal axis displays the structure of the sample in the infinitesimal volume element [‒dx, dx] as seen by the magnetic moments from their initial position. The right side of the sample has less restrictive properties as on the boundary of tissue with a liquid, such as CC and cerebrospinal fluid (CSF). The noiseless isosurfaces consist of two points shown as dots on the graph. Low c-values correspond to the magnetic moments with longer travel paths providing more structural information than high c-values. However, too low values create noisy and disconnected isosurfaces represented with more than two points on the drawing. On the right, isosurfaces from different pixels in the baboon brain marked in Figure 3 demonstrate the effect of c-value on the information content from top = insufficient to bottom = noisy (see the text for the CSF column).

Figure 6 presents different isosurfaces that elucidate the tissue structure on the pons and the mid-brain of the fixed baboon brain sample described in section 2.4. The tracts on the left and right side of the mid-brain are visible with the ellipsoidal looking isosurfaces. Five isosurfaces from the same row of pixels marked on Figure 6 are displayed on Figure 6 corresponding to two different c-values. The green isosurfaces with larger level values are smoother and less informative than the red ones with smaller c-values. Different viewpoints at each row of Figure 6 emphasize that the isosurfaces are 3-D objects. The figure demonstrates one of the challenges of presentation: the displacement integral values must be considered in ${ℝ}^3$ to grasp the complete information offered by CFD–MRI.

FIGURE 6

Figure 6. The isosurfaces (${\overline{P}}_{cfd}=0.14$) on the pons and the mid-brain. Each of the boxes indicate an isosurface presented, respectively on the right. Each column presents the same isosurface from different view angles. The dots on top of the frames are placed for orientative purposes. The surfaces are not necessarily ellipsoids and they have mostly an asymmetric structure. The outer red surface is the set ${\overline{P}}_{\text{cfd}}=0.14$ and the green surface is ${\overline{P}}_{cfd}=0.21$. The red surface envelops the green one.

Overall, the isosurfaces are not constrained to given forms like Gaussians, spherical harmonics or to any expansions. In fact, they are typically not even symmetric. They are structureless, general and direct.

Isosurface visualizations constitute only one method to present the high dimensional information obtained from CFD–MRI. Another example is the dimension reduction by means of computing the so called orientation distribution function (ODF) (Wedeen et al., 2005) obtained from radial integrals. However, for CFD–MRI the ODF raises the concern of inadequately presenting the “outskirts” of the P_cfd because the dispersion of the outgoing rays shown in Figure 7 jeopardizes the inclusion of the values further away from the origin (see also Figure 6).

FIGURE 7

Figure 7. An example of a set of rays from the origin along which the distribution function is integrated to obtain orientation distribution function. As the rays disperse with increasing distance from the origin the points describing the displacement of a smaller number of particles contributes less and less to the numerical integration due to sparse sampling on both surfaces. The encapsulation of the isosurfaces on the right with larger level values by the smaller ones in Figure 6 shows the information that would be missed with numerical radial integration. The utilization of isosurfaces is more informative as discussed in Figure 5.

New methods, additional elaborate schemes such as color coding for representation of three dimensional functions aimed at displaying relevant microstructural information of CFD are left for future studies.

4. Discussion

4.1. Comparison with Existing Methods

From a fundamental point of view, guided by the microstructure that surrounds them, the molecules are displaced due to thermal energy whether they are in the scanner or not. All existing DW–MR methods are designed with the same goal in mind: the reconstruction of the propagator⁴ that describes the displacement of the magnetic moment from the DW–MR signal.

However, as CFD–MRI demonstrates, from a systems science perspective the MRI scanner acts as a time-delay linear system with the input w_i [displacement in Equation (4)], and the output W_i [displacement integral in Equation (8)], in Figure 8. The parameters Δ and δ define the delay and filter parameters. Special attention is paid in CFD–MRI to isolate the Fourier variable, k_D = G_D from these parameters in contrast to the q-space variable (Callaghan, 1991): q = (2π)⁻¹ γ δ G_D and the b-value of DTI: b = γ² ‖G_D‖²₂ δ² (Δ − δ/3) (see Appendix A).

FIGURE 8

Figure 8. The DW–MRI signal from a signals and systems perspective (Özcan, 2010b) describing the MRI observables of the diffusion phenomenon. The input is the magnetic moment displacement w_i and the output is the displacement integral, W^d_i, defined in Equation (8).

The inverse problem of obtaining the propagator from the distribution of the displacement integrals is singular because of the one-to-many relationship between the displacement and its integral:

$\begin{array}{cc}{W}_i^{\text{d}}={\displaystyle \underset{t_{d3}}{\overset{t_{d4}}{\int }}{w}_i}(\tau )d\tau -{\displaystyle \underset{t_{d1}}{\overset{t_{d2}}{\int }}{w}_i(\tau )d\tau }\to w_i(t),& (26)\end{array}$

because all of the displacements with the same low frequency content in time are mapped to the same displacement integral value. This statistical accumulation prevents the determination of the propagator in a general environment from the distribution of displacement integrals⁵).

Existing methods' attempts to estimate the propagator relies on the narrow pulse approximation by assuming negligible pulse duration (δ = t_d2 − t_d1 = t_d4 − t_d3 in Figure 1) compared to pulse separation time, Δ, i.e., (δ << Δ) ⇒ (δ → 0) specifically in (Δ−δ/3) (Callaghan, 1991, p. 342). Under this short integration time assumption, in Wedeen et al. (2005) it is further argued that the approximation, W^d_i ≈ (x_i(t_d3) − x_i(t_d1)) δ = (w_i(t_d3) − w_i(t_d1)) δ, is plausible. Although by the intermediate value theorem and the sample path continuity of the Brownian motion (Shiryaev, 1995), the values of each integral in W^d_i are attained at a time point within the respective integration intervals, [t_d1, t_d2] and [t_d3, t_d4], the nowhere-differentiability of the sample paths (Shiryaev, 1995) implies that the intermediate time points satisfying the equality are not fixed as t_d1 and t_d3, but are themselves random variables. In consequence, without the inference of the displacements at fixed time points the propagator cannot be reconstructed.

Moreover, elaborate derivations carried in Wedeen et al. (2005) to model the propagator as a conditional probability, p(x(t_d3)|x(t_d1)), describing a Markovian process raise concerns specially in environments such as biological tissue since the particle's past collisions with microstructure guides its future displacements. In fact, while it violates the conditions of Wiener process (see Appendix A), this displacement memory provides the inference of the microstructure by way of affecting the displacements and consequently their displacement integral distributions. Accordingly, in CFD–MRI is indifferent to memory properties by modeling P_cfd as a joint distribution function of random variables instead of the conditional probability of a stochastic process.

A summary of CFD–MRI's detailed comparison with existing methods (Özcan, 2011; Özcan et al., 2012) is presented below for completeness of exposure. Namely, there exits two avenues for the path from DW–NMR to DW–MRI in the literature:

1. Model matching methods initiated by diffusion tensor imaging (DTI) (Basser et al., 1994; Mattiello et al., 1994) and further expanded with high angular resolution DW imaging (HARDI) (Frank, 2001), composite hindered and restricted model of diffusion (CHARMED) (Assaf and Basser, 2005), diffusion orientation transform (DOT) (Özarslan et al., 2006), two versions of the generalized DTI (GDTI) (Özarslan and Mareci, 2003; Liu et al., 2004), and diffusional kurtosis imaging (DKI) (Jensen et al., 2005).

2. Spectral methods originating from Callaghan's q-space (Callaghan et al., 1988) followed by the diffusion spectrum imaging (DSI) (Wedeen et al., 2005) and Q-ball imaging (Tuch, 2004).

With the exception of the GDTI presented in Liu et al. (2004) [see also the discussion in Özarslan et al., (2009)], all of the DW–MRI methods estimate symmetric quantities. The model matching methods project the data onto symmetric structures, such as ellipsoids in DTI or spherical harmonics in HARDI. The spectral methods use the magnitude of the signal in the Fourier transform resulting in symmetric functions (see section 2.2). It is difficult to imagine that molecular motion in a biological environment populated with different types of fluids, barriers and tissue would be symmetric at any given location, e.g., at the fiber junctions. Symmetry or lack of it ought to be determined by the data free of any constraints imposed by the model as in the implementation of CFD–MRI's unconstrained structure.

The Fourier relationship between signal and joint distribution function provides a complete understanding of model matching methods. The methods start by applying DFT to the data in the first l_mr (imaging) dimensions. Thus, in the analysis of model matching methods the first l_mr independent variables are the physical location. The three remaining untransformed variables are the independent variables of the Fourier reciprocal of displacement integral space, i.e., they are in the Fourier domain. The goal of model matching methods is to fit the preferred model to the displacement distribution function's Fourier transform, sampled at the (vector) values defined by the diffusion sensitizing gradients. In the case of Brownian motion this mixed variable (physical and frequency variables) approach is applicable because the diffusion coefficient D can be directly estimated from the Fourier domain by Equations (24) and (25). The mixed space, which works well for DW–NMR signal peak attenuation, is translated to DW–MRI at each position x by

$\begin{array}{cc}|I_{k_{\text{mr}}}^{\text{complex}}(x,k_{\text{D}})|=|I_{k_{\text{mr}}}^{\text{complex}}(x,0)|\exp (-H(k_{\text{D}}))& (27)\end{array}$

where I_k^complexmr is given in Equation (18) and the function H defines the model, e.g., the quadratic form of DTI, spherical harmonics of HARDI or higher tensor expansions of GDTI [see Özcan et al. (2012) for a detailed exposition]. In a sense, these methods' aim could be summarized as expanding the Fourier portion of the mixed signal space. The basic example with a single term in the expansion is DTI for which:

$\begin{array}{cc}H(k_{\text{D}})={\gamma }^2{b}_t{k}_{\text{D}}^{T}D{k}_{\text{D}}& (28)\end{array}$

where the calculation of b_t from the PDE approach in Özcan (2010a) and covariance of displacement integrals in Appendix A resulted in the same value: δ² (Δ − δ/3). In CFD parlance, by the Fourier relationship between Gaussians in Equation (25), the diffusion quadratic form, D, is estimated in k_D-space, without recourse to a Fourier transform because of its direct appearance in Equation (28) rather than its inverse, D⁻¹, in the Gaussian of motion space. The coefficient matrix defined by carefully selected vectors in k_D-space that satisfy the invertibility conditions (Özcan, 2005) is used for the estimation of D in the linear algebraic framework of symmetric matrices (Papadakis et al., 1999; Özcan 2010a) [also refer to Özcan (2010a) and Özcan (2012) for the correspondence with the B-matrix formulation of Basser et al., (1994)].

The magnitude-based Fourier relationship presented in q-space methodology (Callaghan et al., 1988) is the origin of spectral methods. In Callaghan's book (Callaghan, 1991), parallel to the historical development of DW models, the theory is first developed for NMR experiments [see Callaghan, (1991, Chap. 6)] using polarized neutron scattering analogy. However, the translation from NMR to MRI is presented [see asserting without proof that the imaging and displacement portions of the signal are separable [see Callaghan, (1991, Chap. 8, pp. 440)]. The derivations of section 2.1 demonstrate that this is not the case.

In addition, by the affine dependence of k_rw on k_mr CFD derivations show that the inseparability partially accounts for the non-Hermitian nature of the S_cfd. Taking the magnitude of the DW–MRI signal, as in the case of DSI (Wedeen et al., 2005), does not count as a phase correction. Figure 9 demonstrates that by preserving Hermitian property, CFD–MRI captures correctly the crossing fibers at the junction of the CC and EC.

FIGURE 9

Figure 9. The comparison of P_cfd (Özcan, 2011) (first row) with the diffusion spectrum imaging orientation distribution function (DSI–ODF) (second row). Both functions are calculated from the same data on the right junction of the corpus callosum (CC) and the external capsule (EC), specifically from the pixels marked on the Figure 3. The isosurface (${\overline{P}}_{\text{cfd}}=0.17$) captures the asymmetric structure of the fiber crossings while the ODF is constrained to be symmetric for all of the pixels. Note that in CSF, ODF detects structure which is not present in reality as indicated by CFD.

4.2. Conclusion and Future Studies

In the biomedical imaging modalities' grand aim of biomarker capability establishment, the discovery path for CFD–MRI passes through the distribution function:

$S_{\text{cfd}}\to P_{\text{cfd}}\to \text{Biologicalproperties}.$

With P_cfd in the middle, both sides of the path present themselves with important challenges.

First and foremost, in DW–MRI, the displacements without reference to initial positions [see Equation (6)] prevent the inference of microstructure position. For example, the distribution function of the biological phantom (Özcan et al., 2011) constructed with two crossing rat trigeminal nerve fibers is always in the form of two crossing bars across the origin regardless of the nerves' position as long as their relative angle is kept the same. Also in the same phantom, the agar gel (isotropic component) appears as a sphere around the origin of P_cfd domain without the possibility of identifying its location. As the distributions from various types of microstructural components accumulate around the origin, the discrimination level of overlaps, more prominent with increasing biological tissue complexity, directly defines the sensitivity and specificity for microstructure changing pathologies. The important goal is the assessment of the theoretical aspects of the distribution function in order to understand whether it can detect in a timely manner, e.g., before significant disease progression, those changes. The determination of biophysical conditions behind the asymmetry (see also Özarslan et al., 2008; Özarslan, 2009) in the distribution functions is also part of the same goal.

However, the absence of analytical descriptions for P_cfd even in simple environments requires the investigations to be conducted with numerical simulations (Özcan et al., 2011) of particle motion within carefully designed geometries (Landman et al., 2010) and locally variable diffusivities. Along with numerical phantoms mimicking biological ones (Özcan et al., 2011), histopathological information is also being used for interpretation and validation (Budde and Frank, 2012; Budde and Annese, 2013). Additionally varying the time parameters δ and Δ will exploit the filtering effects caused by the displacement integral that will determine whether further information extraction is possible by expanding data acquisition with an appropriate set of parameter values.

On the other hand, on the discovery path's initiation by CFD–MRI signal formation, the re-establishment of Hermitian symmetry requires, in addition to the theoretical reasons presented herein, the analysis and quantification of Hermitian disruptive artifacts and systematic conditions in real data. Constructed by initially experimenting with elementary phantoms (e.g., water and mineral oil), this signal model expansion is necessary for the development of more elaborate systematic phase corrections, possibly by utilizing complex analysis theory. Specifically, accurate estimation of the pertinent Fourier transform of P_cfd from real data points in a clinical setup is targeted by the adaptation of CFD–MRI to fast sequences⁶, such as echo planar imaging (EPI) prone to Eddy current artifacts. The model will be expanded up to the point of reaching only minimal incremental improvements with new phase correction algorithms. Thereafter, relying on the residuals' content, which is free from displacement effects consequent to the application of system-wide uniform phase corrections⁷, more effective biomarker construction would be possible by the inclusion of extra information such as tissue susceptibility (Liu et al., 2013). Likewise, on a larger scale CFD–MRI's general aim is to improve outcomes of multimodal imaging, e.g., in prostate cancer strategies (Turkbey and Choyke, 2012).

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This study was supported, in part, by the National Cancer Institute funded Washington University Small Animal Imaging Resource (U24-CA83060), the NIH/NINDS grant Biomarkers and Pathogenesis of MS (P01-NS059560), and the US Army grant, NeuroPerformance Imaging (W23RYX1089N603). Thanks to J. Quirk, T. Conturo, S.-K. Song, K. Uḡurbil, G. Sapiro, K. Wong, N. Tsekos, and S. K. Mun for valuable discussions, the reviewers for their valuable comments and C. D. Kroenke for the sample preparation.

The manuscript is dedicated to Ayla and Halil Özcan.

References

Assaf, Y., and Basser, P. (2005). Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. Neuroimage 27, 48–58.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Baird, A. E., and Warach, S. (1998). Magnetic resonance imaging of acute stroke. J. Cereb. Blood Flow Metab. 18, 583–609.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Basser, P., Mattiello, J., and Lebihan, D. (1994). Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Reson. Ser. B 103, 247–254.

Pubmed Abstract | Pubmed Full Text

Bernstein, M. A., King, K. F., and Zhou, X. J. (2004). Handbook of MRI pulse sequences. Amsterdam; Boston: Academic Press.

Bracewell, R. N. (2000). The Fourier Transform and its Applications, 3rd Edn. Boston, MA: McGraw Hill.

Brigham, E. O. (1988). The Fast Fourier Transform and its Applications. Englewood Cliffs, NJ: Prentice Hall.

Budde, M. D., and Annese, J. (2013). Quantifying anisotropy and fiber orientation in human brain histological sections. Front. Integr. Neurosci. 7:3. doi: 10.3389/fnint.2013.00003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Budde, M. D., and Frank, J. A. (2012). Examining brain microstructure using structure tensor analysis of histological sections. Neuroimage 63, 1–10.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Budde, M. D., Kim, J. H., Liang, H.-F., Schmidt, R. E., Russell, J. H., Cross, A. H., et al. (2007). Toward accurate diagnosis of white matter pathology using diffusion tensor imaging. Magn. Reson. Med. 57, 688–695.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Callaghan, P. T. (1991). Principles of Nuclear Magnetic Resonance Microscopy. Oxford; New York: Oxford University Press.

Callaghan, P. T., Eccles, C. D., and Xia, Y. (1988). NMR microscopy of dynamic displacements: k-space and q-space imaging. J. Phys. E Sci. Instrum. 21, 820.

Carr, H. Y., and Purcell, E. M. (1954). Effects of diffusion on free precession in nuclear magnetic resonance experiments. Phys. Rev. 94, 630–638.

Ciccarelli, O., Catani, M., Johansen-Berg, H., Clark, C., and Thompson, A. (2008). Diffusion-based tractography in neurological disorders: concepts, applications, and future developments. Lancet Neurol. 7, 715–727.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Conturo, T. E., Lori, N. F., Cull, T. S., Akbudak, E., Snyder, A. Z., Shimony, J. S., et al. (1999). Tracking neuronal fiber pathways in the living human brain. Proc. Natl. Acad. Sci. U.S.A. 96, 10422–10427.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Correia, S., Lee, S. Y., Voorn, T., Tate, D. F., Paul, R. H., Zhang, S., et al. (2008). Quantitative tractography metrics of white matter integrity in diffusion-tensor MRI. Neuroimage 42, 568–581.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cory, D. G., and Garroway, A. N. (1990). Measurement of translational displacement probabilities by NMR: an indicator of compartmentation. Magn. Reson. Med. 14, 435–444.

Pubmed Abstract | Pubmed Full Text

Doob, J. L. (1990). Stochastic Processes. New York, NY: Wiley-Interscience.

Frank, L. R. (2001). Anisotropy in high angular resolution diffusion-weighted MRI. Magn. Reson. Med. 45, 935–1141.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haacke, E. M., Brown, R. W., Thompson, M. R., and Venkatesan, R. (1999). Magnetic Resonance Imaging: Physical Principles and Sequence Design. New York, NY: John Wiley and Sons.

Hahn, E. L. (1950). Spin echoes. Phys. Rev. 80, 580–594.

Hinshaw, W. S., and Lent, A. H. (1983). An introduction to NMR imaging: from the bloch equation to the imaging equation. Proc. IEEE 71, 338–350.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jensen, J., Helpern, J., Ramani, A., Lu, H., and Kaczynski, K. (2005). Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging. Magn. Reson. Med. 53, 1432–1440.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kroenke, C. D., Bretthorst, G. L, Inder, T. E., and Neil, J. J. (2005). Diffusion MR imaging characteristics of the developing primate brain. Neuroimage 25, 1205–1213.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Landman, B. A., Bogovic, J. A., Wan, H., ElShahaby, F. E. Z., Bazin, P.-L., and Prince, J. L. (2012). Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri. Neuroimage 59, 2175–2186.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Landman, B. A., Farrell, J. A. D., Smith, S. A., Reich, D. S., Calabresi, P. A., and van Zijl, P. C. M. (2010). Complex geometric models of diffusion and relaxation in healthy and damaged white matter. NMR Biomed. 23, 152–162.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, C., Bammer, R., Acar, B., and Moseley, M. E. (2004). Characterizing non-gaussian diffusion by using generalized diffusion tensors. Magn. Reson. Med. 51, 924–937.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, C., Murphy, N., and Li, W. (2013). Probing white-matter microstructure with higher-order diffusion tensors and susceptibility tensor mri. Front. Integr. Neurosci. 7:11. doi: 10.3389/fnint.2013.00011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mac Donald, C. L., Johnson, A. M., Cooper, D., Nelson, E. C., Werner, N. J., Shimony, J. S., et al. (2011). Detection of blast-related traumatic brain injury in U.S. military personnel. N. Eng. J. Med. 364, 2091–2100.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mattiello, J., Basser, P., and Lebihan, D. (1994). Analytical expressions for the b matrix in NMR diffusion imaging and spectroscopy. J. Magn. Reson. Ser. A 108, 131–141.

McCall, D. W., Douglass, D. C., and Anderson, E. W. (1963). Self-diffusion studies by means of nuclear magnetic resonance spin-echo techniques. Berichte der Bunsengesellschaft für Physikalische Chemie 67, 336–340.

McKinstry, R. C., Mathur, A., Miller, J. H., Özcan, A., Snyder, A. Z., Schefft, G. L., et al. (2002). Radial organization of developing preterm human cerebral cortex revealed by non-invasive water diffusion anisotropy mri. Cereb. Cortex 12, 1237–1243.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mielke, M., Kozauer, N., Chan, K., George, M., Toroney, J., Zerrate, M., et al. (2009). Regionally-specific diffusion tensor imaging in mild cognitive impairment and Alzheimer's disease. Neuroimage 46, 47–55.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mitra, P. P., and Halperin, B. I. (1995). Effects of finite gradient-pulse widths in pulsed-field-gradient diffusion measurements. J. Magn. Reson. Ser. A 113, 94–101.

Mori, S., and van Zijl, P. C. M. (2002). Fiber tracking: principles and strategies – a technical review. NMR Biomed. 15, 468–480.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Moseley, M. E., Cohen, Y., Mintorovitch, J., Chileuitt, L., Shimizu, H., Kucharczyk, J., et al. (1990). Early detection of regional cerebral ischemia in cats: comparison of diffusion- and T2-weighted MRI and spectroscopy. Magn. Reson. Med. 14, 330–346.

Pubmed Abstract | Pubmed Full Text

Özarslan, E. (2009). Compartment shape anisotropy (CSA) revealed by double pulsed field gradient MR. J. Magn. Reson. 199, 56–67.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., Koay, C. G., and Basser, P. J. (2009). Remarks on q-space MR propagator in partially restricted, axially-symmetric, and isotropic environments. Magn. Reson. Imaging 27, 834–844.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., and Mareci, T. H. (2003). Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging. Magn. Reson. Med. 50, 955–965.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., Nevo, U., and Basser, P. J. (2008). Anisotropy induced by macroscopic boundaries: surface-normal mapping using diffusion-weighted imaging. Biophys. J. 94, 2809–2818.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., Shepherd, T. M., Vemuri, B. C., Blackband, S. J., and Mareci, T. H. (2006). Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT). Neuroimage 31, 1086–1103.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2005). (Mathematical) necessary conditions for the selection of gradient vectors in DTI. J. Magn. Reson. 172, 238–241.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2009). “Theoretical and experimental analysis of imaging gradients in DTI,” in Proceedings of the 31st Annual International Conference of the IEEE EMB Society (Minneapolis, MN), 2703–2706.

Özcan, A. (2010a). Characterization of imaging gradients in diffusion tensor imaging. J. Magn. Reson. 207, 24–33.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2010b). “A new model for diffusion weighted MRI: complete fourier direct MRI,” in Proceedings of the 32st Annual International Conference of the IEEE EMB Society (Buenos Aires), 2710–2713.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2011). “Comparison of the complete fourier direct MRI with existing diffusion weighted MRI methods,” in Proceedings of the 2011 IEEE International Symposium on Biomedical Imaging (Chicago, IL), 931–934.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2012). “Complete fourier direct magnetic resonance imaging (CFD–MRI) equations for diffusion MRI signal,” in Proceedings of the 15th International Conference on Medical Image Computing and Computer Assisted Intervention (Nice).

Özcan, A., Quirk, J., Wang, Y., Wang, Q., Sun, P., Spees, W., et al. (2011). “The validation of complete fourier direct MR method for diffusion MRI via biological and numerical phantoms,” in Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE (Boston, MA), 3756–3759.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A., Wong, K. H., Larson-Prior, L., Cho, Z.-H., and Mun, S. K. (2012). Background and mathematical analysis of diffusion mri methods. Int. J. Imaging Syst. Technol. 22, 44–52.

Papadakis, N. A., Xing, D., Huang, C. L.-H., Hall, L. D., and Carpenter, T. A. (1999). A comparative study of acquisition schemes for diffusion tensor imaging using MRI. J. Magn. Reson. 137, 67–82.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rugh, W. J. (1995). Linear System Theory, 2nd Edn. New York, NY: Prentice Hall.

Seal, M. L., Yücel, M., Fornito, A., Wood, S. J., Harrison, B. J., Walterfang, M., et al. (2008). Abnormal white matter microstructure in schizophrenia: a voxelwise analysis of axial and radial diffusivity. Schizophr. Res. 101, 106–110.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shiryaev, A. N. (1995). Probability, 2nd Edn. New York, NY: Springer Verlag.

Song, S.-K., Qu, Z., Garabedian, E. M., Gordon, J. I., Milbrandt, J., and Ackerman, J. J. H. (2002). Improved magnetic resonance imaging detection of prostate cancer in a transgenic mouse model. Cancer Res. 62, 1555–1558.

Pubmed Abstract | Pubmed Full Text

Song, S.-K., Sun, S.-W., Ramsbottom, M. J., Chang, C., Russell, J., and Cross, A. H. (2002). Dysmyelination revealed through MRI as increased radial (but unchanged axial) diffusion of water. Neuroimage 17, 1429–1436.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stejskal, E. O., and Tanner, J. (1965). Spin diffusion measurements: spin echoes in the presence of a time-dependent field. J. Chem. Phys. 42, 288–292.

Torrey, H. C. (1956). Bloch equations with diffusion terms. Phys. Rev. 104, 563–565.

Tuch, D. S. (2004). Q-ball imaging. Magn. Reson. Med. 52, 1358–1372.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turkbey, B., Albert, P. S., Kurdziel, K., and Choyke, P. L. (2009). Imaging localized prostate cancer: current approaches and new developments. Am. J. Roentgenol. 192, 1471–1480.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turkbey, B., and Choyke, P. L. (2012). Multiparametric mri and prostate cancer diagnosis and risk stratification. Curr. Opin. Urol. 22, 310–315.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Voineskos, A. N., Lobaugh, N. J., Bouix, S., Rajji, T. K., Miranda, D., Kennedy, J. L., et al. (2010). Diffusion tensor tractography findings in schizophrenia across the adult lifespan. Brain 133, 1494–1504.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wedeen, V. J., Hagmann, P., Tseng, W.-Y. I., Reese, T. G., and Weisskoff, R. M. (2005). Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magn. Reson. Med. 54, 1377–1385.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xu, J., Humphrey, P. A., Kibel, A. S., Snyder, A. Z., Narra, V. R., Ackerman, J. J., et al. (2009). Magnetic resonance diffusion characteristics of histologically defined prostate cancer in humans. Magn. Reson. Med. 61, 842–850.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Glossary

w_i: ith magnetic moment's displacement from its initial position x_i(t₀).

W^d: The difference of magnetic moment displacement integrals during the periods diffusion (d) gradients are turned on. Its distribution is the main quantity of interest providing information about microstructure. Joined with the other displacement integrals during acquisition and rewind times, W^acq, W^rw, respectively, it forms the total displacement integral: $W=(W^{\text{d}},W^{\text{rw}},W_1^{\text{acq}})\in {ℝ}^7$.

P_cfd(x, W^d): CFD joint distribution function of the number of magnetic moments with initial position x possessing displacement integral values of W^d. It is obtained by marginalizing P^total_cfd(x, W).

S_cfd(k_mr, k_D, k_rw): DW–MRI signal that comes out of the scanner which is equal to the Fourier transform of P^total_cfd. Theoretically it must be Hermitian symmetric since P^total_cfd is real valued.

(k_mr, k_D): k_cfd-space variables corresponding to imaging gradients (mr) and diffusion (D) gradients, respectively.

k_rw: Rewinding (rw) frequency variable affine dependent on k_mr. It cannot be sampled at 0 because of its dependence. For P^total_cfd marginalization, phase corrections are applied to estimate S_cfd at k_rw = 0.

I_k^complexmr(x, k_D): Complex valued images obtained by taking the Fourier transform of S_cfd with respect to imaging frequency k_mr.

I_kmr(x, k_D): Image domain phase corrected I_k^complexmr(x, k_D).

Appendices

A. Second order Statistics for Displacement Integrals of Self Diffusion

In this section, the covariance of displacement integrals is derived for the special case of particles executing Brownian motion. The following preliminaries are necessary to calculate the second order statistics for the displacement integrals of Equation (9). It is assumed for ease of notation below that t_m ≤ t_{m + 1}.

The Wiener process (Shiryaev, 1995), w, which describes the diffusion in an isotropic homogenous medium, is sample path continuous, has w(0) = 0 and independent increments (i.e., Markovian property) with a normal distribution, meaning $w(t)-w(s)~\mathbb{N}(0,(t-s)D)$ where t > s ≥ 0 and D > 0. Using these properties, the covariance of the displacement is:

$\begin{array}{cc}\begin{array}{l}E\left[w(t)w^T(s)\right]=E\left[(w(t)-w(s)+w(s))w^T(s)\right]\\ =E\left[w(s)w^T(s)\right]=\min (t,s)D.\end{array}& (29)\end{array}$

It is straightforward to show that the mean for the displacement integral processes is 0. The covariance in the case of a single interval is calculated as

$\begin{array}{cc}\begin{array}{c}E\left[\left({\displaystyle \underset{t_1}{\overset{t_2}{\int }}{w}_i}(\tau )d\tau \right)\left({\displaystyle \underset{t_1}{\overset{t_2}{\int }}{w}_i^T}(s)ds\right)\right]={\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}E}}\left[w_i(\tau )w_i^T(s)\right]d\tau ds\\ ={\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\min }}(\tau ,s)Dd\tau ds=\frac{1}{3}{\left(t_2-t_1\right)}^2(2t_1+t_2)D,\end{array}& (30)\end{array}$

and for a non-overlapping interval:

$\begin{array}{cc}\begin{array}{l}{\displaystyle \underset{t_3}{\overset{t_4}{\int }}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}E\left[w_i(\tau )w_i^T(s)\right]d\tau ds}}={\displaystyle \underset{t_3}{\overset{t_4}{\int }}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\min (\tau ,s)Dd\tau ds}}\\ =\frac{1}{2}(t_4-t_3)(t_2^2-t_1^2)D.\end{array}& (31)\end{array}$

Finally, using the formulas above, the time points of Figure 1, [t₁ t₂ t₃ t₄] = [t_d1 t_d2 t_d3 t_d4], and these substitutions,

$\begin{array}{l}\delta =t_2-t_1=t_{d2}-t_{d1}=t_4-t_3=t_{d4}-t_{d3},\\ \Delta =t_3-t_1=t_{d4}-t_{d2}=t_4-t_2=t_{d3}-t_{d1},\end{array}$

the following is obtained for W^d_i:

$\begin{array}{cc}\begin{array}{l}E\left[W_i^{\text{d}}{\left(W_i^{\text{d}}\right)}^T\right]=[({(t_2-t_1)}^2(2t_1+t_2)+{(t_4-t_3)}^2\\ (2t_3+t_4))/3-(t_4-t_3)(t_2^2-t_1^2)]D\\ ={\delta }^2(\Delta -\delta /3)D.\end{array}& (32)\end{array}$

The scalar factor, δ² (Δ − δ/3), defined as b_t-value in Özcan (2009, 2010a), is a function of time directly related to the measured quantity, the displacement integrals. It is completely detached from the measurement parameters such as diffusion gradient strength, in contrast to the b-value used in the literature, b = γ² ‖G_D‖²₂ b_t. The dependence of the derivations on Markovian property restricts the validity of the calculations to isotropic samples where the diffusion is characterized by the diffusion coefficient, D, in the distribution of magnetic moment displacement w. For an isotropic homogenous medium, D is estimated from the displacement integral's (W^d) covariance b_t D = δ² (Δ − δ/3)D using DW–MRI acquisitions.

The rigorous treatment of the theory of the stochastic processes in Doob, (1990, p. 62) demonstrates that the sample paths of stochastic processes are Lebesgue integrable under certain conditions and these integrals are well defined random variables. Although a rigorous mathematical analysis, which proves that these conditions are met, is out of the scope of this manuscript, it is important to note that Equation (23) does not prove that the displacement integrals are Gaussian random variables. In this case, the central limit theorem does not apply because the displacement integrals of Equations (8, 9), are not the sums of independent variables since in the integral approximating sum

$\begin{array}{cc}\begin{array}{l}{\displaystyle {\int }_{t_0}^{t_1}{w}_i(\tau )d\tau }\approx \underset{N\to \infty }{\lim }{\displaystyle \sum _{k=1}^N{w}_i(t_0+kd\tau )d\tau }\\ (N\to \infty ,d\tau \to 0,Nd\tau =t_1-t_0),\end{array}& (33)\end{array}$

w_i(τ), rather than independent increments, w_i(t) − w_i(s), is present. A variant of central limit theorem for sums of dependent variables satisfying certain conditions (Shiryaev, 1995, p. 541) might provide the theoretical validation of this highly theoretical open problem.

B. Derivation of the Phase Equation

B.1. Bloch equations and their solutions

The time evolution of the magnetization, ${\overline{m}}_i(t)\in {ℝ}^3$, that generates the MR signal is modeled by the Bloch equations:

$\begin{array}{cc}\begin{array}{l}\frac{d{\overline{m}}_i(t)}{dt}=\gamma {\overline{m}}_i(t)\times B(x_i,t)-\left[\begin{array}{ccc}\frac{1}{T_{2i}}& 0& 0\\ 0& \frac{1}{T_{2i}}& 0\\ 0& 0& \frac{1}{T_{1i}}\end{array}\right]{\overline{m}}_i(t)\\ +\frac{1}{T_{1i}}\left[\begin{array}{c}0\\ 0\\ m_{3i}\end{array}\right]\end{array}& (34)\end{array}$

with T_1i and T_2i denoting spin–lattice and spin–spin relaxation constants, respectively. In the MR scanner the magnetic field as a function of time and position is given by

$\begin{array}{cc}B(x,t)=\left[\begin{array}{c}0\\ 0\\ B_0\end{array}\right]+B_1^m(t)+\left[\begin{array}{c}0\\ 0\\ G(x,t)·x(t)\end{array}\right]& (35)\end{array}$

where B₀ is the static (strong) magnetic field, B^m₁(t) is the radio frequency (RF) pulse modulated at the precession frequency of B₀ and $G(x,t)\in {ℝ}^3$ describes the magnetic field gradients. The expression of the magnetic field simplifies further in the rotating frame:

$\begin{array}{cc}B(x,t)=B_1(t)+\left[\begin{array}{c}0\\ 0\\ G\left(x,\text{ t}\right)\cdot x\left(t\right)\end{array}\right].& (36)\end{array}$

When the magnetic field gradients are turned on without the RF pulse the system matrix becomes

$\begin{array}{cc}{\overline{m}}_i\times \left[\begin{array}{c}0\\ 0\\ \omega (x_i(t))\end{array}\right]=\left[\begin{array}{ccc}0& \omega (x_i(t))& 0\\ -\omega (x_i(t))& 0& 0\\ 0& 0& 0\end{array}\right]{\overline{m}}_i=A(t){\overline{m}}_i.& (37)\end{array}$

with ω(x_i(t)) = G(x_i, t) · x_i(t). As the time-dependent system matrix A(t) commutes with its integral ∫ A(τ)dτ, when the relaxation terms are negligible the solution of Equation (34) is obtained [see Rugh, (1995, pg. 59)] as:

$\begin{array}{cc}{\overline{m}}_i(t)=\exp \left(\gamma {\displaystyle {\int }_{t_0}^{t}A(\tau )d\tau }\right){\overline{m}}_i(t_0)& (38)\end{array}$

by calculating the matrix exponential using Ω(x_i(t), t₀) = ∫^t_t0 ω(x_i(τ)) dτ = ∫^t_t0 G(x_i,τ)· x_i(τ) dτ:

$\exp \left(\text{}\gamma {\displaystyle {\int }_{t_0}^{t}A(\tau )d\tau }\text{}\right)=\left[\begin{array}{ccc}\cos (\gamma \Omega (x_i(t),t_0))& \sin (\gamma \Omega (x_i(t),t_0))& 0\\ -\sin (\gamma \Omega (x_i(t),t_0))& \cos (\gamma \Omega (x_i(t),t_0))& 0\\ 0& 0& 1\end{array}\right].$

Equation (38) defines a system with a rotating transverse magnetization component that can be written in a complex form . The corresponding first order differential equation is

with the solution

For the portion of the pulse sequence involving the RF pulses, in general (Bernstein et al., 2004) the magnetic moments are assumed to be immobile [i.e., w_i(t) ≡ 0 in Equation (4)] thereby defining a constant system matrix that has a standard matrix exponential that solves the differential equation. Herein, however, the effect of the displacements specifically during the π-pulse results in a time varying system matrix. Without loss of generality, the π-pulse is applied in the direction⁸. [1, 0, 0]^T as in Hinshaw and Lent (1983), along with the slice select gradient, G_π = [0, 0, g_π3], active only during the RF pulse. The RF pulse modulated at the same frequency corresponding to the center of the slice, B₀ + g_π3 x_c3, becomes a constant vector, B₁(t) = [B_1a, 0, 0]^T, in Equation (36). The system of equations in this rotating frame is

$\begin{array}{cc}{m}_i\times \left[\begin{array}{c}{B}_{1a}\\ 0\\ {\omega }_{\pi }(x_i(t))\end{array}\right]=\left[\begin{array}{ccc}0& {\omega }_{\pi }(x_i(t))& 0\\ -{\omega }_{\pi }(x_i(t))& 0& B_{1a}\\ 0& -B_{1a}& 0\end{array}\right]m_i=A(t)m_i& (41)\end{array}$

with⁹ ω_π(x_i(t)) = G_π · (x_i(t₀) + w_i(t) − x_c) = g_π3 (x_3i(t₀) + w_3i(t) − x_c3). The procedure would be to solve Equation (41) and then express the solution in the initial rotating frame to obtain a phase component, Φ_πi, that would be a function of the displacement integral and the initial position from the slice center.

However, the new system matrix, A(t), does not commute with its integral preventing the calculation of a matrix exponential formulating a convenient analytical solution as in Equation (40). With the rigorous analysis describing the exact manifestation of the RF pulse left to future studies, herein due to the large magnitude of the RF field in comparison to the magnetic field gradient, its effect will be approximated as a sign reversal of the phase. In this work, the phase change induced by the magnetic field gradient is considered as a factor that is automatically corrected by the phase correction algorithm of section 2.2.

On the other hand, the effect of the displacement induced by the RF slice select gradient (SS–EX in Figure 1) during the excitation period is systematic. As both the RF pulse and the magnetic field gradient remain fixed, statistically speaking, the phase of the total transverse magnetization is affected by the displacements in the same amount at each acquisition prior to the tipping. The position and displacement encoding occurs thereafter during the remainder of the pulse sequence (see also footnote 1 for the phase correction along the slice select direction within the slab) leaving the π/2-pulse out of the derivation of the phase equations below.

B2. Phase equation

For clarity of exposition, the calculations are carried out with the assumption of ideal gradient amplifiers creating rectangular shaped gradient pulses in the time domain and perfect linearity in physical space. In short, during the time interval when the amplifiers are powered on $G(x,t)=G\in {ℝ}^3$.

The derivation of Equations (5– 7) provided in Özcan (2012) is briefly presented in this section for reference purposes. Using the definitions of the variables in Figure 1, the evolution of the phase is described as:

1. First, imaging gradients for read out rewinding, $G_{\text{ro}}\in {ℝ}^3$, phase encoding, $G_{\text{pe}}\in {ℝ}^3$ and slice select $G_{\text{ss}}\in {ℝ}^3$ on the time interval [t₀, t_rw] are turned on. After these gradients are applied the phase and the magnetization become

$\begin{array}{cc}\begin{array}{l}{\Omega }_{\text{rw}}={\displaystyle \underset{t_0}{\overset{t_{\text{rw}}}{\int }}(G_{\text{ro}}(\tau )+G_{\text{pe}}(\tau )+G_{\text{ss}}(\tau ))·x_i(\tau )d\tau }\\ =(G_{\text{ro}}+G_{\text{pe}}+G_{\text{ss}})·\left(\Delta t_{\text{rw}}{x}_i(t_0)+{\displaystyle \underset{t_0}{\overset{t_{\text{rw}}}{\int }}{w}_i(\tau )d\tau }\right)\\ =\Delta t_{\text{rw}}(G_{\text{ro}}+G_{\text{pe}}+G_{\text{ss}})·x_i(t_0)\\ +(G_{\text{ro}}+G_{\text{pe}}+G_{\text{ss}})·{\displaystyle \underset{t_0}{\overset{t_{\text{rw}}}{\int }}{w}_i(\tau )d\tau }.\end{array}& (42)\end{array}$

2. The RF π-pulse between the diffusion gradient pulses, $G_{\text{D}}\in {ℝ}^3$, and the accompanying slice select gradient provide theoretical sign reversal of the phase. The slice select gradient's encoding of magnetic moment motion into the signal is expressed by Φ_πi (see Appendix B1). Since x_i(t_dk) = x_i(t₀) + w_i(t_dk), k = 1, …, 4; at t = t_d4

$\begin{array}{l}\begin{array}{cc}{\Omega }_{\text{D}}={\Phi }_{\pi i}+{\displaystyle \underset{t_{d3}}{\overset{t_{d4}}{\int }}{G}_{\text{D}}·x_i(\tau )d\tau }-{\displaystyle \underset{t_{d1}}{\overset{t_{d2}}{\int }}{G}_{\text{D}}·x_i(\tau )d\tau }& (44)\end{array}\\ ={\Phi }_{\pi i}+G_{\text{D}}·\left((t_{d4}-t_{d3})x_i(t_0)+{\displaystyle \underset{t_{d3}}{\overset{t_{d4}}{\int }}{w}_i(\tau )d\tau }\right)\\ -G_{\text{D}}·\left((t_{d2}-t_{d1})x_i(t_0)+{\displaystyle \underset{t_{d1}}{\overset{t_{d2}}{\int }}{w}_i(\tau )d\tau }\right)\\ ={\Phi }_{\pi i}+G_{\text{D}}·\left({\displaystyle \underset{t_{d3}}{\overset{t_{d4}}{\int }}{w}_i(\tau )d\tau }-{\displaystyle \underset{t_{d1}}{\overset{t_{d2}}{\int }}{w}_i(\tau )d\tau }\right)\\ \begin{array}{cc}+((t_{d4}-t_{d3})-(t_{d2}-t_{d1}))G_{\text{D}}·x_i(t_0).& (45)\end{array}\end{array}$

If the diffusion gradient times are equal, i.e., t_d4 − t_d3 = t_d2 − t_d1 = δ, then the last term in Equation (6) is equal to zero, erasing the influence of initial position from motion encoding part of the signal. This is the insight gained by using the formulation with displacement integrals, ∫ w_i(τ) dτ, rather than the center of mass (COM) of random walk, ∫ x(τ) dτ, introduced in Mitra and Halperin (1995).

The sign reversal affects also previously accumulated phase:

3. The last part of the sequence is where the data are collected:

$\begin{array}{cc}\begin{array}{l}{\Omega }_{\text{ro}}(t)={\displaystyle \underset{t_{\text{acq}}}{\overset{t}{\int }}{G}_{\text{ro}}·(x_i(t_0)+w_i(\tau ))d\tau }=(t-t_{\text{acq}})G_{\text{ro}}\\ ·x_i(t_0)+{\displaystyle \underset{t_{\text{acq}}}{\overset{t}{\int }}{G}_{\text{ro}}·w_i(\tau )d\tau },\end{array}& (47)\end{array}$

leading to

and Equations (5– 7).

C. Total number of Particles from CFD-MRI

The fundamental Fourier relationship of Equation (16), S_cfd(k_mr, k_D, k_rw) = ℱ{P^total_cfd}(k_mr, k_D, k_rw) establishes the relationship between the standard MR image space and the higher dimensional CFD space. By Equation (16)

and by Equation (18), I_k^complexmr(x, k_D, k_rw) = ℱ⁻¹_kmr{S_cfd(k_mr, k_D,k_rw)}, the image obtained without the diffusion encoding magnetic field gradients is

$\begin{array}{cc}\begin{array}{c}={\displaystyle \int P_{\text{cfd}}^{\text{total}}}\left(x^{\prime },W\right)\delta \left(\frac{\gamma }{2\pi }(x^{\prime }-x)\right)d{x}^{\prime }dW\\ =\frac{2\pi }{\gamma }{\displaystyle \int P_{\text{cfd}}^{\text{total}}}(x,W)dW\\ ~\text{total number of particles at }x.\end{array}& \begin{array}{c}(52)\\ (53)\\ (54)\end{array}\end{array}$

Footnotes

1. ^The refocusing lobe of the slice select gradient after the RF pulse denoted by SS in Figure 1 adjusts ϕ_slice to be constant throughout the slice in the slice select direction.
2. ^The unit for k_mr is $\frac{\text{magnetic field}\times \text{time}}{\text{distance}}$ but for k_D and k_rw it is $\frac{\text{magnetic field}}{\text{distance}}$. They are both in accordance with the units of the position and the displacement integrals in Equation (13). After k_cfd is multiplied by the gyromagnetic ratio γ, with unit (magnetic field × time)⁻¹, the product becomes unitless in Equation (14).
3. ^Another approach in the literature is to present the value of the function over a sphere (Tuch, 2004; Wedeen et al., 2005). However, uniqueness is lost in this representation as demonstrated by these functions: f(x, y, z) = x² + y² + z², h(x, y, z) = (f(x, y, z))², which have the same value on the unit sphere but different isosurfaces.
4. ^Also named nuclear spin self-correlation function (Callaghan et al., 1988) and translation probability (Cory and Garroway, 1990).
5. ^One exception is the basic case of isotropic diffusion analyzed in Appendix A where the diffusion coefficient that describes the Gaussian propagator can be recovered using the adjustment factor for the displacement integral covariance namely b_t in Equation (24) and the b-value in DTI formulation.
6. ^Acquisition time is shortened with reduced sampling schemes aimed at specific goals, e.g., compressed sensing for tractography (Landman et al., 2012) at the expense of some information loss.
7. ^Instead of pixel by pixel corrections that would completely eliminate the imaginary part in Figure 3.
8. ^The choice of direction does not matter (Haacke et al., 1999, p. 387).
9. ^The absence of a rewinding slice select gradient, in contrast to π/2-pulse in Figure 1, requires the addition of the slice center, x_c, in the formulation as there will be a phase shift along the slice select direction in the slab [see Equation (11) and footnote 1].