Generalised fractional diffusion equations for subdiffusion on arbitrarily growing domains

TL;DR

This paper derives new fractional partial differential equations to model subdiffusion on arbitrarily growing domains, incorporating domain growth dynamics and a novel comoving fractional derivative, based on a Continuous Time Random Walk framework.

Contribution

It introduces a new fractional derivative and formulates generalized equations for subdiffusion on growing domains, addressing challenges of history dependence.

Findings

Derived fractional PDEs for subdiffusion on growing domains

Introduced a new comoving fractional derivative

Framework applicable to various physical phenomena involving domain growth

Abstract

Many physical phenomena occur on domains that grow in time. When the timescales of the phenomena and domain growth are comparable, models must include the dynamics of the domain. A widespread intrinsically slow transport process is subdiffusion. Many models of subdiffusion include a history dependence. This greatly confounds efforts to incorporate domain growth. Here we derive the fractional partial differential equations that govern subdiffusion on a growing domain, based on a Continuous Time Random Walk. This requires the introduction of a new, comoving, fractional derivative.

Equations

$$y=\underset{n\to \infty }{\lim }\sum _{i=1}^n\delta y_i={\int }_0^x\exp \left({\int }_0^t\mu (z,s)ds\right)dz=g(x,t).$$

$${\nu }^{\ast }(x,t)=\frac{\partial g(x,t)}{\partial x}=e^{{\int }_0^t\mu (x,s)ds},$$

$${\eta }^{\ast }(x,t)=\frac{\partial g(x,t)}{\partial t}={\int }_0^x\mu (z,t)e^{{\int }_0^t\mu (z,s)ds}dz.$$

$$\nu (y,t)={\nu }^{\ast }(g^{-1}(y,t),t),$$

$$\eta (y,t)={\eta }^{\ast }(g^{-1}(y,t),t).$$

$$\delta y(x,t)=e^{{\int }_0^t\mu (x,s)ds}\delta x.$$

$$\rho (g(x,t),t)\delta y(x,t)={\int }_0^t\Phi (t-t^{\prime })q(g(x,t^{\prime }),t^{\prime })\delta y(x,t^{\prime })d{t}^{\prime },$$

$$\begin{array}{rl}\rho (g(x,t),t)e^{{\int }_0^t\mu (x,s)ds}& ={\int }_0^t\Phi (t-t^{\prime })q(g(x,t^{\prime }),t^{\prime })e^{{\int }_0^{t^{\prime }}\mu (x,s)ds}d{t}^{\prime }.\end{array}$$

$$\rho (y,t)=\rho (g(x,t),t)={\rho }^{\ast }(x,t),q(g(x,t),t)=q^{\ast }(x,t).$$

$${\rho }^{\ast }(x,t)e^{{\int }_0^t\mu (x,s)ds}={\int }_0^t\Phi (t-t^{\prime })q^{\ast }(x,t^{\prime })e^{{\int }_0^{t^{\prime }}\mu (x,s)ds}d{t}^{\prime }.$$

$$\frac{\partial {\rho }^{\ast }(x,t)}{\partial t}=q^{\ast }(x,t)-{\int }_0^t\psi (t-t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}{q}^{\ast }(x,t^{\prime })d{t}^{\prime }-\mu (x,t){\rho }^{\ast }(x,t).$$

$$i^{\ast }(x,t)={\int }_0^t\psi (t-t^{\prime })q^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }.$$

$$\begin{array}{rl}{q}^{\ast }(x,t)& ={\int }_0^{L(0)}\lambda (x|x^{\prime },t)i^{\ast }(x^{\prime },t)d{x}^{\prime }.\end{array}$$

$$\begin{array}{rl}\frac{\partial {\rho }^{\ast }(x,t)}{\partial t}& ={\int }_0^{L(0)}\lambda (x|x^{\prime },t){\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x^{\prime },t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x^{\prime },s)ds}d{t}^{\prime }d{x}^{\prime }\\ & -{\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }-\mu (x,t){\rho }^{\ast }(x,t).\end{array}$$

