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.


INTRODUCTION
Since the conception of mathematical models for the effect of the magnetic moment diffusion in nuclear magnetic resonance (NMR) experiments by Hahn (1950), Carr and 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 noninvasive, 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 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  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 reexamining 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).

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): (1) 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 0 ) denoting the initial magnetization tipped to the transverse plane, describes the phase (when multiplied by γ) as a function of the magnetic field gradients G(x, t) ∈ R 3 , and the position of the magnetic moment is x i ∈ R 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 where w i (t) ∈ R 3 represents the displacement of the magnetic moment from its initial position, 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 * ∈ R 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 0 , is the end time SS−PI t 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.
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): The second term in Equation (6) removes the injection of the initial position into the DW signal because of equal pulse duration times, 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) : 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, k mr ∈ R 2 : with the additional term, which is constant because the slice select axis component of x i (t 0 ) is the slice position 1 . 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 mr1 /g ro1 + t acq + t rw and g pe2 = −k mr2 / t 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 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. affine function of k mr : Consequently, the phase, i , of Equation (3) is expressed concisely by defining k D = G D : reflecting the effect of initial position and displacement integrals on the phase 2 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 0 ) = m 0 , 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 positiondisplacement properties and then to count the numbers elements in the groups. The signal in Equation (14) is calculated by integrating over the whole position-displacement space (absorbing m 0 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 1 . Among the elements of W, the focus is on the most descriptive MRI observable, W d . Its distribution is obtained by marginalizing 2 The unit for k mr is magnetic field×time distance but for k D and k rw it is magnetic field 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 (magneticfield × time) −1 , the product becomes unitless in Equation (14).
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 2 1 + k 2 2 ) , with the "undesirable" variable k 2 sampled on a line k 2 = a k 1 + b : The result is complex valued in comparison to the real valued Fourier transform of exp(−k 2 1 ) which can be obtained by setting a = b = 0.

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) F (Hermitian symmetric function).
Accordingly CFD-MRI reconstruction is based on the following: Since by definition P total cfd is real valued, S cfd is Hermitian symmetric.
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 reestablishment 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 phaseencode (e.g., due to sample shaking) and read-out directions are corrected by first determining from the Fourier transform in k mr , the angle ( I complex k mr ) from the signal regions along the center lines of each physical direction at k D = 0, The phase corrections are then applied systematically at each value of k D (see Figure 3): (20) The Fourier transform in the remaining variables, is evaluated sequentially in each k D -dimension with the aim of reestablishing the Hermitian property, , 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 , respectively for l = 1, 2, 3), e.g., I k mr (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 k mr (x c , 0) remains unchanged. For example, for the first direction, at selected values of k D1 the polynomial looks like 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 k mr to be equal to ) for all k D2 , k D3 and x. 6. 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 phase shift around 0 frequency, indicative of coherent motion. After the phase corrections obtained using the polynomial 0.266 k 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.
without diffusion gradients, I k mr (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.
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 : (see Appendix C). CFD phase correction algorithm outperformed the fitting of the phase values up to the fourth degree multinomials in R 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.

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) (see Özcan, 2009, 2010a) because the Fourier transform of a Gaussian with a covariance matrixD is also a Gaussian with covarianceD −1 : 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.

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 2 and 0.5 mm slice thickness. The k Dspace 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 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.

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 3 of normalized P cfd , 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), 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 2 + y 2 + z 2 , h(x, y, z) = (f (x, y, z)) 2 , which have the same value on the unit sphere but different isosurfaces.  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 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).
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 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 cvalues. 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 R 3 to grasp the complete information offered by CFD-MRI.
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).
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.

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 4 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) The inverse problem of obtaining the propagator from the distribution of the displacement integrals is singular because of FIGURE 6 | The isosurfaces (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 P cfd = 0.14 and the green surface isP cfd = 0.21. The red surface envelops the green one.

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.
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 5 . 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)  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, ) δ, 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.
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 where I complex k mr 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 t 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 k D -space, without recourse to a Fourier transform because of its direct appearance in Equation (28)  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 et al. (2012) for the correspondence with the B-matrix formulation of Basser et al. (1994)].
The magnitude-based Fourier relationship presented in qspace 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 Callaghan (1991, Chap. 8)] 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.

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 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 (45)] 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 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 (P 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.

Frontiers in Integrative Neuroscience
www.frontiersin.org April 2013 | Volume 7 | Article 18 | 11 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 6 , 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 7 , more effective biomarker construction would 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.
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).
Finally, using the formulas above, the time points of Figure 1, [t 1 t 2 t 3 t 4 ] = [t d1 t d2 t d3 t d4 ], and these substitutions, the following is obtained for W d i : The scalar factor, δ 2 ( − δ/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 = γ 2 G D 2 2 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 = δ 2 ( − δ/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 (32) 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 Frontiers in Integrative Neuroscience www.frontiersin.org April 2013 | Volume 7 | Article 18 | 14 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.