Front Propagation of Exponentially Truncated Fractional-Order Epidemics

Abstract

The existence of landscape constraints in the home range of living organisms that adopt Lévy-flight movement patterns, prevents them from making arbitrarily large displacements. Their random movements indeed occur in a finite space with an upper bound. In order to make realistic models, by introducing exponentially truncated Lévy flights, such an upper bound can thus be taken into account in the reaction-diffusion models. In this work, we have investigated the influence of the $\lambda$-truncated fractional-order diffusion operator on the spatial propagation of the epidemics caused by infectious diseases, where $\lambda$ is the truncation parameter. Analytical and numerical simulations show that depending on the value of $\lambda$, different asymptotic behaviours of the travelling-wave solutions can be identified. For small values of $\lambda$ ($\lambda \gtrsim 0$), the tails of the infective waves can decay algebraically leading to an exponential growth of the epidemic speed. In that case, the truncation has no impact on the superdiffusive epidemics. By increasing the value of $\lambda$, the algebraic decaying tails can be tamed leading to either an upper bound on the epidemic speed representing the maximum speed value or the generation of the infective waves of a constant shape propagating at a minimum constant speed as observed in the classical models (second-order diffusion epidemic models). Our findings suggest that the truncated fractional-order diffusion equations have the potential to model the epidemics of animals performing Lévy flights, as the animal diseases can spread more smoothly than the exponential acceleration of the human disease epidemics.

1. Introduction

Infectious pathogens, such as viruses are the microorganisms that can cause infectious diseases following their presence and growth in humans and animals (hosts) [1]. One of the characteristics of infectious diseases is the ability of pathogens to be transmitted between individuals by either close contact between infected and non-infected (susceptible) individuals or a secondary host that carries the pathogens. Infectious diseases are thus also called communicable or transmissible diseases. Influenza, novel coronavirus (COVID-19), malaria, dengue fever, and West Nile virus are examples of infectious diseases. The propagation of infectious pathogens into non-infected regions leads to the spatial expansion of transmissible diseases and as a result, a large number of susceptible individuals can be infected so that an epidemic occurs. In some cases, a human epidemic can propagate over the entire globe, which is known as a pandemic. The 1918 Spanish influenza and COVID-19 are examples of classic and novel pandemics, respectively [2,3].

In epidemiology, mathematical models have been widely used to study the temporal and spatial dynamics of infectious diseases. Assuming the spatial homogeneity of the environment leads to the models that encompass the temporal dynamics of infectious diseases [4,5,6,7]. These models vary based on the characteristics of the particular disease, such as the transmission mechanism by either close contact or secondary hosts (see [8] for an extensive review of various mathematical models of animal and human infectious diseases). One of the purposes of such models is to derive the basic reproductive number that indicates the potential for the propagation of infectious diseases within a population. When its value is greater than unity, it shows that the infection is capable of invading the susceptible population, leading to the occurrence of an epidemic. Since this number is dependent on the model parameters, some strategies can be taken into account to control the propagation of the epidemics [9,10]. However, considering the random mobility of hosts and also the spatial heterogeneity of the environment lead to the reaction-diffusion models that describe the spatial dynamics of infectious diseases [11,12,13]. As an original model of this type, Källén et al. [14] studied the spatial spread of rabies among foxes. The main advantage of the reaction-diffusion epidemic models is to understand how fast an infectious disease can propagate. In that case, in order to control the spatial spread of the diseases, the parameters that govern the epidemic speed can be distinguished [15].

The type of host random mobility is a key factor that affects the epidemic speed and the patterns of the epidemic spread. Brownian motion and Lévy flight are two different patterns of random movements. Brownian motion is a normal dispersion process in which particles (walkers) make short movements such that the probability distribution function (pdf) of the displacement is obtained by the normal (Gaussian) distribution [16]. Considering Brownian motion pattern of susceptible and infective individuals leads to a second-order reaction-diffusion epidemic model. In that case, the epidemic waves propagate at a constant speed and the tails of the epidemic’s fronts decay exponentially. Such classical models have been applied to estimate the propagation speed of the plague in the 14th century and rabies in Europe [15,17]. However, Lévy-flight mobility pattern represents an anomalous diffusion process in which the walkers make occasionally large jumps between many short movements such that the pdf of the displacement is obtained by a power-law distribution with an asymptotic behaviour as $l^{-(\alpha +1)}$, where l is the length of the displacement and $0<\alpha <2$ [18]. Different studies have shown that the random movement of many living organisms can be described by Lévy flights, including large marine predators [19,20], fruit flies [21], bumblebees [22], honeybees [23,24], albatrosses [25] and even humans [26]. Considering a population of the walkers performing a Lévy-flight mobility pattern leads to a fractional-order reaction-diffusion equation whose solution represents the walker densities [18,27]. From an applied point of view, fractional-order models have been used in a broad range of problems in ecology [28,29], biology [30], plasma turbulence [31,32], finance [33] and also recently numerous studies have been conducted on the subjects of bifurcation analysis, stability, and optimal control of fractional-order systems, such as predator-prey models [34,35], diffusive mussel-algae models [36], and neural networks [37,38,39]. As an application of space fractional-order diffusion equations in epidemiology, Hanert et al. [40] have shown that such equations can be applied to model the modern epidemics, such as avian influenza, and even SARS and also that the fractional-order diffusion operators make the epidemic waves travel at an exponential speed and the tails of the wave solutions decay algebraically as a power-law.

The existence of landscape and physiological constraints in the home range of living organisms following a pure Lévy-flight mobility pattern, can prevent arbitrarily large movements. To make realistic displacements, the tails of the power-law distribution must be truncated. In that case, the pure Lévy distribution is replaced by a truncated power-law distribution whose tail decays more quickly than it does in the pure Lévy case. Mantegna and Stanley [41] and Koponen [42] originally introduced truncated Lévy processes and after that Rosiński [43] proposed an exponentially truncated (tempered) Lévy process that yields a smoother decay than the abrupt cutoff. The pure power-law distributions result in diverging measurable quantities, such as variance and moments that cannot be observed in physically realizable systems. However, the truncated Lévy flight leads to the finite moments. The exponentially tempered Lévy distribution exhibits an asymptotic behaviour as ∼$e^{-\lambda l}{l}^{-(\alpha +1)}$, where l is the length of the displacement, $0<\alpha <2$, and $\lambda >0$ is the truncation parameter. When $\lambda =0$, it leads to a pure Lévy flight distribution. Considering a large number of random walkers performing a truncated Lévy flight leads to a truncated fractional-order diffusion equation whose solution represents the walker densities [44,45]. As an application of such equations, Vallaeys et al. [29] estimated isolation distances between genetically modified (GM) and non-GM crops by considering the fact that the honeybees mobility pattern can be well described by truncated Lévy flights. del Castillo Negrete [46] examined the influence of the truncation on the travelling-wave solutions of a truncated fractional-order Fisher–Kolmogorov equation. It would be of interest to investigate the results of such an equation in the epidemic models, as modern epidemics do not follow Brownian dynamics and can lead to accelerated waves of infectious diseases [47,48]. Hence, in this study, we apply the random mobility of infective individuals following truncated Lévy flights to an epidemic model and explore the epidemic speed based on the different values of the truncation parameter.

