Fractional Diffusion with Geometric Constraints: Application to Signal Decay in Magnetic Resonance Imaging (MRI)

We investigate diffusion in three dimensions on a comb-like structure in which the particles move freely in a plane, but, out of this plane, are constrained to move only in the perpendicular direction. This model is an extension of the two-dimensional version of the comb model, which allows diffusion along the backbone when the particles are not in the branches. We also consider memory effects, which may be handled with different fractional derivative operators involving singular and non-singular kernels. We find exact solutions for the particle distributions in this model that display normal and anomalous diffusion regimes when the mean-squared displacement is determined. As an application, we use this model to fit the anisotropic diffusion of water along and across the axons in the optic nerve using magnetic resonance imaging. The results for the observed diffusion times (8 to 30 milliseconds) show an anomalous diffusion both along and across the fibers.


Introduction
Mathematical models of diffusion are needed to interpet experimental measurements designed to probe the micro-structure of heterogeneous materials [1][2][3]. For example, identifying the role of diffusion in the electrical response of electrolytic cells, in electrochromic charge transfer in tungsten oxide films [4], in the dispersion of gold nanoparticles in a polymer melt [5], in the neuronal growth on surfaces with controlled geometries [6], and in the movement of proteins in living cells [7,8]. Depending on the interaction between the media and the diffusing particles, these processes manifest different features that promote normal or anomalous diffusion regimes. In the normal case, the mean square displacement exhibits a linear time dependence, i.e., (r − r ) 2 ∼ t, where r is the position coordinate and t the time, which is typical of Markovian processes. In the anomalous case, the anomalous regime is characterized by a nonlinear time dependence of the mean-square displacement, e.g., (r − r ) 2 ∼ t γ , where γ < 1 corresponds to subdiffusion, i.e., a diffusion process slower than the normal one, whereas γ > 1, to superdiffusion-a diffusion process faster than the normal one [9]. This power-law behavior of the mean-square displacement is widely found in disordered structures [10], fractals [11], and percolation clusters [12], with α = 2/d ω , where d ω is the fractal dimension. It is also possible to find behaviors such as (r − r ) 2 ∼ ln γ t, which are related to ultraslow diffusion [13]. Motivated by these anomalous regimes, a comb-like structure has been proposed as the main ingredient of a mathematical diffusive model to investigate anomalous diffusion in percolation clusters with topological bias [14,15]. In this model, according to Ref. [10], the branches of the comb structure play the same role as the dangling ends of the percolation cluster; the backbone of the comb is analogous to the quasilinear structure of the backbone of the cluster. This model describes a diffusive process on a comb-like structure consisting of a "backbone" (a single infinite line in the x-direction) and "branches" (parallel lines in the y-direction that intersect the x-axis). Here the comb model provides a simplified description of the fractal geometry of percolation clusters, where the backbone represents the large bond and the branches are the remaining bonds or "dangling ends" of percolation clusters. It retains the essential properties of diffusion on fractals, with the advantage of providing exact results for such complicated systems such as spiny dendrites (see, for example, the Refs. [16,17]). Moreover, random walks on comb-like structures establish the sojourn times of walkers in the teeth as the underlying mechanism of anomalous diffusion in the backbone.
These features can be observed in biological tissues using magnetic resonance imaging (MRI). The MRI scanner is a device capable of forming both images and spectra from a sample via the application of strong magnetic fields, RF pulses, and gradients, to display the internal structure and composition of the human body [18]. MRI offers various tools that encode water mobility in terms of the magnetic resonance signal and the decay or the relaxation of its diffusion components [19]. Diffusion MRI is used to analyze the complexity the complexity of the orthotropic biological tissues, such as skeletal muscle, tendons, ligaments, and the optic nerve, at the mesoscale. In terms of diffusion, the Bloch-Torrey equation was the first to be considered with which to analyze MRI experiments, followed by several extensions, accomplishing different effects [20].
This approach may also be applied to fractional diffusion equations with geometric constraints, i.e., an anisotropic backbond-like structure, such as the white-matter connections in the brain by bundles of nerve axons [21]. In this paper, we investigate the diffusion process of a system on a backbone structure, such as the one illustrated in Figure 1, where the particles at z = 0 may diffuse on the xy-plane, and, at z = 0 the particles only diffuse along the z-direction. We consider the diffusion processes connected to the following equation: in which ρ(r, t) represents the distribution of particles, with D xy (t) and D z (t) being the kernels related to memory effects, l z is a characteristic thickness, r = (x, y, z), and r xy = (x, y). This equation extends the equations in Refs. [22][23][24][25] used to investigate the diffusion process on a comb-like structure by incorporating fractional derivatives in space and time in both diffusive terms. The kernels may be connected with the waiting-time distributions of random walks and represent a coarse-grained description of the environment's randomness. In particular, the kernels of the time-convoluted operator represent a density memory and not a trajectory memory [26]. The spatial fractional operators ∇ with ψ(r xy , k xy ) = J 0 (k xy r xy ) (J ν (x), the Bessel function of the first kind and order ν [27]) and Both definitions recover the standard differential of integer order for µ xy = 2 and µ z = 2. It is worth mentioning that the first definition of the spatial fractional differential operator, Equation (2), is essentially an extension of the one-dimension spatial fractional operator to the two-dimension case, where the radial symmetry is considered [28]. The second one defined by Equation (3) is essentially the Riesz-Weyl fractional operator. The kernels D xy (t) and D z (t) can be connected with different integro-differential operators, such as the Riemann-Liouville fractional derivative or the others integro-differential operators, depending on the mathematical expression they assume.In our analysis, after performing a general development in the Laplace space, we consider to evaluate the inverse of Laplace transform and analyze the mean square displacement as well as the experimental data obtained for the bovine optical nerve with the technique of MRI [29]. It is worth mentioning that this choice for D z (t) corresponds to the Riemann-Liouville fractional time derivative with 0 < η z < 1. In addition, the choice for the kernels may be related to structure of the media where the diffusion proceeds with barriers, traps, and tortuosity, among others, which may be connected to random walk with a long-tailed distribution for the waiting time distribution. These developments are performed in the next section, Section 2. In Section 3, we present our discussion and conclusions.  (1), diffusing on a comb-like structure.

