Time Fractional Fisher–KPP and Fitzhugh–Nagumo Equations

Abstract

A standard reaction–diffusion equation consists of two additive terms, a diffusion term and a reaction rate term. The latter term is obtained directly from a reaction rate equation which is itself derived from known reaction kinetics, together with modelling assumptions such as the law of mass action for well-mixed systems. In formulating a reaction–subdiffusion equation, it is not sufficient to know the reaction rate equation. It is also necessary to know details of the reaction kinetics, even in well-mixed systems where reactions are not diffusion limited. This is because, at a fundamental level, birth and death processes need to be dealt with differently in subdiffusive environments. While there has been some discussion of this in the published literature, few examples have been provided, and there are still very many papers being published with Caputo fractional time derivatives simply replacing first order time derivatives in reaction–diffusion equations. In this paper, we formulate clear examples of reaction–subdiffusion systems, based on; equal birth and death rate dynamics, Fisher–Kolmogorov, Petrovsky and Piskunov (Fisher–KPP) equation dynamics, and Fitzhugh–Nagumo equation dynamics. These examples illustrate how to incorporate considerations of reaction kinetics into fractional reaction–diffusion equations. We also show how the dynamics of a system with birth rates and death rates cancelling, in an otherwise subdiffusive environment, are governed by a mass-conserving tempered time fractional diffusion equation that is subdiffusive for short times but standard diffusion for long times.

1. Introduction

Reaction–diffusion partial differential equations are among the most widely used equations in applied mathematics modelling. These equations govern the time evolution of concentrations, or population densities, of species, at different spatial locations, that are diffusing and reacting. Applications include the spatio-temporal spread of epidemics, the spatial spread of invasive species and the development of animal coat patterns [1,2,3]. In these modelling equations, diffusion is represented by a spatial Laplacian operating on the population densities, and reactions are included as additive terms representing changes per unit time in population densities through reaction rates. In well-mixed systems the reaction rate equations can often be derived from the law of mass-action [4]. A famous example of a reaction–diffusion equation is the Fisher–KPP equation named after Fisher [5] and Kolmogorov, Petrovsky and Piskunov [6]. The standard reaction–diffusion representation of this equation is

$$\frac{\partial u(x,t)}{\partial t}=D\frac{{\partial }^{2}u(x,t)}{\partial x^2}+ru(x,t)(1-u(x,t)),D>0,r>0.$$

(1)

Here, $u(x,t)$ represents the population density of a species, $D\frac{{\partial }^{2}u(x,t)}{\partial x^2}$ represents the diffusion of the species and $ru(x,t)(1-u(x,t\left)\right)$ represents the reactions of the species. In the absence of diffusion, the time rate of change in the population density is the same at all points in space and is given by

$$\frac{\partial u(x,t)}{\partial t}=ru(x,t)(1-u(x,t)).$$

(2)

In this example and in the following, for simplicity, we have considered systems in one spatial dimension. Extensions to higher spatial dimensions are possible.

Over the past two decades, there has been a growing awareness of fractional diffusion, where diffusion cannot be modelled using a standard Laplacian and the mean square displacement of diffusing species does not grow linearly in time, as anticipated by Einstein’s famous modelling of Brownian motion [7]. In particular, following widespread observations in biological systems, there has been a great deal of attention focussed on fractional subdiffusion, characterized by the mean square displacement of a population spreading as a sublinear power law in time. It is now generally accepted that if subdiffusion arises from particles being trapped for arbitrarily long periods of time, the appropriate equation to model subdiffusion is the time fractional diffusion equation [8]

$$\frac{\partial u(x,t)}{\partial t}={}_0{D}_t^{1-\gamma }\frac{{\partial }^{2}u(x,t)}{\partial x^2},0<\gamma <1,$$

(3)

which can be derived [9,10] from a continuous time random walk (CTRW) [11] with a power law waiting time density. In this equation,

$${}_0{D}_t^{1-\gamma }y(x,t)=\frac{1}{\Gamma \left(\gamma \right)}\frac{\partial }{\partial t}{\int }_0^t\frac{y(x,t^{\prime })}{{(t-t^{\prime })}^{1-\gamma }}d{t}^{\prime }$$