$${\mathcal{L}}_t\{K(t)\}=\frac{{\mathcal{L}}_t\{\psi (t)\}}{{\mathcal{L}}_t\{\Phi (t)\}},$$

$$\begin{array}{rl}\lambda (x|x^{\prime },t)& =\frac{1}{2}\left(\delta (x-x^{\prime }-{ϵ}^{+})+\delta (x-x^{\prime }+{ϵ}^{-})\right),\end{array}$$

$$\Delta y$$

$$\begin{array}{rl}& \frac{\partial {\rho }^{\ast }(x,t)}{\partial t}=\frac{\Delta y^2{e}^{-2{\int }_0^t\mu (x,s)ds}}{2}(\left(\frac{{\partial }^2}{\partial x^2}{\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }\right)\\ & -\left({\int }_0^t\frac{\partial \mu (x,s)}{\partial x}ds\right)\left(\frac{\partial }{\partial x}{\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }\right))-\mu (x,t){\rho }^{\ast }(x,t)+O(\Delta y^3).\end{array}$$

$$\psi (t)=\frac{t^{\alpha -1}}{{\tau }^{\alpha }}{E}_{\alpha ,\alpha }\left(-{\left(\frac{t}{\tau }\right)}^{\alpha }\right),$$

$$E_{\alpha ,\beta }(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )}.$$

$$K(t)={\mathcal{L}}_s^{-1}\left\{\frac{s^{1-\alpha }}{{\tau }^{\alpha }}\right\}.$$

$${}_0{\mathcal{D}}_t^{1-\alpha }\left(f(t)\right)=\frac{1}{\Gamma (\alpha )}\frac{d}{dt}{\int }_0^{t}f(t^{\prime })(t-t^{\prime }{)}^{\alpha -1}d{t}^{\prime }.$$

$${\mathcal{L}}_t\left\{{}_0{\mathcal{D}}_t^{1-\alpha }\left(f(t)\right)\right\}=s^{1-\alpha }{\mathcal{L}}_t\left\{f(t)\right\}.$$

$$\begin{array}{rl}\frac{\partial {\rho }^{\ast }(x,t)}{\partial t}=& \frac{\Delta y^2{e}^{-2{\int }_0^t\mu (x,s)ds}}{2{\tau }^{\alpha }}(\frac{{\partial }^2}{\partial x^2}\left(\frac{{}_0{\mathcal{D}}_t^{1-\alpha }\left({\rho }^{\ast }(x,t)e^{{\int }_0^t\mu (x,s)ds}\right)}{e^{{\int }_0^t\mu (x,s)ds}}\right)\\ & -\left({\int }_0^t\frac{\partial \mu (x,s)}{\partial x}ds\right)\frac{\partial }{\partial x}\left(\frac{{}_0{\mathcal{D}}_t^{1-\alpha }\left({\rho }^{\ast }(x,t)e^{{\int }_0^t\mu (x,s)ds}\right)}{e^{{\int }_0^t\mu (x,s)ds}}\right))-\mu (x,t){\rho }^{\ast }(x,t)+\mathcal{O}(\Delta y^3).\end{array}$$

$$D_{\alpha }=\underset{\Delta y,\tau \to 0}{\lim }\frac{\Delta y^2}{2{\tau }^{\alpha }},$$

