Fractional Order Complexity Model of the Di ﬀ usion Signal Decay in MRI

: Fractional calculus models are steadily being incorporated into descriptions of di ﬀ usion in complex, heterogeneous materials. Biological tissues, when viewed using di ﬀ usion-weighted, magnetic resonance imaging (MRI), hinder and restrict the di ﬀ usion of water at the molecular, sub-cellular, and cellular scales. Thus, tissue features can be encoded in the attenuation of the observed MRI signal through the fractional order of the time- and space-derivatives. Speciﬁcally, in solving the Bloch-Torrey equation, fractional order imaging biomarkers are identiﬁed that connect the continuous time random walk model of Brownian motion to the structure and composition of cells, cell membranes, proteins, and lipids. In this way, the decay of the induced magnetization is inﬂuenced by the micro- and meso-structure of tissues, such as the white and gray matter of the brain or the cortex and medulla of the kidney. Fractional calculus provides new functions (Mittag-Le ﬄ er and Kilbas-Saigo) that characterize tissue in a concise way. In this paper, we describe the exponential, stretched exponential, and fractional order models that have been proposed and applied in MRI, examine the connection between the model parameters and the underlying tissue structure, and explore the potential for using di ﬀ usion-weighted MRI to extract biomarkers associated with normal growth, aging, and the onset of disease.


Introduction
"One of the principal objects of theoretical research in my department of knowledge is to find the point of view from which the subject appears in the greatest simplicity". J. Willard Gibbs, Letter in Proc. Amer. Acad. Arts & Sci. (1881), pp. 420-421.
Mathematical models are a prism through which we can view the complexity of nature [1]. Just as a prism separates sunlight into the colors of its optical spectrum-portraying hidden features (frequency, intensity, and polarization)-the formulation of a model identifies features not displayed in the raw data. Parameter extraction and estimation mimic a spectrograph by selecting individual spectral components for analysis. In both measurement and modeling, we seek to isolate specific aspects of a physical phenomenon for further study. Success is measured by the degree to which the spectra (or model) captures particular features of the image or visual scene. In paint, the pigments can be identified by a UV-Vis spectrophotometer, but in painting color is described by hue, value, tone, tint, shade, and saturation, all of which are processed by the eye as characteristics of a picture. In modeling biomedical images, we use mathematics to help the brain recognize patterns by identifying new measures of image Two research paradigms are commonly explored in DW-MRI. First, one can build simplified models of tissue structure by selecting two, three, or four tissue compartments or components, each with a specific, and to be determined, parameter set (e.g., intra-and extra-cellular diffusion and membrane permeability). Diffusion Tensor Imaging (DTI) applies this approach to each image voxel via gradient pulses applied in multiple (at least six) directions. Other approaches (CHARMED, NOODI, AxCaliper, etc.) [15], extend the data collection to multiple b-value shells sampled in hundreds of directions. Alternatively, one can select a heuristic model for S(b), such as the kurtosis, S(b) = S0exp(−bD)exp[(K/6)(bD) 2 ] [16], or the stretched exponential model, S(b) = S0 exp[−(bD) α ] [17], and fit the available clinical data to the extended set of model parameters {S0,D,K} or {S0,D,α}, respectively. Such models, since they have an extra parameter, will provide an improved fit to S(b) data that extends over multiple b-values. The sword of Damocles in the first approach is the added imaging time, while for the second, it is the dichotomy between the precision of the parameter fit and the ambiguity in the connection between the tissue composition and fitting parameter. Another factor, most often overlooked for both models is the likelihood that the fitting process is degenerate, that is, more than one set of parameters can fit the data with the same least squares fit error [18].
In this paper, we will combine the two modeling approaches described above by using a varying diffusion coefficient D(b) to account for the complexity of the tissue as an inverse power law for higher values of diffusion weighting (b-values), and by generalizing the stretched exponential and kurtosis diffusion decay models by extending the governing diffusion decay to fractional order. The fractional order model naturally accounts for both multiple compartments and non-Gaussian diffusion. After a few preliminary definitions, we will outline an extended form for the fractional diffusion coefficient model, and Two research paradigms are commonly explored in DW-MRI. First, one can build simplified models of tissue structure by selecting two, three, or four tissue compartments or components, each with a specific, and to be determined, parameter set (e.g., intra-and extra-cellular diffusion and membrane permeability). Diffusion Tensor Imaging (DTI) applies this approach to each image voxel via gradient pulses applied in multiple (at least six) directions. Other approaches (CHARMED, NOODI, AxCaliper, etc.) [15], extend the data collection to multiple b-value shells sampled in hundreds of directions. Alternatively, one can select a heuristic model for S(b), such as the kurtosis, S(b) = S 0 exp(−bD)exp[(K/6)(bD) 2 ] [16], or the stretched exponential model, S(b) = S 0 exp[−(bD) α ] [17], and fit the available clinical data to the extended set of model parameters {S 0 ,D,K} or {S 0 ,D,α}, respectively. Such models, since they have an extra parameter, will provide an improved fit to S(b) data that extends over multiple b-values. The sword of Damocles in the first approach is the added imaging time, while for the second, it is the dichotomy between the precision of the parameter fit and the ambiguity in the connection between the tissue composition and fitting parameter. Another factor, most often overlooked for both models is the likelihood that the fitting process is degenerate, that is, more than one set of parameters can fit the data with the same least squares fit error [18].
In this paper, we will combine the two modeling approaches described above by using a varying diffusion coefficient D(b) to account for the complexity of the tissue as an inverse power law for higher values of diffusion weighting (b-values), and by generalizing the stretched exponential and kurtosis diffusion decay models by extending the governing diffusion decay to fractional order. The fractional order model naturally accounts for both multiple compartments and non-Gaussian diffusion. After a few preliminary definitions, we will outline an extended form for the fractional diffusion coefficient model, and the corresponding integer and fractional order diffusion decay equations. Examples of the expected functional behavior of the diffusion signal decay as a function of b-value will be plotted and fit to ex vivo bovine optic nerve tissue. In the Discussion we will describe how these models expand the available modeling tools for DW-MRI data and describe how they can be used to fit animal and clinical data. Overall, this approach provides a way to extend the heuristic stretched exponential approach toward more complete multiple-compartment models. This bridge uses fractional order derivatives and varying diffusion coefficients as connecting links.