(4)

is the Riemann–Liouville fractional derivative of order $1-\gamma$, see, for example, reference [12]. It might be anticipated that the appropriate evolution equation to model subdiffusion, with reactions governed by the reaction rate equation,

$$\frac{\partial u(x,t)}{\partial t}=f\left(u(x,t)\right),$$

(5)

would be

$$\frac{\partial u(x,t)}{\partial t}={}_0{D}_t^{1-\gamma }\frac{{\partial }^{2}u(x,t)}{\partial x^2}+f\left(u(x,t)\right).$$

(6)

Indeed, such an equation had been derived from an underyling CTRW model, under certain assumptions, [13], however it is not valid in general. For example, the simple model equation

$$\frac{\partial u(x,t)}{\partial t}={}_0{D}_t^{1-\gamma }\frac{{\partial }^{2}u(x,t)}{\partial x^2}-u(x,t),$$

(7)

can have unphysical negative solutions [14].

The time fractional subdiffusion equation is also often written as [15]

$$\frac{{\partial }^{\gamma }u(x,t)}{\partial t^{\gamma }}=\frac{{\partial }^{2}u(x,t)}{\partial x^2},0<\gamma <1,$$

(8)

where

$$\frac{{\partial }^{\gamma }}{\partial t^{\gamma }}y(x,t)=\frac{1}{\Gamma (1-\gamma )}{\int }_0^t\frac{\frac{\partial }{\partial t^{\prime }}y(x,t^{\prime })}{{(t-t^{\prime })}^{\gamma }}d{t}^{\prime }$$

(9)

denotes a Caputo fractional derivative, see, for example, reference [12]. There has been quite a bit written in the published literature on the greater physical practicality of the Caputo derivative over the Riemann–Liouville derivative, but this is largely unfounded [12]. Note, however, that if one takes Equation (8) as the starting evolution equation for subdiffusion then this is suggestive of the following reaction–subdiffusion equation,

$$\frac{{\partial }^{\gamma }u(x,t)}{\partial t^{\gamma }}=\frac{{\partial }^{2}u(x,t)}{\partial x^2}+f\left(u(x,t)\right).$$

(10)

Equations along the lines of Equation (10) are particularly widespread in the literature with the motivation that fractional derivatives incorporate a history dependence, and solutions of Equation (10) remain positive. Equation (10) can be derived from a CTRW where particles are being removed or added instantaneously at the start of the waiting times between jumps, but only under the contrived constraint that $\frac{{\partial }^{1-\gamma }f\left(u(x,t)\right)}{\partial t^{1-\gamma }}$ represents the cumulative total of additions and removals to the arrival density of particles at position x and time t [14].

The derivation of reaction–subdiffusion equations from physically consistent CTRWs has been carried out in a series of papers [14,16,17,18,19,20,21,22,23,24,25,26]. The main lessons from this body of work are: (i) The governing equations are different depending on whether or not new born particles inherit the waiting times of their parents. (ii) Birth terms and death terms must be treated differently. (iii) In the case where particles are removed, but not instantaneously at the start of the waiting time between jumps, the reaction and subdiffusion terms are not additive. The following equation [21,24],

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& =& D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-{\int }_0^{t}a(u(x,t^{\prime }),x,t^{\prime })d{t}^{\prime }}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{{\int }_0^{t}a(u(x,t^{\prime }),x,t^{\prime })d{t}^{\prime }}u(x,t)\right)\right]\hfill \\ & & +c\left(u\right(x,t),x,t)-a\left(u\right(x,t),x,t)u(x,t),\hfill \end{array}$$

(11)

which was derived from a continuous time random walk model, provides the evolution equation for particles undergoing subdiffusion with particles annihilated at a per capita rate, $a\left(u\right(x,t),x,t)$ and created at a rate $c\left(u\right(x,t),x,t)$. In the derivation of this equation it was assumed that newborn particles do not inherit the waiting times of their parents.

