The Mathematics of Quasi-Diffusion Magnetic Resonance Imaging

: Quasi-diffusion imaging (QDI) is a novel quantitative diffusion magnetic resonance imaging (dMRI) technique that enables high quality tissue microstructural imaging in a clinically feasible acquisition time. QDI is derived from a special case of the continuous time random walk (CTRW) model of diffusion dynamics and assumes water diffusion is locally Gaussian within tissue microstructure. By assuming a Gaussian scaling relationship between temporal ( 𝛼 ) and spatial ( 𝛽 ) fractional exponents, the dMRI signal attenuation is expressed according to a diffusion coefficient, 𝐷 (in mm 2 s −1 ), and a fractional exponent, 𝛼 . Here we investigate the mathematical properties of the QDI signal and its interpretation within the quasi-diffusion model. Firstly, the QDI equation is derived and its power law behaviour described. Secondly, we derive a probability distribution of underlying Fickian diffusion coefficients via the inverse Laplace transform. We then describe the functional form of the quasi-diffusion propagator, and apply this to dMRI of the human brain to perform mean apparent propagator imaging. QDI is currently unique in tissue microstructural imaging as it provides a simple form for the inverse Laplace transform and diffusion propagator directly from its representation of the dMRI signal. This study shows the potential of QDI as a promising new model-based dMRI technique with significant scope for further development.


Introduction
Quasi-diffusion imaging [1] is a novel diffusion magnetic resonance imaging (dMRI) technique based on a special case of the continuous time random walk (CTRW) diffusion model [2][3][4][5].The diffusion dynamics are represented by an effective normal diffusion where the mean squared-displacement of the diffusing particles is proportional to time, 〈 2 〉~.This corresponds to the fractional exponents representing the probability density functions of the waiting times between steps, , and the step lengths, , having the Gaussian scaling relationship, Consequently,  and  exponents are not independent and both processes are described by related inverse power dependencies for step length ( 2(−1) ) and waiting time ( −1 ).
Diffusion MRI provides a probe of tissue microstructure by application of diffusion encoding gradients that sensitise signals to water diffusion within an image voxel [6][7][8][9][10].In vivo dMRI voxel sizes are typically on a mm scale, whereas tissue microstructural properties are measured in microns (μm).Hence, the dMRI signal includes components from all cells within a voxel, and contributions from the extra-cellular space.When there is restriction of water diffusion at surface boundaries, the pulsed gradient spin echo (PGSE) Nuclear Magnetic Resonance (NMR) diffusion measurement provides a probe of the length scales of anisotropic microstructures [11].
The observed dMRI signal in tissue is non-Gaussian.This has led to the development of numerous techniques to infer attributes of the diffusion environment and to several theories regarding the origin of the non-Gaussian signal [8,9,12,13].These theories relate to the effect on diffusion by cell membranes, their internal structures, and the complex extracellular spaces in tissues being imaged.Early measurements using a Gaussian model of water diffusion identified that the apparent diffusion coefficient (ADC) in tissue was lower than for free water diffusion, and was related to cell density, cell size and water concentration within the tissue.Modern clinical MRI systems have revealed non-Gaussian behaviour of the dMRI signal creating a need for the development and application of more complex diffusion models [8,9,12,13].
In the CTRW diffusion model, the underlying diffusion process can be space-fractional super-diffusion for  = 1, 0 <  ≤ 2 [14][15][16][17], time fractional sub-diffusion for 0 <  ≤ 1 [16][17][18], or governed by the general CTRW model, allowing the characteristics of the underlying diffusion process to be determined [19][20][21][22][23]. Outside the CTRW model, techniques have been developed that assume the underlying diffusion process to be Gaussian.These include representing the diffusion signal attenuation as a second moment expansion known as Diffusional Kurtosis Imaging (DKI) [24,25], or as a signal in the experimentally acquired "q-space" from which, via application of the inverse Fourier transform leads to a Mean apparent diffusion propagator (MAP) of molecular displacement [26][27][28] allowing measurement of the length scale of the diffusion environment.This is of interest in clinical imaging as models of the diffusion process that enable quantitative assessment of tissue compartment or cell size can be used to assess pathological change, with potential applications in aiding diagnosis [29,30], predicting disease outcome [31] and monitoring treatment effects [32].
The Random Permeable Barriers Model (RPBM) is a prominent theory in dMRI based on the time-dependence of diffusion coefficients observed in a voxel [33][34][35][36].This model assumes that diffusion is locally Gaussian and that water can exchange freely between pores within a disordered medium.The RPBM formulation leads to a power-law time dependence for the observed signal and a time-varying diffusion coefficient.A limitation of this technique is that it does not allow derivation of an explicit form in the time domain, or an overall form for the diffusion signal attenuation.
Other techniques consider non-exchanging [37][38][39][40], or exchanging [41,42] diffusion compartments constructed of different mixtures of geometrical shapes assumed a priori to represent tissue microstructure.These formulations do not constitute a model of the underlying diffusion dynamics and rely instead on developing a highly idealised geometrical model of the structure of the diffusion environment.
Here we consider the mathematics of quasi-diffusion and quasi-diffusion imaging.The fundamental solution of the space-time diffusion-wave equation within the CTRW diffusion model has been extensively studied [43][44][45][46][47], as have the properties of the quasidiffusion model [47][48][49][50].In particular, quasi-diffusion is referred to as -fractional diffusion by [48,49] who suggest this special case of the CTRW diffusion model represents a natural fractionalisation of the diffusion process [49].
The quasi-diffusion equation is given by the stretched Mittag-Leffler function [51][52][53].We describe its asymptotic properties and Laplace transform, and consider the properties of the quasi-diffusion propagator.Finally, we show how the quasi-diffusion propagator can be used in quasi-diffusion MAP imaging to calculate sub-micron measurements of cell radii.Examples of quasi-diffusion MAP imaging are given in which we show the technique can be performed using a clinical MR scanner and that dMRI data for analysis of quasi-diffusion can be acquired in a clinically feasible scan time.