Definitions and Properties
The signals analyzed in DW-MRI, f (b) = S(b)/S 0 , are expressed in terms of b, the diffusion-weighting variable defined above for a pair of rectangular diffusion pulses [2][3][4]. Typically, f (b) is a single valued function that is positive, real, and monotonically decreasing for b-values in the range from 0 to 5000 s/mm 2 . For each pixel in an MR image, f (b) is fit to data acquired by varying b in magnitude and gradient direction. The default function used in most DW-MRI is the single exponential, f (b) = exp(−bD). This model reflects the underlying assumptions of a linear, first order relaxation model appropriate for free or hindered Gaussian diffusion. In biological tissues, a sum of exponentials or the stretched exponential, f (b) = exp[(−bD) α ], where 0 < α < 1, is often needed to describe the data, particularly at high b-values. When neither of these approaches are satisfactory, non-Gaussian behavior is considered, and in this paper we will describe it by using a fractional-order relaxation model. Our model involves the Caputo fractional derivative, C D α [f (b)], so below we have defined this derivative operator and its properties. In addition, we also define the Mittag-Leffler function because this function frequently appears in various forms when solving fractional-order differential equations. We will show that the Mittag-Leffler function is also a generalization of the exponential.
Applying the Caputo fractional derivative term by term to the power series expansion of the single exponential f(b) = exp(b), yields: Rather than tabulate all the changes incurred when applying the Caputo fractional derivative to familiar functions, mathematicians have developed particular special functions with simpler behavior under the fractional derivative operation. One such function was defined by the Swedish mathematician, Magnus Gustaf Mittag-Leffer; it is a generalization of the exponential.

Mittag-Leffer Function
The Mittag-Leffer function, like the Bessel function, comes in different kinds (first, second, modified, etc.), and occurs in both integer and fractional order. It was developed out of the need to simplify the solution of certain classes of ordinary differential equations. Also, like the Bessel function, the Mittag-Leffer function can be defined in terms of a convergent power series. For more information on the Mittag-Leffer function, please consult the monograph by Mainardi and colleagues [20]. There are many versions of the Mittag-Leffler function which subsume almost all of the basic functions of mathematical physics. In this article, we only need three kinds of Mittag-Leffler functions, which we will define in terms of their corresponding power series.