In the remainder of this paper we explore examples related to Equation (11). These examples have been selected to emphasize the importance of considering the details of the reaction kinetics when dealing with reaction–subdiffusion problems. Whilst there have been many papers published on various methods of solution for variants of Equation (10) (see, for example, [27,28,29,30,31]), there have been very few papers published considering algebraic or numerical solution methods for variants of Equation (11). We hope that the examples below will stimulate further activity in this area, where the physical motivation for the modelling equation is stronger.

2. Examples

2.1. Birth and Death Balance

As a first example, we consider a population density of $u(x,t)$ particles per unit volume that are diffusing with a per capita death rate $\alpha$ and a birth rate $\alpha u(x,t)$. The reaction rate equation reflecting this balance between births and deaths, in a well-mixed population, at a location x is

$$\frac{\partial u(x,t)}{dt}=0,$$

(12)

and thus the standard reaction–diffusion equation describing this system is

$$\frac{\partial u(x,t)}{dt}=D\frac{{\partial }^{2}u(x,t)}{\partial x^2}.$$

(13)

The simple generalization of this equation for subdiffusive transport is

$$\begin{array}{ccc}\frac{\partial u(x,t)}{dt}& =& D_{\gamma }{}_0{\mathcal{D}}_t^{1-\gamma }\frac{{\partial }^{2}u(x,t)}{\partial x^2}\hfill \\ & =& D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[{}_0{\mathcal{D}}_t^{1-\gamma }u(x,t)\right].\hfill \end{array}$$

(14)

Indeed, if there were no births or deaths then the reaction rate equation would still be given by Equation (12); and Equation (14) is the appropriate equation to describe subdiffusion without births or deaths. However, the reaction–subdiffusion equation, following Equation (11), and using the reaction rate kinetics $a\left(u\right(x,t),x,t)=\alpha$ and $c\left(u\right(x,t),x,t)=\alpha u(x,t)$, which are also consistent with the rate equation, Equation (12), is remarkably different;

$$\frac{\partial u(x,t)}{\partial t}=D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-\alpha t}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{\alpha t}u(x,t)\right)\right],\alpha >0.$$

(15)

The fundamental difference between Equations (14) and (15) is that in the former equation the Laplacian operates on a time fractional derivative and in the latter the Laplacian operates on a tempered time fractional derivative [32,33]. In the more general time fractional reaction–diffusion equation, Equation (11), the term in brackets following the Laplacian defines a generalized tempered time fractional derivative. The physical interpretation of the tempering is that if particles are being annihilated at a given rate while they wait then they cannot wait an arbitrarily long time at a given location. Note that both Equations (14) and (15) are mass conserving and thus Equation (15) then defines a mass conserving, tempered, time fractional diffusion equation.

The mean square displacement of the diffusing particles, $\langle x^2\left(t\right)\rangle$, provides a clear measurable difference between particles following Equation (14) or Equation (15). In the former case, identified as $\langle x_I^2\left(t\right)\rangle$, we have [8],

$$\langle x_I^2\left(t\right)\rangle =\frac{2D_{\gamma }}{\Gamma (1+\gamma )}{t}^{\gamma },$$

(16)

and in the latter case, identified as $\langle x_{II}^2\left(t\right)\rangle$, we have (Appendix A)

$$\langle x_{II}^2\left(t\right)\rangle =2D_{\gamma }{e}^{-\alpha t}{t}^{\gamma }{E}_{1,\gamma }^{\left(1\right)}\left(\alpha t\right),$$

(17)

where

$$E_{1,\gamma }^{\left(1\right)}\left(z\right)=\frac{d}{dz}\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\gamma +k)}$$

(18)

is the derivative of a generalized Mittag–Leffler function [34]. Note that at short times,

$$\langle x_{II}^2\left(t\right)\rangle \sim \frac{2D_{\gamma }}{\Gamma (1+\gamma )}{t}^{\gamma },$$

(19)

but at large times, using the asymptotic expansion of the generalized Mittag–Leffler function (Equation (6) in [35]),

$$\langle x_{II}^2\left(t\right)\rangle \sim 2D_{\gamma }{\alpha }^{1-\gamma }t.$$