$$\begin{array}{rl}\frac{\partial {\rho }^{\ast }(x,t)}{\partial t}& =D_{\alpha }{e}^{-2{\int }_0^t\mu (x,s)ds}(\frac{{\partial }^2}{\partial x^2}\left(\frac{{}_0{\mathcal{D}}_t^{1-\alpha }\left({\rho }^{\ast }(x,t)e^{{\int }_0^t\mu (x,s)ds}\right)}{e^{{\int }_0^t\mu (x,s)ds}}\right)\\ & -\left({\int }_0^t\frac{\partial \mu (x,s)}{\partial x}ds\right)\frac{\partial }{\partial x}\left(\frac{{}_0{\mathcal{D}}_t^{1-\alpha }\left({\rho }^{\ast }(x,t)e^{{\int }_0^t\mu (x,s)ds}\right)}{e^{{\int }_0^t\mu (x,s)ds}}\right))-\mu (x,t){\rho }^{\ast }(x,t).\end{array}$$

$$\lambda (x|x^{\prime },t)=\delta (x-x^{\prime }+{ϵ}^{-}),$$

$$\lambda (x|x^{\prime },t)=\delta (x-x^{\prime }-{ϵ}^{+}),$$

$$\begin{array}{rl}\frac{\partial {\rho }^{\ast }(x,t)}{\partial t}& =\Delta y{e}^{-{\int }_0^t\mu (x,s)ds}\frac{\partial }{\partial x}\left({\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }\right)\\ & +\frac{\Delta y^2{e}^{-2{\int }_0^t\mu (x,s)ds}}{2}(\left(\frac{{\partial }^2}{\partial x^2}{\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }\right)\\ & -\left({\int }_0^t\frac{\partial \mu (x,s)}{\partial x}ds\right)\left(\frac{\partial }{\partial x}{\int }_0^{t}K(t-t^{\prime }){\rho }^{\ast }(x,t^{\prime })e^{-{\int }_{t^{\prime }}^t\mu (x,s)ds}d{t}^{\prime }\right))-\mu (x,t){\rho }^{\ast }(x,t)+O(\Delta y^3).\end{array}$$

$${\frac{\partial }{\partial x}\left(\frac{{}_0{\mathcal{D}}_t^{1-\alpha }\left({\rho }^{\ast }(x,t)e^{{\int }_0^t\mu (x,s)ds}\right)}{e^{{\int }_0^t\mu (x,s)ds}}\right)|}_{x=0}=0.$$

Full text

A wide range of important physical phenomena involve transport on expanding, and contracting, domains. Fundamental examples include, the diffusion of proteins within growing cells, the interactions of cells in a growing organism, and diffusion in an expanding universe. The governing equations for reaction diffusion on growing domains and related studies of pattern formation have been considered in a series of publications, see for example, CGK1999 ; M2001 ; CM2001 ; BYE2010 ; WBGM2011 ; YBEM2012 ; SSMB2015 ; YAE2016 . Domain growth has been shown to be fundamentally important to the development of patterns KA1995 . Here we consider the problem of subdiffusive transport on a growing domain by constructing a continuous time random walk (CTRW) and limiting to a fractional order partial differential equation (PDE).

Subdiffusion, which is characterised by a sub-linear power-law scaling in time of the mean squared displacement, is common in biological systems with traps and obstacles S2007 , such as diffusion of molecules in spiny nerve cells SWDA2006 , diffusion across potassium channels in membranes MNHMT2010 ; WSTK2011 , and diffusion of HIV virions in cervical mucous SACSA2013 . Subdiffusion is also present in other physical systems such as cosmic rays S2005 , porous media LB2003 , and volcanic earthquakes AS2017 . The generalisation of canonical mathematical diffusion models to incorporate subdiffusive transport, such as, reaction-diffusion PDEs HLW2006 ; SSS2006 ; F2010 ; AYL2010 ; ADH2013mmnp , and Fokker-Planck PDEs BMK2000 ; SK2006 ; HLS2010 ; ADH2013mmnp , has proven non-trivial. In the work below we show that this is also true for subdiffusion on a growing domain.