Exponential Function
Since Γ(n + 1) = n!, the simple decaying exponential has the following power series representation:

One Parameter Mittag-Leffer Function
The one parameter Mittag-Leffler function [20] is a generalization of the exponential with the parameter, α, inserted a multiplying the integer n. It has the following power series representation: Note, that E 1 (−b) = exp(−b) and E 2 (−b 2 ) = cos(b).

Two Parameter Mittag-Leffer Function
The two parameter Mittag-Leffler function [20] includes a second parameter β in the gamma function of the power series as: Note that here,

Three Parameter Mittag-Leffer Function
There are several ways to extend the Mittag-Leffler function to three parameters-one credited to Prabhakar [21], and another to Kilbas and Saigo [22]. We will use the Kilbas-Saigo definition, which in power law form can be written as: where c 0 = 1 and: Note that E 1,1, Two useful results can be established by applying the Caputo derivative term by term to these power series expressions. First, we find the following for the one parameter Mittag-Leffler function: Thus, the stretched Mittag-Leffler function solves the homogeneous fractional order differential And, for α = 1, the above expression reduces to the simple relationship d[exp(−bD)]/db = -D exp(−bD), and gives the solution to the corresponding first order differential equation.

Fractional Order D(b) Models
A key feature of the anomalous diffusion exhibited in DW-MRI is the fall-off or reduction of the decay rate at b-values above 1000 s/mm 2 . A recent paper [14] surveyed integer order models of this phenomena, where the apparent diffusion coefficient D 0 was expressed as a decaying function of b with a separate decay rate, D 1 , for example, Here, we focus on fractional order models for systems with both a constant and an inverse power law decay rate. A summary of our approach is given in Figure 2, which illustrates the different integer and fractional order models that will be described. For the integer order case, with a constant diffusion coefficient, D(b) = D 0 , the linear differential equation model, yields the classical exponential decay, A common extension of this model [23,24] assumes an inverse power law decay rate, This integer order, linear differential equation has a stretched exponential solution, , which has often been used to fit DW-MRI signals (see, [17,24,25]). The fractional β parameter is zero for the cases of normal and hindered diffusion, as in CSF and gray matter, but less than one for situations where anomalous diffusion appears-such as in brain white matter. Nevertheless, this signal decay model is still exponential, and will decay faster than a power law at high b-values. Fractional order models asymptotically decay as a power law, so another way to capture the power law tail of the anomalous diffusion signal decay is to consider the decay process to be governed by a Caputo fractional order derivative (order α, with 0 < α < 1) and a constant diffusion coefficient, D0. The solution for this case is the single parameter stretched Mittag-Leffer function, [23]. It has been fit to DW-MRI data by Magin, Ingo, and others [12,26]. For completeness, we also show the case of a fractional order Caputo derivative model (order α) with an inverse power law decay rate (order β). The solution for this case is a three parameter Mittag-Leffler function (Kilbas-Saigo form, see Equation (15) and Table 1), but for the inverse power law decay rate is expressed using only α and β as S(b)/S0 = Eα,1+β/α,β/α[-(bD0) α+β ], [23]. Two things should be noted about this function. First, the conflation of the separate fractional parameters-the α connected with the order of the Caputo fractional derivative, while the β is connected with the presumed inverse power law decay of the diffusion decay rate, D(b). The second thing to note is that, the α and β in this solution are not the α, β of the second order Mittag-Leffler function, but reflect a particular condition on the α, m, and l parameters of the three parameter Mittag-Leffler function. In addition, from its construction, we note that the fractional order derivative with the inverse power law decay rate case includes all of the other models as special cases. All four are displayed in Table 1 for reference and ease of comparison.
In order to display the behavior of the four signal decay models, we have prepared plots of the different diffusion signal decay functions for a selected set of the α and β parameters. The functions were evaluated using MATLAB ® (MathWorks, Natick, MA 01760) and the code is publicly available at [27]. In the MRI field, the stretched exponential was first applied by Bennett et al. [17] to model DW-MRI data.  Fractional order models asymptotically decay as a power law, so another way to capture the power law tail of the anomalous diffusion signal decay is to consider the decay process to be governed by a Caputo fractional order derivative (order α, with 0 < α < 1) and a constant diffusion coefficient, D 0 . The solution for this case is the single parameter stretched Mittag-Leffer function, This model is attractive because at low b-values it decays as a stretched exponential (approx. 1 − b α /Γ(1 + α) + . . . ) and at high b-values decays as a power law of the form (b −α /Γ(1 − α)) [23]. It has been fit to DW-MRI data by Magin, Ingo, and others [12,26]. For completeness, we also show the case of a fractional order Caputo derivative model (order α) with an inverse power law decay rate (order β). The solution for this case is a three parameter Mittag-Leffler function (Kilbas-Saigo form, see Equation (15) and Table 1), but for the inverse power law decay rate is expressed using only α and β as [23]. Two things should be noted about this function. First, the conflation of the separate fractional parameters-the α connected with the order of the Caputo fractional derivative, while the β is connected with the presumed inverse power law decay of the diffusion decay rate, D(b). The second thing to note is that, the α and β in this solution are not the α, β of the second order Mittag-Leffler function, but reflect a particular condition on the α, m, and l parameters of the three parameter Mittag-Leffler function. In addition, from its construction, we note that the fractional order derivative with the inverse power law decay rate case includes all of the other models as special cases. All four are displayed in Table 1 for reference and ease of comparison.
In order to display the behavior of the four signal decay models, we have prepared plots of the different diffusion signal decay functions for a selected set of the α and β parameters. The functions were evaluated using MATLAB ® (MathWorks, Natick, MA 01760) and the code is publicly available at [27]. In the MRI field, the stretched exponential was first applied by Bennett et al. [17] to model DW-MRI data. 0.37 when bD 0 = 1, and for smaller b-values, all the curves fall below the exponential, while for bD 0 > 1, the decay curves are all higher. Thus, the stretched exponential function mimics the combination of a fast plus a slow decay. The deviation from the exponential decreases as the value of (1 + β) increases. Second, slope of the stretched exponential becomes increasingly steep near b = 0, and in fact diverges when (1 + β) < 1. This behavior can be eliminated by shifting the stretched exponential, but that changes the character of the decay [28]. Third, for high b-values, the stretched exponential decays as a power law on a semilogarithmic graph, but not on a linear scale. The multiexponential character of the stretched exponential can be extracted using the inverse Laplace transform (see the work of Berberan-Santos et al., [29]) and it also finds applications in the cumulative function of the Weibull probability distribution [30]. values, all the curves fall below the exponential, while for bD0 > 1, the decay curves are all higher. Thus, the stretched exponential function mimics the combination of a fast plus a slow decay. The deviation from the exponential decreases as the value of (1 + β) increases. Second, slope of the stretched exponential becomes increasingly steep near b = 0, and in fact diverges when (1 + β) < 1. This behavior can be eliminated by shifting the stretched exponential, but that changes the character of the decay [28]. Third, for high b-values, the stretched exponential decays as a power law on a semilogarithmic graph, but not on a linear scale. The multiexponential character of the stretched exponential can be extracted using the inverse Laplace transform (see the work of Berberan-Santos et al., [29]) and it also finds applications in the cumulative function of the Weibull probability distribution [30]. The last value in this sequence gives the single exponential function. As the value of (1 + β) moves toward 1, the shape of the decay curve transitions from a multiple exponential with fast and slow components to a single exponential decay.
The single parameter Mittag-Leffer function is well known as a relaxation model for dielectric and viscoelastic materials [31,32], and more recently in MRI [12,26]. The last value in this sequence gives the single exponential function. As the value of (1 + β) moves toward 1, the shape of the decay curve transitions from a multiple exponential with fast and slow components to a single exponential decay.
The single parameter Mittag-Leffer function is well known as a relaxation model for dielectric and viscoelastic materials [31,32], and more recently in MRI [12,26]. is a plot S(b)/S 0 = E α [-(bD 0 ) α ] versus b for b-values between 0 and 4,000 s/mm 2 , and D 0 = 1 × 10 −3 mm 2 /s. We selected five values for α = 0.2, 0.4, 0.6, 0.8, and 1.0. The last value in this sequence, corresponding to α = 1, gives the single exponential function. The one parameter Mittag-Leffler function behaves, overall, in a manner similar to that of the stretched exponential-falling below the exponential at low b-values, and above it at high b-values. However, as shown in Figure 4, all the curves do not pass through bD 0 = 1, but most go through S(b)/S 0 = 0.5 near bD 0 = 0.7. And, for high b-values, on linear scales, the curves flatten out as the α values decrease, reflecting a significant reduction in the apparent diffusion constant-a hallmark of restricted diffusion. On logarithmic scales at high b-values, the one parameter Mittag-Leffler will appear as a straight line, which is characteristic of a pure power decay (see, Carpinteri and Mainardi [31] for these plots). The one parameter Mittag-Leffler has been used to fit DW-MRI data by several groups [12,26] and there is a growing recognition that the fractional order, α, is a homogeneity measure of the sub-voxel tissue environment-as the cellular heterogeneity increases (homogeneity decreases), the α values fall below 1.0. This view is consistent with the spectral distribution model of the one parameter Mittag-Leffler presented by Carpinteri and Mainardi [31], and the inverse Laplace transform analysis of Berberan-Santos [33]. Both models coalesce to a Dirac delta function when α = 1.0, and broaden significantly as α < 1.0, indicating a growing population of compartments with low diffusion rates (coefficients). Hence, one can interpret the one parameter Mittag-Leffler function as a concise multiexponential model of complex, heterogeneous material, and a natural candidate to use when searching for imaging biomarkers of tissue changes or pathology.
Mittag-Leffler has been used to fit DW-MRI data by several groups [12,26] and there is a growing recognition that the fractional order, α, is a homogeneity measure of the sub-voxel tissue environment-as the cellular heterogeneity increases (homogeneity decreases), the α values fall below 1.0. This view is consistent with the spectral distribution model of the one parameter Mittag-Leffler presented by Carpinteri and Mainardi [31], and the inverse Laplace transform analysis of Berberan-Santos [33]. Both models coalesce to a Dirac delta function when α = 1.0, and broaden significantly as α < 1.0, indicating a growing population of compartments with low diffusion rates (coefficients). Hence, one can interpret the one parameter Mittag-Leffler function as a concise multiexponential model of complex, heterogeneous material, and a natural candidate to use when searching for imaging biomarkers of tissue changes or pathology. The three parameter Mittag-Leffler function (Kilbas-Saigo function), when expressed using only α and β was derived and first applied in MRI by Hanyga and Seredynska [34], but plotted only over a narrow range. Here, we follow the recent work of Capelas de Oliveira, Mainardi, and Vaz [22] and, in addition, suggest that the α and β parameters of the Kilbas-Saigo function might provide separate tissue characterization parameters. We selected different values for α and β subject to the conditions, 0 < α < 1 and 0 < α + β < 1 to ensure complete monotonicity and non-negativity of the spectral distribution [22].   As the value of α decreases, the shape of the decay curve transitions from stretched exponential to an asymptotic power law decay.
The three parameter Mittag-Leffler function (Kilbas-Saigo function), when expressed using only α and β was derived and first applied in MRI by Hanyga and Seredynska [34], but plotted only over a narrow range. Here, we follow the recent work of Capelas de Oliveira, Mainardi, and Vaz [22] and, in addition, suggest that the α and β parameters of the Kilbas-Saigo function might provide separate tissue characterization parameters. We selected different values for α and β subject to the conditions, 0 < α < 1 and 0 < α + β < 1 to ensure complete monotonicity and non-negativity of the spectral distribution [22].   [22]), decreasing β only broadens the distribution about its mean value, which attenuates the diffusion signal in a more uniform manner over the entire range of b-values. and 1.0. In the figure, we see for a fixed power law decay of D(b) the shift from lower to higher decay rates as α approaches the exponential function. The separate effects of the two fractional order parameters are also portrayed in the spectral distribution (inverse Laplace transform) of the Kilbas-Saigo function, where for β = 0 (Figure 4 in [22]) decreasing α shifts the spectrum to lower frequencies (smaller D values), which is the basis for the power law tail in the signal decay curve. On the other hand, for fixed α (Figures 5 and 6 in [22]), decreasing β only broadens the distribution about its mean value, which attenuates the diffusion signal in a more uniform manner over the entire range of b-values.