(20)

Thus, mass conserving tempered time fractional diffusion is not anomalous at long times.

We can also write down explicit expressions for solutions to Equations (14) and (15), labelled as $u_I(x,t)$ and $u_{II}(x,t)$, respectively. For simplicity we consider the infinite domain Greens function solutions with initial condition $u(x,0)=\delta \left(x\right)$.

The Greens function solution of the fractional diffusion equation Equation (14) can be written as [8]

$$u_I(x,t)=\frac{1}{\sqrt{4\pi D_{\gamma }{t}^{\gamma }}}{H}_{1,2}^{2,0}\left[\frac{x^2}{4D_{\gamma }{t}^{\gamma }}\left|\begin{array}{cc}(1-\frac{\gamma }{2},\gamma )\hfill & \\ (0,1)\hfill & (\frac{1}{2},1)\hfill \end{array}\right.\right],$$

(21)

where H denotes a Fox H-function [36], see Equation (A11).

To find the Greens function solution $u_{II}(x,t)$ we first note that Equation (15) can be re-written as

$$\frac{\partial v(x,t)}{\partial t}=D_{\gamma }\frac{{\partial }^2}{\partial x^2}{}_0{\mathcal{D}}_t^{1-\gamma }v(x,t)+\alpha v(x,t),$$

(22)

where

$$v(x,t)=e^{\alpha t}{u}_{II}(x,t).$$

(23)

The Greens function solution of Equation (22) can be obtained as a special case of the more general results in Appendix B of [14], yielding

$$v(x,t)=\frac{1}{\sqrt{4\pi D_{\gamma }{t}^{\gamma }}}\sum _{j=0}^{\infty }\frac{{\left(\alpha t\right)}^j}{j!}{H}_{1,2}^{2,0}\left[\frac{x^2}{4D_{\gamma }{t}^{\gamma }}\left|\begin{array}{cc}(1-\frac{\gamma }{2}+j,\gamma )\hfill & \\ (0,1)\hfill & (\frac{1}{2}+j,1)\hfill \end{array}\right.\right],$$

(24)

and then using Equation (23) we have

$$u_{II}(x,t)=e^{-\alpha t}\frac{1}{\sqrt{4\pi D_{\gamma }{t}^{\gamma }}}\sum _{j=0}^{\infty }\frac{{\left(\alpha t\right)}^j}{j!}{H}_{1,2}^{2,0}\left[\frac{x^2}{4D_{\gamma }{t}^{\gamma }}\left|\begin{array}{cc}(1-\frac{\gamma }{2}+j,\gamma )\hfill & \\ (0,1)\hfill & (\frac{1}{2}+j,1)\hfill \end{array}\right.\right].$$

(25)

In Appendix B, we show that the Fox functions in Equations (21) and (25) can be simplified for $\gamma =\frac{1}{2}$ in terms of Miejer G-Functions [37], see Equation (A12), which have the advantage that they can readily be evaluated using computer algebra packages such as MATHEMATICA and MAPLE. Using the result of Equation (A19) from the Appendix B, we can write (see also [8] in the case of $u_I(x,t)$)

$$u_I(x,t)=\frac{1}{\sqrt{8{\pi }^{3}D{t}^{\frac{1}{2}}}}{G}_{0,3}^{3,0}\left[{\left(\frac{x^2}{16D{t}^{\frac{1}{2}}}\right)}^2\left|\begin{array}{c}-\hfill \\ 0,\frac{1}{4},\frac{1}{2}\hfill \end{array}\right.\right],$$

(26)

and

$$u_{II}(x,t)=e^{-\alpha t}\frac{1}{\sqrt{8{\pi }^{3}D{t}^{\frac{1}{2}}}}\sum _{j=0}^{\infty }\frac{{\left(2\alpha t\right)}^j}{j!}{G}_{1,4}^{4,0}\left[{\left(\frac{x^2}{16D{t}^{\frac{1}{2}}}\right)}^2\left|\begin{array}{c}\frac{3}{4}+j\hfill \\ 0,\frac{1}{2},\frac{1}{4}+\frac{j}{2},\frac{3}{4}+\frac{j}{2}\hfill \end{array}\right.\right].$$