The Mathematical Problem and Experimental Data
We start our discussion with the model, represented by Equation (1), which will be used to analyze the experimental data analyzed in Refs. [29,30] for the bovine optical nerve. The model is an extension of the standard comb model formed by incorporating the memory effects in space and time through the derivatives of fractional order via convolution kernels. One of the consequences obtained by incorporating these fractional operator is the ability to describe different scenarios related to anomalous diffusion. Equation (1) may be solved by using the Green's function approach [27]. Thus, the Green's function equation to be solved is subjected to the boundary conditions lim r xy →∞ G(r, r , t) = 0, lim r xy →0 ∂ r xy G(r, r , t) = 0, lim z→±∞ G(r, r , t) = 0, and G(r, r , t) = 0 for t ≤ 0. We consider the Green's function expressed in terms of a superposition of the eigenfunction ψ(r xy , k xy ) as follows: with the inverse integral transform given by Note that the choice of the eigenfunction is related to the spatial operator present in Equation (6) for the spatial variable r xy . By applying Equations (8) and (9) in Equation (2), we have where the Laplace on the t variable and Fourier transform on the z variable ae also used to simplify Equation (2). By performing some calculations, it is possible to show that and G(k xy , 0, r , t) = G z,µ z (z, s) By substituting Equation (12) in Equation (9), we obtain By using Equation (13), it is possible to obtain the solution for Equation (1) in the Fourier space when the initial condition ϕ(r) = 1/r xy δ(r xy − r xy )δ(z) is considered. The solution for this case allows us to use the model to analyze the experimental data obtained for the bovine optical nerve in Refs. [29,30] with the technique of MRI. It is given by Note in Equation (14) that the relaxation process on the xy-plane is connected to the behavior of the particles in z-direction. Now, we consider the reduced distribution in each direction in the Fourier space to perform the comparison with the experimental data. For the z-direction, we have and after performing the inverse Laplace transform with D z (s) = D z /s η z , we obtain Equation (16) is expressed in terms of the Mittag-Leffler function, i.e., which recovers the exponential form when α = 1. The behavior of the Mittag-Leffler function for x → −∞ is which implies a power-law behavior in the asymptotic limit. This feature shows that the Equation (16) has an unconventional relaxation process directly connected to the choice of th-kernel D z (t). The case µ z = 2 in Equation (16) may be linked to the Lévy like distributions. The random walk connected to Equation (16) exhibits a mixing between two different regimes for η z = 1 and µ z = 2. For the distribution on the plane, we have ρ(k xy , 0, s) = 1 and, consequently, with Equation (20) has the same form as Equation (16), i.e., it is expressed in terms of the Mittag-Leffler function. However, the indexes present in Equation (20) are different from Equation (16) and, in particular, also depend on η z and µ z , which are connected with the processes occurring along the z-direction. The reduced distribution related to this case, i.e., ρ xy (k xy , t), may also be connected to a fractional diffusion equation in space and time. This feature suggests that the stochastic processes behind these distributions imply a random walks with long-tailed distributions for the waiting time and the jumping probabilities and mixing between the different regimes related to each case. In particular, the inverse Fourier transform of Equations (19) and (20) yields (see Figure 2). For the xyplane, the distribution is given by ρ xy (r xy , t) = 1 µ xy r xy D xy t (see Figure 3). The mean-square displacement may be obtained from the previous equations by using the scaling arguments, as done in Ref. [9] for different fractional diffusion equations. For Equations (21) and (22), it is possible to show that and ρ xy (r xy , t) = 1 Thus, for each case represented by Equations (21) and (22), for the previous initial condition, we have σ 2 z ∝ t 2η z /µ z and σ 2 xy ∝ t 2η xy /µ xy . ρ z (z, t) = (D z t η z ) 1/µ z ρ z (z, t) and z = z/(D z t η z ) 1/µ z . The black dashed-dotted line corresponds to the case η z = 1 and µ z = 3/2. The red solid line corresponds to the case η z = 4/5 and µ z = 2 and the blue dotted line considers the case η z = 4/5 and µ z = 3/2.  where ρ xy (r xy , t) = (D xy t η xy ) 2/µ xy ρ xy (r xy , t) and r xy = r xy /(D xy t η xy ) 1/µ xy . The black solid line corresponds to the case η z = 1 and µ z = 3/2. The red dashed-dotted line corresponds to the case η z = 4/5 and µ z = 2 and the blue dotted line considers the case η z = 4/5 and µ z = 3/2.
Let us apply these results, with µ z = 2 and η xy = 0, to the experimental data published in Refs. [29,30] for the bovine optic nerve for different values of ∆ and δ = 3 ms. We may associate the experimental scenario described in Ref. [29] with a comb-like structure by considering that the z-direction corresponds to the parallel direction along the optic nerve and the xy-plane with the perpendicular direction in a slice of the fiber bundle. Hence, the reduced distributions in the Fourier space for each scenario are connected with the attenuation functions S(b)/S(0), in terms of b and an apparent diffusion coefficient. In addition, we also adapt the previous results to analyze the experimental data obtained with the MRI. We consider the relation between the parameters as follows for the z-direction: and b xy = |k xy | 2 ∆ (where k xy → γGδ and ∆ = ∆ − δ/3). In these equations, γ is the gyromagnetic ratio, G is the amplitude of diffusion gradient, ∆ is the observation time between gradient pulses, and δ is the pulse length. Thus, the attenuation for each case is given by and Equations (25) and (26) extend the standard scenario characterized by an exponential relaxation, by incorporating different aspects related to anomalous diffusion and the geometric constraint imposed by the comb model. To proceed, we first apply this model to the experimental data of the bovine optic nerve with ∆ = 8 ms. Figure 4 shows the curve fits, represented by Equations (25) and (26), and the experimental data.
They are in a good agreement for this case, which demonstrates an anomalous relaxation for the attenuation functions different from the exponential behavior obtained with the standard approach. Next, we analyzed the experimental data for the case in which ∆ = 30 ms. The experimental data and the model are shown in Figure 5 and, as in the case of Figure 4, a good agreement between the experimental data and the model is observed. We also consider the standard form of the attenuation function to illustrate that the experimental data have aspects that are suitable described in terms of the usual diffusion. This feature implies that the system has characteristics in addition to normal diffusion and, consequently, opens the possibility of considering different approaches. In Tables 1 and 2, we show the results for the parameters in the vertical and parallel directions and the estimate variance (equivalent to the squared sum of fit residuals divided by the degrees of freedom n − p where n is the length of the data set and p is the number of parameters), when different values of ∆ are considered.  . Behavior of the model (lines, green, and red) and the experimental data (symbols, square and circles), when ∆ = 8 ms and δ = 3 ms. We observe that the model captures the experimental behavior for both. The parameters values obtained in the fitting process by using the Mathematica function NonlinearModelFit for each case are µ xy /2 = 0.96, η xy = 0.60, l z = 6.2 × 10 −2 mm; D xy = 2.90 × 10 −4 mm 2 /s, error = 0.61 × 10 −4 (estimated variance) for S xy and ηz = 0.76; and D z = 7.30 × 10 −4 mm 2 /s, error = 1.7 × 10 −4 (estimated variance) for S z . For each case, by using the scaling argument as discussed above, the mean square displacement for each case is σ 2 xy ∝ t 0.28 and σ 2 z ∝ t 0.76 . We also consider the standard form to the attenuation function given by in terms of an exponential, i.e., S z(xy) (b) = e −D z(xy) b . It corresponds to the dashed-dotted-dotted cyan lines. For S z , we have that D z = 6.32 × 10 −4 mm 2 /s and error = 2.22 × 10 −3 (estimated variance). For the other direction, S xy , we have that D xy = 7.30 × 10 −4 mm 2 /s and error = 3.79 × 10 −3 (estimated variance).