There are different theoretical approaches that have been used to model subdiffusive transport. One of the more rigorous approaches is to derive the governing equations from the stochastic process of a CTRW MW1965 . The CTRW describes transport of particles on a mesoscopic scale in which particles wait for a time, governed by waiting time probability density, before randomly jumping, governed by a jump length probability density, to another location. If the jump length density is symmetric with a finite variance and the expected waiting time is convergent, then the CTRW limits to the standard diffusion PDE HA1995 ; MK2000 . If the waiting time density is replaced with a heavy tailed power-law waiting time density, then the CTRW limits to a time subdiffusion fractional diffusion PDE HA1995 ; MK2000 .

In the following we start with the underlying stochastic process of a CTRW to derive master equations for subdiffusive transport on a growing domain. In our derivation we first consider a mapping between the position $y$ on the growing domain at any time $t$ and a corresponding position $x$ on the original domain at time $t=0$. With this mapping we then transform the CTRW from the coordinates on the growing domain to a non-growing fixed domain. An auxiliary master equation for the evolution of the density on the fixed domain is derived. The auxiliary master equation is constructed so that the value of the density at a given $x$ and $t$ equates to the probability density on the growing domain for $y$ and $t$. The diffusion limit of the master equation is taken to produce a fractional diffusion equations on both the fixed and growing domains.

Our approach enables us to model subdiffusive transport of particles on arbitrarily growing domains, and the solution of the auxiliary master equation on the fixed domain could be used as the basis for numerical simulations of subdiffusive transport on growing domains. The equations we derive on the growing domain can be interpreted phenomenologically as a reaction sub-diffusion process with an additional advective term. In this context, the reaction represents the dilution of the concentration due to the growing domain.

We wish to construct a mapping between a location on the initial fixed domain, $x\in[0,L_{0}]$, to the corresponding location at some later time $t$, $y\in[0,L(t)]$. To characterise how the domain is changing in time we begin by partitioning the domain $[0,L_{0}]$ into $m$ cells of width $\delta x=\frac{1}{m}$. We will denote the boundary positions of this partition such that $x_{i}=i\delta x$. As the domain grows, the width of the partitions, now denoted by $\delta y_{i}(t)$ will have grown with the domain and formed a partition of $[0,L(t)]$. Note that whilst the initial cell widths were constant this is no longer the case in the growing domain, i.e. $\delta y_{i}$ is a function of both the initial position $x_{i}$ and time. The mapping is defined through a growth function, $\mu(x_{i},t)$, which defines the growth rate of the interval at $x_{i}$ at a time $t$. Explicitly it can be shown that the mapping $g(x,t)$ from a position in the fixed domain, $x$, to a corresponding position on the growing domain, $y$, is given by,

[TABLE]

This is illustrated schematically in Figure 1. Note that, $g(0,t)=0$ and the initial condition, $y=g(x,0)=x$ for all $x\in[0,L_{0}]$, places a physical restriction on the mapping between $y$ and $x$. For future notational convenience we will denote the spatial derivative of $g(x,t)$ as $\nu^{*}(x,t)$, so that,

[TABLE]

and the time derivative as,

[TABLE]

As the mapping is invertible, so that $x=g^{-1}(y,t)$, these can be expressed on the growing domain, giving,

[TABLE]

and

[TABLE]

It should also be noted that if we consider the growth of a small interval in the initial domain $(x,x+\delta x)$, then the width of the interval at some later time, in the limit of small $\delta x$ can be written as,

[TABLE]

We now consider a CTRW on a growing domain, such that a particle will jump to a location, wait for some time, and then jump to a new location. We will assume that the waiting time and jump length densities are independent. The waiting time probability density for a particle that arrived at a location at time $t^{\prime}$ to jump at time $t$ will be denoted by $\psi(t-t^{\prime})$, where $t-t^{\prime}$ is the amount of time that the particle waited. The jump length density for a particle that is at a location $z^{\prime}$ to jump to location $z$ at time $t$ is denoted by $\lambda(z|z^{\prime},t)$. In the following we consider a CTRW on the growing domain $z=y$ and an auxiliary CTRW on the fixed domain $z=x$. In taking the diffusion limit we will restrict ourselves to fixed length jumps on the growing domain, $\Delta y$, where the particle may jump either left or right. The corresponding jumps in the auxiliary CTRW on the initial fixed domain will therefore have lengths that change in both time and space as the domain grows.