Fitting Fractional Order D(b) Models to DW-MRI Data
As an example application of these DW-MRI signal decay models, we fit the stretched Kilbas-Saigo function to published ex vivo data acquired from an intact bovine optic nerve [35]. Fits to the exponential, stretched exponential, stretched one parameter Mittag-Leffer were also performed and the results are provided as Supplementary Material. Since the optic nerve consists of bundles of axons, the diffusion signals were measured both along (parallel) and across (perpendicular) the axons. Given the cylindrical structure of the individual nerve axons, we would expect the diffusion to be hindered in the direction along the nerve, and restricted in the direction across the nerve. The Stejskal-Tanner diffusion pulse sequence [5] consists of a pair of rectangular gradient pulses of height g, duration δ, and separation ∆. It was applied by varying the gradient strength, g, and direction for a series of increasing diffusion times (∆ = 8, 10, 20, and 30 ms) with δ set to 3 ms in all cases. The increase in the diffusion time between the gradient pulses can be expected to increase the restricted diffusion contribution for both gradient directions (i.e., when the diffusion length, L D = (2D∆) 1/2 exceeds the radius of the axons, typically 5 microns in the parallel case and less than one micron in the perpendicular case). The data was extracted from Figure 1 in [35], and fit using nonlinear least squares optimization in MATLAB ® . The mean squared error was calculated for the Kilbas-Saigo function in both the parallel and perpendicular gradient directions. These results are presented in Table 2. As expected, the parallel gradient direction exhibited larger values of the diffusion coefficient than the perpendicular gradient direction (0.72 × 10 −3 mm 2 /s versus 0.2 × 10 −3 mm 2 /s, when ∆ = 30 ms). The stretched Kilbas-Saigo function was able to fit both the parallel and the perpendicular data for all diffusion times. Figure 8 is a plot of the parallel and perpendicular cases when ∆ = 30 ms. In the parallel case the α value (ca., 0.75) did not change with increasing diffusion time, while the β value decreased from 0.33 to 0.06. In the perpendicular case, the α value fell from 0.62 to 0.38 with increasing diffusion time, while the β value decreased only from 0.31 to 0.20. These results suggest that the stretched Kilbas-Saigo function is sensitive to both gradient direction and diffusion time. Our goal in this proof of concept paper was only to introduce the Kilbas-Saigo function as a candidate for modeling DW-MRI data. Nevertheless, the nature of the α and β in the model offers the potential of distinguishing between multiexponential (compartmental) contributions via the α parameter and tissue complexity (surface to volume ratio and membrane permeability) contributions through the β parameter. Further analysis of whole imaging slices is needed to establish whether or not this, or another model provides better contrast.   [35] and are used with permission from the journal Magnetic Resonance in Medicine.