(27)

Note that the expression for $u_{II}(x,t)$ simplifies to the expression for $u_I(x,t)$ if $\alpha =0$. If $\left|x\right|\gg 4\sqrt{D{t}^{\frac{1}{2}}}$ then we can use asymptotic expansions for $G_{0,3}^{3,0}\left(z\right)$ and $G_{1,4}^{4,0}\left(z\right)$ with $z\gg 1$ (see Appendix B) to write

$$u_I(x,t)\sim \frac{1}{\sqrt{8{\pi }^{3}D{t}^{\frac{1}{2}}}}exp(-3{\left(\frac{x^2}{16D{t}^{\frac{1}{2}}}\right)}^{\frac{2}{3}}){\left(\frac{x^2}{16D{t}^{\frac{1}{2}}}\right)}^{\frac{1}{2}}\frac{M_0}{{\left(\frac{x^2}{16D{t}^{\frac{1}{2}}}\right)}^{\frac{2}{3}}-1},$$

(28)

and

$$u_{II}(x,t)\sim M{e}^{\alpha t}{u}_I(x,t),$$

(29)

where $M_0$ and M are constant terms. The solutions $u_I\left(t\right)$, Equation (26), and $u_{II}\left(t\right)$, Equation (27) are plotted in Figure 1, with $\alpha =1$ and $D=1$, at times $t=0.1$, $t=1.0$ and $t=10.0$. The solutions are very similar at early times but the corner at the origin, which is characteristic of subdiffusion, is less sharp at longer times in the solution of Equation (27).

The lesson from this simple example is that reaction dynamics equations do not contain sufficient information on their own to provide model equations for reaction–subdiffusion systems even in well-mixed systems. In the case of standard diffusion, the evolution of the population density is only affected by the overall reaction rates, in a well-mixed system, but not the details of the reaction kinetics. In a standard reaction–diffusion system, the dynamics with no births and no deaths is the same as if there were births and deaths but the rates cancelled out. The reaction–diffusion equation with these reaction kinetics has no memory of the birth and death processes. This is very different in the case of subdiffusion where the details of the reaction kinetics are important to the overall dynamics of the system. The subdiffusive system retains a memory that there were particles that were created and annihilated. Moreover, the particle deaths temper the fractional diffusion. The example in the next section further highlights the significance of the reaction kinetics in reaction–subdiffusion systems.

2.2. Fractional Fisher–KPP Equation

The reaction rate equation for the Fisher–KPP Equation (1) is given in Equation (2). There are many different reaction kinetics that could be considered that are consistent with Equation (2). For example, the term $(1-u(x,t\left)\right)$ in its entirety could represent a per capita birth rate if is is strictly positive, or a per capita death rate if it is strictly negative. This term could also be regarded as being composed of two terms, a constant per capita birth term and a linear per capita death term. These three possibilities are highlighted for illustrative purposes below to show how different subdiffusion–reaction equations apply depending on the reaction kinetics.

(i)

Constant per capita birth rate, $c\left(u\right(x,t),x,t)=ru(x,t)$, linear per capita death rate, $a\left(u\right(x,t),x,t)=ru(x,t)$,

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& =& D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-{\int }_0^{t}ru(x,t^{\prime })d{t}^{\prime }}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{{\int }_0^{t}ru(x,t^{\prime })d{t}^{\prime }}u(x,t)\right)\right]\hfill \\ & & +ru(x,t)(1-u(x,t\left)\right),u(x,t)\ge 0\hfill \end{array}$$

(30)

(ii)

No births, $c\left(u\right(x,t),x,t)=0$, linear per capita death rate, $a\left(u\right(x,t),x,t)=r(1-u(x,t\left)\right)$,

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& =& D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-{\int }_0^{t}r(1-u(x,t^{\prime }))d{t}^{\prime }}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{{\int }_0^{t}r(1-u(x,t^{\prime }))d{t}^{\prime }}u(x,t)\right)\right]\hfill \\ & & +r(1-u(x,t\left)\right)u(x,t),u(x,t)\ge 1.\hfill \end{array}$$