Quasi-Diffusion Imaging
Application of the CTRW diffusion model to dMRI was first described by Ingo et al. [19,54].In the conventional case of particles characterised by Brownian motion in a homogeneous and isotropic geometry, their unrestricted self-diffusion is described by the second-order partial differential equation, where the solution in -dimensional space is the Gaussian propagator,   (, ), which defines the probability density of finding a molecule at position, , at time, , as, ) ,  ≥ 0, and  is the diffusion coefficient in mm 2 s −1 .The extension to the CTRW diffusion model is provided by the fractional partial diffusion equation, where    0  is the Caputo fractional derivative (which is the αth fractional order time derivative for 0 <  < 1), and   ||  ⁄ is the Reisz fractional derivative (which is the βth fractional order space derivative for 0 <  < 2) and  , is the effective diffusion coefficient in units of mm β s −α .It is assumed that at time,  = 0, all the material is located at the origin and is given by the Dirac Delta function (, 0) = ().In the special case of quasidiffusion, the Gaussian scaling relationship, where  ,2 is an effective normal diffusion coefficient in units of mm 2α s −α .In the quasidiffusion case, a Gaussian diffusion coefficient in mm 2 s −1 can be recovered from the effective normal diffusion coefficient, The fractional partial differential equation in (5) can be solved using the Laplace-Fourier transform [5].Firstly, the transform   (, ) →   (, ) is given by, then by application of the inverse Laplace transform to (7) we perform the transform   (, ) →   (, ) to obtain the unique characteristic equation, which is a fractional relaxation curve for each wavenumber, , given by, where   () is the one-parameter Mittag-Leffler function which we make use of via its series representation, for  > 0 and  ∈ ℂ.For application to dMRI, we make the change of variable from  to  in (8) to indicate data acquisition by experimentally controlled PGSE parameters in space.Experimental parameters are set to provide diffusion encoding in different orientations in -space where  = 1

2𝜋
(assumed to be in mm −1 ), γ is the gyromagnetic ratio of hydrogen (for quantification of water diffusion), and  is the diffusion encoding gradient strength (in mTm −1 ).The effective diffusion time of the pulse sequence is denoted as (in s) for a given diffusion gradient pulse duration, δ, and separation, .In clinical applications, dMRI are typically acquired by keeping ∆ ̅ constant while altering  by changing the diffusion encoding gradient strength, .
Mathematical analysis of diffusion imaging in -space makes the assumption that dMRI data are acquired using an infinitely short diffusion gradient pulse duration, , with  ≪ ∆.This is known as the short pulse approximation.The Fourier relationship between q-space (in mm −1 ) and r-space (in mm) is only exact in this limit.In practice, this assumption is routinely violated on clinical MRI systems due to technical and safety limitations of in vivo MRI, meaning that gradient pulses have finite duration usually in the range 20 ms <  <  < 70 ms.
Equation ( 8) can also be written with respect to an overall diffusion sensitisation,  =  2 ∆ ̅ .This overall sensitivity is due to the cumulative effect of the diffusion encoding gradient and the effective diffusion time, ∆ ̅ .In a radial direction through -space starting at the origin, the characteristic equation for quasi-diffusion is given by, where () is the diffusion signal intensity at a given diffusion-weighting  (in s mm −2 ), and (0) is the signal intensity at  = 0 s mm −2 .Equation ( 10) is a stretched Mittag-Leffler function that describes Gaussian diffusion when  = 1 and an effective normal diffusion for 0 <  < 1.The diffusion coefficient,  1,2 , and fractional exponent, , are independent parameters that together parameterise a family of decay curves according to the rate of diffusion signal decay ( 1,2 ) and the shape of the power law tail, .
A key application of QDI is in imaging of the human body where diffusion of free water at body temperature is   = 3 × 10 −3 mm 2 s −1 and 0.5 <  < 1 in typical healthy brain tissue [1]. Figure 1 illustrates the family of signal decay curves described by (10). Figure 1a shows the quasi-diffusion signal attenuation parameterised by  for an arbitrary diffusion coefficient,  1,2 = 1.5 × 10 −3 mm 2 s −1 for 0.1 ≤  ≤ 0.99 with Figure 1b showing the quasi-diffusion signal attenuation parameterised by .The quasi-diffusion model can be demonstrated to provide good quality fits to experimental data.Figure 2 shows QDI maps of  1,2 and  in the brain of a young (age 28 years) healthy subject for which a comprehensive data acquisition was acquired with 29 b-values over the range 0 ≤  ≤ 5000 s mm −2 in 6 non-collinear gradient encoding directions evenly distributed across the sphere.Experimental diffusion times were  = 23.5 ms and Δ = 43.7 ms giving an effective diffusion time of ∆ ̅ = 35.9ms.Data were acquired on a clinical 3T MR scanner at St George's, University of London (SGUL) with a voxel size of 1.5 mm × 1.5 mm × 5 mm in 35 min 12 s.Full image acquisition parameters are given in Appendix A. Data analysis were performed using the technique described in [1] to estimate  1,2 and  values in each diffusion encoding direction and their mean values within each image voxel.The top row of Figure 2 shows the exceptional quality of fit of the quasi-diffusion model to observed data in individual grey (Figure 2a) and white matter voxels (Figure 2b) across the full range of -factors.The  1,2 and  maps exhibit similar mean apparent diffusion coefficients in brain tissue, with the bright signal pertaining to cerebrospinal fluid-filled (CSF) spaces where diffusion is Gaussian with   = 3 × 10 −3 mm 2 s −1 .The bottom row of Figure 2 shows the quality of data reconstruction from three points,  = {0,1080,5000} s mm −2 within the dMRI data acquisition, which were chosen as their modelled quasi-diffusion signal attenuation closest to that of the full 29 -value dataset across the entire image [55].Overall, Figure 2 highlights how well the quasi-diffusion model fits acquired data and demonstrates that high quality images can be acquired in a clinically feasible time of 120 s without the need for an extensive set of different b-value images to accurately define the signal decay curve.We will now explore some properties of the quasi-diffusion model that are relevant to its application in diffusion magnetic resonance imaging.