The remainder of the paper is organized as follows. In Section 2, we introduce the space fractional-order operators. In Section 3, we first review the classical and untruncated fractional-order epidemic models, then formulate a truncated fractional-order epidemic model and finally, discuss the epidemic speed of the new model. In Section 4, we explore our analytical results by numerically solving the model and conclusion is given in Section 5.

2. Preliminaries for Fractional-Order Operators

In this section, we shall give some definitions and notations that are required for this work. In the field of fractional models, the classical operators in time and/or space are replaced by the fractional derivative operators that are defined as follows.

For a given function f, the untruncated space fractional-order derivative operator of order $\alpha$ with a shorthand notation ${\mathcal{D}}_x^{\alpha }$ or $d^{\alpha }/d{\left|x\right|}^{\alpha }$ is defined as [44]:

$${\mathcal{D}}_x^{\alpha }f\left(x\right)=l{}_{-\infty }{D}_x^{\alpha }f\left(x\right)+r{}_x{D}_{\infty }^{\alpha }f\left(x\right),$$

(1)

where ${}_{-\infty }{D}_x^{\alpha }$ and ${}_x{D}_{\infty }^{\alpha }$ are the positive (left) and negative (right) space-fractional Riemann–Liouville derivatives, respectively, defined as [49,50]:

$${}_{-\infty }{D}_x^{\alpha }f\left(x\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\displaystyle \frac{{\partial }^n}{\partial x^n}}{\int }_{-\infty }^x{\displaystyle \frac{f\left(\xi \right)}{{(x-\xi )}^{\alpha -n+1}}}d\xi ,$$

(2)

$${}_x{D}_{\infty }^{\alpha }f\left(x\right)={\displaystyle \frac{{(-1)}^n}{\mathsf{\Gamma }(n-\alpha )}}{\displaystyle \frac{{\partial }^n}{\partial x^n}}{\int }_x^{\infty }{\displaystyle \frac{f\left(\xi \right)}{{(\xi -x)}^{\alpha -n+1}}}d\xi ,$$

(3)

where $\mathsf{\Gamma }(.)$ denotes Euler’s gamma function and $n=1+\left[\alpha \right]$ such that $\left[\alpha \right]=\max \{m\in \mathbb{Z}|m\le \alpha \}$ and $\mathbb{Z}$ is the set of integers. For instance, for $1<\alpha \le 2$, we get $n=2$, and the parameters l and r are the weighting factors defined as:

$$l=-{\displaystyle \frac{1-\beta }{2cos(\alpha \pi /2)}},r=-{\displaystyle \frac{1+\beta }{2cos(\alpha \pi /2)}},$$

(4)

where the parameter $\beta \in [-1,1]$ is a skewness parameter that shows a preferred direction of displacements that can be seen in heterogeneous systems. When $\beta =0$, the distribution is symmetric and space derivative represents a symmetric Riesz derivative.

Since the Riemann–Liouville derivatives are singular on the boundaries of a bounded domain $[0,L]$, where $L>0$, instead, Caputo derivatives can be used. The left and right Caputo derivatives of order $\alpha$ are defined as follow [49,50]:

$${}_0^c{D}_x^{\alpha }f\left(x\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\int }_0^x{\displaystyle \frac{{\partial }^{n}f\left(\xi \right)/\partial {\xi }^n}{{(x-\xi )}^{\alpha -n+1}}}d\xi ,$$

(5)

$${}_x^c{D}_L^{\alpha }f\left(x\right)={\displaystyle \frac{{(-1)}^n}{\mathsf{\Gamma }(n-\alpha )}}{\int }_x^L{\displaystyle \frac{{\partial }^{n}f\left(\xi \right)/\partial {\xi }^n}{{(\xi -x)}^{\alpha -n+1}}}d\xi .$$

(6)

The $\lambda -$truncated (tempered) fractional-order derivative operator of order $\alpha$ with a shorthand notation ${\mathcal{D}}_x^{\alpha ,\lambda }$ is defined as [44].

For $0<\alpha <1$,

$${\mathcal{D}}_x^{\alpha ,\lambda }f\left(x\right)=l{e}^{-\lambda x}{}_{-\infty }{D}_x^{\alpha }\left(e^{\lambda x}f\left(x\right)\right)+r{e}^{\lambda x}{}_x{D}_{\infty }^{\alpha }\left(e^{-\lambda x}f\left(x\right)\right)+{\displaystyle \frac{{\lambda }^{\alpha }}{cos(\alpha \pi /2)}}f\left(x\right),$$

(7)

For $1<\alpha <2$,

$$\begin{array}{ccc}\hfill {\mathcal{D}}_x^{\alpha ,\lambda }f\left(x\right)& =& l{e}^{-\lambda x}{}_{-\infty }{D}_x^{\alpha }\left(e^{\lambda x}f\left(x\right)\right)+r{e}^{\lambda x}{}_x{D}_{\infty }^{\alpha }\left(e^{-\lambda x}f\left(x\right)\right)+{\displaystyle \frac{{\lambda }^{\alpha }}{cos(\alpha \pi /2)}}f\left(x\right)\hfill \\ & +& {\displaystyle \frac{\alpha \beta {\lambda }^{\alpha -1}}{|cos(\alpha \pi /2\left)\right|}}{\displaystyle \frac{\partial f}{\partial x}},\hfill \end{array}$$

(8)

where the operators ${}_{-\infty }{D}_x^{\alpha }$ and ${}_x{D}_{\infty }^{\alpha }$ and the parameters l and r are defined by Equations (2)–(4), respectively, and $\lambda \ge 0$ is a truncation parameter with SI units m${}^{-1}$. Now, for $1<\alpha <2$, we can define the left and right-sided truncated fractional-order operators ${}_{-\infty }{\mathcal{D}}_x^{\alpha ,\lambda }$ and ${}_x{\mathcal{D}}_{\infty }^{\alpha ,\lambda }$, for $\beta =-1$ and $\beta =1$, respectively, as follows:

$${}_{-\infty }{\mathcal{D}}_x^{\alpha ,\lambda }f\left(x\right)={\displaystyle \frac{-1}{cos(\alpha \pi /2)}}\left(e^{-\lambda x}{}_{-\infty }{D}_x^{\alpha }\left(e^{\lambda x}f\left(x\right)\right)-{\lambda }^{\alpha }f\left(x\right)-\alpha {\lambda }^{\alpha -1}{\displaystyle \frac{\partial f}{\partial x}}\right),$$

(9)

$${}_x{\mathcal{D}}_{\infty }^{\alpha ,\lambda }f\left(x\right)={\displaystyle \frac{-1}{cos(\alpha \pi /2)}}\left(e^{\lambda x}{}_x{D}_{\infty }^{\alpha }\left(e^{-\lambda x}f\left(x\right)\right)-{\lambda }^{\alpha }f\left(x\right)+\alpha {\lambda }^{\alpha -1}{\displaystyle \frac{\partial f}{\partial x}}\right),$$

(10)

3. Spatial Propagation of an Epidemic

