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

The foundation for an accurate and unifying Fourier-based theory of diffusion weighted magnetic resonance imaging (DW-MRI) is constructed by carefully re-examining the first principles of DW-MRI signal formation and deriving its mathematical model from scratch. The derivations are specifically obtained for DW-MRI signal by including all of its elements (e.g., imaging gradients) using complex values. Particle methods are utilized in contrast to conventional partial differential equations approach. The signal is shown to be the Fourier transform of the joint distribution of number of the magnetic moments (at a given location at the initial time) and magnetic moment displacement integrals. In effect, the k-space is augmented by three more dimensions, corresponding to the frequency variables dual to displacement integral vectors. The joint distribution function is recovered by applying the Fourier transform to the complete high-dimensional data set. In the process, to obtain a physically meaningful real valued distribution function, phase corrections are applied for the re-establishment of Hermitian symmetry in the signal. Consequently, the method is fully unconstrained and directly presents the distribution of displacement integrals without any assumptions such as symmetry or Markovian property. The joint distribution function is visualized with isosurfaces, which describe the displacement integrals, overlaid on the distribution map of the number of magnetic moments with low mobility. The model provides an accurate description of the molecular motion measurements via DW-MRI. The improvement of the characterization of tissue microstructure leads to a better localization, detection and assessment of biological properties such as white matter integrity. The results are demonstrated on the experimental data obtained from an ex vivo baboon brain.

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 lmr-dimensional MRI slice, the DW–MRI complex valued signal that comes out of the scanner is the (lmr + 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 lmr 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 (lmr + 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):

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, is written in complex number form with m(t0) denoting the initial magnetization tipped to the transverse plane, describes the phase (when multiplied by γ) as a function of the magnetic field gradients , and the position of the magnetic moment is . The time-dependent position, x, in Equation (3) affects the phase, Ω, 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 where represents the displacement of the magnetic moment from its initial position, x(t0), [i.e., w(t0) = 0]. The only physical requirement is the continuity of w(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 (lmr = 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 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, t0, 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):

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 − t = t − t. The term Φπ 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(t): (Wd, Wacq, Wrw) corresponding to the displacement integrals for diffusion (d), analog to digital conversion acquisition (acq), and initial rewind (rw) gradient time periods, respectively

First term in Equation (5) is the definition of the k-space in regular MRI, : with the additional term, which is constant because the slice select axis component of x(t0) 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, Gro = [gro1, 0, 0], Gpe = [0, gpe2, 0], Gss = [0, 0, gss3], time and gpe2 become functions of kmr = [kmr1, kmr2, kmr3] using Equation (10):

Accordingly, Wacq(t) becomes a function of kmr1 and the coefficients of Wrw in Equation (7) are written as a vector which is an affine function of kmr:

Consequently, the phase, Ω, of Equation (3) is expressed concisely by defining kD = GD: 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(t0) = m0, 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,

The signal in Equation (14) is calculated by integrating over the whole position-displacement space (absorbing m0 into Ptotalcfd for ease of notation):

Equation (15) is by definition the Fourier transform of Ptotalcfd with non-linearities added by Wacq1.

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

However, the affine dependence of krw in Equation (12) makes it impossible to fix krw = 0 and to sample in (kmr, kD, 0) subspace. The following example demonstrates how the affine dependence affects the measurements by using a two dimensional Gaussian function, exp(−(k21 + k22)), with the “undesirable” variable k2 sampled on a line k2 = a k1 + b:

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

2.2. Phase corrections for the estimation of Pcfd

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:

Since by definition Ptotalcfd is real valued, Scfd is Hermitian symmetric.

Furthermore, an immediate implication of the transform property and the Hermitian symmetry of Scfd 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 Scfd 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 Pcfd, 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 kmr-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 kmr, the angle (∠ I) from the signal regions along the center lines of each physical direction at kD = 0,

The phase corrections are then applied systematically at each value of kD (see 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(x, 0) and on the bottom row, Pcfd(x, 0) are displayed. The imaginary part is approximately 10% of the real part in both cases. On the right, full representation of Pcfd with isosurfaces (, 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).

The Fourier transform in the remaining variables, is evaluated sequentially in each kD-dimension with the aim of re-establishing the Hermitian property, I(x, −kD) = I(x, kD), using the following steps.

CFD Phase correction algorithm:

Pick a pixel at location x, preferably near the center of the image where tissue or a good signal area is located.

Starting from the first direction, l = 1 of the kD space calculate the phase on the line passing through the origin (i.e., [kD1, 0, 0], [0, kD2, 0], [0, 0, kD3], respectively for l = 1, 2, 3), e.g., ∠I(x,(kD1, 0, 0)).

Investigate the plot of the phase versus kD. Pick as many as possible consecutive values of kD near 0 without sudden changes to assure high signal to noise value.

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(x, 0) remains unchanged. For example, for the first direction, at selected values of kD1 the polynomial looks like as demonstrated in Figure 2.

Figure 2

Özcan ( The plots show data acquired at each diffusion gradient value kD1 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.266kD1 estimated from the points kD1 = ‒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.

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, kD1, update I to be equal to for all kD2, kD3 and x.

Repeat steps 2–5 for the remaining directions.

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

Figure 5

On the left, one dimensional graphical representation of the choice of 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).

In Pcfd(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 Pcfd:

(see Appendix C).

CFD phase correction algorithm outperformed the fitting of the phase values up to the fourth degree multinomials in . 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 kmr-space, the infinite bandwidth in kD-space due to Pcfd's finite support in Wd-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 kD 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, Pcfd's bandwidth (i.e., Scfd's support) is shrunk by increasing δ and Δ causing the dispersion (covariance) of the displacement integral Wd (and therefore Pcfd's support) to increase. This is directly visible in the special case of Brownian Motion characterized by the diffusion tensor D in Figure 4. Pcfd and Scfd are zero mean Gaussians with covariances, respectively equal to (see Appendix A)

Figure 4

Effects of fitting the bandwidth within the field of view (FOV) in the Fourier space on the 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.

(see Özcan (2009, 2010a) because the Fourier transform of a Gaussian with a covariance matrix is also a Gaussian with covariance :

The procedure is graphically displayed in Figure 4 also emphasizing the effect of ripples on the small values of Pcfd 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 kD-space to prevent aliasing artifacts on Pcfd. This is constrained by the time available for acquisitions as each point in kD-space requires the scan time of the entire kmr-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 kmr-space was sampled to result in images of 128 × 128 pixels with a FOV 64 × 64 mm2 and 0.5 mm slice thickness. The kD-space was sampled in a uniformly spaced Cartesian grid in a cube [−30, 30 Gauss/cm]3 with 6 Gauss/cm sampling intervals at each dimension resulting in 11× 11 × 11 voxels. The repetition time T = 1 s, echo time T = 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 (lmr = 2), plus three more for displacement integrals] creates a visualization challenge which is addressed herein by using Pcfd(x, 0) as the background image. Furthermore, the isosurface of normalized Pcfd, 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), over all locations. The key point is the choice of an appropriate c-value that will reveal the outskirts of Pcfd corresponding to the small number of “scout” magnetic moments that travel further away thereby portraying the microstructure of the environment. In summary,

Too high values do not provide enough structural information (see first rows in Figure 5).

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).

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 Pcfd which implies a larger support for Scfd 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 ∫ Pcfd(x, Wd) dWd = 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 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 to grasp the complete information offered by CFD–MRI.

Figure 6

The isosurfaces ( 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 and the green surface is . 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 Pcfd 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

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 [displacement in Equation (4)], and the output W [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, kD = GD from these parameters in contrast to the q-space variable (Callaghan, 1991): q = (2π)−1 γ δ GD and the b-value of DTI: b = γ2 ‖ GD ‖22 δ2 (Δ − δ/3) (see Appendix A).

Figure 8

The DW–MRI signal from a signals and systems perspective (Özcan, The input is the magnetic moment displacement w and the output is the displacement integral, W, 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: 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 − t = t − t 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, Wd ≈ (x(t) − x(t)) δ = (w(t) − w(t)) δ, 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 Wd are attained at a time point within the respective integration intervals, [t, t] and [t, t], the nowhere-differentiability of the sample paths (Shiryaev, 1995) implies that the intermediate time points satisfying the equality are not fixed as t and t, 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)|x(t)), 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 Pcfd 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:

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).

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 lmr (imaging) dimensions. Thus, in the analysis of model matching methods the first lmr 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 where I 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: where the calculation of b from the PDE approach in Özcan (2010a) and covariance of displacement integrals in Appendix A resulted in the same value: δ2 (Δ − δ/3). In CFD parlance, by the Fourier relationship between Gaussians in Equation (25), the diffusion quadratic form, D, is estimated in kD-space, without recourse to a Fourier transform because of its direct appearance in Equation (28) rather than its inverse, D−1, in the Gaussian of motion space. The coefficient matrix defined by carefully selected vectors in kD-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 krw on kmr CFD derivations show that the inseparability partially accounts for the non-Hermitian nature of the Scfd. 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

The comparison of 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 () 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:

With Pcfd 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 Pcfd 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 Pcfd 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 Pcfd 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.