For a particle undergoing a CTRW on the growing domain, we let $\rho(y,t)\delta y(x,t)$ denote the probability of finding the particle in the region $(y,y+\delta y(x,t))$, in the time $(t,t+\delta t)$ for a small $\delta y(x,t)$. Thus $\rho(y,t)$ is the probability density of finding the particle, which we can express as follows,

[TABLE]

where $\Phi(t-t^{\prime})$ is the survival function associated with the waiting time density $\psi(t-t^{\prime})$. The inbound flux, $q(g(x,t),t)$, is defined such that the probability of the particle entering the region $(y,y+\delta y(x,t))$ in the time $(t,t+\delta t)$, given $y=g(x,t)$, is $q(g(x,t),t)\delta y(x,t)\delta t$. This equation states that for a particle to be in the region, it must have previously arrived in the region and not jumped away.

Equation (7) can be simplified by using Eq. (6),

[TABLE]

To transform the evolution equation to a master equation it is necessary to replace the explicit dependence on $q(g(x,t),t)$ with a dependence on $\rho(g(x,t,),t)$. The growth of the domain requires us to utilise non-standard techniques to achieve this. As the region is moving and growing this is most easily expressed by mapping the required functions back to the fixed $x$ domain. The formulation of the CTRW on the fixed domain will be referred to as an auxiliary CTRW.

To formulate the auxiliary CTRW on the fixed domain, we relate the associated densities to densities on the growing domain, such that,

[TABLE]

Here we use a star to denote a function associated with the auxiliary process on the fixed domain. Hence we can write the auxiliary form of Eq. (8) as,

[TABLE]

Note that this left hand side, $\rho^{*}(x,t)e^{\int_{0}^{t}\mu(x,s)ds}$, is a conserved probability density. Differentiating Eq. (10) with respect to time and simplifying, we arrive at an evolution equation for the probability density,

[TABLE]

In this equation the second term on the right hand side is the flux out of the neighbourhood around $x$ in the time interval around $t$, while the third term is the reduction in concentration of particles, around $x$ around $t$, due to the growth of the domain. Explicitly we define the flux out as,

[TABLE]

In this equation the incoming flux, $q^{*}(x,t)$, can itself be expressed in terms of the flux out resulting in the relation,

[TABLE]

Using Eq. (13), noting the semi-group property of the exponential function, we can rewrite Eq. (11) and using Laplace transform methods, we can express the evolution equation for the auxiliary CTRW as the auxiliary master equation,

[TABLE]

In this equation, the memory kernel, $K(t)$, is defined by,

[TABLE]

where $\mathcal{L}_{t}$ denotes a Laplace transform with respect to time.

The master equation, Eq. (14), has been derived for arbitrary waiting time and jump densities. To obtain a diffusion limit of the master equation we will require specific forms for these densities. We wish to consider the case of a fixed jump length on the growing domain, where the particle will jump either right or left with equal probability. In this case the jump length for the auxiliary master equation will change with both space and time. The jump probability density can therefore be written as,

[TABLE]

where $\delta(x)$ is the Dirac delta function, and $\epsilon^{+}$ and $\epsilon^{-}$ are time and space dependent. To relate the $\epsilon$’s to the fixed jump length, $\Delta y$, we note that from Eq. (1) we have,

[TABLE]

Using the relations from Eqs. (17) and (18), we perform a Taylor expansion of Eq. (14) with the jump distribution given by, Eq. (16) around $\Delta y=0$ to arrive at,

[TABLE]

We now take a Mittag-Leffler waiting time density, given by,

[TABLE]