(31)

(iii)

Linear per capita birth rate, $c\left(u\right(x,t),x,t)=ru(x,t)(1-u(x,t\left)\right)$, no deaths, $a\left(u\right(x,t),x,t)=0$,

$$\frac{\partial u(x,t)}{\partial t}=D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[{}_0{\mathcal{D}}_t^{1-\gamma }u(x,t)\right]+ru(x,t)(1-u(x,t)),0\le u(x,t)<\le 1.$$

(32)

Note that none of the factional Fisher–KPP reaction–diffusion equations can be expressed in the form

$$\frac{{\partial }^{\gamma }u(x,t)}{\partial t^{\gamma }}=D_{\gamma }\frac{{\partial }^{2}u(x,t)}{\partial x^2}+ru(x,t)(1-u(x,t)),$$

(33)

which results from simply replacing the integer order time derivative with a fractional order Caputo derivative. As noted above, an equation of this form could only be obtained from a CTRW if $\frac{{\partial }^{1-\gamma }}{\partial t^{1-\gamma }}\left(ru(x,t)(1-u(x,t)\right)$ is contrived as the cumulative instantaneous creation and annihilation of particles at the start of the waiting time between particle jumps at position x and time t [14].

The Greens function solutions for the nonlinear fractional reaction–diffusion equations, Equations (30)–(32), cannot be obtained simply using Fourier–Laplace transform methods. However, it is possible to find numerical solutions using the discrete time random walk methods described in [38].

The Fisher–KPP reaction rate equation, Equation (2) can be motivated by different chemical reactions consistent with the law of mass action [4]. One possibility is that of a single species A which undergoes coalescence reactions $A+A\stackrel{r}{\to }A$, and degradation reactions $A\stackrel{r}{\to }A+A$; also referrred to as reversible coagulation dynamics [39]. In this scenario the creation term, $ru(x,t)$, arises from degradation and the annihilation term, $-r\left(u^2(x,t)\right)$, arises from coalescence. Another possibility is a branching–coalescence scheme [17], $B+X\rightleftharpoons X+X$, with the concentration of B maintained at a constant level. Equation (30) is a fractional Fisher–KPP reaction–diffusion equation consistent with each of the reaction schemes described here and it was obtained earlier for the branching–coalescence reaction scheme in [17].

2.3. Fractional Fitzhugh–Nagumo Equation

A widely studied reaction–diffusion system used to model wave propagation and pattern formation in excitable media is the Fitzhugh–Nagumo system of equations [40,41]

$$\begin{array}{ccc}\hfill \frac{\partial v(x,t)}{\partial t}& =& D_v\frac{{\partial }^{2}v(x,t)}{\partial x^2}+v(x,t)(v(x,t)-a)(1-v(x,t))-w(x,t),D_v\ge 0,a\ge 0\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \frac{\partial w(x,t)}{\partial t}& =& D_w\frac{{\partial }^{2}w(x,t)}{\partial x^2}+ϵ\left(v(x,t)-bw(x,t)\right)D_w\ge 0,ϵ\ge 0,b\ge 0,\hfill \end{array}$$

(35)

named after Fizthugh [42] and Nagumo [43]. In recent years, the single component fractional equation

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=D_u\frac{{\partial }^{2}u(x,t)}{\partial x^2}+u(x,t)(u(x,t)-a)(1-u(x,t))$$

(36)

has been studied as a test equation for various methods of solution of time fractional reaction–diffusion equations (see, for example, [27,30] and references there-in).

A time fractional Fitzhugh–Nagumo system of equations consistent with Equation (11), derived from a CTRW formalism, can be obtained by identifying per capita annihilation rates, $a_v$ and $a_w$, and creation rates, $c_v$ and $c_w$, as follows:

$$\begin{array}{ccc}\hfill a_v(v(x,t),w(x,t))& =& a+v^2(x,t)+\frac{w(x,t)}{v(x,t)},\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill c_v(v(x,t),w(x,t))& =& (1+a)v^2(x,t),\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill a_w(v(x,t),w(x,t))& =& ϵb,\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill c_w(v(x,t),w(x,t))& =& ϵv(x,t).\hfill \end{array}$$

(40)

The corresponding time fractional Fitzhugh–Nagumo system is given by

$$\begin{array}{ccc}\hfill \frac{\partial v(x,t)}{\partial t}& =& D_{v,\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-{\int }_0^t(v^2(x,t^{\prime })+a+w(x,t^{\prime }))d{t}^{\prime }}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{{\int }_0^t(v^2(x,t^{\prime })+a+w(x,t^{\prime }))d{t}^{\prime }}v(x,t)\right)\right]\hfill \end{array}$$