In this section, we investigate the effect of different random movements of individuals on the spatial spread and speed of an epidemic. To do so, we consider a simple version of the epidemic models. The model takes into account the densities of two populations at location x at time t, the susceptibles $S(x,t)$ and the infectives $I(x,t)$. We assume that the susceptibles catch the disease from the infectives at a constant transmission-efficiency rate $\delta$ and that the infectives disappear either by recovery from the disease or disease-induced mortality, at a constant rate $\gamma$. We model the random movement of susceptibles and infectives by considering diffusion operators in the model equations. With these assumptions, the densities of susceptiples and infectives vary according to the following equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial S}{\partial t}}& =& -\delta S(x,t)I(x,t)+\underset{\underset{\mathrm{random}\mathrm{motion}}{⏟}}{\mathrm{diffusion}\mathrm{operator}},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial I}{\partial t}}& =& \delta S(x,t)I(x,t)-\gamma I(x,t)+\underset{\underset{\mathrm{random}\mathrm{motion}}{⏟}}{\mathrm{diffusion}\mathrm{operator}},\hfill \end{array}$$

(12)

In what follows, we shall formulate “diffusion operator” by considering three patterns of random mobility of individuals, namely Brownian motion, pure (untruncated) Lévy flight and truncated Lévy flight, and then discuss the epidemic speed in each case.

3.1. Brownian Motion

Brownian motion is a normal dispersion process that the probability distribution function (pdf) for the displacements is a Gaussian (Normal) distribution. In that case, if we assume that the mobility of susceptible and infective individuals is described by Brownian motion, the diffusion operators in Equations (11) and (12) are the second-order diffusion denoted $K{\partial }^{2}S(x,t)/\partial x^2$ and $K{\partial }^{2}I(x,t)/\partial x^2$, respectively, where K represents the diffusion coefficient and has units of m${}^2$s${}^{-1}$. Therefore, we get

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial S}{\partial t}}& =& -\delta S(x,t)I(x,t)+K{\displaystyle \frac{{\partial }^{2}S(x,t)}{\partial x^2}},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial I}{\partial t}}& =& \delta S(x,t)I(x,t)-\gamma I(x,t)+K{\displaystyle \frac{{\partial }^{2}I(x,t)}{\partial x^2}},\hfill \end{array}$$

(14)

To obtain the epidemic speed, we consider a situation in which the wave of infectives propagates into the susceptible population whose density is uniformly set to be $S(x,0)=S_0$ throughout the entire domain. In that case, the constant epidemic speed at the leading edge, i.e., where $S\approx S_0$ and $I\approx 0$ is given by $c=2\sqrt{K(\delta S_0-\gamma )}$ for $K(\delta S_0-\gamma )>0$. For more details, refer to [15,17].

3.2. Pure (Untruncated) Lévy Flights

For pure Lévy flights, the power-law asymptotic behaviour of the pdf for the displacements leads to the space fractional-order diffusion operator ${\mathcal{D}}_x^{\alpha }$ defined by Equation (1). If we assume that the susceptible and infective populations perform Lévy flights, the fractional-order epidemic model reads

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial S}{\partial t}}& =& -\delta S(x,t)I(x,t)+K_{\alpha }{\mathcal{D}}_x^{\alpha },\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial I}{\partial t}}& =& \delta S(x,t)I(x,t)-\gamma I(x,t)+K_{\alpha }{\mathcal{D}}_x^{\alpha },\hfill \end{array}$$

(16)

where $K_{\alpha }$ is the fractional-order diffusivity with units of m${}^{\alpha }$s${}^{-1}$. Hanert et al. [40] investigated the effect of the left-sided fractional diffusion operator on the left-and right-moving fronts. It was shown that the epidemic waves travel to the right side at an unbounded and exponential speed given by $c\left(t\right)\sim (1-\lambda )e^{\left(\right(1-\lambda )/(\alpha +1\left)\right)t}$ for large values of t and to the left side at a constant speed given by $c=\alpha {((1-\lambda )/(\alpha -1))}^{(\alpha -1)/\alpha }$, where $\lambda =\gamma /\delta S_0$. It should be noted that in the case of the symmetric fractional operator $(\beta =0)$, the infective waves propagating to the left and right sides have the same speed, i.e., both the left- and right-moving waves accelerate exponentially.

3.3. Truncated Lévy Flights

In this subsection, we investigate the effect of the truncated Lévy flights on the epidemic speed. To do so, we consider the left-sided truncated fractional-order diffusion operator ${}_{-\infty }{\mathcal{D}}_x^{\alpha ,\lambda }$ defined by Equation (9). Since the movement of susceptible individuals does not change the front speed, for the sake of simplicity, we suppose that there is no diffusion term in the S equation. However, infectious individuals move into the susceptible population. The model equations then read as follows:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial S(x,t)}{\partial t}}& =& -\delta S(x,t)I(x,t),\hfill \\ \hfill {\displaystyle \frac{\partial I(x,t)}{\partial t}}& =& \delta S(x,t)I(x,t)-\gamma I(x,t)\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& +& D_{\alpha }\left(e^{-\lambda x}{}_{-\infty }{D}_x^{\alpha }\left(e^{\lambda x}I(x,t)\right)-{\lambda }^{\alpha }I(x,t)\right),\hfill \end{array}$$

(18)

where $D_{\alpha }=K_{\alpha }/|cos(\alpha \pi /2)|$. It should be noted that similar to the assumption used by del Castilo Negrete [46] for the truncated fractional-order Fisher–Kolmogorov equation, we do not consider the term $\partial I(x,t)/\partial x$ in the I equation. In order to analyse the solutions of Equations (17) and (18), we first make them dimensionless by introducing the following variables:

$$\overline{S}={\displaystyle \frac{S}{S_0}},\overline{I}={\displaystyle \frac{I}{S_0}},\overline{x}={\displaystyle \frac{x}{x_c}},\overline{t}={\displaystyle \frac{t}{t_c}},$$

where $x_c$ and $t_c$ are characteristic values of x and t, respectively. In what follows, we estimate these values. The next task is to insert the new dimensionless variables into Equations (17) and (18). That is, we replace S by $S_0\overline{S}$, I by $S_0\overline{I}$, t by $t_c\overline{t}$ and x by $x_c\overline{x}$. The derivative with respect to $\overline{t}$ is derived through the chain rule as

$${\displaystyle \frac{\partial I}{\partial t}}={\displaystyle \frac{\partial \left(S_0\overline{I}\right)}{\partial \overline{t}}}{\displaystyle \frac{d\overline{t}}{dt}}=S_0{\displaystyle \frac{\partial \overline{I}}{\partial \overline{t}}}{\displaystyle \frac{1}{t_c}}={\displaystyle \frac{S_0}{t_c}}{\displaystyle \frac{\partial \overline{I}}{\partial \overline{t}}},{\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{S_0}{t_c}}{\displaystyle \frac{\partial \overline{S}}{\partial \overline{t}}},$$

and the left fractional-order derivative with respect to $\overline{x}$, ${}_{-\infty }{D}_{\overline{x}}^{\alpha }$ is also derived as follows:

$${}_{-\infty }{D}_x^{\alpha }\left(e^{\lambda x}I\right)=S_0{}_{-\infty }{D}_{\overline{x}}^{\alpha }\left(e^{\lambda x_c\overline{x}}\overline{I}\right){\left({\displaystyle \frac{d\overline{x}}{dx}}\right)}^{\alpha }={\displaystyle \frac{S_0}{x_c^{\alpha }}}{}_{-\infty }{D}_{\overline{x}}^{\alpha }\left(e^{\lambda x_c\overline{x}}\overline{I}\right).$$

Therefore, Equations (17) and (18) now become

