G-Subdiffusion Equation as an Anomalous Diffusion Equation Determined by the Time Evolution of the Mean Square Displacement of a Diffusing Molecule

1. Introduction

Anomalous diffusion models are often based on the assumption of constant parameters, which means that the structure of the medium remains constant over time. A frequently used model to derive the normal and anomalous diffusion equations is the continuous time random walk (CTRW) model [1,2,3,4,5]. Within this model, random walk of a single molecule is considered. The process is described by the probability density (Green’s function) $P(x,t|x_0)$ of finding the molecule at point x at time t, $x_0$ being the initial position of the molecule. In a one-dimensional unbounded homogeneous system (such a system is considered in this paper) the mean square displacement (MSD) of a molecule ${\sigma }^2$ is calculated as

$${\sigma }^2\left(t\right)={\int }_{-\infty }^{\infty }{(x-x_0)}^{2}P(x,t|x_0)dx.$$

(1)

The relation

$${\sigma }^2\left(t\right)=\beta t^{\alpha }$$

(2)

is often used to define the type of diffusion. For $0<\alpha <1$ we have subdiffusion, then

$$\beta ={\displaystyle \frac{2D_{\alpha }}{\Gamma (1+\alpha )}},$$

(3)

$D_{\alpha }$ is the subdiffusion coefficient measured in the units of $m^2/s^{\alpha }$. Subdiffusion occurs in a system in which the movement of diffusing molecules is very hindered, as occurs in gels, porous media, and bacterial biofilms [2,6,7,8,9,10]. For $\alpha =1$ we have normal diffusion with $\beta =2D_1$. When $\alpha >1$ we have superdiffusion (facilitated diffusion), the examples are diffusion in turbulent media and in random velocity fields [11,12,13,14], cell migration in some biological processes [15], movement of endogeneous intracellular particles in some pathogens [16], and mussels movement [17]. However, the CTRW model gives $\beta =\infty$. The superdiffusion model, based on the g-subdiffusion equation, which provides $\beta <\infty$ is shown in Refs. [18,19]. The relation (2) is complemented by the following relation defining ultraslow diffusion (slow subdiffusion)

$${\sigma }^2\left(t\right)\sim v\left(t\right),$$

(4)

where v is a slowly varying function that satisfies the condition $v\left(at\right)/v\left(t\right)\to 1$ when $t\to \infty$ for any $a>0$. In practice, v is a combination of logarithmic functions. Ultraslow diffusion is an extremely slow process, qualitatively different from ordinary subdiffusion. This process has been observed in diffusion of water in aqueous sucrose glasses [20] and languages dynamics [21]. Superdiffusion is usually described by an equation with a fractional derivative with respect to the spatial variable, while subdiffusion is described by an equation with a fractional derivative with respect to time. Ultraslow diffusion is described by integro–differential equations with integral operators which are not usually identified as fractional time derivatives [22,23,24,25,26,27,28].

Differential equations mentioned above with constant parameters are used to describe diffusion in a medium which properties does not change with time. However, diffusion parameters depend on the interaction of diffusing molecules with the environment and on a structure of the medium, both can change over time. The single molecule tracking method allows for the experimental determination of ${\sigma }^2\left(t\right)$. We mention that there have been used other power functions with respect to time, experimentally measurable, from which subdiffusion parameters can be determined. An example of this is the time evolution of the so-called thickness of the membrane layer [8]. Processes with a time-varying diffusion exponent have been observed in bacterial motion on small beads in a freely suspended soap film [29], in transport of colloidal particles between two parallel plates [30], microspheres in a living eukaryotic cell [31], endogenous lipid granules in living yeast cells [32], and in the diffusion of passive molecules in the active bath with moving particles [33,34]. According to the Stokes-Einstein formula, $D_1\sim T$, the change of the temperature T of a liquid generates the change in the diffusion coefficient. It was found that the diffusion coefficient of chloride ions in concrete shows a dependence on time [35,36]. Other examples are $D_{\alpha }\left(t\right)\sim log\left(t\right)$ which is caused by the aging process of a complex system in which anomalous diffusion occurs [37] and $D_{\alpha }\left(t\right)=a{\mathrm{e}}^{\pm 2bt}$, where a and b are constant parameters [38]. Experimental study provided that the diffusion model with power-law $D_{\alpha }\left(t\right)$ well describes water diffusion in brain tissues [39,40]. The function $D_{\alpha }\left(t\right)$ can have a more complicated form, e.g. it may contain an oscillatory component. Such a dependency can occur, for example, in diffusion of antibiotics in bacterial biofilms. Bacteria activate various defense mechanisms against the action of antibiotics [41,42]. One such mechanism is the thickening of the biofilm, which significantly impedes diffusion of the antibiotic and reduces the diffusion coefficient. Slowing down diffusion of the antibiotic reduces the risk of the antibiotic having an effective effect on the bacteria. This process can cause a relaxation of the defense mechanisms of bacteria and increase diffusion of the antibiotic. Bacteria, feeling a greater threat, intensify their defense mechanisms again, and so on. The subdiffusion coefficient of the antibiotic in the biofilm may then undergo periodic changes. Diffusion coefficients that oscillate with time, reflecting complex memory or frictional effects in the system, have been considered in some fractional diffusion models [43].