General Properties of the Mittag-Leffler Function
This section reviews several general properties of the Mittag-Leffler function.For properties beyond those stated below see [51,[56][57][58].The Mittag−Leffler function ( 9) is an entire function that is convergent in the whole complex plane with a singularity as  → 0 + such that, for || < 1.The Mittag-Leffler function is considered to be a generalisation of the exponential function and simplifies to mono-exponential decay for the case  = 1, as  1 () = exp().

Asymptotic Properties of the Quasi-Diffusion Characteristic Equation
The asymptotic behaviour of the quasi-diffusion characteristic Equation ( 8) as  → 0 + and  → ∞ can be determined from its Mittag-Leffler function series representation (9).In the low -value approximation, the behaviour of the quasi-diffusion equation as  → 0 is derived from the convergent power series representation of ( 9), and indicates that as  → 0, the function is asymptotic to a stretched exponential function.
In the high -value approximation, the behavior of the quasi-diffusion equation as  → ∞ + is derived from the asymptotic power series representation [58], and is given by, with first order approximation, In practice, for in vivo tissue microstructural imaging, the diffusion coefficient of water is finite and unlikely to be less than  1,2 = 1 × 10 −5 mm 2 s −1 and has an upper limit of   for free water at human body temperature.In this case, the behaviour of ( 12) and (15) will be dominated by -value as both  → 0 and  → ∞.
These asymptotic properties show that the quasi-diffusion imaging equation interpolates between a stretched exponential at low -values and a negative power law at high -values.This is consistent with the fast diffusion signal decay at low -values and slow signal attenuation at high -values which are observed in Figure 2. The asymptotic representations of the quasi-diffusion equation in tissue microstructural imaging can be stated as, after substitution of . The behaviour of the quasi-diffusion equation is in contrast to the Mittag-Leffler function ( 9) which tends to an exponential function as  → 0 + .The quasi-diffusion equation is a stretched Mittag-Leffler function that has the same form and behaviour as the fractional relaxation equation which has been extensively studied in application to several physical systems including fractional linear viscoelasticity [53,59,60] and dielectric models [60].The negative power law behaviour of the quasidiffusion characteristic equation is highlighted in Figure 1c over the range 10 1 ≤  ≤ 10 5 for an arbitrary diffusion coefficient,  1,2 = 1.5 × 10 −3 mm 2 s −1 and 0.1 ≤  ≤ 0.99.

The Laplace Transform of the Quasi-Diffusion Characteristic Equation
The Laplace transform and inverse Laplace transform of the quasi-diffusion characteristic equation are of interest in quasi-diffusion imaging.For instance, the inverse Laplace transform of (8) provides a decomposition of the quasi-diffusion signal into a probability density function of Fickian diffusion coefficients, as will be seen later.The mathematics of the Laplace transform of the stretched Mittag-Leffler function have been described in detail [53,[58][59][60].The Laplace transform of the Mittag-Leffler function ( 9) is given by e.g., [56], which is fundamental in the evaluation of the Laplace transform of the fractional relaxation curve.In our case, the form of the fractional relaxation curve is given by gives, As  1,2 is positive, then   (−( 1,2 )  ) is completely monotonic for  > 0 when 0 <  ≤ 1.The inverse Laplace transform can be formulated as a Laplace inverse integral along a Hankel path, Ha, as shown by, e.g., [59,60], where  > 0 and  ∈ ℂ.The Hankel path is the integration around a loop that starts and ends at −∞ and encircles the disc || ≤ | 1,2 | 1/ in the positive sense: − ≤   ≤ , e.g., [53,59,60].The result of this integration is the normalised spectrum of quasi-diffusion signal decay in , which is a non-negative probability density function that is locally integrable in ℝ + given by [53,59,60], and is analogous to the normalised spectrum of relaxation in time [53,59].The spectrum of apparent diffusion coefficients is given when  = / 1,2 such that   (/ 1,2 ) =   () 1,2 [53,59], and represents the quasi-diffusion signal as a distribution of Gaussian diffusion signals,  (in mm 2 s −1 ).This provides further support for quasi-diffusion being a natural generalisation and fractionalisation of diffusion dynamics when the process occurs within complex structures that hinder or restrict free diffusion.The quasi-diffusion spectral function is illustrated in Figure 3a for unit  1,2 = 1 mm 2 s −1 and 0.1 <  < 0.99.At  close to unity (e.g.,  = 0.99), the probability density function tends towards the Dirac delta function with a characteristic ADC equivalent to the resultant mean of an unrestricted Gaussian diffusion process.This indicates that diffusion dynamics within the structure being studied are Gaussian with no boundaries that hinder or restrict diffusing particles.As  decreases (from 0.9 to 0.7), the probability density function broadens to include significant contributions at low ADCs representing characteristics of diffusion in a restricting environment and increasing tissue heterogeneity.For  < 0.65, the probability density function smoothly deforms to a hyperbolic shape that becomes dominated by large contributions from restriction and greater tissue heterogeneity.It should be noted that in practice for dMRI in brain tissue, we are typically in the range 0.5 <  ≤ 1 [1].
Spectra of ADCs are shown for healthy brain tissue in CSF (Figure 3b), grey matter (Figure 3c) and white matter (Figure 3d).The spectra were calculated based on measurements of mean, axial and radial  1,2 and  from Barrick et al. [1].Axial measurements are calculated along the direction parallel to tissue microstructure, with radial measures being perpendicular to the tissue microstructural axis.Figure 3b shows that diffusing water molecules in CSF-filled spaces exhibit isotropic diffusion with an ADC of 3.0 × 10 −3 mm 2 s −1 .The slight apparent anisotropy in measurements is likely related to partial volume effects with brain tissue that are not fully removed in the CSF segmentation, combined with the effects of noise which have broadened the distribution from being a Dirac-delta function.In grey matter (Figure 3c), which contains layers of neurons, the probability density functions represent an almost isotropic diffusion within a heterogeneous tissue microstructure.In white matter (Figure 3d), which consists of axons that provide the wiring of the brain, there is anisotropic diffusion such that diffusion parallel to tissue microstructure is through a more homogeneous diffusion environment (i.e., along axons surrounded by myelin sheaths) than diffusion perpendicular to axons which is more hindered and/or restricted by tissue microstructural boundaries.The white matter spectra indicate considerable anisotropy in the axial and radial ADCs, with large contributions of restricted diffusion to the radial quasi-diffusion coefficient,  1,2 , representing a low .

The Quasi-Diffusion Propagator
Here we describe mathematical results for deriving and understanding the behaviour of the  -dimensional quasi-diffusion propagator.We use notation whereby the Green's function in -dimensional space for the solution of ( 5) is given by   (, ).Equations for the propagator in one-, two-and three-dimensions are presented in integral form, or closed form where possible.The following is given in greater detail by [48].
The inverse Fourier transform of the quasi-diffusion characteristic Equation ( 9) is given by, As this function belongs to the functional space  1 (ℝ  ) with respect to , the inverse Fourier transform can be written as [61], The Mittag-Leffler function ( 9) is axially symmetric in  and is a radial function, so by the formula [62], where   denotes the Bessel function of index , Equation ( 24) can be written as the dimensional quasi-diffusion propagator, whenever the integral converges absolutely or conditionally.For || ≠ 0, by use of the asymptotic formulae for the Mittag-Leffler function (13) and for the Bessel function, Equation ( 25) is absolutely convergent for  < 4 − 1 and conditionally convergent for  < 4 + 1.In the one-dimensional case, it is absolutely convergent for 0.25 <  ≤ 1 and conditionally convergent for 0 <  ≤ 1.In two-dimensions, ( 25) is absolutely convergent for 0.75 <  ≤ 1 and conditionally convergent for 0.25 <  ≤ 1. Conditional convergence is guaranteed for the one-, two-, and three-dimensional cases for 0.5 <  ≤ 1.
We now consider the case of zero net displacement where || = 0. From ( 22) we have, which by the formula [62], can be expressed as, Using the asymptotic formula for the Mittag-Leffler function (13) it can be seen that Equation ( 28) is absolutely convergent if  < 2 and, consequently, only in the one-dimensional case when 0.5 <  ≤ 1.For the one-dimensional case we have, with a closed form for  1 (0, ) given by [47,48], which simplifies to the Gaussian solution,  1 (0, ) = 1 √4 1,2  ⁄ , when  = 1.For || ≠ 0, the quasi-diffusion propagator,   (, ), given by ( 25) cannot be written in simple closed form for the one-dimensional and three-dimensional cases, however, a closed form has been derived in the two-dimensional case.In one-dimension, the quasidiffusion propagator (25) is given by the integral, e.g., [47,48], as the Bessel function of order −1/2 is, Although (31) cannot be written in a simple closed form, it can be written in closed form using the Fox H-function [19,63].In two-dimensions, the closed form of ( 25) is given by e.g., [46][47][48], where  , () is the two parameter Mittag-Leffler function (see [61] for properties of the two-parameter Mittag-Leffler function), In the Gaussian case, when  = 1, Equation ( 34) takes the conventional form, Finally, in three dimensions, by substituting the Bessel function of order 1/2, into Equation ( 25), we have the integral form of the three-dimensional quasi-diffusion propagator for || ≠ 0 which is given by, e.g., [47,48], The quasi-diffusion propagator has also been investigated using subordination principles [45][46][47] (see [5] for a general description of subordination processes in the CTRW model).Subordination principles in stochastic processes involve the definition of a stochastic process in time (the subordinating function) that is within another stochastic process (the subordinated stochastic process).In the case of the diffusion propagator, this allows subordination formulae to be constructed that include, for example, Gaussian [47] or Poisson [45] distributions.Here we present the subordination equations for the Gaussian distribution [47], where the fundamental solution to the conventional diffusion equation is, and the function   () is the inverse Laplace transform of the stretched Mittag-Leffler function, As   (−  ) is a stretched Mittag-Leffler function, the inverse Laplace transform is a spectral probability density function, for this case, in time given by, Equations ( 38) to (41) show that the quasi-diffusion propagator is the integral of a multiplication of the conventional Gaussian probability density function and a spectral probability density function in time.The quasi-diffusion propagator is unimodal.Figure 4 illustrates the one-dimensional quasi-diffusion propagator in space, , in microns for an effective diffusion time of 35.9 ms (Figure 4a) and time, , in seconds for a distance of 5 μm (Figure 4b) for a diffusion coefficient of  1,2 = 1.5 × 10 −3 mm 2 s −1 and fractional exponents in the range 0.5 <  ≤ 0.99, which cover the range of  typically found in the human brain.The propagator was calculated based on Equation (38).The shape of the propagator in space (Figure 4a) is such that the probability density function becomes more kurtotic as  decreases from the Gaussian case ( = 1) corresponding to a higher probability of both smaller and larger step lengths, in correspondence with greater heterogeneity of the diffusion environment.This is consistent with assumptions made in DKI [24,25].The shape of the propagator in time (Figure 4b) is also consistent with quasi-diffusion representing an effective normal diffusion process as it is unimodal and Gaussian-like.Several interesting mathematical results have been derived from the Green's function in the special case of the CTRW model for 2/ = 1.It has been shown that the one-and two-dimensional quasi-diffusion propagators are probability density functions that evolve spatially in time [48][49][50].Furthermore, the entropy production rate of the quasidiffusion propagator in the one-dimensional and two-dimensional cases is the same as for Gaussian diffusion process [48][49][50].However, in the one-dimensional case, the second spatial moment of the probability density function,  1 (, ), does not exist and the mean squared displacement of the diffusing particles is not finite, indicating the diffusion process is anomalous [48][49][50].Similarly, for the two-dimensional case, the second moment of  2 (, ), does not exist for 0.5 <  ≤ 1 and the diffusion is anomalous with long-tailed waiting time and jump lengths.In the three-dimensional case, the interpretation of the process governed by ( 5) is different.Instead of a spatial probability density function for an anomalous diffusion process evolving in time, it has been suggested to represent an anomalous damped wave propagating with a time-dependent phase velocity [48,50].

Application of Quasi-Diffusion MRI to Mean Apparent Propagator Imaging
In this section, we consider how the quasi-diffusion propagator can be applied to estimate the mean apparent propagator and an effective pore size within a dMRI voxel.Several techniques have been previously described for computing effective pore size measurements from dMRI data, e.g., [11,28,[37][38][39]64], but here we consider the Mean apparent propagator (MAP) imaging technique [28] as this uses the inverse Fourier transform for dMRI data acquired in q-space and provides molecular displacements in mm.Our formulation differs from [28] by applying the functional form of the quasi-diffusion propagator rather than using Hermite functions to represent the diffusion propagator in -space.
MAP imaging provides measures of the three-dimensional return-to-the-origin probability (), the two-dimensional return-to-the-axis probability () and the onedimensional return-to-the-plane probability ().These maps represent the zero net displacement probabilities for molecular displacement and can be used to provide estimations of effective measurements of mean pore volume, 〈〉, mean pore cross-sectional area, 〈〉, and mean pore axial length scale, 〈〉.
In tissue microstructural dMRI, the "pores" under measurement are cells.In MAP imaging, the tissue microstructure is assumed to consist of isolated pores with a distribution of arbitrary shapes, sizes and orientations.In the brain, for example, the cortex and deep grey matter consist of neurons, whereas the connections between cortical regions are provided by axons in white matter.Grey matter may be considered to contain spheres with a distribution of radii (i.e., neurons) [26].White matter consists of bundles of axons that are surrounded by a myelin sheath and may be considered as cylinders with a distribution of radii (i.e., axons) [26].Diffusion parallel to axons is not free, but is greater than in axonal cross-section.Along axons, the presence of the nodes of Ranvier (where the myelin sheath periodically narrows), and astrocytes (which hold axonal bundles together) are considered to hinder diffusing water molecules.This may be considered as consecutive slabs with a distribution of spacings [26].
In dMRI, under the short pulse approximation, the zero net displacement probability can be computed to provide estimates of  in units of mm −3 ,  in mm −2 and  in mm −1 from which pore characteristics may be calculated.For the one-dimension quasi-diffusion case,  is given by ( 29) as, Although the quasi-diffusion propagator for || = 0 is only convergent in the onedimensional case, it is possible to provide estimations of the probability of zero net displacement in two-and three-dimensions.Such a calculation is of interest as estimation of  and  maps are sensitive to tissue microstructure and its alteration with disease [29][30][31].To ensure consistency in estimation of the two-and three-dimensional estimations we perform integration over a predefined range of -space.Here we perform integration up to a pre-defined limit of   = 5000 mm −1 , far higher than the capability of any MR scanner or NMR system. and  are then given by ( 28) as, The details of the relationship between zero net displacement probability and pore size in dMRI have been previously described in detail [11,[26][27][28].These amount to the following simple relationships for estimating effective mean pore volume, 〈〉, cross-sectional area 〈〉, and axial length scale, 〈〉, 〈〉 =  −1 , (45) To enable calculation of quasi-diffusion MAP images, there is a requirement to firstly estimate  1,2 and  values parallel to the tissue microstructure (i.e., axial  1,2 and ), perpendicular to the tissue microstructure (i.e., radial  1,2 and ) and obtain their mean isotropic characteristics (i.e., mean  1,2 and ) within the voxel.In practice, the axial, radial and mean  1,2 and  values are calculated using techniques described by Barrick et al. [1].In addition to the effective pore size estimates given by ( 46) to (48), the effective mean spherical pore radius, 〈〉, and effective mean cylindrical pore radius, 〈 ⊥ 〉, may be obtained from ( 45) and ( 46) as, In MAP imaging, the effective diffusion time of the dMRI acquisition is required to be long enough to allow the pores to be fully explored by the diffusing water molecules [28].In practice, the effective diffusion times are between 20 and 40 ms on clinical and research MR scanners and are sufficient for this purpose.However, according to Equation (30), quantification of the zero net displacement quasi-diffusion propagators at the effective diffusion time of the dMRI experiment will overestimate pore sizes.To avoid this problem, the quasi-diffusion propagators may be calculated at the time at which diffusing molecules first interact with the tissue microstructural environment.This is referred to as the short time limit,   .In quasi-diffusion imaging, we define this as the time at which the ensemble of diffusing molecules within a voxel ceases to be Gaussian, but is instead represented by an effective normal diffusion.This corresponds to the time at which the meansquared displacement, 〈 2 〉  = 2 1,2 ∆ ̅ , is equivalent to the mean-squared displacement for free water at body temperature, i.e.,   = 3 × 10 −3 mm 2 s −1 .Consequently, the short time limit is, and can be substituted into the quasi-diffusion propagators to provide a more accurate estimate of apparent mean pore size characteristics.A further observation is that, as in the case of Gaussian diffusion, the ratio,   , of the length scales at the effective diffusion time over the short time limit is given by both the diffusion times and diffusion coefficients as,

Quasi-Diffusion Mean Apparent Propagator Imaging of the Corpus Callosum
The corpus callosum is a large commissural white matter bundle that connects cortical regions in the left and right cerebral hemispheres of the brain.The structure of the corpus callosum is such that it does not intersect with other white matter bundles through the interhemispheric fissure (which is in the mid-sagittal plane of the brain).This means that dMRI measurements obtained where the corpus callosum crosses the interhemispheric fissure are only indicative of white matter tissue microstructure in the corpus callosum.In this example, we apply quasi-diffusion MAP imaging to estimate axon radii within a mid-sagittal section of the corpus callosum.
Diffusion MRI data were acquired on a healthy 27-year-old subject using a 3T clinical MR scanner at SGUL.Images were acquired with isotropic voxel size of 2 mm in 10 min 48 s.Data were acquired at  = {0, 1100,4000} s mm −2 in 32 non-collinear diffusion encoding directions equally spaced on the sphere with  = 28.7 ms and Δ = 43.9ms, giving an effective diffusion time of ∆ ̅ = 34.33 • ms.Full data acquisition details are provided in Appendix A.
Figure 5 shows  (Figure 5a),  (Figure 5b) and  (Figure 5c) maps in an axial slice through the corpus callosum.A sagittal cross-section through the midplane of the corpus callosum is shown in Figure 5d.The pseudo colour map indicates estimated axon radii in microns.A histogram of axon radii through a 10 mm thick sagittal section of the corpus callosum is shown in Figure 5e.This indicates axon radii are between 0.25 μm and 3.5 μm (mean 1.03 μm, standard deviation 0.5 μm, median 0.91 μm) with larger radii in the body than the genu and splenium.These measurements are in accordance with histological studies [67].

Quasi-Diffusion Imaging of Brain Tumour
As quasi-diffusion imaging can be acquired in 120 s with high signal to noise ratios [1], it could be feasibly acquired in the clinic.Here we consider an application of quasidiffusion imaging to a high-grade glial tumour, the glioblastoma multiforme (GBM).Clinical assessment of glial tumour grade and genetic subtype is relevant for predicting patient survival time and treatment strategy which includes surgical resection, chemo-and radiotherapy, and the need for accurate delineation of tumour tissue boundaries [68,69].
Diffusion MRI data were acquired using a 3T clinical MR scanner at SGUL.Images were acquired with a voxel size of 1.5 mm × 1.5 mm × 5 mm in 120 s.Data were acquired at  = {0, 1100, 5000} s mm −2 in 6 non-collinear diffusion encoding directions equally spaced on the sphere with δ = 23.5 ms, and Δ = 43.7 ms, giving an effective diffusion time of ∆ ̅ = 35.9ms.Full patient and data acquisition details are provided in Appendix A. Image analysis was performed using the same techniques as Section 3.1 with the exception that effective mean spherical pore radius was calculated.
In Figure 6, we compare the standard clinical Gadolinium contrast-enhanced T1weighted (T1wGd) image with parameters derived from quasi-diffusion imaging.In the T1wGd image (Figure 6a), the hyperintense region is typical for GBM and indicates the high-grade tumour core region (where there is angiogenesis and breakdown of the bloodbrain barrier) that would be the target for resection and/or highest radiotherapy dose; the central dark region is necrotic tissue.Surrounding this region is the area of tumour infiltration and oedema that is hypointense on the T1wGd image, and is more clearly seen as elevated D1,2 values (Figure 6b).This region plus a 2 cm margin would be the target for radiotherapy.The quasi-diffusion imaging maps provide different image contrasts with potential to distinguish these regions from a single image modality.Elevated α is evident across the whole lesion (Figure 6c) indicating that the oedema and tumour core are less structured than normal white matter.In addition, α and RTOP 1/3 (Figure 6d) allow distinction between grey and white matter for which contrast is present in the T1wGd image but not the D1,2 map.Within the RTOP 1/3 hypointense region that delineates tumour oedema/infiltration, there is an area of elevated signal relating to T1wGd enhancement, which is a result of the presence of a greater density of tumour microstructure.Within the T1wGd enhancing region, the effective pore size is approximately 2.5 µm (Figure 6e).and (e) effective mean pore size, 〈〉, (which in the absence of oedema or necrosis will relate to cell size).Arrows indicate the following regions: white-tumour core, red-necrosis, blue-oedema, yellow-grey matter, greenwhite matter.

Discussion
We have described the mathematical properties of quasi-diffusion, a special case of the CTRW diffusion model and presented the form of its Laplace and Fourier transforms.These functional forms provide insight into the behaviour and interpretation of quasidiffusion as a model of effective normal diffusion dynamics within a dMRI voxel.The quasi-diffusion equations provide a functional form describing diffusion in MRI.Together, the diffusion coefficient and fractional exponent of the quasi-diffusion characteristic equation define the form of the inverse Laplace transform which provides a spectrum of Fickian diffusivities, and an explicit form for the propagator.
In dMRI, there have been several attempts to describe the functional form of the signal attenuation [8,9,12,13].Unlike quasi-diffusion, these applied representations and models have not previously provided a simple functional form from which both the inverse Laplace transform and inverse Fourier transform of the signal can be derived.The quasi-diffusion model is compatible with previous frameworks and models that consider the underlying diffusion process in tissue microstructure to be locally Fickian [33][34][35] and lead to a signal that decays with a negative power law [13,19,[33][34][35]67].In the CTRW model, the assumption of a Gaussian scaling relationship between the fractional exponents of time, , and space, , leads to a mean-squared displacement that is linearly proportional to time.This leads to the quasi-diffusion model, that is valid in an image voxel beyond the short diffusion time limit, and prior to the long diffusion time limit where tortuosity effects will dominate signal behaviour.For accurate measurement of quasi-diffusion, there is a requirement that diffusing molecules have had sufficient time to explore the microstructural environment.
We consider quasi-diffusion to be the effect of a locally Gaussian diffusion that, within an image voxel, leads to observation of non-Gaussian characteristics due to the ensemble of diffusing molecules being within a heterogeneous environment containing pores with a distribution of structures and radii.The quasi-diffusion model makes no assumptions regarding the shape and distribution of the pores or that diffusing molecules are restricted to locations in the microstructural environment during the diffusion time.
The quasi-diffusion model predicts the existence of dMRI signal at extremely high bvalues due to the negative power law behaviour of the stretched Mittag-Leffler function as  → ∞ + .This prediction is consistent with experimental results that observe signal and power law decay at high diffusion sensitisations of  = 10,000 s mm −2 in human white matter in vivo [70], and  = 100,000 s mm −2 in rat and human ex vivo corpus callosum [67].Furthermore, the quasi-diffusion model predicts a stretched exponential form of the signal as  → 0 + consistent with low -values ( < 300 s mm −2 ) containing a significant signal proportion attributable to blood perfusion and a potentially super-diffusive dynamic [71].
The quasi-diffusion signal can be decomposed into a spectrum of Fickian apparent diffusion coefficients via the inverse Laplace transform.Such a transform has been suggested to relate to specific length scales within the diffusion environment [72,73].Visualisation of the dMRI signal as a diffusion spectrum was first introduced by Yablonskiy et al. [72], where a Gaussian probability density function was assumed.More recently, the assumption of a gamma variate distribution [74] has indicated spectral regions related to pathological changes in disease [74][75][76].However, each of these studies make assumptions as to the functional form of the probability distribution without supporting evidence.In contrast, the quasi-diffusion model provides a functional form for the spectrum with the  1,2 and  parameters together providing an indication of heterogeneity and restriction within the diffusion environment.
The long diffusion times afforded by clinical MR scanners and the complexity of the tissue environment ensure that diffraction patterns from groups of similar pores [11,77] cannot be observed from the dMRI signal [8,78].Instead, the signal is a relaxation curve with few inherent features.The quasi-diffusion model provides a functional form for the signal which offers an opportunity for developing new quantitative model-based dMRI analysis techniques.We have shown here that the assumption of quasi-diffusion leads to model-based power law and 'kurtosis' information, and that the propagator can be applied to obtain MAP images.Possible further applications include use of the diffusion spectrum to provide model-based Multi-dimensional Diffusion MRI [79,80], and use of the propagator to obtain model-based three-dimensional orientation distribution functions of tissue microstructure (similar to those described by [28,[81][82][83]) that could be employed within multi-tissue constrained spherical deconvolution techniques for identifying white matter fibre crossing [84].
Advanced diffusion MRI techniques, such as DKI, have been shown to be a sensitive probe of the tissue microstructural environment in healthy and diseased tissue [21,[85][86][87][88].We expect  1,2 ,  and quasi-diffusion MAP images to show a similar sensitivity to disease to previous clinical studies of DKI and MAP imaging.Our examples highlight the potential for quasi-diffusion imaging in clinical research and demonstrate that effective mean pore size measurements may be obtained.Our initial pore size estimations appear to be more accurate in the white matter of the corpus callosum, where the high density of myelinated axons can be considered as isolated pores, than in grey matter and tumour tissue, where more complex distributions of cell structures and shapes are present with permeable membranes, potentially leading to a breakdown of the assumption of isolated pores.Further research is required to determine the efficacy of the quasi-diffusion MAP approach and its clinical potential.
A significant advantage of quasi-diffusion imaging over current dMRI techniques is its ability to provide robust and accurate fits to minimal imaging data [1,55].This allows acquisition of  1,2 and  maps with high signal to noise ratios in a clinically feasible time.New dMRI analysis techniques developed using the quasi-diffusion model will benefit from this advantage, which will potentially enable rapid translation of new techniques from large scale clinical feasibility studies into clinical practice.
There are limitations to the quasi-diffusion imaging technique.For instance, assumption of the short pulse approximation in -space imaging is generally violated in practice, and the effect of this on measurements is currently unclear.However, the quasi-diffusion model is unique in anomalous dMRI techniques as it does not undermine the Gaussian phase approximation (GPA) [8,11].Further studies are required to understand the conditions under which quasi-diffusion dynamics emerge within a voxel.In particular, the relationship of the signal to the underlying microstructure, and the effect of the PGSE  and  on the signal are not well-understood.It is also not yet clear whether the quasidiffusion model provides imaging measures that are invariant to diffusion time.Nevertheless, the quasi-diffusion model provides a well-defined mathematical basis for further research to elucidate characteristics of the diffusion environment.

Conclusions
We have described the mathematics and initial interpretation of the quasi-diffusion model in application to dMRI.This new model of the diffusion dynamics within MRI voxels provides a functional form for the signal from which a diffusion spectrum and propagator can be derived.The quasi-diffusion model is consistent with current understanding of diffusion in tissue microstructure and can be used to provide high quality images of parameters relating to non-Gaussian diffusion in clinically feasible time.The functional form of the quasi-diffusion model represents a novel focus for future mathematical and imaging research.This study shows the potential of quasi-diffusion imaging as a promising new dMRI technique with significant scope for further development.

Patents
The quasi-diffusion imaging technique is covered by the patent, "A method and apparatus for quasi-diffusion magnetic resonance imaging", by inventors Dr. T.R. Barrick without diffusion-sensitisation ( = 0 s mm −2 ) were acquired 16 times.Diffusion encoding gradients were applied in 6 non-collinear directions.Data were acquired in 35 min 12 s.
Acquisition of imaging data illustrated in Section 3.1 and Figure 5: Quasi-diffusion tensor imaging data were acquired with:  = 90 ms,  = 9000 ms,  = 28.7 ms,  = 43.9ms, field of view 224 mm × 224 mm with 50 2 mm thick slices acquired at 3 mm × 3 mm × 2 mm resolution that were zero-filled (by use of the Fourier transform) to provide 2 mm isotropic voxels.dMRI were acquired once in 3  -values with  = {0, 1100, 4000} s mm −2 .Images without diffusion-sensitisation were acquired 8 times.Diffusion encoding gradients were applied in 32 non-collinear directions.Data were acquired in 10 min 48 s.
Acquisition of imaging data illustrated in Section 3.2 and Figure 6: Quasi-diffusion tensor imaging data were acquired with:  = 90 ms,  = 6000 ms, δ = 23.5 ms,  = 43.9ms, field of view 210 mm × 210 mm with 22 5 mm thick slices acquired at 2.3 mm × 2.3 mm × 5 mm resolution that were zero-filled (by use of the Fourier transform) to provide 1.5 mm × 1.5 mm × 5 mm voxels.dMRI were acquired once in 3 -values with  = {0, 1100, 5000} s mm −2 .Images without diffusion-sensitisation were acquired 8 times.Diffusion encoding gradients were applied in 6 non-collinear directions.Data were acquired in 120 s.

Figure 1 .
Figure 1.The family of quasi-diffusion imaging signal attenuation curves for an arbitrary diffusion coefficient of  1,2 = 1.5 × 10 −3 mm 2 s −1 and a range of fractional exponents, 0.1 ≤  ≤ 0.99.Normalised signal is shown parameterised by (a) diffusion-sensitisation, , and (b) .Graph (c) shows the power law behaviour of the signal decay parameterised by  on a logarithmic scale.

Figure 2 .
Figure 2. Quasi-diffusion model fits to experimental diffusion magnetic resonance imaging data acquired in the brain of a young, healthy subject (age 28 years).The top row shows the results of fitting the quasi-diffusion model to an acquisition of 29 -values over the range 0 <  ≤ 5000 s mm −2 , with the bottom row showing the results of fitting the quasidiffusion model to 3 -values,  = {0,1080,5000} s mm −2 .All imaging data were acquired in 6 diffusion encoding directions at an effective diffusion time of Δ ̅ = 35.9ms.Normalised signal attenuation is shown for a grey matter voxel indicated by the blue arrow (graphs (a,e)), and a white matter voxel indicated by the red arrow (graphs (b,f)).Axial slices are shown for maps of mean  1,2 (images (c,g)) and mean α (images (c,g)).

Figure 3 .
Figure 3. Decomposition of the signal into a spectrum of Fickian apparent diffusion coefficients via the inverse Laplace transform.Graph (a) shows the probability density functions for a unit diffusion coefficient,  1,2 = 1 mm 2 s −1 for different fractional exponents, 0.1 ≤  ≤ 0.99.Spectra are shown for axial, radial and mean  1,2 and  values in (b) cerebrospinal fluid, (c) grey matter and (d) white matter. 1,2 and  values used to calculate the spectra are from Barrick et al. [1].

Figure 5 .
Figure 5. Quasi-diffusion mean apparent propagator imaging results for a young, healthy subject (age 27 years).The top row shows axial slices through (a) √ 3 , (b) √, and (c)  maps.The bottom row shows (d) effective pore radius estimates in a cross-section through the mid-sagittal plane of the corpus callosum, and (e) a probability density function of micron radii for the analysed 10 mm thick section of the corpus callosum (graph (e)).The location of the corpus callosum is identified by red arrows.The red and blue lines in graph (e) indicate the mean and median, respectively.

Figure 6 .
Figure 6.Quasi-diffusion imaging results for axial slices through a high-grade glial tumour (WHO Grade IV glioblastoma).Diffusion magnetic resonance imaging data were acquired at  = {0, 1100,5000} s mm −2 in 6 diffusion encoding directions at an effective diffusion time of Δ = 35.9ms.Imaging data were acquired in a clinically feasible time of 120 s.Image (a) shows a T1-weighted image acquired after injection of gadolinium contrast agent, with maps of (b) mean  1,2 , (c) mean , (d) √