$$\begin{array}{ccc}& & +v(x,t)\left(v\right(x,t)-a)(1-v(x,t)-w(x,t),\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill \frac{\partial w(x,t)}{\partial t}& =& D_{w,\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-ϵbt}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{ϵbt}w(x,t)\right)\right]+ϵv(x,t)-ϵbw(x,t).\hfill \end{array}$$

(42)

If $w(x,t)=0$ this identifies a single component time fractional equation

$$\begin{array}{ccc}\hfill \frac{\partial u(x,t)}{\partial t}& =& D_{\gamma }\frac{{\partial }^2}{\partial x^2}\left[e^{-{\int }_0^t(u^2(x,t^{\prime })+a)d{t}^{\prime }}{}_0{\mathcal{D}}_t^{1-\gamma }\left(e^{{\int }_0^t(u^2(x,t^{\prime })+a)d{t}^{\prime }}u(x,t)\right)\right]\hfill \\ & & +u(x,t)\left(u\right(x,t)-a)(1-u(x,t),\hfill \end{array}$$

(43)

which could be called a time fractional Fitzhugh–Nagumo equation, although the nomenclature could be misleading because a single component equation, without external sources or sinks, could not display Fitzhugh–Nagumo dynamics. Equation (43) is, however, well posed as a nonlinear time fractional reaction–diffusion equation that can be derived from a physically consistent CTRW, and thus it should be preferred for testing numerical methods of solution over the single component model equation, Equation (36), obtained by replacing an integer order time derivative with a Caputo fractional order derivative.

3. Discussion

Over the past two decades there have been large numbers of papers published on numerical methods for nonlinear fractional reaction–diffusion equations. The original motivation for including time fractional derivatives in reaction–diffusion equations was based on a CTRW description of diffusion with traps and reactions [13]. This description was refined and improved in a series of papers [14,16,17,18,19,20,21,22,23,24], leading to the formulation of time fractional reaction–diffusion equations along the lines of Equation (11). However, many investigations of time fractional reaction–diffusion equations have been carried out on systems obtained by simply replacing integer order time derivatives with Caputo fractional order derivatives. These studies may be interesting from a mathematical analysis point of view but they may not be directly relevant to mathematical modelling applications.

In this paper we have illustrated, through examples, how different time fractional reaction–diffusion equations can be formulated, consistent with an underlying CTRW formalism, taking into account the reaction kinetics. There are three points worth noting in this context: (i) The fractional reaction–diffusion systems considered in this approach are relevant to well-mixed reactions that are not diffusion limited. The reaction dynamics can often be formulated using the law of mass action in these systems. (ii) Different time fractional reaction–diffusion systems can be formulated that are consistent with the same equation for the reaction dynamics. It is important to know the reaction kinetics. (iii) Reaction–subdiffusion equations typically involve a spatial Laplacian operating on a generalized tempered time fractional derivative. The solution of these types of equations would typically require very different numerical approaches than those proposed for reaction–diffusion systems with a fractional order Caputo time derivative replacing the integer order time derivative.

It is hoped that the physically motivated time fractional reaction–diffusion equations, such as Equations (30) and (43), will become more widely used, replacing the simpler ad-hoc equations, such as Equations (33) and (36), as a test for different methods of solution of nonlinear fractional reaction–diffusion systems. Beyond this, there is a real need for physical experiments to be devised and carried out to validate and calibrate time fractional reaction–diffusion models.