When the diffusion parameters are not constant, various equations have been used to describe the diffusion processes such as subdiffusion equations with a fractional time derivative of the order depending on time (and on a spatial variable) [44,45,46,47,48,49,50,51,52,53,54,55], and the equation with a linear combination of fractional time derivatives of different orders [56].

If the function ${\sigma }^2\left(t\right)$ has a complicated form, the question arises as to what equation describes the subdiffusion process and whether there are methods for solving such an equation. We will show that by using the g-subdiffusion equation with an appropriately chosen function g one can describe the diffusion process defined by a time-increasing function ${\sigma }^2\left(t\right)$. The Green’s function for this equation generates the assumed function ${\sigma }^2\left(t\right)$. The g-subdiffusion equation contains the Caputo derivative with respect to the function g [57,58]. This equation can be solved using the Laplace transform with respect to the function g (the g-Laplace transform).

2. G-Subdiffusion Equation

The CTRW model provides the following ordinary subdiffusion equation

$${\displaystyle \frac{C{\partial }^{\alpha }P(x,t|x_0)}{\partial t^{\alpha }}}=D_{\alpha }{\displaystyle \frac{{\partial }^{2}P(x,t|x_0)}{\partial x^2}},$$

(5)

where $0<\alpha <1$, the Caputo fractional derivative is defined for $0<\alpha <1$ as