Discussion
"One cannot collect all the beautiful shells on the beach" Anne Morrow Lindbergh, Gift from the Sea, 1955.
Fractional calculus interpolates between the integer order operations of calculus just as the real numbers are interspersed between integers. Investigators can employ various fractional derivatives (e.g., Grünwald-Letnikov, Riemann-Liouville, Caputo, Riesz, Hadamard, and Erdélyi-Kober) to stitch together earlier or surrounding data in an extended space-time model of a physical process such as diffusion or dielectric relaxation [19]. The fractional derivatives threads are woven into the model via generalized transcendental functions (Mittag-Leffler, Kilbas-Saigo, Wright, and the H-function of Fox) which are functions that often fit the experimental data with greater fidelity and precision [30][31][32]. In this paper, addressing the decay of diffusion-weighted magnetic resonance signals, we view the problem from a simple perspective-that of a fractional order relaxation with either a fixed or a power-law varying rate constant. The resulting fractional order Kilbas-Saigo model encompasses the simple exponential, the stretched exponential and the stretched Mittag-Leffler function as special cases [21][22][23].
In previous papers, we investigated the challenge of extracting specific complexity measures from DW-MRI data [12,13]. The choices are myriad, while the data is often sparse. Nevertheless, by tuning the applied MR pulse sequence (e.g., by changing the gradient direction or by increasing the diffusion time between gradient pulses) to match a particular aspect of the tissue under study, savings in time and an