$$\begin{array}{ccc}\hfill {\displaystyle \frac{S_0}{t_c}}{\displaystyle \frac{\partial \overline{S}}{\partial \overline{t}}}& =& -\delta S_0^2\overline{S}\overline{I},\hfill \\ \hfill {\displaystyle \frac{S_0}{t_c}}{\displaystyle \frac{\partial \overline{I}}{\partial \overline{t}}}& =& \delta S_0^2\overline{S}\overline{I}-\gamma S_0\overline{I}+{\displaystyle \frac{D_{\alpha }{S}_0}{x_c^{\alpha }}}{e}^{-\lambda x_c\overline{x}}{}_{-\infty }{D}_{\overline{x}}^{\alpha }\left(e^{\lambda x_c\overline{x}}\overline{I}\right)-D_{\alpha }{\lambda }^{\alpha }{S}_0\overline{I},\hfill \end{array}$$

By choosing $t_c={\displaystyle \frac{1}{\delta S_0}}$ and $x_c={\left({\displaystyle \frac{D_{\alpha }}{\delta S_0}}\right)}^{1/\alpha }$, we get the following model in a dimensionless form:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial S(x,t)}{\partial t}}& =& -S(x,t)I(x,t),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial I(x,t)}{\partial t}}& =& S(x,t)I(x,t)-\theta I(x,t)+e^{-\mu x}{}_{-\infty }{D}_x^{\alpha }\left(e^{\mu x}I(x,t)\right),\hfill \end{array}$$

(20)

where we have dropped the overbar “−”, $\theta ={\displaystyle \frac{\gamma +D_{\alpha }{\lambda }^{\alpha }}{\delta S_0}}$, and $\mu =\lambda {\left({\displaystyle \frac{D_{\alpha }}{\delta S_0}}\right)}^{1/\alpha }$ are dimensionless parameters.

3.3.1. Theoretical Analysis (Right-Propagating Front)

Here, we look for the speed of the infective waves propagating to the right side. Similar to the theoretical analysis of the model studied in [40], we replace $S(x,t)$ by $1-s(x,t)$, where $s(x,t)$ is the deviation with respect to the initial density of the susceptible population. Since the epidemic speed can be obtained by the speed of the right-moving fronts in the leading edge region, i.e., where $s\approx 0$ and $I\approx 0$, Equations (19) and (20) in that region can be expressed as

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial s(x,t)}{\partial t}}& =& I(x,t),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial I(x,t)}{\partial t}}& =& (1-\theta )I(x,t)+e^{-\mu x}{}_{-\infty }{D}_x^{\alpha }\left(e^{\mu x}I(x,t)\right).\hfill \end{array}$$

(22)

Here, we assume that $1-\theta >0$, which guarantees the existence of the right-moving fronts.

Now, we find the analytical solution of Equation (22), by defining the following localized initial condition:

$$I(x,t=0)=\left\{\begin{array}{cc}A\hfill & \mathrm{if}x<0\hfill \\ A{e}^{-\nu x}\hfill & \mathrm{if}x\ge 0,\hfill \end{array}\right.$$

(23)

where A and $\nu$ are positive and non-zero constants. We next follow the same method as the one used by del Castilo Negrete [46] for the truncated fractional-order Fisher–Kolmogorov equation. First, we assume a solution of Equation (22) in the form of $I(x,t)=e^{-\mu x+(1-\theta )t}\mathsf{\Phi }(x,t)$. By substituting this solution into (22) and initial condition (23), one gets the following fractional-order diffusion equation

$${\displaystyle \frac{\partial \mathsf{\Phi }(x,t)}{\partial t}}{=}_{-\infty }{D}_x^{\alpha }\mathsf{\Phi }(x,t),$$

(24)

with the initial condition

$${\mathsf{\Phi }}_0\left(x\right)=\mathsf{\Phi }(x,t=0)=\left\{\begin{array}{cc}A{e}^{\mu x}\hfill & \mathrm{if}x<0\hfill \\ A{e}^{(\mu -\nu )x}\hfill & \mathrm{if}x\ge 0.\hfill \end{array}\right.$$

(25)

Equation (24) with initial condition (25) has the general solution of the form

$$\mathsf{\Phi }(x,t)={\int }_{-\infty }^{\infty }G(x-y,t){\mathsf{\Phi }}_0\left(y\right)dy,$$

(26)

where $G(x,t)=t^{-1/\alpha }{p}_{\alpha }(x/t^{1/\alpha })$ is the Green function of Equation (24) and $p_{\alpha }\left(\eta \right)=(1/2\pi ){\int }_{-\infty }^{\infty }{e}^{{\left(ik\right)}^{\alpha }+ik\eta }dk$ is a skewed Lévy distribution with exponent $\alpha$. By substituting initial condition (25) into Equation (26), we get

$$\begin{array}{ccc}\hfill \mathsf{\Phi }(x,t)& =& {\int }_{-\infty }^{x{t}^{-1/\alpha }}A{e}^{(\mu -\nu )(x-t^{1/\alpha }\eta )}{p}_{\alpha }\left(\eta \right)d\eta +{\int }_{x{t}^{-1/\alpha }}^{\infty }A{e}^{\mu (x-t^{1/\alpha }\eta )}{p}_{\alpha }\left(\eta \right)d\eta \hfill \\ & =& A{e}^{(\mu -\nu )x}{\int }_{-\infty }^{x/\tau }{e}^{(\nu -\mu )\tau \eta }{p}_{\alpha }\left(\eta \right)d\eta +A{e}^{\mu x}{\int }_{x/\tau }^{\infty }{e}^{-\mu \tau \eta }{p}_{\alpha }\left(\eta \right)d\eta ,\hfill \end{array}$$

(27)

where $\tau =t^{1/\alpha }$. Now, one can find the solution of Equation (22) in the following form

$$I(x,t)=A{e}^{-\nu x+(1-\theta )t}{I}_1+A{e}^{(1-\theta )t}{I}_2,$$

(28)

where ${\displaystyle I_1={\int }_{-\infty }^{x/\tau }{e}^{(\nu -\mu )\tau \eta }{p}_{\alpha }\left(\eta \right)d\eta }$ and ${\displaystyle I_2={\int }_{x/\tau }^{\infty }{e}^{-\mu \tau \eta }{p}_{\alpha }\left(\eta \right)d\eta }$.

The next task is to look for the asymptotic behaviour of the integrals $I_1$ and $I_2$ for $x/\tau \to \infty$ with a fixed $\tau$. We first consider the asymptotic behaviour of $I_2$. By considering the fact that for $\eta >0$ the Lévy distribution $p_{\alpha }\left(\eta \right)$ exhibits fat tails as ∼${\eta }^{-(\alpha +1)}$ [46,51], and integration by parts, we get

$$\begin{array}{ccc}\hfill I_2\sim {\int }_{x/\tau }^{\infty }{e}^{-\mu \tau \eta }{\eta }^{-(\alpha +1)}d\eta & =& {\left[{\displaystyle \frac{-1}{\mu \tau }}{e}^{-\mu \tau \eta }{\eta }^{-(\alpha +1)}\right]}_{\eta =x/\tau }^{\eta =\infty }-{\displaystyle \frac{\alpha +1}{\mu \tau }}{\int }_{x/\tau }^{\infty }{e}^{-\mu \tau \eta }{\eta }^{-(\alpha +2)}d\eta \hfill \\ & =& {\displaystyle \frac{{\tau }^{\alpha }}{\mu }}{\displaystyle \frac{e^{-\mu x}}{x^{\alpha +1}}}-{\displaystyle \frac{\alpha +1}{\mu \tau }}{\int }_{x/\tau }^{\infty }{e}^{-\mu \tau \eta }{\eta }^{-(\alpha +2)}d\eta \sim {\displaystyle \frac{{\tau }^{\alpha }}{\mu }}{\displaystyle \frac{e^{-\mu x}}{x^{\alpha +1}}}.\hfill \end{array}$$

(29)

For the first integral $I_1$, if we define a cutoff $\mathsf{\Omega }$ with $1\ll \mathsf{\Omega }<x\tau$ so that $p_{\alpha }\left(\eta \right)\sim {\eta }^{-(\alpha +1)}$, we can then write the integral $I_1$ as

$$I_1\sim {\int }_{-\infty }^{\mathsf{\Omega }}{e}^{(\nu -\mu )\tau \eta }{p}_{\alpha }\left(\eta \right)d\eta +{\int }_{\mathsf{\Omega }}^{x/\tau }{e}^{(\nu -\mu )\tau \eta }{\eta }^{-(\alpha +1)}d\eta .$$

(30)

Because of the exponential decay of $p_{\alpha }\left(\eta \right)$ at minus infinity (for more details, refer to [46,51]), the first integral converges to a finite value C. By applying an integration by parts to the second integral of Equation (30), we get

$$\begin{array}{ccc}\hfill {\int }_{\mathsf{\Omega }}^{x/\tau }{e}^{(\nu -\mu )\tau \eta }{\eta }^{-(\alpha +1)}d\eta & =& {\left[{\displaystyle \frac{1}{(\nu -\mu )\tau }}{e}^{(\nu -\mu )\tau \eta }{\eta }^{-(\alpha +1)}\right]}_{\eta =\mathsf{\Omega }}^{\eta =x/\tau }\hfill \\ & +& {\displaystyle \frac{\alpha +1}{(\nu -\mu )\tau }}{\int }_{\mathsf{\Omega }}^{x/\tau }{e}^{(\nu -\mu )\tau \eta }{\eta }^{-(\alpha +2)}d\eta \sim {\displaystyle \frac{{\tau }^{\alpha }}{\nu -\mu }}{\displaystyle \frac{e^{(\nu -\mu )x}}{x^{\alpha +1}}}\hfill \end{array}$$

Therefore,

$$I_1\sim C+{\displaystyle \frac{{\tau }^{\alpha }}{\nu -\mu }}{\displaystyle \frac{e^{(\nu -\mu )x}}{x^{\alpha +1}}}.$$

(31)

By substituting Equations (29) and (31) into Equation (28), we get

$$I(x,t)\sim AC{e}^{-\nu x+(1-\theta )t}+A\left({\displaystyle \frac{1}{\nu -\mu }}+{\displaystyle \frac{1}{\mu }}\right){\displaystyle \frac{t}{x^{\alpha +1}}}{e}^{-\mu x+(1-\theta )t},\nu \ne \mu ,\mu \ne 0.$$

(32)

In order to find the asymptotic behaviour of $I(x,t)$ for large x, one can find the dominant term in Equation (32) by considering the relative values of $\nu$ and $\mu$.

If $\mu <\nu$, the truncation decays slower than the initial condition. In that case, the second term of Equation (32) is the dominant term for large x, i.e.,

$$I(x,t)\sim A\left({\displaystyle \frac{1}{\nu -\mu }}+{\displaystyle \frac{1}{\mu }}\right){\displaystyle \frac{t}{x^{\alpha +1}}}{e}^{-\mu x+(1-\theta )t}.$$

(33)

By considering the expansion of the exponential function $e^{-\mu x}$, we get

$$I(x,t)\sim A\left({\displaystyle \frac{1}{\nu -\mu }}+{\displaystyle \frac{1}{\mu }}\right)t{e}^{(1-\theta )t}{\displaystyle \frac{1}{x^{\alpha +1}}}\left(1-\mu x+{\displaystyle \frac{{\mu }^2}{2}}{x}^2-\cdots \right).$$

(34)

For $\mu x\ll 1$, the solution of $I(x,t)$ thus decays algebraically as

$$I(x,t)\sim t{e}^{(1-\theta )t}{x}^{-(\alpha +1)},$$

(35)

However, for $\mu x\gg 1$, the algebraic decaying tail of the solution is tempered with the exponential factor $e^{-\mu x}$ as

$$I(x,t)\sim t{e}^{-\mu x+(1-\theta )t}{x}^{-(\alpha +1)}.$$

(36)

If $\mu >\nu$, the truncation decays faster than the initial condition. In that case, the first term of Equation (32) is the dominant term for large x, i.e.,

$$I(x,t)\sim e^{-\nu x+(1-\theta )t},$$

(37)

which shows an exponentially decaying tail of the solution.

Until now, we have assumed that $\nu \ne \mu$. To explore the case $\nu =\mu$, we cannot use Equation (32). Instead, we need to go back to Equation (30) and write the integral $I_1$ as

$$I_1\sim {\int }_{-\infty }^{\mathsf{\Omega }}{p}_{\alpha }\left(\eta \right)d\eta +{\int }_{\mathsf{\Omega }}^{x/\tau }{\eta }^{-(\alpha +1)}d\eta .$$

(38)

By defining the finite value ${\displaystyle C={\int }_{-\infty }^{\mathsf{\Omega }}{p}_{\alpha }\left(\eta \right)d\eta }$ and solving the second integral of Equation (38), we have

$$I_1\sim C-{\displaystyle \frac{{\tau }^{\alpha }}{\alpha }}{\displaystyle \frac{1}{x^{\alpha }}}.$$

(39)

By substituting Equations (29) and (39) into (28), we obtain the solution $I(x,t)$ as

$$I(x,t)\sim AC{e}^{-\nu x+(1-\theta )t}+A\left({\displaystyle \frac{1}{\mu x^{\alpha +1}}}-{\displaystyle \frac{1}{\alpha x^{\alpha }}}\right)t{e}^{-\nu x+(1-\theta )t}.$$

(40)

In this case, the first term of Equation (40) decays slower than the other term. Therefore, it is the dominant term for large x. We conclude that for the case $\nu =\mu$, the solution tail exhibits an exponentially decaying as

$$I(x,t)\sim e^{-\nu x+(1-\theta )t}$$

(41)

Obviously, the solution of Equation (21), the susceptibles equation, is obtained from the solution of the infectives equation. Thus, the asymptotic behaviour of the susceptible front is the same as the infective front.

The above results are based on the assumption $\mu \ne 0$. However, in the absence of the truncation parameter, i.e., $\mu =0$, we have to explore the asymptotic behaviour of the following integrals

$$I_1={\int }_{-\infty }^{x/\tau }{e}^{\nu \tau \eta }{p}_{\alpha }\left(\eta \right)d\eta ,I_2={\int }_{x/\tau }^{\infty }{p}_{\alpha }\left(\eta \right)d\eta .$$

(42)

The analysis of the asymptotic behaviour of these integrals has been studied in [40] and it was shown that the solution $I(x,t)$ has the same algebraically decaying tail as Equation (35).

In order to derive the epidemic speed, we first compute a Lagrangian trajectory $x_L\left(t\right)$ of infective fronts such that $I(x_L\left(t\right),t)=I_0$, where $I_0\approx 0$. Given $x_L\left(t\right)$, the epidemic speed is given by $c\left(t\right)=d{x}_L\left(t\right)/dt$.

In the case of a power-law decaying tail, we find $x_L\left(t\right)$ from Equation (35) such that $I(x_L\left(t\right),t)=t{e}^{(1-\theta )t}{x}_L{\left(t\right)}^{-(\alpha +1)}=I_0$. After some calculation, $x_L\left(t\right)$ and $c\left(t\right)$ can be computed as

$$x_L\left(t\right)=I_0^{-1/(\alpha +1)}{t}^{1/(\alpha +1)}{e}^{\left(\right(1-\theta )/(\alpha +1\left)\right)t},c\left(t\right)\sim (1-\theta )e^{\left(\right(1-\theta )/(\alpha +1\left)\right)t}.$$

(43)

For the truncated algebraically decaying tail, we use Equation (36) to find $I(x_L\left(t\right),t)=t{e}^{-\mu x_L\left(t\right)+(1-\theta )t}{x}_L{\left(t\right)}^{-(\alpha +1)}=I_0$. Thus, $x_L\left(t\right)$ satisfies the following equation

$$lnt+(1-\theta )t-\mu x_L\left(t\right)-(\alpha +1)ln{x}_L\left(t\right)=ln{I}_0,$$

(44)

and $c\left(t\right)$ can be expressed as

$$c\left(t\right)={\displaystyle \frac{(1-\theta )+{\displaystyle \frac{1}{t}}}{\mu +{\displaystyle \frac{\alpha +1}{x_L\left(t\right)}}}}.$$

(45)

When $t\to \infty$, the velocity of the fronts converges towards a constant value representing the maximum value of the epidemic speed given by $c_{max}=(1-\theta )/\mu$.

Finally, in the case of exponentially decaying tail, we can use Equation (37) to find $I(x_L\left(t\right),t)=e^{-\nu x_L\left(t\right)+(1-\theta )t}=I_0$. Thus, $x_L\left(t\right)=(-ln{I}_0+(1-\theta )t)/\nu$ and the front moves at a constant velocity $c\left(t\right)=\overline{c}=(1-\theta )/\nu$.

3.3.2. Theoretical Analysis (Left-Propagating Front)

For exploring the impact of the left-sided truncated operator defined by Equation (9) on the left-moving fronts, we consider a travelling-wave solution of a constant shape as follows:

$$I(x,t)=I\left(z\right),s(x,t)=s\left(z\right),z=x+ct,$$

where c represents the speed of the epidemic waves moving to the left side. Substituting this solution into Equations (21) and (22), one can find the following system:

$$\begin{array}{ccc}\hfill c{\displaystyle \frac{\partial s}{\partial z}}& =& I,\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill c{\displaystyle \frac{\partial I}{\partial z}}& =& (1-\theta )I+e^{-\mu z}{}_{-\infty }{D}_z^{\alpha }\left(e^{\mu z}I\right).\hfill \end{array}$$

(47)

Now we define a solution of the following form

$$I\left(z\right)=\widehat{I}{e}^{kz},s\left(z\right)=\widehat{s}{e}^{kz},$$

(48)

where k is a parameter depending on the value of c. Considering the fact that ${}_{-\infty }{D}_z^{\alpha }{e}^{kz}=k^{\alpha }{e}^{kz}$ and substituting the solutions (48) into Equations (46) and (47), we get the following system:

$$\left[\begin{array}{cc}ck& -1\\ 0& ck-{(k+\mu )}^{\alpha }-(1-\theta )\end{array}\right]\left[\begin{array}{c}\widehat{s}\\ \widehat{I}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$

(49)

If $ck(ck-{(k+\mu )}^{\alpha }-(1-\theta ))=0$, we can find a non-trivial solution. In that case, we get the following relation representing the value of the epidemic speed as a function of the parameter k:

$$c\left(k\right)={\displaystyle \frac{1-\theta }{k}}+{\displaystyle \frac{{(k+\mu )}^{\alpha }}{k}}.$$

(50)

In order to obtain the minimum epidemic speed, one can minimize Equation (50) with respect to k. For $0<\alpha <2$ and $\mu \ne 0$, we cannot analytically calculate the minimum speed $c_{min}$. In the next section, we thus find it numerically. However, for $0<\alpha <2$ and $\mu =0$, as it is discussed in [40], $c_{min}=\alpha {\left((1-\theta )/(\alpha -1)\right)}^{(\alpha -1)/\alpha }$ and the corresponding rate of the exponential decay is equal to $k_{min}={\left((1-\theta )/(\alpha -1)\right)}^{1/\alpha }$. For $\alpha =2$, the minimum speed of the classical model is given by $c_{min}=2\sqrt{1-\theta }$.

4. Numerical Examples

In order to illustrate the analytical results obtained for the asymptotic behaviour of the infective waves, we numerically solve the model equations on a finite domain $[0,L]$, where $L>0$. Here, we have used the same method as the one performed by Vallaeys et al. [29] for solving a symmetric truncated fractional-order diffusion equation. To study more numerical methods for solving the space fractional-order diffusion equations, refer to [52,53,54]. We discretize Equations (19) and (20) by finite-element (FE) method based on a Galerkin formulation. To illustrate this method more precisely, we discretize the following simple diffusion equation:

$${\displaystyle \frac{\partial I(x,t)}{\partial t}}=e^{-\lambda x}{}_0{D}_x^{\alpha }\left(e^{\lambda x}I(x,t)\right)-{\lambda }^{\alpha }I(x,t).$$

(51)

In this method, the exact solution (unknown variable) $I(x,t)$ is expressed as a sum of the unknown coefficients $e_j\left(t\right)$ and basis functions ${\varphi }_j\left(x\right)$ as follows:

$$I(x,t)\approx \tilde{I}(x,t)=\sum _{j=1}^N{e}_j\left(t\right){\varphi }_j\left(x\right).$$

(52)

By introducing a partition of the domain $[0,L]$ into $N-1$ subintervals $[x_j,x_{j+1}]$ with a constant length h, i.e., $x_1=0,x_N=L,$ and $x_{j+1}-x_j=h$ for $j=1,\cdots ,N-1$, we can consider the piecewise linear basis functions ${\varphi }_j\left(x\right)$ for $j=1,2,\cdots ,N$ as follows:

$${\varphi }_1\left(x\right)=\left\{\begin{array}{cc}{\displaystyle \frac{x_2-x}{x_2-x_1}}\hfill & :x_1\le x\le x_2,\hfill \\ 0\hfill & :x\notin [x_1,x_2].\hfill \end{array}\right.$$

for $j=2,\dots ,N-1$,

$${\varphi }_j\left(x\right)=\left\{\begin{array}{cc}{\displaystyle \frac{x-x_{j-1}}{x_j-x_{j-1}}}\hfill & :x_{j-1}\le x\le x_j,\hfill \\ {\displaystyle \frac{x_{j+1}-x}{x_{j+1}-x_j}}\hfill & :x_j\le x\le x_{j+1},\hfill \\ 0\hfill & :x\notin [x_{j-1},x_{j+1}].\hfill \end{array}\right.$$

and for $j=N$,

$${\varphi }_N\left(x\right)=\left\{\begin{array}{cc}{\displaystyle \frac{x-x_{N-1}}{x_N-x_{N-1}}}\hfill & :x_{N-1}\le x\le x_N,\hfill \\ 0\hfill & :x\notin [x_{N-1},x_N].\hfill \end{array}\right.$$

We also discretize the exponentially truncated unknown solution, i.e., $e^{\lambda x}I(x,t)$ as follows

$$e^{\lambda x}I(x,t)\approx {\tilde{I}}^l(x,t)=\sum _{j=1}^N{e}_j^l\left(t\right){\varphi }_j\left(x\right).$$

(53)

where $e_j^l\left(t\right)$ is the unknown coefficient corresponding to the left truncated solution.

By using a Galerkin formulation, we replace $\tilde{I}(x,t)$ and ${\tilde{I}}^l(x,t)$ in Equation (51) and then by orthogonalizing the discrete equation with respect to all ${\varphi }_j$, we get the following equation:

$$\begin{array}{ccc}\hfill \left({\int }_0^L{\varphi }_i{\varphi }_{j}dx\right){\displaystyle \frac{d{e}_j}{dt}}\left(t\right)& =& \left({\int }_0^L{\varphi }_i{e}^{-\lambda x}{}_0{D}_x^{\alpha }{\varphi }_{j}dx\right)e_j^l\left(t\right)-{\lambda }^{\alpha }\left({\int }_0^L{\varphi }_i{\varphi }_{j}dx\right)e_j\left(t\right),\hfill \end{array}$$

for $i,j=1,\dots ,N$.

By introducing the matrices $\mathbf{M}:={\int }_0^L{\varphi }_i{\varphi }_{j}dx$ and $\mathbf{D}:={\int }_0^L{\varphi }_i{e}^{-\lambda x}{}_0{D}_x^{\alpha }{\varphi }_{j}dx$, we get the following semi-discrete equation in a matrix form

$$\mathbf{M}{\displaystyle \frac{d\mathbf{e}}{dt}}\left(t\right)=\mathbf{D}{\mathbf{e}}^l\left(t\right)-{\lambda }^{\alpha }\mathbf{M}\mathbf{e}\left(t\right),$$

(54)

where $\mathbf{e}\left(t\right)={\left[e_1\left(t\right)\cdots e_N\left(t\right)\right]}^T$ and ${\mathbf{e}}^l\left(t\right)={\left[e_1^l\left(t\right)\cdots e_N^l\left(t\right)\right]}^T$ are the vector of unknown coefficients at time t. It should be noted that by considering the same method in [55], the diffusion matrix $\mathbf{D}$ can be expressed as: $\mathbf{D}=-{\int }_0^L{e}^{-\lambda x}({\displaystyle \frac{d{\varphi }_i}{dx}}-\lambda {\varphi }_i){}_0{D}_x^{\alpha -1}{\varphi }_{j}dx$. Here we use the left Caputo derivative of order $\alpha -1$ of ${\varphi }_j$ defined by Equation (5) as Caputo derivatives are easier to handle and also for many applications could prevent mass-balance errors on bounded domains, while such errors can be made by Riemann–Liouville derivatives (for more details, see [56]).

In order to solve Equation (54), the vector ${\mathbf{e}}^l\left(t\right)$ needs to be computed in terms of the vector $\mathbf{e}\left(t\right)$. To do so, we use Galerkin formulation for Equations (52) and (53) as follows:

$${\int }_0^L{\varphi }_i{e}^{\lambda x}\tilde{I}dx={\int }_0^L{\varphi }_i{\tilde{I}}^{l}dx\to \left({\int }_0^L{\varphi }_i{e}^{\lambda x}{\varphi }_{j}dx\right)e_j\left(t\right)=\left({\int }_0^L{\varphi }_i{\varphi }_{j}dx\right)e_j^l\left(t\right).$$

By defining the matrix $\mathbf{W}:={\int }_0^L{\varphi }_i{e}^{\lambda x}{\varphi }_{j}dx$, we get ${\mathbf{e}}^l\left(t\right)={\mathbf{M}}^{-1}\mathbf{W}\mathbf{e}\left(t\right)$. we can thus express Equation (54) as follows

$$\mathbf{M}{\displaystyle \frac{d\mathbf{e}}{dt}}\left(t\right)=\left(\mathbf{D}{\mathbf{M}}^{-1}\mathbf{W}-{\lambda }^{\alpha }\mathbf{M}\right)\mathbf{e}\left(t\right),$$

(55)

Finally, in order to discretize Equation (55), we use a third-order Adams–Bashforth method.

In order to solve the model Equations (19) and (20), we use the following initial conditions:

$$S(x,t=0)=S_0\left(x\right)=1-h_S\left(1\pm tanh\left({\displaystyle \frac{x-x_0}{w_{s}L}}\right)\right),$$

$$I(x,t=0)=I_0\left(x\right)=\left\{\begin{array}{cc}{h}_I{e}^{-\nu (x-x_0)}\hfill & :x\ge x_0,\hfill \\ h_I{e}^{\nu (x-x_0)}\hfill & :x\le x_0,\hfill \end{array}\right.$$

where $h_S=0.37,w_S=0.003,h_I=0.14,$ and $\nu =0.001$. we consider the “+” sign and $x_0=9L/10$ for the left-moving front and the “−” sign and $x_0=L/30$ for the right-moving front. The dimensionless length of the domain is set to $L={10}^5$ and $L=3\times {10}^5$ for the left- and right-moving fronts, respectively. For all simulations, $\theta =0.5$, and the value of the fractional order $\alpha$ is set to $1.2$. According to the analytical results, we can consider three different asymptotic behaviours for the right-moving fronts, i.e., $y_1=x^{-(\alpha +1)}$ (power-law decay), $y_2=e^{-\mu x}{x}^{-(\alpha +1)}$ (truncated power-law decay), and $y_3=e^{-\mu x}$ (exponential decay). In what follows, we shall show how different values of the truncation parameter $\mu$ lead to these different asymptotic behaviours.

The algebraic decaying tail of the right-moving fronts are shown in Figure 1. Here, we take $\mu =5\times {10}^{-6}<\nu$. Figure 1a shows the densities of both the susceptibles and infectives at different times. Figure 1b,c show that the right-moving fronts for both the susceptibles and infectives exhibit a power-law decaying tail as $y_1=x^{-(\alpha +1)}$. The Lagrangian trajectory $x_L\left(t\right)$ and the epidemic speed $c\left(t\right)$ corresponding to this case are shown in Figure 2a,b, respectively. Figure 2a shows the time evolution of the Lagrangian trajectory $x_L\left(t\right)$ such that $I(x_L\left(t\right),t)=I_0$, where $I_0=0.01$. Figure 2b shows the time evolution of the numerical estimation of the instantaneous speed of the points with density $I_0=0.01$. Both Figure 2a,b also show that the results obtained numerically agree well with the asymptotic expansions obtained by Equation (43). It should be noted that the truncation parameter here has no influence on the rapid propagation of the epidemic and the fronts move at an exponential speed similar to the fractional-order epidemic model based on pure Lévy flights [40].

By increasing the value of the truncation parameter $\mu$, the algebraically decaying tail of the solutions can be tamed. If we take $\mu ={10}^{-4}<\nu$, both the susceptible and infective fronts propagate so that their tails exhibit a tempered algebraic decaying, i.e., $1-S$ and $I\sim y_2=e^{-\mu x}{x}^{-(\alpha +1)}$ (see Figure 3). Figure 4 shows that the Lagrangian trajectory $x_L\left(t\right)$ and the velocity $c\left(t\right)$ obtained by numerically solving Equations (19) and (20) are in good agreement with the results obtained analytically by Equations (44) and (45). As expected, the infective fronts initially propagate at the minimum average speed $c_{min}=5\times {10}^2$ and then the front speed gradually increases so that after some time, it converges towards the maximum front speed value $c_{max}=5\times {10}^3$.

Figure 5 shows the dynamics of the right-moving susceptible and infective fronts for $\mu ={10}^{-3}$. As expected, for $\mu \ge \nu$, the tails of the solutions for the susceptible and infective fronts exhibit an exponential decay, i.e, $1-S$ and $I\sim y_3=e^{-\nu x}$ (see Figure 5b,c). Figure 6 shows that the Lagrangian trajectory $x_L\left(t\right)$ and the velocity $c\left(t\right)$ obtained by numerically solving the model equations are in good agreement with the analytical results. In this case, similar to the classical epidemic models, the infective fronts propagate at a constant speed $c=5\times {10}^2$. It should be noted that since in the first case $(\mu =5\times {10}^{-6})$, the epidemic fronts propagate over the entire domain during a shorter time compared to the other cases $(\mu ={10}^{-4}$ and $\mu ={10}^{-3})$, for the first case, we consider the simulation duration equals 30. However, for the others, it is equal to 90.

Figure 7 shows the time evolution of the susceptible and infective waves propagating to the left side. Figure 7a shows the densities of the susceptibles and infectives at different times. Figure 7b,c show that the tails of the solutions decay exponentially as ∼$y=e^{\nu x}$. Here, the simulation duration is equal to 90. In this case, for all values of the truncation parameter, the fronts move at a constant velocity that is equal to $c=5\times {10}^2$.

It should be noted that when we consider both the left-and right-sided truncated fractional-order diffusion operators by choosing the strict values of the skewness parameter $\beta$, i.e., $-1<\beta <1$, ($\beta =0$ leads to a symmetric operator), both the left- and right-sided derivatives affect the epidemic speed and the asymptotic behaviour of the front’s tail. In that case, both the left- and right-moving front tails exhibit the same asymptotic behaviour based on the values of the truncation parameter. For instance, if the truncation parameter $\mu$ is set to ${10}^{-4}$, the tails of both the left- and right-moving fronts decay as ∼$e^{-\mu x}{x}^{-(\alpha +1)}$.

Figure 8a shows the time evolution of the Lagrangian trajectory $x_L\left(t\right)$ for the constant value of the fractional-order derivative $\alpha =1.2$ and different values of the truncation parameter $\mu$. As we see, in all cases, the values of $x_L\left(t\right)$ initially increase linearly at a constant speed on average equals $5\times {10}^2$ and then, for $\mu =5\times {10}^{-6}$, $x_L\left(t\right)$ increases exponentially leading to an exponential and unbounded velocity of the infective fronts (see Figure 2b), but for $\mu ={10}^{-4}$, $x_L\left(t\right)$ increases gradually leading to the convergence of the epidemic speed to the maximum velocity value $c_{max}\approx 5\times {10}^3$ (see Figure 4b), and finally for $\mu ={10}^{-3}$, $x_L\left(t\right)$ increases linearly leading to the constant speed $c=5\times {10}^2$ (see Figure 6b).

Figure 8b shows the time evolution of the Lagrangian trajectory $x_L\left(t\right)$ for different values of the fractional-order derivative $\alpha$ and a constant value of the truncation parameter $\mu ={10}^{-4}$. Since, the value of $\mu$ is constant, in all cases the Lagrangian velocity converges towards a plateau showing the maximum average speed, but the time to reach that plateau is not the same for different values of $\alpha$. As we see, the smaller the value of $\alpha$, i.e., the closer it is to 1, the faster the Lagrangian trajectory accelerates leading to the faster convergence of the front speed towards the maximum speed compared to values of $\alpha$ that are closer to 2.

5. Conclusions

In this work, we have investigated the spatial propagation of the epidemics caused by infectious diseases. We have considered a simple version of the epidemic models consisting of susceptible and infective populations. The diffusion process resulting from the random mobility of the infective individuals is taken into account in our model. Based on the previous studies, in the case of the Brownian motion, the infective waves propagate into the susceptible population at a constant speed and the tails of the travelling-wave solutions exhibit an exponential decay, while in the case of pure (untruncated) Lévy flights, the left-sided fractional derivative leads to an exponential and unbounded speed of the infective waves moving to the right side, but generates a constant speed at the left side. In that case, right-moving fronts have an algebraic tail. However, left-moving fronts exhibit an exponential tail.

As a new study, we have applied the truncated Lévy flight to the proposed model. Similar to the results obtained for the $\lambda$-truncated fractional-order Fisher–Kolmogorov Equation [46], for the right-moving fronts, we have shown that the epidemic speed is dependent on the level of the truncation parameter. In this case, we have considered a left-sided truncated fractional-order diffusion operator. For small values of $\lambda$ ($\lambda \gtrsim 0$), the tails of the infective waves can decay algebraically leading to an exponential growth of the epidemics. In that case, the truncation has no impact on the superdiffusive epidemics. By increasing the value of $\lambda$, the algebraic decaying tails can be tamed leading to either an upper bound on the epidemic speed representing the maximum speed or the generation of the infective waves of a constant shape propagating at a minimum constant speed similar to the classical models (second-order diffusion epidemic models). For the left-moving fronts, our numerical results show that for different values of the truncation parameter, the truncated left-sided fractional derivative generates the waves of a constant shape moving at a constant speed. Obviously, considering the strict values of the skewness parameter $\beta$, i.e., $-1<\beta <1$, leads to the effect of both the left-and right-sided truncated fractional-order diffusion operators on the epidemic speed. In that case, based on the value of the truncation parameter, the tails of both the left- and right-moving fronts have the same asymptotic behaviour.

From an applied perspective, the untruncated power laws (pure Lévy flights) have some drawbacks. As discussed in [57], in nature, due to the finite space, landscape and physiological limitations, the occurrence of arbitrarily large displacements by the individuals following a pure Lévy flight is not realistic. In that case, one individual cannot make displacements beyond an upper bound. By introducing truncated power-laws, such an upper bound can thus be taken into account in the fractional-order diffusion models. As mentioned by Hanert et al. [40], since modern epidemics caused by human infectious diseases can propagate over the entire globe very quickly, the untruncated fractional-order diffusion epidemic models can better represent such epidemics. However, concerning the infectious diseases of animals doing a Lévy flight, the propagation of the epidemics is more smooth, our findings thus suggest that truncated fractional-order diffusion models are more appropriate for modelling the animal diseases.

The space fractional-order epidemic models could be further improved by considering the non-Markovian diffusion processes in which the pdf for waiting times between displacements have a power-law asymptotic behaviour as ∼${\tau }^{-(\gamma +1)}$, where $0<\gamma <1$. In that case, the time derivative of order one is replaced by a time fractional derivative of order $\gamma$ [18]. As an example of such diffusion processes, Brockmann et al. [26] observed that the dispersion of bank notes has memory effects. In order to numerically solve time-space fractional diffusion equations, Hanert [52] has proposed an efficient and flexible scheme.