$${\displaystyle \frac{C{d}^{\alpha }f\left(t\right)}{d{t}^{\alpha }}}={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{(t-u)}^{-\alpha }{f}^{\left(1\right)}\left(u\right)du,$$

(6)

$f^{\left(1\right)}\left(u\right)=df\left(u\right)/du$. Formally, the normal diffusion equation

$${\displaystyle \frac{\partial P(x,t|x_0)}{\partial t}}=D_1{\displaystyle \frac{{\partial }^{2}P(x,t|x_0)}{\partial x^2}}$$

(7)

can be treated as a special case of the subdiffusion equation for $\alpha =1$. We mention that Eq. (5) can be transformed to its equivalent form with the fractional Riemann–Liouville time derivative of the order $1-\alpha$ [2,3,4,7],

The g–subdiffusion equation can be interpreted as a modified form of Eq. (5). The modification consists of changing the time variable t to a function $g\left(t\right)$, $t\mapsto g\left(t\right)$, where $g\left(t\right)$ is given in units of time and meets the conditions $g\left(0\right)=0$, $g(\infty )=\infty$, and $g^{\left(1\right)}\left(t\right)>0$. The g–subdiffusion equation is [57,58]

$${\displaystyle \frac{C{\partial }_g^{\alpha }P(x,t|x_0)}{\partial t^{\alpha }}}=D_{\alpha }{\displaystyle \frac{{\partial }^{2}P(x,t|x_0)}{\partial x^2}},$$

(8)

where the g-Caputo fractional derivative of the order $\alpha$ with respect to the function g is defined for $0<\alpha <1$ as [59]

$${\displaystyle \frac{C{d}_g^{\alpha }f\left(t\right)}{d{t}^{\alpha }}}={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{(g\left(t\right)-g\left(u\right))}^{-\alpha }{f}^{\left(1\right)}\left(u\right)du.$$

(9)

When $g\left(t\right)=t$, the g-Caputo fractional derivative takes the form of the ordinary Caputo derivative. The Green’s function $P(x,t|x_0)$ is a solution to a subdiffusion equation for the initial conditions $P(x,0|x_0)=\delta (x-x_0)$, where $\delta$ is the delta–Dirac function. In an unbounded region the boundary conditions are assumed to be $P(\pm \infty ,t|x_0)=0$.

The g–subdiffusion equation can be solved by means of the g–Laplace transform method. The g-Laplace transform is defined as [60,61]

$${\mathcal{L}}_g\left[f\left(t\right)\right]\left(s\right)={\int }_0^{\infty }{\mathrm{e}}^{-sg\left(t\right)}f\left(t\right)g^{\prime }\left(t\right)dt.$$

(10)

The g–Laplace transform is related to the ordinary Laplace transform $\mathcal{L}\left[f\left(t\right)\right]\left(s\right)={\int }_0^{\infty }{\mathrm{e}}^{-st}f\left(t\right)dt$,

$${\mathcal{L}}_g\left[f\left(t\right)\right]\left(s\right)=\mathcal{L}\left[f\left(g^{-1}\left(t\right)\right)\right]\left(s\right).$$

(11)

Eq. (11) provides the relation

$${\mathcal{L}}_g\left[f\left(t\right)\right]\left(s\right)=\mathcal{L}\left[h\left(t\right)\right]\left(s\right)\iff f\left(t\right)=h\left(g\left(t\right)\right).$$

(12)

The above formula is helpful in calculating the inverse g–Laplace transform if the inverse ordinary Laplace transform is known. For example, since ${\mathcal{L}}^{-1}[1/s^{1+\nu }]\left(t\right)=t^{\nu }/\Gamma (1+\nu )$, $\nu >-1$, and ${\mathcal{L}}^{-1}\left[s^{\nu }{\mathrm{e}}^{-a{s}^{\mu }}\right]\left(t\right)=t^{-1-\nu }{\sum }_{k=0}^{\infty }(1/k!\Gamma (-\nu -\mu k)){(-a/t^{\mu })}^k$, $a,\mu >0$, we get [62,63]

$${\mathcal{L}}_g^{-1}\left[{\displaystyle \frac{1}{s^{1+\nu }}}\right]\left(t\right)={\displaystyle \frac{g^{\nu }\left(t\right)}{\Gamma (1+\nu )}},\nu >-1,$$

(13)

and

$$\begin{array}{c}\hfill {\mathcal{L}}_g^{-1}\left[s^{\nu }{\mathrm{e}}^{-a{s}^{\mu }}\right]\left(t\right)={\displaystyle \frac{1}{g^{1+\nu }\left(t\right)}}\sum _{k=0}^{\infty }{\displaystyle \frac{1}{k!\Gamma (-\nu -\mu k)}}{\left(-{\displaystyle \frac{a}{g^{\mu }\left(t\right)}}\right)}^k\equiv f_{\nu ,\mu }(g\left(t\right);a),\end{array}$$

(14)

$a,\mu >0$. The function $f_{\nu ,\mu }$ is a special case of the Wright function and the H-Fox function.

The calculations for solving Eq. (8) by means of the g–Laplace transform method are similar to those for solving Eq. (5) using the ordinary Laplace transform. Due to the relation [60,61]

$${\mathcal{L}}_g\left[{\displaystyle \frac{C{d}_g^{\alpha }f\left(t\right)}{d{t}^{\alpha }}}\right]\left(s\right)=s^{\alpha }{\mathcal{L}}_g\left[f\left(t\right)\right]\left(s\right)-s^{\alpha -1}f\left(0\right),$$

(15)

where $0<\alpha \le 1$, the g–Laplace transform of Eq. (8) reads

$$\begin{array}{c}\hfill s^{\alpha }{\mathcal{L}}_g\left[P(x,t|x_0)\right]\left(s\right)-s^{\alpha -1}P(x,0|x_0)=D_{\alpha }{\displaystyle \frac{{\partial }^2{\mathcal{L}}_g\left[P(x,t|x_0)\right]\left(s\right)}{\partial x^2}}.\end{array}$$

(16)

The g–Laplace transform of Green’s function, that is the following solution to Eq. (16) for the boundary conditions ${\mathcal{L}}_g\left[P(\pm \infty ,t|x_0)\right]\left(s\right)=0$, is

$${\mathcal{L}}_g\left[P(x,t|x_0)\right]\left(s\right)={\displaystyle \frac{1}{2\sqrt{D_{\alpha }}{s}^{1-\alpha /2}}}{\mathrm{e}}^{-{\displaystyle \frac{|x-x_0|}{\sqrt{D_{\alpha }}}}{s}^{\alpha /2}}.$$

(17)

From Eqs. (14) and (17) we obtain

$$P(x,t|x_0)={\displaystyle \frac{1}{2\sqrt{D_{\alpha }}}}{f}_{-1+\alpha /2,\alpha /2}\left(g\left(t\right);{\displaystyle \frac{|x-x_0|}{\sqrt{D_{\alpha }}}}\right).$$

(18)

The g–Laplace transform of Eq. (1), and Eqs. (13), (17), provide ${\mathcal{L}}_g\left[{\sigma }^2\left(t\right)\right]\left(s\right)=2D_{\alpha }/s^{1+\alpha }$. Finally, we get

$${\sigma }^2\left(t\right)={\displaystyle \frac{2D_{\alpha }{g}^{\alpha }\left(t\right)}{\Gamma (1+\alpha )}}.$$

(19)

3. Time Evolution of ${\sigma }^2$ as a Function Defining the Diffusion Process

As mentioned in Sec. Section 1, the function ${\sigma }^2\left(t\right)$ is usually used as a definition of the type of diffusion. The function is experimentally measurable. The single particle tracking method is used when a random walk of a single molecule is observed [64,65,66,67]. The function ${\sigma }^2\left(t\right)$ can also be determined in another way, e.g. by studying the release of a substance from one vessel to another through a thin membrane. At the initial moment, the vessel A contains a homogeneous solution of the diffusing substance with an initial concentration $C_0$, the subdiffusion parameters in the vessel are $\alpha$ and $D_{\alpha }$, and the vessel B contains a pure solvent. When the membrane allows free passage from the vessel A to B, and the molecule return passage is practically impossible, then the total amount of the substance N in the vessel B evolves in time as $N\left(t\right)=\lambda t^{\alpha /2}$, where $\lambda =C_0\sqrt{D_{\alpha }}/\Gamma (1+\alpha /2)$[68]. Combining these equations with Eqs. (2) and (3) we get ${\sigma }^2\left(t\right)=2{\Gamma }^2(1+\alpha /2)N^2\left(t\right)/\left[C_0^2\Gamma (1+\alpha )\right]$. Another method is to measure the temporal evolution of the thickness of near-membrane layer $\rho \left(t\right)$. It is defined as the distance from the membrane to the point, where the substance concentration drops k times with respect to the membrane surface in the vessel B. We get ${\sigma }^2\left(t\right)=\lambda {\rho }^2\left(t\right)$, where $\lambda$ is controlled by $\alpha$, $D_{\alpha }$, and k, see [8].

In the case where ${\sigma }^2$ evolves in time according to Eq. (2), the equations describing the diffusion process are known. However, as mentioned in Sec. Section 1, more complicated forms of ${\sigma }^2\left(t\right)$ are possible for which the equation describing the process may not be known. Based on our considerations in Sec. Section 2, we conclude that such a process can be described by the g-subdiffusion equation Eq. (8), in which

$$g\left(t\right)={\left({\displaystyle \frac{\Gamma (1+\alpha ){\sigma }^2\left(t\right)}{2D_{\alpha }}}\right)}^{1/\alpha },$$

(20)

see Eq. (19). The Green’s function is given by Eq. (18). When the molecules diffuse independently of each other, the concentration C of the diffusing molecules can be calculated using the formula

$$C(x,t)={\int }_{-\infty }^{\infty }P(x,t|x_0)C(x_0,0)d{x}_0.$$

(21)

Then, the g-subdiffusion equation is also satisfied by the function $C(x,t)$.

As an example, we consider four cases of the function ${\sigma }^2\left(t\right)$, two of which contain an oscillating component.

• For

  $${\sigma }_1^2\left(t\right)={\displaystyle \frac{2D_{\alpha }{t}^{\alpha }{(1+at)}^{\kappa }}{\Gamma (1+\alpha )}},$$

  (22)

  we have

  $$g_1\left(t\right)=t{(1+at)}^{\kappa /\alpha }.$$

  (23)
• When

  $${\sigma }_2^2\left(t\right)={\displaystyle \frac{2D_{\alpha }{t}^{\alpha }}{\Gamma (1+\alpha ){(1+at)}^{\kappa }}},$$

  (24)

  we get

  $$g_2\left(t\right)={\displaystyle \frac{t}{{(1+at)}^{\kappa /\alpha }}}.$$

  (25)
• When

  $${\sigma }_3^2\left(t\right)={\displaystyle \frac{2D_{\alpha }{t}^{\alpha }{(1+at+b\sin \left(\omega t\right))}^{\kappa }}{\Gamma (1+\alpha )}},$$

  (26)

  there is

  $$g_3\left(t\right)=t{(1+at+b\sin \left(\omega t\right))}^{\kappa /\alpha }.$$

  (27)
• For

  $${\sigma }_4^2\left(t\right)={\displaystyle \frac{2D_{\alpha }{t}^{\alpha }}{\Gamma (1+\alpha ){(1+at+b\sin \left(\omega t\right))}^{\kappa }}},$$

  (28)

  we get

  $$g_4\left(t\right)={\displaystyle \frac{t}{{(1+at+b\sin \left(\omega t\right))}^{\kappa /\alpha }}}.$$

  (29)

In the above equations it is assumed that $D_{\alpha },\kappa ,\omega >0$, $\alpha \in (0,1)$, and $a>b\omega >0$.

The time evolutions of the mean square displacement of diffusing particle ${\sigma }^2\left(t\right)$ are given in Figs. Figure 1 and Figure 2. Figs. Figure 3 and Figure 4 show the Green’s function plots for times $t=1$ and $t=10$, respectively. The plots of the functions are compared with the functions obtained for ordinary subdiffusion with constant parameters $\alpha$ and $D_{\alpha }$, for which $g\left(t\right)\equiv t$ and ${\sigma }^2\left(t\right)=2D_{\alpha }{t}^{\alpha }/\Gamma (1+\alpha )$, that are marked with a thick solid line without symbols. All plots are made for the function $g_i$, $i=1,2,3,4$, defined by Eqs. (23), (25), (27), and (29), respectively; the numbering of other functions is consistent with the numbering of the function $g_i$. The plots are made for $D_{\alpha }=10$, $\alpha =0.7$, $a=2$, $b=1$, $\omega =1$, and $x_0=0$, the values of all parameters are given in arbitrarily chosen units. In each case two values of $\kappa$, namely $\kappa =0.2$ and $\kappa =0.5$, are considered. Comparing ${\sigma }_1^2$ and ${\sigma }_2^2$ with ${\sigma }_3^2$ and ${\sigma }_4^2$, respectively, in Fig. Figure 2 we see how the oscillatory effect changes ${\sigma }^2\left(t\right)$. The effect, involved in the functions ${\sigma }_3^2$ and ${\sigma }_4^2$, is visible for relatively short times. This fact can also be seen by analyzing the plots of Green’s functions. In Fig. Figure 3, for $t=1$, the plots of the functions generated by ${\sigma }_1^2$ and ${\sigma }_3^2$ differ from each other, and the same applies to the functions generated by ${\sigma }_2^2$ and ${\sigma }_4^2$. These differences almost disappear in Fig. Figure 4 for $t=10$.

4. Final Remarks

We have used the g-subdiffusion equation to describe diffusion process defined by the function ${\sigma }^2\left(t\right)$. This equation contains the time fractional Caputo derivative with respect to an increasing function $g\left(t\right)$. As $g\left(t\right)\equiv t$, the g-subdiffusion equation becomes the ordinary subdiffusion equation. The function $g\left(t\right)$ given by Eq. (20) provides an equation describing the process characterized by a time-increasing function ${\sigma }^2\left(t\right)$ that can be determined experimentally. The presented model confirms the usefulness of the g-subdiffusion equation in modeling anomalous diffusion processes.

The function ${\sigma }^2\left(t\right)$ is criticized as not defining the type of diffusion unambiguously. As shown in Ref. [69], the appropriate combination of superdiffusion and subdiffusion effects provides the relation ${\sigma }^2\left(t\right)\sim t$, but such a process cannot be interpreted as normal diffusion. However, ${\sigma }^2\left(t\right)$ is an important characteristic of the diffusion process. In our considerations, we use this function as defining the diffusion process, without defining the type of diffusion based on Eq. (2).

For diffusion processes defined by Eq. (2), the diffusion equations are derived from the stochastic CTRW model. The g-subdiffusion equation also has its stochastic interpretation [58]. It can be derived from the modified CTRW model. The interpretation of this equation is based on the assumption that the time variable t is replaced by a time-increasing function $g\left(t\right)$ with the additional assumptions $g\left(0\right)=0$ and $g(\infty )=\infty$. If the function ${\sigma }^2\left(t\right)$ has the same properties, then, together with the parameters $\alpha$ and $D_{\alpha }$, the function defines $g\left(t\right)$. Then, the process can be described by the g-subdiffusion equation. The parameters $\alpha$ and $D_{\alpha }$ are determined from additional considerations. For example, the drug transport in a system containing densely packed gel beads is described by the g-subdiffusion equation, which is confirmed by empirical studies [68]. However, in the initial time interval, this process is well described by the ordinary subdiffusion equation Eq. (5), for which the parameters can be determined.