Discussion
"One cannot collect all the beautiful shells on the beach" Anne Morrow Lindbergh, Gift from the Sea, 1955.
Fractional calculus interpolates between the integer order operations of calculus just as the real numbers are interspersed between integers. Investigators can employ various fractional derivatives (e.g., Grünwald-Letnikov, Riemann-Liouville, Caputo, Riesz, Hadamard, and Erdélyi-Kober) to stitch together earlier or surrounding data in an extended space-time model of a physical process such as diffusion or dielectric relaxation [19]. The fractional derivatives threads are woven into the model via generalized transcendental functions (Mittag-Leffler, Kilbas-Saigo, Wright, and the H-function of Fox) which are functions that often fit the experimental data with greater fidelity and precision [30][31][32]. In this paper, addressing the decay of diffusion-weighted magnetic resonance signals, we view the problem from a simple perspective-that of a fractional order relaxation with either a fixed or a power-law varying rate constant. The resulting fractional order Kilbas-Saigo model encompasses the simple exponential, the stretched exponential and the stretched Mittag-Leffler function as special cases [21][22][23].
In previous papers, we investigated the challenge of extracting specific complexity measures from DW-MRI data [12,13]. The choices are myriad, while the data is often sparse. Nevertheless, by tuning the applied MR pulse sequence (e.g., by changing the gradient direction or by increasing the diffusion time between gradient pulses) to match a particular aspect of the tissue under study, savings in time and an economy of interpretation are possible. For example, in diffusion tensor imaging (DTI), one uses a simple exponential decay function, but applies the diffusion pulse gradients in six or more directions to capture the anisotropy of skeletal muscle or brain white matter [10]. In a similar manner, in diffusion kurtosis imaging (DKI), by applying gradients in multiple directions and at multiple b-values, the signal decay encodes a measure of restricted diffusion [16]. In some situations, we acquire diffusion trace data to reduce the contributions of tissue anisotropy [36], but extend the range of b-values to 4000 s/mm 2 or greater to capture aspects of sub-voxel tissue complexity in the fractional order parameters of the Mittag-Leffler function [37,38].
The connection between tissue complexity and the fractional order diffusion decay process is explicit when extremely short diffusion gradient pulses are employed [12,26]. For this situation, the time-space fractional order diffusion equation is known to express-via the Continuous Time Random Walk (CTRW) generalization of Brownian motion-the underlying traps and jumps in the movement of restricted and hindered tissue water. And, while the short gradient pulse duration sequence (∆ > δ) is often applied in laboratory NMR and animal MRI, the larger gradient coils used in human scanners are not currently able to meet this condition. Also, complete mapping of the diffusion signal {g, ∆, δ} carries a time penalty that precludes exploration of a wide range of values in a typical clinical scanning protocol. Hence, an underlying theme in our work is to optimize imaging time, model selection and full data space characterization.
What are the trade-offs between the respective fractional order models (decay functions) considered in the paper? In general, as the b-value increases beyond 1000 s/mm 2 (e.g., varying g, for fixed ∆ and δ) the diffusion regime moves from the domain of hindered diffusion (where diffusion is Gaussian, and the decay is exponential) to one of restricted diffusion (where the diffusion is non-Gaussian, and the decay is often multi-exponential). In the Gaussian case, S(b) = S 0 exp[−(bD 0 )], we note the reduction in the apparent diffusion coefficient from its value of 2.4 × 10 −3 mm 2 /s in cerebral spinal fluid to 0.8 × 10 −3 mm 2 /s in brain gray matter by a tortuosity factor of 3, so the apparent diffusion coefficient, ADC = D CSF /3. For intermediate b-values (1000-3000 s/mm 2 ), the stretched exponential decay function, S(b) = S 0 exp[−(bD 0 ) α ], can be used, but so can a bi-exponential model, a kurtosis model or the varying diffusion coefficient (VDC) model where D(b) = D 0 exp(−bD 1 ) [14,39]. When the b-values increase above 3000 s/mm 2 , the exponential models all rapidly decay while the signal often persists due to trapped tissue water, which exhibits an inverse power law decrease in signal intensity. This behavior can be captured by the VDC model, but it is also naturally encoded in the Mittag-Leffler decay, S(b) = S 0 E α [−(bD 0 ) α ], which for low b-values approximates the stretched exponential, and for high b-values decays asymptotically as (bD 0 ) −α . The Kilbas-Saigo generalization of the Mittag-Leffer function is the mother function (model) for all three cases, and has been applied before by Hanyga and Seredynska [34] to diffusion in MRI and by Lin to problems in NMR [40,41]. The behavior of the Kilbas-Saigo modeling approach for the case of varying ∆ (for fixed g and δ) has not yet been fully examined. Here, we were successful in capturing the anisotropy of the bovine optic nerve and the emergence of anomalous diffusion as the diffusion time increases. These results are similar to the behavior of the q-space model, S(b) = S 0 E α [−D α,β q β (∆−δ/3) α ], where q = (γgδ)/2π is a measure of the strength of the diffusion gradient phase encoding [26]. Future such studies would involve normal and tumor tissues and could employ oscillating diffusion gradients for short times and stimulated echo imaging for long times.
The motivation for fractional calculus models of diffusion signal decay in MRI has both a practical and a theoretical basis [42]. The practical emerges from advances in diffusion−weighted imaging sequences that use stronger gradients (b-values, above 2000 s/mm 2 ) to probe sub-voxel tissue compartments. The signal from regions with relatively unhindered diffusion (e.g., cerebral spinal fluid) decays very quickly, while the signal from regions exhibiting pronounced hindered and restricted diffusion persists. Researchers and clinicians fit these signals to a variety of models (e.g., multi-exponential, stretched exponential, and kurtosis) to account for the change in the apparent diffusion coefficient as the b-value increases-fits that provide additional imaging biomarkers of sub-voxel structure [43][44][45]. In this ongoing effort, the medical community seeks to link the model parameters to changes in tissue due to local and diffuse disease. The theoretical basis for considering fractional calculus models emerges from the growing recognition in the biophysics community of the prevalence of anomalous diffusion-due to molecular crowding-in cells and cell membranes [8,9]. Mathematical models of anomalous diffusion include Gaussian models (Brownian and fractional Brownian motion), Langevin equations, the non-Gaussian CTRW, and Lorentz obstructed transport models [46,47]. Within this group, there is a recurrence of power law behavior in the statistics of stochastic models, the correlation coefficients of dynamic models and the waiting time and jump increment probabilities in non-Gaussian models. Since the integral and differential operators of fractional calculus naturally accommodate power laws, via convolution, it is not surprising that anomalous diffusion can be cast in the mathematical language of fractional calculus. Alternatively, the diffusion process can be described by generalizing the diffusion constant into a diffusion coefficient that is a function of either b-value or diffusion time.

Conclusions
In this paper, we described DW-MRI signal decay using sub-diffusion models that can be expressed simply in terms of b-value: exponential, stretched exponential, Mittag-Leffler and Kilbas-Saigo. The stretched Kilbas-Saigo function encompasses this set, and depending on the gradient direction and range of b-values acquired in the experiment, can fit the data with three parameters: (i) D 0 , the apparent diffusion coefficient; (ii) α, the order of the fractional derivative; and (iii) β, the fractional exponent of the inverse power law decay of D(b). The goal was to illustrate how this model can capture the anomalous decay of S(b) observed when data collected over a range of diffusion times. The resulting model parameters encode aspects of the sub-voxel tissue structure, hence offer the potential for describing changes that are brought about by the onset of disease or the success of treatment.
Author Contributions: This work involved all coauthors. R.L.M. conceived the idea, wrote the original draft and contributed to the software, the investigation and the editing. H.K. analyzed the functions, wrote the plotting software and contributed to the visualization. S.W. analyzed the data and investigated the results. Y.L. investigated the functional behavior, supervised the data analysis, and edited the manuscript.
Funding: This research received no external funding.