with $0<\alpha<1$ and $\tau>0$ HA1995 , where $E_{\alpha,\beta}$ is a two parameter Mittag-Leffler function defined by,

[TABLE]

The Mittag-Leffler probability density is heavy tailed, which is asymptotically $\psi(t)\sim t^{-1-\alpha}$ for long times. The memory kernel of a Mittag-Leffler probability density can be calculated from the inverse Laplace transform of Eq. (15),

[TABLE]

The Riemann-Liouville fractional derivative of order $1-\alpha$ is defined as,

[TABLE]

As we are considering smooth real valued functions, the initial condition term in the Laplace transform of the Riemann-Liouville fractional derivative will be zero LD2007 , so that,

[TABLE]

Using Mittag-Leffler distributed waiting times the auxiliary master equation on the fixed domain, Eq. (19), becomes,

[TABLE]

The fractional diffusion limit is one in which the length and time scales are taken to zero, such that,

[TABLE]

exists. The fractional diffusion limit of Eq. (25) is,

[TABLE]

This is the auxiliary fractional diffusion equation defined on the fixed domain. Note that, apart from the advective type term, this is the same form as a fractional reaction subdiffusion equation ADH2013mmnp , with the additional feature of a space and time dependent diffusivity. In writing the equation in terms of the growing domain coordinates the diffusivity will be constant.

Boundary conditions may be implemented by considering different jump length densities near the boundary. Explicitly, a zero flux boundary will be implemented by taking,

[TABLE]

for $x\in[L(0)-\epsilon^{-},L(0)]$ and,

[TABLE]

for $x\in[0,\epsilon^{+}]$. This jump density guarantees that there is no flux across the boundary, and in the diffusive limit the master equation at the boundary point will be consistent with the master equation in the bulk.

Using the jump length density for the left boundary, Eq. (29), and taking a Taylor expansion around $\Delta y=0$, the master equation, Eq. (14), becomes,

[TABLE]

for $x\in[0,\epsilon^{+}]$. The difference between this equation and the bulk result is the occurrence of a first order spatial derivative. With the Mittag-Leffler waiting time density in order for the diffusion limit, Eq. (26) to exist, we require the first order spatial derivative term to be,

[TABLE]

Only holding at the boundary point as $\Delta y\to 0$. This zero flux boundary condition is equivalent to the zero flux boundary derived for fractional reaction subdiffusion equations ADHJL2016 . The derivation for the right hand side of the boundary results in an equivalent condition.

The fractional diffusion equation can be found by mapping the auxiliary equation, Eq. (27), to the growing domain. Using the mapping $y=g(x,t)$, with Eqs. (4) and (5), we perform a change of variables and find,

[TABLE]

Here we have defined a new comoving fractional derivative, $\,{}_{0}^{g}\mathcal{C}_{t}^{1-\alpha}$, which operates along the curve, $y=g(x,t)$, for a fixed $x$. Formally this is defined as,

[TABLE]

Informally, the history of the function is not integrated over a fixed value of $y$ but rather along the trajectory of the point in the domain as it grows. We note that,

[TABLE]

The boundary condition, Eq. (31), on the growing domain is,

[TABLE]

We note that when the rate of the domain growth has no spatial dependence, i.e. $\nu(y,t)=f(t)$, the boundary conditions for the Mittag-Leffler waiting time density is simplified to,

[TABLE]

on the growing domain.

In this work we have derived evolution equations that describe subdiffusive transport on a growing domain. Equation (27) describes the transport on a rescaled fixed domain whilst Eq. (32) describes the same process on the growing domain. The evolution equation on the growing domain required the definition of a new fractional order differential operator that follows the domain growth, Eq. (33). Our work provides a sound basis for further explorations of models of subdiffusive transport on a growing domain, including reactions, pattern formation, and morphogenesis.

This work was supported by the Australian Commonwealth Government (ARC DP140101193). We gratefully acknowledge useful discussions with S. B. Yuste.