Discussion and Conclusions
We have analyzed fractional anisotropic diffusion subject to the comb-like structure defined in Figure 1. The particles were restricted to move in the xy−plane for z = 0. Above and below this plane, the particles move in the z-direction. These conditions have a direct influence on the diffusion of the particles by coupling particle displacement in the z− direction with the movement in the xy−-lane. Another aspect of this model is the presence of fractional order convolution kernels in the diffusive terms. These introduce memory effects associated with the long-tailed behavior in the waiting time distribution and the spatial correlations connected to the long tailed behavior of the jumping distributions. These features yield fractional order relaxation and diffusion processes that are different from the standard case.
The models can be applied to the attenuation of the function S(b)/S(0), used to model the experimental data obtained with the technique of MRI [31]. Thus, we first connect the experimental scenario presented in [29] for the optical nerve with structure defined in Figure 1. In this case, the z direction represents the direction along the optic nerve and the xy plane with the vertical direction or the perpendicular direction. The attenuation function for each direction may be obtained by considering the reduced distributions. We have used them to fit the experimental data as shown in Figures 2 and 3. Tables 1 and 2 show the parameter values for different experimental data obtained for the bovine optical nerve [29]. In this case, we observe a subdiffusive behavior for all scenarios, which may be related to the structural organization of the optic nerve. We underline that different models have been used to model the experimental data obtained from the MRI technique. For example, in Refs. [29,32] the authors consider attenuation of the function S(b)/S(0) in terms of the Kilbas-Sago function, which is the solution of the anomalous diffusion equation expressed in terms of fractional derivatives in time and space. In this case, however, proceed independently in along and across the bundles of axons in the optic nerve. On the other hand, the spatial structure presented in the comb model couples the different directions, which influences the diffusion of the particles. It is also possible to consider a varying diffusion coefficient, D(b), or to consider a bi-exponential model to fit S(b)/S(0). Application of these models [33,34] to optic nerve data (see Supplemental Material) have also provided good fits to this data.
The main purpose of these models is to capture the experimental behavior observed with the MRI technique and to use these results to explain the behavior of the particles in the intra-and extracellular space. Thus, the comb model uses the known structure of the optic nerve and establishes a dynamic connection between a preferential direction, which is coupled to the other directions. This structure permits us to model scenarios in which the diffusion in one direction influences diffusion in another direction. For the experimental scenario described here, the fiber orientation may be identified with the z direction and the fiber cross sections with the xy plane, where the reduced distribution for each case was identified with the attenuation function S(b)/S(0). The agreement between the theory and the experimental data was satisfactory when the comb model was considered. These features suggest that the comb model (or extension in this direction) brings us the possibility of investigating how the structural aspects of the system influence the appearance of anomalous diffusion.