Dynamical Analysis of an SIS Epidemic Integrodifference Model

Miaomiao Zhang, Hui Cao^*

School of Mathematics & Data Science, Shaanxi University of Science & Technology, Xi’an, China.

DOI: 10.4236/jamp.2023.1111220                

Abstract

In order to analyze the impact of dispersal on disease transmission, we establish an SIS epidemic integrodifference model with a nonlinear incidence function. Firstly, the discrete-time SIS epidemic model is established and studied, including the existence and stability of equilibria, the existence of a flip bifurcation, and chaos. Secondly, the SIS epidemic integrodifference model is built based on the discrete-time SIS epidemic model with dispersal. The dynamic analysis of the model includes the existence and stability of equilibria, the existence of a traveling wave solution, and a minus-one bifurcation. Finally, the results suggest that dispersal causes the system to become more unstable and accelerates the spread of the disease when the equilibrium is unstable. Numerical examples are provided to demonstrate the theoretical results.

Keywords

Integrodifference Model, Flip Bifurcation, Travelling Wave Solution, Minus-One Bifurcation

1. Introduction

Discrete-time epidemic models have attracted the attention of many scholars [1] - [13] based on the following facts: 1) Since epidemiological data is often measured in discrete units of time, for example, weekly, or monthly, it is natural to use discrete-time models to describe the spread of epidemics within a population; 2) Discrete-time model exhibit rich dynamic behaviors compared to the corresponding continuous models [14] [15] [16] [17] . As with continuous-time models, a discrete-time model can describe dispersal using continuous space. It has led to the emergence of epidemic models using the integrodifference equations.

The integrodifference epidemic model considers the influence of time and spatial location on disease transmission. We add diffusion to the discrete-time model to obtain an integrodifference model, analyze the dynamic behavior of the integrodifference model, and then obtain the transmission law of infectious diseases. Finally, measures to control the spread of infectious diseases are obtained.

Therefore, the research results of integrodifference epidemic model gradually increase [18] [19] [20] [21] [22] , including the existence of a traveling wave solution, minimum wave speed, bifurcation, and so on.

At present, relatively little research has been done on the SIS epidemic integrodifference model, not considering that the incidence rate is nonlinear. At the same time, the classical discrete infectious disease model does not consider the influence of individual diffusion on disease transmission, but in real life, the spatial movement of individuals is inevitable, and the dynamic behavior of spatial models also deserves further study.

In this paper, we gain further insights into the impact of dispersal on the spread of the disease by analyzing the dynamic behavior of the epidemic integrodifference model with nonlinear incidence. Our interest is on the dynamic behavior of the discrete-time epidemic model with and without dispersal. To this end, we first establish the discrete-time SIS epidemic model based on the difference equations.

Let $N_t$ be the total number of the population at time t, which is divided into two compartments, namely the susceptible and the infectious compartment, respectively. We use $S_t$ and $I_t$ to represent the number of the individuals in the susceptible compartment and infectious compartment at time t, respectively, $\Lambda >0$ to represent the constant recruitment of the population, $0<p<1$ to represent the probability of survival of an individual, and $0<\gamma <1$ to represent the recovery probability of an infectious individual.

Based on the idea in [5] , a discrete-time SIS epidemic model with a nonlinear incidence function is defined as follows:

$\{\begin{array}{l}{S}_{t+1}=\Lambda +p{S}_t\left(1-\beta {\text{e}}^{-m{I}_t}{I}_t\right)+p\gamma I_t,\hfill \\ I_{t+1}=p\beta {\text{e}}^{-m{I}_t}{S}_t{I}_t+p\left(1-\gamma \right)I_t,\hfill \end{array}$ (1)

where $\beta {\text{e}}^{-m{I}_t}{S}_t{I}_t$ is the number of newly infected individuals in the time interval $\left[t\mathrm{,}t+1\right)$ , $0<\beta <1$ is the probability of infection, $m>0$ weighs the infectious capacity of the infectious individuals, and ${\text{e}}^{-m{I}_t}$ is a control factor, which indicates that the implementation of certain effective control measures can reduce the spread of the disease. The impact of control measures is mainly based on the number of infectious individuals: Stricter controls will be applied to more infected individuals, which will reduce new infections. The nonlinear incidence function $\beta {\text{e}}^{-m{I}_t}{S}_t{I}_t$ has been discussed in many papers to characterize the decrease in infection rate caused by the control measures [6] [10] [23] [24] , and satisfies that $\beta {\text{e}}^{-m{I}_t}{S}_t{I}_t\le S_t$ , which represents the number of susceptible individuals who are infected cannot exceed the total number of susceptible individuals during per unit of time.

Since $\beta {\text{e}}^{-m{I}_t}{S}_t{I}_t\le S_t$ , we have that $1-\beta {\text{e}}^{-m{I}_t}{I}_t>0$ . Therefore, the non negativity of the solution of model (1) can be demonstrated.

By using of $N_t=S_t+I_t$ , it can be obtained from model (1)

$N_{t+1}=\Lambda +p{N}_t,$

which implies that

$N_t=\frac{\Lambda }{1-p}+p^t\left(N_0-\frac{\Lambda }{1-p}\right)\mathrm{.}$

It is easy to obtain ${\lim }_{t\to \infty }{N}_t=\frac{\Lambda }{1-p}\triangleq N^{*}$ since $0<p<1$ .

Let $X_t=\frac{S_t}{N^{*}}$ , $Y_t=\frac{I_t}{N^{*}}$ , then ${\lim }_{t\to \infty }\left(X_t+Y_t\right)=1$ , and model (1) can be written as

$\{\begin{array}{l}{X}_{t+1}=1-p+p{X}_t-a{\text{e}}^{-b{Y}_t}{X}_t{Y}_t+\left(p-c\right)Y_t,\hfill \\ Y_{t+1}=a{\text{e}}^{-b{Y}_t}{X}_t{Y}_t+c{Y}_t,\hfill \end{array}$ (2)

where $a=p\beta N^{\mathrm{<mo>\; *\; </mo>}}$ , $b=m{N}^{\mathrm{<mo>\; *\; </mo>}}$ , and $c=p\left(1-\gamma \right)$ . It is clear that the dynamic behaviors of model (1) and (2) are equivalent.

Let $X_t\left(x\right)$ and $Y_t\left(x\right)$ be the proportion of susceptible and infectious individuals in the population at time t and at position x, respectively, and $K\left(x\mathrm{,}y\right)$ denote the dispersal redistribution kernel. If an individual is located at position y before diffusion, then $K\left(x\mathrm{,}y\right)$ gives the probability that the individual arrives at position x after diffusion. Suppose that,

1) $K\left(x\mathrm{,}y\right)$ is related to the distance of spread, that is, $K\left(x\mathrm{,}y\right)=K\left(x-y\right)$ ;

2) there is no death in diffusion, ${\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\text{d}y=1$ , which means that growth and dispersal phase of population are independent of each other;

3) $K\left(x\mathrm{,}y\right)$ is exponentially bounded, continuous and non-negative.

Based on model (2) and the above assumptions, we constructed an integrodifference SIS epidemic model:

$\{\begin{array}{l}{X}_{t+1}\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left[\frac{\Lambda }{N^{*}}+p{X}_t\left(y\right)-a{\text{e}}^{-b{Y}_t\left(y\right)}{X}_t\left(y\right)Y_t\left(y\right)+\left(p-c\right)Y_t\left(y\right)\right]\text{d}y,\hfill \\ Y_{t+1}\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left[a{\text{e}}^{-b{Y}_t\left(y\right)}{X}_t\left(y\right)Y_t\left(y\right)+c{Y}_t\left(y\right)\right]\text{d}y.\hfill \end{array}$ (3)

In the remaining parts, we will focus on the dynamic behavior of model (2) and model (3). The structure of this paper as follows: In Section 2, we discuss the dynamic behavior of model (2), which include the basic reproductive number, the existence and the stability of equilibria, and the existence of a Flip bifurcation. In Section 3, we study the dynamic behavior of model (3), which include the basic reproductive number, the existence and the stability of equilibria, the existence of a traveling wave solution, and the existence of a minus-one bifurcation. In Section 4, we give concluding observations and discussions.

2. The Discrete-Time Epidemic Model

Here, we analyze the dynamic behavior of the discrete-time model (2), including the basic reproduction number ${\mathcal{R}}_0$ , the existence and the stability of equilibria, and the existence of a flip bifurcation.

Following the next generation matrix approach in [25] , the basic reproduction number ${\mathcal{R}}_0$ of model (2) can be defined as

${\mathcal{R}}_0=\frac{a}{1-c}\mathrm{,}$

where $\frac{1}{1-c}=\frac{1}{1-p\left(1-\gamma \right)}$ represents the average infection period, and ${\mathcal{R}}_0=\frac{a}{1-c}=p\times \beta \times \frac{1}{1-p\left(1-\gamma \right)}\times N^{*}$ represents the average number of new cases produced by an infected member throughout the period of infection.

Since ${\lim }_{t\to \infty }\left(X_t+Y_t\right)=1$ , substituting $X_t=1-Y_t$ into the second equation of model (2) leads to the following limiting system

$Y_{t+1}=a{\text{e}}^{-b{Y}_t}\left(1-Y_t\right)Y_t+c{Y}_t.$ (4)

The equivalence of a dynamical system and its limiting system has been established by many scholars [26] [27] . Although there is no rigorous theoretical conclusions for discrete dynamical system, numerical simulation implies that the dynamic behaviors of model (2) and (4) are equivalent under some assumptions [28] . In fact, there are many scholars have used this method to study the dynamic behavior of discrete-time dynamic models [5] [10] [27] .

By direct calculation, it can be concluded that model (4) always has zero equilibrium $Y=0$ . In order to discuss the existence of the positive equilibrium $Y=Y_{*}$ , we denote

$f\left(Y\right)=a{\text{e}}^{-bY}\left(1-Y\right),Y\in \left[0,1\right].$

Obviously, $f\left(Y\right)\ge 0$ for any $Y\in \left[\mathrm{0,1}\right]$ , $f\left(0\right)=a$ , and $f\left(1\right)=0$ . Because

$f^{\prime }\left(Y\right)=-a{\text{e}}^{-bY}\left(b-bY+1\right)\mathrm{,}Y\in \left[\mathrm{0,1}\right]\mathrm{,}$

and $b-bY+1>0$ , we can get $f^{\prime }\left(Y\right)<0$ when $Y\in \left[\mathrm{0,1}\right]$ . This implies that $f\left(Y\right)$ is a monotonically decreasing function of Y in $Y\in \left[\mathrm{0,1}\right]$ . Therefore, when $0<1-c<a$ , that is, ${\mathcal{R}}_0>1$ , there is a unique $0<Y_{*}<1$, that satisfies $f\left(Y_{*}\right)=1-c$ . Namely, model (4) has a unique positive equilibrium $Y_{\mathrm{<mo>\; *\; </mo>}}$ in $\left(\mathrm{0,1}\right)$ .

Theorem 1. Model (2) always has the disease free equilibrium $E_0=\left(\mathrm{1,0}\right)$ , and when ${\mathcal{R}}_0>1$ , model (2) also has a unique endemic equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}\left(X_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}_{\mathrm{<mo>\; *\; </mo>}}\right)$ , where $X_{*}=1-Y_{*}$ , $0<Y_{*}<1$ , and satisfies $a{\text{e}}^{-b{Y}_{\mathrm{<mo>\; *\; </mo>}}}\left(1-Y_{\mathrm{<mo>\; *\; </mo>}}\right)=1-c$ .

2.1. The Stability of the Equilibria

Theorem 2. If ${\mathcal{R}}_0<1$ , then the disease free equilibrium E₀ of model (2) is globally asymptotically stable; while if ${\mathcal{R}}_0>1$ , $E_0$ is unstable.

Proof. Since the Jacobian matrix of model (2) at E₀ is

$J_0=\left(\begin{array}{cc}p& p-c-a\\ 0& c+a\end{array}\right).$

We have $\left|\lambda E-J_0\right|=\left(\lambda -p\right)\left(\lambda -c-a\right)$ , where E is the identity matrix. It is easy to see that the characteristic roots are ${\lambda }_1=p$ , and ${\lambda }_2=c+a$ , respectively. Obviously, $\left|{\lambda }_1\right|<1$ . When ${\mathcal{R}}_0<1$ , we have $0<c+a<1$, that is ${\lambda }_2<1$ , which implies that $E_0$ is locally asymptotically stable. When ${\mathcal{R}}_0>1$ , we have ${\lambda }_2=c+a>1$ , namely, $E_0$ is unstable. To prove the global stability of $E_0$ , we consider the linear difference equation

${\hat{Y}}_{t+1}=\left(a+c\right){\hat{Y}}_t\mathrm{.}$ (5)

It is clear that Equation (5) has unique zero equilibrium $\hat{Y}=0$ , and ${\lim }_{t\to +\infty }{\hat{Y}}_t=0$ when $0<a+c<1$ . Then, by model (4), we obtain

$\begin{array}{c}{Y}_{t+1}=a{\text{e}}^{-b{Y}_t}\left(1-Y_t\right)Y_t+c{Y}_t\\ \le a\left(1-Y_t\right)Y_t+c{Y}_t\\ \le a{Y}_t+c{Y}_t\\ =\left(a+c\right)Y_t.\end{array}$

The comparison principle means that $0\le Y_t\le {\hat{Y}}_t$ when the difference equation and Equation (5) have the same initial condition. Namely, ${\lim }_{t\to +\infty }{Y}_t\le {\lim }_{t\to +\infty }{\hat{Y}}_t=0$ when ${\mathcal{R}}_0<1$ . Therefore, when ${\mathcal{R}}_0<1$ , $E_0$ is globally asymptotically stable.

In order to analyze the stability of the endemic equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ , we give the Jacobian matrix of the model (2) at $E_{\mathrm{<mo>\; *\; </mo>}}$ :

$J_{\mathrm{<mo>\; *\; </mo>}}=\left(\begin{array}{cc}p-\frac{Y_{\mathrm{<mo>\; *\; </mo>}}\left(1-c\right)}{1-Y_{\mathrm{<mo>\; *\; </mo>}}}& p-1+b{Y}_{\mathrm{<mo>\; *\; </mo>}}\left(1-c\right)\\ \frac{Y_{\mathrm{<mo>\; *\; </mo>}}\left(1-c\right)}{1-Y_{\mathrm{<mo>\; *\; </mo>}}}& 1-b{Y}_{\mathrm{<mo>\; *\; </mo>}}\left(1-c\right)\end{array}\right).$ (6)

Then, the characteristic equation is

$F_{\mathrm{<mo>\; *\; </mo>}}\left(\lambda \right)={\lambda }^2-A\lambda +p\left(A-p\right)\mathrm{,}$

where $A=p-\frac{\left(1-c\right)Y_{\mathrm{<mo>\; *\; </mo>}}}{1-Y_{\mathrm{<mo>\; *\; </mo>}}}+1-b{Y}_{\mathrm{<mo>\; *\; </mo>}}\left(1-c\right)$ . Denote

$\begin{array}{l}{g}_1\left(Y_{*}\right)=\left(1+p\right)\left[b\left(1-c\right)Y_{*}^2-\left(2+\left(1+b\right)\left(1-c\right)\right)Y_{*}+2\right],\\ g_2\left(Y_{*}\right)=\left(1-p\right)\left(1-c\right)\left[-b{Y}_{*}^2+\left(1+b\right)Y_{*}\right].\end{array}$

We have

$F_{*}\left(-1\right)=\frac{1}{1-Y_{*}}{g}_1\left(Y_{*}\right),F_{*}\left(1\right)=\frac{1}{1-Y_{*}}{g}_2\left(Y_{*}\right).$

Since $g_1\left(0\right)=2\left(1+p\right)>0$ , $g_1\left(1\right)=\left(1+p\right)\left(c-1\right)<0$ , the continuity of function $g_1$ about $Y_{\mathrm{<mo>\; *\; </mo>}}$ implies there exists $Y_1\in \left(\mathrm{0,1}\right)$ , such that $g_1\left(Y_1\right)=0$ . Therefore, we have $g_1\left(Y_{*}\right)>0$ for $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ , and $g_1\left(Y_{*}\right)<0$ for $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(Y_1\mathrm{,1}\right)$ . In addition, both $g_2\left(0\right)=0$ and $g_2\left(1\right)=\left(1-p\right)\left(1-c\right)>0$ means that $g_2\left(Y_{*}\right)>0$ for $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,1}\right)$ .

Theorem 3. Suppose ${\mathcal{R}}_0>1$ .

1) $E_{\mathrm{<mo>\; *\; </mo>}}$ is global asymptotic stability if $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ ;

2) $E_{\mathrm{<mo>\; *\; </mo>}}$ is non-hyperbolic if $Y_{*}=Y_1$ ;

3) $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable if $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(Y_1\mathrm{,1}\right)$ .

Proof. (i) When $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ , we have $g_1\left(Y_{*}\right)>0$ , and $g_2\left(Y_{*}\right)>0$ . It implies that $F_{*}\left(1\right)>0$ , and $F_{*}\left(-1\right)>0$ . In addition, $0<p<1$ and $0<c<1$ implies that $1-p\left(A-p\right)>0$ . Therefore, the modulus of both roots of $F_{\mathrm{<mo>\; *\; </mo>}}\left(\lambda \right)=0$ is less than 1. It implies that $E_{\mathrm{<mo>\; *\; </mo>}}$ is locally asymptotically stable when $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ .

Let us prove the global stability of $E_{\mathrm{<mo>\; *\; </mo>}}$ as following. Denote

$g\left(Y\right)=a{\text{e}}^{-bY}\left(1-Y\right)Y+cY,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}Y\in \left[\mathrm{0,1}\right]\mathrm{.}$

By direct calculation, we obtain

$\begin{array}{c}1+g^{\prime }\left(Y_{*}\right)=1+c+a{\text{e}}^{-b{Y}_{*}}\left[b{Y}_{*}^2-\left(b+2\right)Y_{*}+1\right]\\ =1+c+\frac{1-c}{1-Y_{*}}\left[b{Y}_{*}^2-\left(b+2\right)Y_{*}+1\right]\\ =\frac{1}{1-Y_{*}}\left[b\left(1-c\right)Y_{*}^2-\left(2+\left(1+b\right)\left(1-c\right)\right)Y_{*}+2\right]\\ =\frac{1}{\left(1-Y_{*}\right)\left(1+p\right)}{g}_1\left(Y_{*}\right)\\ >0.\end{array}$

It implies that $1+g^{\prime }\left(Y_{\mathrm{<mo>\; *\; </mo>}}\right)\ne 0$ for $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ . Suppose g has a nontrivial 2-period solution $\left\{q_1\mathrm{,}{q}_2\right\}$ , $q_1\mathrm{,}{q}_2\in \left(\mathrm{0,}{Y}_1\right)$ , and satisfies $g\left(q_1\right)=q_2$ and $g\left(q_2\right)=q_1$ , where $q_1\ne q_2$ . Therefore, it is known from the mean theorem that

there exists $l\in \left(q_1\mathrm{,}{q}_2\right)$ , such that $g^{\prime }\left(l\right)=\frac{g\left(q_1\right)-g\left(q_2\right)}{q_1-q_2}=-1$ , that is, $1+g^{\prime }\left(l\right)=0$ , which contradicts the above conclusion. Therefore, $g\left(Y_{\mathrm{<mo>\; *\; </mo>}}\right)$ has no non-trivial 2-period solution on $\left(\mathrm{0,}{Y}_1\right)$ , namely, $E_{\mathrm{<mo>\; *\; </mo>}}$ is global asymptotically stable when $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ .

(ii) When $Y_{*}=Y_1$ , we have $g_1\left(Y_{\mathrm{<mo>\; *\; </mo>}}\right)=0$ , which leads to $F_{\mathrm{<mo>\; *\; </mo>}}\left(-1\right)=0$ . Because of $p\left(A-p\right)=p\left(1-\frac{\left(1-c\right)Y_{\mathrm{<mo>\; *\; </mo>}}}{1-Y_{\mathrm{<mo>\; *\; </mo>}}}-b{Y}_{\mathrm{<mo>\; *\; </mo>}}\left(1-c\right)\right)$ , it is easy to get that $p\left(A-p\right)\ne 1$ , thus $-A\ne 2$ . And then, the characteristic equation $F_{\mathrm{<mo>\; *\; </mo>}}\left(\lambda \right)=0$ of model (2) at $E_{\mathrm{<mo>\; *\; </mo>}}$ has a eigenvalue of -1, i.e. $E_{\mathrm{<mo>\; *\; </mo>}}$ is non-hyperbolic.

(iii) When $Y_{\mathrm{<mo>\; *\; </mo>}}\in \left(Y_1\mathrm{,1}\right)$ , we have $g_1\left(Y_{*}\right)<0$ , namely, $F_{*}\left(-1\right)<0$ . Based on the Jury criterion, $F_{\mathrm{<mo>\; *\; </mo>}}\left(\lambda \right)=0$ at least has a characteristic root, whose modulus is greater than 1. Therefore, $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable.

2.2. Flip Bifurcation

From the stability conclusion of Theorem 3, it can be seen that $F_{\mathrm{<mo>\; *\; </mo>}}\left(\lambda \right)=0$ has a characteristic root-1 when $Y_{*}=Y_1$ , that is, the equilibrium is non-hyperbolic. Therefore, when the parameter passes the critical point, a flip bifurcation may generate near $E_{\mathrm{<mo>\; *\; </mo>}}$ , and the existence of this bifurcation will be illustrated by the following example.

There are difficulties with the calculation process of the flip bifurcation, so we take $a=2\left(1-c\right){\text{e}}^{\frac{b}{2}}$ to simplify it. In this case, we can get $E_{*}=E_{\frac{1}{2}}=\left(\frac{1}{2},\frac{1}{2}\right)$ , then

$J_{\frac{1}{2}}=\left(\begin{array}{cc}p-1+c& p-c+\frac{b}{2}\left(1-c\right)\\ 1-c& 1-\frac{b}{2}\left(1-c\right)\end{array}\right)\mathrm{,}$

and

$F_{\frac{1}{2}}\left(\lambda \right)={\lambda }^2-\left[p+c-\frac{b}{2}\left(1-c\right)\right]\lambda +p\left[c-\frac{b}{2}\left(1-c\right)\right]\mathrm{.}$

When $b=\frac{2\left(1+c\right)}{1-c}$ , and $p+c-\frac{b}{2}\left(1-c\right)\ne -2$ , the eigenvalues of $J_{\frac{1}{2}}$ are ${\lambda }_1=-1$ , ${\lambda }_2=p<1$ , respectively. The following will prove the existence of the Flip bifurcation.

Theorem 4. When $b=\frac{2\left(1+c\right)}{1-c}$ , and $c\in \left(\mathrm{0.42024,1}\right)$ , the model (4) produces a flip bifurcation at the endemic equilibrium $Y_{*}=\frac{1}{2}$ . In this condition, the 2-period solution bifurcating from the endemic equilibrium $Y_{\mathrm{<mo>\; *\; </mo>}}$ is stable.

Proof. Let $\bar{Y}=Y-\frac{1}{2}$ , $\mu =b-\frac{2\left(1+c\right)}{1-c}$ , where $\mu$ is the bifurcation parameter, the model (4) can be rewritten as

${\bar{Y}}_{t+1}=a{\text{e}}^{-\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\left({\bar{Y}}_t+\frac{1}{2}\right)}\left(\frac{1}{2}-{\bar{Y}}_t\right)\left({\bar{Y}}_t+\frac{1}{2}\right)+c\left({\bar{Y}}_t+\frac{1}{2}\right)-\frac{1}{2}\mathrm{.}$

Substituting $a=2\left(1-c\right){\text{e}}^{\frac{b}{2}}$ into the above model, it can be obtained that

${\bar{Y}}_{t+1}=2\left(1-c\right){\text{e}}^{-\left(\mu +\frac{2\left(1+c\right)}{1-c}\right){\bar{Y}}_t}\left(\frac{1}{4}-{\bar{Y}}_t^2\right)+c{\bar{Y}}_t+\frac{1}{2}\left(c-1\right)\mathrm{.}$ (7)

Let $G\left(\bar{Y}\mathrm{,}\mu \right)=2\left(1-c\right){\text{e}}^{-\left(\mu +\frac{2\left(1+c\right)}{1-c}\right){\bar{Y}}_t}\left(\frac{1}{4}-{\bar{Y}}_t^2\right)+c{\bar{Y}}_t+\frac{1}{2}\left(c-1\right)$ . Because $G\left(\mathrm{0,}\mu \right)=0$ , and

$\frac{\partial G\left(\bar{Y},\mu \right)}{\partial \bar{Y}}=2\left(1-c\right){\text{e}}^{-\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\bar{Y}}\left[\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\left({\bar{Y}}^2-\frac{1}{4}\right)-2\bar{Y}\right]+c,$

$\frac{\partial G^2\left(\bar{Y},\mu \right)}{\partial \bar{Y}\partial \mu }=2\left(1-c\right){\text{e}}^{-\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\bar{Y}}\left[\left({\bar{Y}}^2-\frac{1}{4}\right)\left(1-\bar{Y}\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\right)+2{\bar{Y}}^2\right],$

$\begin{array}{c}\frac{\partial G^3\left(\bar{Y},\mu \right)}{\partial {\bar{Y}}^3}=2\left(1-c\right)\left(\mu +\frac{2\left(1+c\right)}{1-c}\right){\text{e}}^{-\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\bar{Y}}[{\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)}^2\left({\bar{Y}}^2-\frac{1}{4}\right)\\ -6\left(\mu +\frac{2\left(1+c\right)}{1-c}\right)\bar{Y}+6].\end{array}$

We have

$\begin{array}{l}\frac{\partial G\left(0,0\right)}{\partial \bar{Y}}=-1,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{\partial G^2\left(0,0\right)}{\partial \bar{Y}\partial \mu }=-\frac{1-c}{2}<0,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\\ \frac{\partial G^3\left(0,0\right)}{\partial {\bar{Y}}^3}=\frac{4\left(1+c\right)\left(5c^2-14c+5\right)}{{\left(1-c\right)}^2}.\end{array}$

Since $5c^2-14c+5<0$ when $c\in \left(\mathrm{0.42024,1}\right)$ (see Figure 1(a)), we can get $\frac{\partial G^3\left(0,0\right)}{\partial {\bar{Y}}^3}<0$ . Hence, the model (4) will experience a stable flip bifurcation at the endemic equilibrium $Y_{*}=\frac{1}{2}$ .

We can also use the numerical simulations to verify the theoretical results. Taking $a=3241.23357$ , $b=18$ , and $c=0.8$ , a flip bifurcation for model (2) at $E_{\mathrm{<mo>\; *\; </mo>}}=\left(\frac{1}{2}\mathrm{,}\frac{1}{2}\right)$ is given in Figure 1(b). The numerical simulation results show that when $Y_{*}=Y_1$ , there is a flip bifurcation and the two period solution is stable.

In fact, let $b=9$ , $c=0.63$ , and $p=0.4$ . We have ${\mathcal{R}}_0>1$ , and $Y_{*}>Y_1$ , the dynamic behavior of model (2) is shown in Figure 2. When $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable, Figure 2(a) and Figure 2(b) displays that the model (4) may undergo the 4-period solution and 8-period solution, respectively.

We let $b=0.0055$ , $c=0.00054$ , and $p=0.27$ , the bifurcation diagram of model (2) show that in Figure 3. Figure 3 displays that $E_0$ is global asymptotic stability when ${\mathcal{R}}_0<1$ . With the increase of ${\mathcal{R}}_0$ , the equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ is gradually unstable. We can see that $E_{\mathrm{<mo>\; *\; </mo>}}$ is global asymptotic stability when $1<{\mathcal{R}}_0<3$ . When ${\mathcal{R}}_0=3$ , the model (2) produces Flip bifurcation at $E_{\mathrm{<mo>\; *\; </mo>}}$ , and the 2-period

Figure 1. Model (2) experience a Flip bifurcation at $E_{\mathrm{<mo>\; *\; </mo>}}=\left(\frac{1}{2}\mathrm{,}\frac{1}{2}\right)$ .

Figure 2. Model (2) may experience periodic oscillations when $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable.

Figure 3. Bifurcation graph of model (2) about ${\mathcal{R}}_0$ .

solution is stable. The model (2) experience 4-period solution when $3.457<{\mathcal{R}}_0<3.55$ . Finally, chaos will occur as ${\mathcal{R}}_0$ increases.

3. The Integrodifference Epidemic Model

Here, we will investigate the impact of the dispersal on the spread of diseases by analyzing the dynamic behavior of the epidemic integrodifference model (3). To this end, based on (4), we give an equivalent system of the integrodifference model (3) as follows:

$Y_{t+1}\left(x\right)=F\left[Y_t\right]\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y_t\left(y\right)\right)\text{d}y\mathrm{,}$ (8)

where $g\left(Y_t\right)=a{\text{e}}^{-b{Y}_t}\left(1-Y_t\right)Y_t+c{Y}_t$ .

In the following, we will use (8) to discuss the dynamic behavior of model (3).

3.1. The Basic Reproduction Number

Obviously, the basic reproduction number of the epidemic integrodifference model (8) is related to $g^{\prime }\left(0\right)$ , and $g^{\prime }\left(0\right)=a+c$ . Let $Y_t\left(y\right)=Y^0+{\tilde{Y}}_t\left(y\right)$ , and substituting it into the model (8), we get

$\begin{array}{c}{Y}^0+{\tilde{Y}}_{t+1}\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y^0+{\tilde{Y}}_t\left(y\right)\right)\text{d}y\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)a{\text{e}}^{-b\left(Y^0+{\tilde{Y}}_t\left(y\right)\right)}\left(1-Y^0-{\tilde{Y}}_t\left(y\right)\right)\left(Y^0+{\tilde{Y}}_t\left(y\right)\right)\\ +c\left(Y^0+{\tilde{Y}}_t\left(y\right)\right)\text{d}y.\end{array}$ (9)

Furthermore, the Taylor series expansion of (9) at $Y^0=0$ is

${\tilde{Y}}_{t+1}\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)B{\tilde{Y}}_t\left(y\right)\text{d}y,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}B=a+c.$

Next, we explore the stability of the equilibrium $Y^0=0$ based on the perturbation system ${\tilde{Y}}_t\left(x\right)={\lambda }^{t}z\left(x\right)$ , where $\lambda \ne 0$ , $z\left(x\right)>0$ , ${\displaystyle {\int }_{-\infty }^{\infty }}z\left(x\right)\text{d}x<\infty$ .

For the system ${\tilde{Y}}_t\left(x\right)={\lambda }^{t}z\left(x\right)>0$ , if $\lambda >1$ , the zero equilibrium is unstable, while if $0<\lambda <1$ , the zero equilibrium is locally asymptotically stable. Because of

$\begin{array}{l}{\lambda }^{t+1}z\left(x\right)=B{\displaystyle {\int }_{-\infty }^{\infty }}{\lambda }^{t}z\left(y\right)K\left(x-y\right)\text{d}y,\\ \lambda z\left(x\right)=B{\displaystyle {\int }_{-\infty }^{\infty }}z\left(y\right)K\left(x-y\right)\text{d}y,\end{array}$

We have

$\begin{array}{c}\lambda {\displaystyle {\int }_{-\infty }^{\infty }}z\left(x\right)\text{d}x=B{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}z\left(y\right)K\left(x-y\right)\text{d}y\text{d}x\\ =B\left({\displaystyle {\int }_{-\infty }^{\infty }}z\left(y\right)\text{d}y\right)\left({\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\text{d}x\right).\end{array}$

The assumption 2) ensures $\lambda =B>0$ , which implies that the basic reproduction number ${\bar{\mathcal{R}}}_0$ for the integrodifference model (8) can be defined as

${\bar{\mathcal{R}}}_0=\frac{a}{1-c}.$

It is clear that ${\bar{\mathcal{R}}}_0={\mathcal{R}}_0$ , that is, dispersal does not affect the basic reproduction number.

Based on the assumption 2), the constant solution of model (8) satisfies $Y_{t+1}=g\left(Y_t\right)$ . Theorem 1 implies that the following conclusions hold.

Theorem 5 The zero equilibrium $Y^0=0$ for the integrodifference model (8) always exists, and if ${\bar{\mathcal{R}}}_0>1$ , model (8) has the unique positive equilibrium $Y^{\mathrm{<mo>\; *\; </mo>}}$ , where $0<Y^{*}<1$ , and satisfies $a{\text{e}}^{-b{Y}^{\mathrm{<mo>\; *\; </mo>}}}\left(1-Y^{\mathrm{<mo>\; *\; </mo>}}\right)=1-c$ .

From Theorem 5, it can be seen that model (3) always has a disease free equilibrium $E^0=\left(X^0,Y^0\right)=\left(1,0\right)$ , and if ${\bar{\mathcal{R}}}_0>1$ , model (3) has the unique endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}=\left(X^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , where

$X^{*}=1-Y^{*},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}a{\text{e}}^{-b{Y}^{*}}\left(1-Y^{*}\right)=1-c,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>0}<Y^{*}<1.$

3.2. The Stability of the Equilibria

Theorem 6. When ${\bar{\mathcal{R}}}_0<1$ , the zero equilibrium $Y^0=0$ of the integrodifference model (8) is globally asymptotically stable.

Proof. Since the kernel function $K\left(x-y\right)$ is nonnegative, for a continuous space $C_{[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}]}$ over $ℝ$ , it is known by the Lipschitz condition, there exists constant L, such that

$\begin{array}{c}\left|Y_{t+1}\left(x\right)-0\right|=\left|{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y_t\left(y\right)\right)\text{d}y-{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(0\right)\text{d}y\right|\\ \le K\left(x-y\right){\displaystyle {\int }_{-\infty }^{\infty }}\left|g\left(Y_t\left(y\right)\right)-g\left(0\right)\right|\text{d}y\\ \le LK\left(x-y\right){\displaystyle {\int }_{-\infty }^{\infty }}\left|Y_t\left(y\right)-0\right|\text{d}y.\end{array}$

Theorem 2 implies that $\left|Y_t-0\right|\to 0$ as $t\to +\infty$ when ${\bar{\mathcal{R}}}_0<1$ , that is, ${\lim }_{t\to +\infty }\left|Y_{t+1}\left(x\right)-0\right|=0$ . Namely, $Y^0$ is globally asymptotically stable if ${\bar{\mathcal{R}}}_0<1$ .

Because of $Y^{\mathrm{<mo>\; *\; </mo>}}=Y_{\mathrm{<mo>\; *\; </mo>}}$ , Theorem 3 implies that $\left|Y_t-Y^{\mathrm{<mo>\; *\; </mo>}}\right|\to 0$ as $t\to +\infty$ when ${\bar{\mathcal{R}}}_0>1$ . Based on a similar proof of the Theorem6, we can conclude as follows:

Theorem 7. If ${\bar{\mathcal{R}}}_0>1$ , and $0<Y^{*}<Y_1$ , the positive equilibrium $Y^{\mathrm{<mo>\; *\; </mo>}}$ of (8) is globally asymptotically stable.

Theorem 6 and Theorem 7 implies that the disease free equilibrium $E^0$ of model (3) is globally asymptotically stable when ${\bar{\mathcal{R}}}_0<1$ ; whereas the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable when ${\bar{\mathcal{R}}}_0>1$ , and $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ .

Next, we will use the fast Fourier transform to numerically demonstrate the dynamic behaviors of model (3) when $E^{\mathrm{<mo>\; *\; </mo>}}$ is unstable. Taking $K\left(x-y\right)$ as a Gaussian kernel, namely, $K\left(x-y\right)=\frac{1}{\sigma \sqrt{2\pi }}{\text{e}}^{-\frac{{\left(x-y\right)}^2}{2{\sigma }^2}}$ , ${\sigma }^2=0.25$ , $b=9$ , $c=0.63$ , $p=0.4$ , and the spatial variable as $-40\le x\le 40$ , the dynamic behavior of model (3) is shown in Figure 4. Figure 4(a) and Figure 4(b) show that model (3) may experience period oscillations when $a=80$ and $a=233$ . Moreover, we see that the addition of diffusion increases the wavelength, which accelerates the spread of the disease.

Figure 4. Model (3) may experience periodic oscillation when $E^{\mathrm{<mo>\; *\; </mo>}}$ is unstable.

3.3. The Existence of a Traveling Wave Solution

If a solution of the integrodifference model (8) can be in the form of $Y_t\left(x\right)=w\left(\varphi \right)$ , where $\varphi =x-ct$ , then the solution $w\left(\varphi \right)$ for this model is called a traveling wave solution with a constant dispersal speed of c. Substituting it into (8) leads to the following system

$w\left(x-ct-c\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-ct-y\right)g\left(w\left(y-ct\right)\right)\text{d}\left(y-ct\right)\mathrm{.}$

Let $\eta =y-ct$ , then the above equation can be simplified to

$w\left(\varphi -c\right)=F\left[w\right]\left(\varphi \right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(\varphi -\eta \right)g\left(w\left(\eta \right)\right)\text{d}\eta \mathrm{.}$ (10)

In this case, the equilibria $Y=Y^{*}$ and $Y=Y^0$ are solutions to model (10) and satisfy the following boundary conditions

$w\left(-\infty \right)=Y^{*},w\left(\infty \right)=0.$ (11)

The asymptotic spreading speed of the integrodifference model (8) is the minimum wave speed of a nontrivial traveling wave solution. It is defined as

$c^{*}=\underset{s>0}{\min }\left\{\frac{1}{s}\ln \left[g^{\prime }\left(0\right)m\left(s\right)\right]\right\},$

where $m\left(s\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x\right){\text{e}}^{sx}\text{d}x$ .

According to the method in [29] , we explore the traveling wave solution of the integrodifference model (8) when ${\bar{\mathcal{R}}}_0>1$ .

Theorem 8. Consider that when $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ , the operator F on the space $C_{[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}]}$ of continuous functions has the following properties.

(i) For all $n\in ℝ$ , there is $F\left[Y\left(\cdot -n\right)\right]\left(x\right)=F\left[Y\right]\left(x-n\right)$ ;

(ii) If $Y\in \left[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right]$ , then $F\left[Y\right]\in \left[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right]$ ;

(iii) $F\left(0\right)=0$ , $F\left(Y^{\mathrm{<mo>\; *\; </mo>}}\right)=Y^{\mathrm{<mo>\; *\; </mo>}}$ , $F\left(n\right)>n$ for $n\in \left(\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ ;

(iv) When $0\le Y_1\le Y_2$ , $F\left(Y_1\right)\le F\left(Y_2\right)$ ;

(v) If $\left\{g_j\right\}\subset C_{[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}]}$ and $g_j$ converges uniformly to $g_1$ on a compact subset of $ℝ$ , then $F\left[g_j\right]$ pointwise converge $F\left[g_1\right]$ as $j\to \infty$ ;

(vi) Every sequence $\left\{g_j\right\}$ exists subsequence $\left\{g_{j_i}\right\}$ in $C_{[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}]}$ , and make $\left\{F\left[g_{j_i}\right]\right\}$ converges uniformly on every bounded subset of $ℝ$ .

If F in (8) holds for (i)-(v), then F exists an asymptotic spreading speed $c^{*}>0$ . Furthermore, if (vi) also holds, then for any $c\ge c^{\mathrm{<mo>\; *\; </mo>}}$ , there is a monotonic traveling wave solution $Y_t\left(x\right)=w\left(x-ct\right)$ of F, and satisfies the boundary conditions (11).

Proof. (i) Since

$\begin{array}{l}F\left[Y\left(\cdot -n\right)\right]\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y_t\left(y-n\right)\right)\text{d}y\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left[a{\text{e}}^{-b{Y}_t\left(y-n\right)}\left(1-Y_t\left(y-n\right)\right)Y_t\left(y-n\right)+c{Y}_t\left(y-n\right)\right]\text{d}y\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(\left(x-n\right)-\left(y-n\right)\right)[a{\text{e}}^{-b{Y}_t\left(y-n\right)}\left(1-Y_t\left(y-n\right)\right)Y_t\left(y-n\right)\\ {}_{}{}^{}+c{Y}_t\left(y-n\right)]\text{d}\left(y-n\right)\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(\left(x-n\right)-y\right)g\left(Y_t\left(y\right)\right)\text{d}y\\ =F\left[Y\right]\left(x-n\right),\end{array}$

we know (i) holds.

(ii) When $0\le Y^{\mathrm{<mo>\; *\; </mo>}}\le Y_1$ , the global stability of $E^{\mathrm{<mo>\; *\; </mo>}}$ implies that $0\le g\left(Y\right)\le Y^{\mathrm{<mo>\; *\; </mo>}}$ with $0\le Y\le Y^{\mathrm{<mo>\; *\; </mo>}}$ . Furthermore, by using of ${\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\text{d}y=1$ , We can get

$0\le F\left[Y\right]={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y_t\right)\text{d}y\le Y^{*}{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\text{d}y=Y^{*}.$

That is, (ii) holds.

(iii) $g\left(0\right)=0$ leads to $F\left(0\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(0\right)\text{d}y=0$ . In addition, $g\left(Y^{\mathrm{<mo>\; *\; </mo>}}\right)=Y^{\mathrm{<mo>\; *\; </mo>}}$ leads to $a{\text{e}}^{-b{Y}^{*}}\left(1-Y^{*}\right)Y^{*}+c{Y}^{*}=Y^{*}$ . Therefore,

$\begin{array}{c}F\left(Y^{*}\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y^{*}\right)\text{d}y\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left[a{\text{e}}^{-b{Y}^{*}}\left(1-Y^{*}\right)Y^{*}+c{Y}^{*}\right]\text{d}y\\ =Y^{*}{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\text{d}y\\ =Y^{*}.\end{array}$

The condition ${\bar{\mathcal{R}}}_0>1$ implies that $g^{\prime }\left(0\right)>1$ . Since $g\left(Y\right)=Y$ has a unique positive equilibrium $Y^{\mathrm{<mo>\; *\; </mo>}}$ on $\left(\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right]$ , for each $n\in \left(\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , there are $g\left(n\right)>n$ . Therefore, there is

$F\left[n\right]={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(n\right)\text{d}y>{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)n\text{d}y=n.$

That is, (iii) holds.

(iv) When ${\bar{\mathcal{R}}}_0>1$ , if ${\left\{Y_t\right\}}_{t\ge 0}$ increases monotonically in $\left(\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , let $Y_1\mathrm{,}{Y}_2\in \left(\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , you might as well set $Y_1\le Y_2$ , then

$Y_2=g\left(Y_1\right)=a{\text{e}}^{-b{Y}_1}\left(1-Y_1\right)Y_1+c{Y}_1\ge Y_1.$

From (iii) it can be seen that for each $Y_2\in \left(\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , there is $g\left(Y_2\right)>Y_2$ , so $0\le g\left(Y_1\right)\le g\left(Y_2\right)$ , that is, $g\left(Y_t\right)$ increases monotonically in this interval. Since $K\left(x-y\right)$ is not negative, so

$\begin{array}{c}F\left(Y_1\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y_1\left(y\right)\right)\text{d}y\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left[a{\text{e}}^{-b{Y}_1\left(y\right)}\left(1-Y_1\left(y\right)\right)Y_1\left(y\right)+c{Y}_1\left(y\right)\right]\text{d}y\\ \le {\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left[a{\text{e}}^{-b{Y}_2\left(y\right)}\left(1-Y_2\left(y\right)\right)Y_2\left(y\right)+c{Y}_2\left(y\right)\right]\text{d}y\\ ={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(Y_2\left(y\right)\right)\text{d}y\\ =F\left(Y_2\right).\end{array}$

That is, (iv) holds.

(v) If $\left\{g_j\right\}\subset C_{[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}]}$ , then for each l and Lipschitz constant $L_2$ , there is

$\begin{array}{l}\left|F\left[g_j\right]\left(x\right)-F\left[g_1\right]\left(x\right)\right|\\ =\left|{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left(g\left(g_j\left(y\right)\right)-g\left(g_1\left(y\right)\right)\right)\text{d}y\right|\\ \le \left|{\displaystyle {\int }_{\left|y\right|>l}}K\left(x-y\right)\left(g\left(g_j\left(y\right)\right)-g\left(g_1\left(y\right)\right)\right)\text{d}y\right|\\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}+\left|{\displaystyle {\int }_{\left|y\right|\le l}}K\left(x-y\right)\left(g\left(g_j\left(y\right)\right)-g\left(g_1\left(y\right)\right)\right)\text{d}y\right|\\ \le 2{\displaystyle {\int }_{\left|y\right|>l}}K\left(x-y\right)\text{d}y+L_2{\displaystyle {\int }_{\left|y\right|\le l}}K\left(x-y\right)\left|g_j\left(y\right)-g_1\left(y\right)\right|\text{d}y.\end{array}$

Since the integral of $K\left(x-y\right)$ is 1, we can take l to be arbitrarily large and make the first term arbitrarily small. And because $g_j$ converges uniformly to $g_1$ , then $F\left[g_j\right]$ pointwise converge $F\left[g_1\right]$ . That is, (v) holds.

(vi) If $\left\{g_j\right\}$ is a sequence in $C_{[\mathrm{0,}{Y}^{\mathrm{<mo>\; *\; </mo>}}]}$ , then for any constant $l>0$ , define the sequence in the compact set $\Theta$ as follows

${\tilde{g}}_j\left(x\right)={\displaystyle {\int }_{-l}^l}K\left(x-y\right)g\left(g_j\left(y\right)\right)\text{d}y\mathrm{,}x\in \Theta \mathrm{.}$

Because g is a bounded function, there is constant $Q>0$ , such that $\left|g\left(g_j\right)\right|\le Q$ . Furthermore, by the assumptions 3), we know there exists $M>0$ , such that $\left|K\left(x-y\right)\right|\le M$ , so for any $x\in \Theta$ , there is

$\begin{array}{c}\left|{\tilde{g}}_j\left(x\right)\right|=\left|{\displaystyle {\int }_{-l}^l}K\left(x-y\right)g\left(g_j\left(y\right)\right)\text{d}y\right|\\ \le {\displaystyle {\int }_{-l}^l}\left|K\left(x-y\right)g\left(g_j\left(y\right)\right)\right|\text{d}y\\ <2lMQ.\end{array}$

Therefore, $\left\{{\tilde{g}}_j\right\}$ is consistently bounded in the compact set $\Theta$ . Since the kernel function $K\left(x-y\right)$ is continuous, then for each $\epsilon >0$ , there is $\delta >0$ , when $x_1,x_2\in \Theta$ and satisfy $\left|x_1-x_2\right|<\delta$ , we have $\left|K\left(x_1-y\right)-K\left(x_2-y\right)\right|<\frac{\epsilon }{2lQ}$ , so the kernel function is consistently continuous on the compact set $\Theta$ . For each $x_1,x_2\in \Theta$ , when $\left|x_1-x_2\right|<\delta$ , there is

$\begin{array}{c}\left|{\tilde{g}}_j\left(x_1\right)-{\tilde{g}}_j\left(x_2\right)\right|=\left|{\displaystyle {\int }_{-l}^l}\left[K\left(x_1-y\right)-K\left(x_2-y\right)\right]g\left(g_j\left(y\right)\right)\text{d}y\right|\\ \le Q{\displaystyle {\int }_{-l}^l}\left|K\left(x_1-y\right)-K\left(x_2-y\right)\right|\text{d}y\\ <2lQ\frac{\epsilon }{2lQ}\\ =\epsilon ,\end{array}$

thus $\left\{{\tilde{g}}_j\right\}$ is isochronous continuous. According to the Arzelà-Ascoli theorem, arbitrary sequence $\left\{g_j\right\}$ exists subsequence $\left\{g_{j_i}\right\}$ and converges at $g_1$ on the compact set $\Theta$ . It is shown below that any bounded subset of $ℝ$ holds.

Because $\left\{g_j\right\}$ is a bounded function and g is continuous, so for any ${\epsilon }_1>0$ , there exists $l_{{\epsilon }_1}>0$ , there is

${\displaystyle {\int }_{\overline{ℝ}\backslash \left(-l_{{\epsilon }_1}\mathrm{,}{l}_{{\epsilon }_1}\right)}}K\left(x-y\right)g\left(g_j\left(y\right)\right)\text{d}y<{\epsilon }_1\mathrm{.}$

And because the subsequence $\left\{g_{ji}\right\}$ converges to $g_1$ on the compact set $\Theta$ , so we have

$\underset{x\in \Theta }{\sup }\left|{\displaystyle {\int }_{-l_{{\epsilon }_1}}^{l_{{\epsilon }_1}}}K\left(x-y\right)g\left(g_j\left(y\right)\right)\text{d}y-F\left[g_1\right]\left(x\right)\right|<{\epsilon }_2\mathrm{.}$

Let $\epsilon ={\epsilon }_1+{\epsilon }_2$ , then

$\begin{array}{l}\underset{x\in \Theta }{\sup }\left|F\left[g_{j_i}\right]\left(x\right)-F\left[g_1\right]\left(x\right)\right|\\ =\underset{x\in \Theta }{\sup }\left|{\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)g\left(g_{j_i}\left(y\right)\right)\text{d}y-F\left[g_1\right]\left(x\right)\right|\\ \le \underset{x\in \Theta }{\sup }[\left|{\displaystyle {\int }_{-l_{{\epsilon }_1}}^{l_{{\epsilon }_1}}}K\left(x-y\right)g\left(g_{j_i}\left(y\right)\right)\text{d}y-F\left[g_1\right]\left(x\right)\right|\\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}+{\displaystyle {\int }_{ℝ\backslash \left(-l_{{\epsilon }_1},l_{{\epsilon }_1}\right)}}K\left(x-y\right)g\left(g_{j_i}\left(y\right)\right)\text{d}y]\\ <{\epsilon }_1+{\epsilon }_2=\epsilon .\end{array}$

Therefore, $\left\{F\left[g_{j_i}\right]\right\}$ is converges consistently on each bounded subset $ℝ$ .

The existence of a traveling wave solution of the integrodifference model (8) is shown in Figure 5. Where the parameters take $a=55$ , $b=9$ , $c=0.63$ . Moreover $K\left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}$ , where ${\sigma }^2=0.0325$ . By calculation, we can get that $c^{*}=0.51109$ , $Y^{*}=0.4825$ , $Y_1=0.4927$ , so $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{Y}_1\right)$ . The numerical simulation results show that when $c\ge c^{\mathrm{<mo>\; *\; </mo>}}$ , the integrodifference model (8) has a monotonic traveling wave solution and satisfies the boundary condition (11).

3.4. Minus-One Bifurcation

This section analyzes a bifurcation of the integrodifference model (3) when ${\bar{\mathcal{R}}}_0>1$ . We first consider the perturbation at $E^{\mathrm{<mo>\; *\; </mo>}}=\left(X^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$

$\{\begin{array}{l}{X}_t\left(x\right)=X^{*}+{\tilde{X}}_t\left(x\right),\hfill \\ Y_t\left(x\right)=Y^{*}+{\tilde{Y}}_t\left(x\right).\hfill \end{array}$

Then the integrodifference model (3) linearizes at the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}=\left(X^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ .

Figure 5. The existence of a travelling wave solution about model (8).

$\{\begin{array}{l}{\tilde{X}}_{t+1}\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left(a_{11}{\tilde{X}}_t\left(y\right)+a_{12}{\tilde{Y}}_t\left(y\right)\right)\text{d}y,\hfill \\ {\tilde{Y}}_{t+1}\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}K\left(x-y\right)\left(a_{21}{\tilde{X}}_t\left(y\right)+a_{22}{\tilde{Y}}_t\left(y\right)\right)\text{d}y.\hfill \end{array}$ (12)

By using of the Fourier transform, the model (12) can be rewritten as

${\left[\begin{array}{c}\hat{X}\left(\omega \right)\\ \hat{Y}\left(\omega \right)\end{array}\right]}_{t+1}=\hat{K}{J}_{*}{\left[\begin{array}{c}\hat{X}\left(\omega \right)\\ \hat{Y}\left(\omega \right)\end{array}\right]}_t,$

where $J_{\mathrm{<mo>\; *\; </mo>}}$ is the Jacobian matrix (6) about the discrete model (2) at the endemic equilibrium $E_{*}=E^{*}$ . $\hat{K}$ is the diagonal matrix of the kernel function $K\left(x-y\right)$ after the Fourier transform. Meanwhile, $0<\hat{K}\le 1$ . If the kernel

function is $K\left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}$ , then $\hat{K}\left(\omega \right)={\text{e}}^{-\frac{{\sigma }^2{\omega }^2}{2}}$ . There is no dispesal when K is the unity matrix I. According to the Jury conditions, we give the following lemma.

Lemma 9. [29] For inequalities

$1-\hat{K}\left(a_{11}+a_{22}\right)+{\hat{K}}^2\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)>0,$ (13)

$1+\hat{K}\left(a_{11}+a_{22}\right)+{\hat{K}}^2\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)>0,$ (14)

$1-{\hat{K}}^2\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)>0.$ (15)

If inequality (13) does not hold for the integrodifference model (3), then a plus-one bifurcation will be generated; If inequality (14) does not hold, then a minus-one bifurcation will be generated and a stable two-cycle solution will be formed; If inequality (15) does not hold, then the integrodifference model (3) has a characteristic value at the endemic equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ as a complex root, i.e. there is a Hopf bifurcation.

Theorem 10. When ${\bar{\mathcal{R}}}_0>1$ , $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left[Y_1\mathrm{,1}\right)$ and satisfies ${\hat{K}}^2\left(1-p\right)>\frac{1}{2}$ , the integrodifference model (3) will experience a minus-one bifurcation and form a stable two-cycle solution.

Proof. When $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,1}\right)$ , then $1-p\left(A-p\right)>0$ , i.e. $1-det\left(J_{\mathrm{<mo>\; *\; </mo>}}\right)>0$ . In addition, $0\le \hat{K}\le 1$ leads to $1-{\hat{K}}^{2}det\left(J_{*}\right)>0$ . Therefore, the integrodifference model (3) does not undergo Hopf bifurcation.

When $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,1}\right)$ , $F\left(1\right)>0$ , i.e. $1-A+p\left(A-p\right)>0$ , thus $\left(1+p\right)\left(p-A+1\right)>0$ . And because $0<p,\hat{K}<1$ , so $p-A+1>0$ . For

$\begin{array}{l}1-\hat{K}A+{\hat{K}}^{2}p\left(A-p\right)\\ =1-\hat{K}A+{\hat{K}}^{2}pA-{\hat{K}}^2{p}^2\\ =\left(1-\hat{K}p\right)\left(\hat{K}p-\hat{K}A+1\right)\\ >\left(1-\hat{K}p\right)\left(1-\hat{K}\right)\\ >0.\end{array}$

Therefore, the integrodifference model (3) has no plus-one bifurcation under this condition from the lemma 9.

When $Y^{\mathrm{<mo>\; *\; </mo>}}\in \left(Y_1\mathrm{,1}\right)$ , then $p\left(A-p\right)<A$ . Since $F\left(-1\right)<0$ , then $1+A+p\left(A-p\right)<0$ , it is easy to know $A<p-1$ . And because ${\hat{K}}^2\left(1-p\right)>\frac{1}{2}$ , thus

$\begin{array}{l}1+\hat{K}A+{\hat{K}}^{2}p\left(A-p\right)\\ <1+\hat{K}A+{\hat{K}}^{2}A\\ <1+2{\hat{K}}^{2}A\\ <1+2{\hat{K}}^2\left(p-1\right)\\ <0.\end{array}$

Therefore, the lemma 9 knows that the integrodifference model (3) has a minus-one bifurcation under this condition.

When $Y^{*}=Y_1$ there is $F\left(-1\right)=0$ , and the same goes for $A-p+1=0$ . It is easy to have that

$\begin{array}{l}1+\hat{K}A+{\hat{K}}^{2}p\left(A-p\right)\\ =\left(1+\hat{K}p\right)\left(\hat{K}A-\hat{K}p+1\right)\\ =\left(1+\hat{K}p\right)\left(1-\hat{K}\right)\\ \ge 0.\end{array}$

A minus-one bifurcation of the integrodifference model (3) is shown in Figure 6. Let $a=3241.23$ , $b=18$ , $c=0.8$ , ${\sigma }^2=2.5$ , and $p=0.4$ . In this case, we have ${\bar{\mathcal{R}}}_0=16206.15>1$ , $Y^{*}=Y_1=\frac{1}{2}$ , ${\hat{K}}^2\left(1-p\right)>\frac{1}{2}$ . The numerical simulation results show that a minus-one bifurcation exists in the integrodifference model (3), and the two-cycle solution generated by the bifurcation is stable.

Figure 6. The minus-one bifurcation of the model (3).

4. Conclusion and Discussion

We considered an SIS epidemic integrodifference model with nonlinear incidence function. Because of the ease of transportation, people will move within a certain range, which has a certain impact on the spread of the disease. Therefore, we aimed to explore the impact of dispersal on the spread of disease.

Firstly, we conducted the theoretical analysis of model (2) using the theory of dynamical systems, mainly including the existence and stability of the disease free equilibrium $E_0$ and the endemic equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ , and the existence conditions of a flip bifurcation. In addition, numerical simulations show that model (2) may exhibit stable 4-cycle, 8-cycle, and eventually chaos will occur as ${\mathcal{R}}_0$ increases.

Secondly, we analyzed the dynamic behavior of model (8) based on the discrete-time model (4) with diffusion. We proved the global asymptotically stability of the equilibria $Y^{\mathrm{<mo>\; *\; </mo>}}$ and $Y^0$ of model (8), gave the conditions for the existence of a traveling wave solution and a minus-one bifurcation. Since the dynamic behavior of model (3) is the same as that of model (8), it indicates that model (3) has a traveling wave solution, and experiences a minus-one bifurcation under some conditions.

Finally, by comparing the dynamic behavior of model (2) and model (3), we concluded that the addition of diffusion does not affect its stability if the equilibrium is stable, at this time, the disease will persist. Besides, there will exist a traveling wave solution under this condition. However, adding diffusion increases the instability of system when the equilibrium is unstable, which makes the spread speed faster and accelerates the spread of disease. Besides, there may be no traveling wave solution.

This paper only analyzes the existence of a traveling wave solution. The conjecture that a traveling wave solution does not exist can be given when $c<c^{*}$ , but no detailed proof is given. In the dispersal process, we only choose the Gaussian kernel function, and the influence of different kernel functions on the stability of the model is not discussed. It is expected that it can be further studied in subsequent work.

Acknowledgements

This research was partially supported by National Natural Science Foundation of China through Grants No. 12071268 and No. 11971281.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Shoyimardonov, S.K. (2021) A Discrete-Time Epidemic SISI Model. Uzbek Mathematical Journal, 65, 153-162.
[2]	Shoyimardonov, S.K. (2021) A Non-Linear Discrete-Time Dynamical System Related to Epidemic SISI Model. Communications in Mathematics, 29, 505-525.
[3]	Eskandari, Z., Ghaziani, K. and Avazzadeh, Z. (2023) Bifurcations of a Discrete-Time SIR Epidemic Model with Logistic Growth of the Susceptible Individuals. International Journal of Biomathematics, 16, Article 2250120. https://doi.org/10.1142/S1793524522501200
[4]	Chaihao, B. and Khomrutai, S. (2023) Extinction and Permanence of a General Non-Autonomous Discrete-Time SIRS Epidemic Model. AIMS Mathematics, 8, 9624-9646. https://doi.org/10.3934/math.2023486
[5]	Castillo-Chavez, C. and Yakubu, A.A. (2002) Discrete-Time SIS Models with Simple and Complex Population Dynamics. In: Castillo-Chavez, C., Blower, S., van den Driessche, P., Kirschner, D. and Yakubu, A.A., Eds., Mathematical Approaches for Emerging and Reemerging Infectious Diseases: A Introduction, Springger-Verlag, New York, 153-163. https://doi.org/10.1007/978-1-4757-3667-0_9
[6]	Yakubu, A.A. (2007) Allee Effects in a Discrete-Time SIS Epidemic Model with Infected Newborns. Journal of Difference Equations and Applications, 13, 341-356. https://doi.org/10.1080/10236190601079076
[7]	Li, L., Sun, G.Q. and Jin, Z. (2010) Bifurcation and Chaos in an Epidemic Model with Nonlinear Incidence Rates. Applied Mathematics and Computation, 216, 1226-1234. https://doi.org/10.1016/j.amc.2010.02.014
[8]	Hu, Z., Teng, Z. and Jiang, H. (2012) Stability Analysis in a Class of Discrete SIRS Epidemic Models. Nonlinear Analysis Real World Applications, 13, 2017-2033. https://doi.org/10.1016/j.nonrwa.2011.12.024
[9]	Hui, C. and Yi, Z. (2011) The Discrete Age-Structured SEIT Model with Application to Tuberculosis Transmission in China. Mathematical and Computer Modelling, 55, 385-395. https://doi.org/10.1016/j.mcm.2011.08.017
[10]	Hui, C. and Yi, Z. (2013) Bifurcation Analysis of a Discrete SIS Model with Bilinear Incidence Depending on New Infection. Mathematical Biosciences and Engineering, 10, 1399-1417. https://doi.org/10.3934/mbe.2013.10.1399
[11]	Eskandari, Z. and Alidousti, J. (2020) Stability and Codimension 2 Bifurcations of a Discrete Time SIR Model. Journal of the Franklin Institute, 357, 10937-10959. https://doi.org/10.1016/j.jfranklin.2020.08.040
[12]	Schreiber, S.J., Huang, S. and Jiang, J. (2021) Extinction and Quasi-Stationarity for Discrete-Time, Endemic SIS and SIR Models. Society for Industrial and Applied Mathematics, 5, 2195-2217. https://doi.org/10.1137/20M1339015
[13]	Al-Basyouni, K.S. and Khan, A.Q. (2022) Discrete-Time COVID-19 Epidemic Model with Chaos, Stability and Bifurcation. Results in Physics, 43, Article ID: 106038. https://doi.org/10.1016/j.rinp.2022.106038
[14]	May, R.M. (1975) Biological Populations Obeying Difference Equations: Stable Points, Stable Cycles, and Chaos. Journal of Theoretical Biology, 51, 511-524. https://doi.org/10.1016/0022-5193(75)90078-8
[15]	May, R.M. (1975) Deterministic Models with Chaotic Dynamics. Nature, 256, 165-166. https://doi.org/10.1038/256165a0
[16]	Li, J., Ma, Z. and Brauer, F. (2017) Global Analysis of Discrete-Time SI and SIS Epidemic Models. Mathematical Biosciences and Engineering, 4, 99-710.
[17]	Watkinson, A.R. (1980) Density-Dependence in Single-Species Populations of Plants. Journal of Theoretical Biology, 83, 345-357. https://doi.org/10.1016/0022-5193(80)90297-0
[18]	Burgin, A.M. (1998) Deterministic and Stochastic Discrete-Time Epidemics with Spatial Considerations. Master’s Thesis, Texas Tech University, Lubbock.
[19]	Karen, R., Ros-Soto, K. and Castillo-Chavez, C. (2006) Epidemic Spread in Populations at Demographic Equilibrium. Contemporary Mathematics, 410, 297-310. https://doi.org/10.1090/conm/410/07733
[20]	Ros-Soto, K. (2008) Dispersal and Disease Dynamics in Populations with and without Demography. Master’s Thesis, Cornell University, Ithaca.
[21]	Pan, S. (2022) Minimal Wave Speed in an HIV-1 Virus Integrodifference System. International Journal of Biomathematics, 15, Article ID: 2150078. https://doi.org/10.1142/S1793524521500789
[22]	Aydogmus, O. (2017) A Discrete Time Model for Epidemic Spread: Traveling Waves and Spreading Speeds. Communications Faculty of Sciences University of Ankara Series A1 Mathematics and Statistics, 66, 243-252. https://doi.org/10.1501/Commua1_0000000815
[23]	Cui, J., Sun, Y. and Zhu, H. (2008) The Impact of Media on the Control of Infectious Diseases. Journal of Dynamics and Differential Equations, 20, 31-53. https://doi.org/10.1007/s10884-007-9075-0
[24]	Song, P. and Xiao, Y. (2018) Global Hopf Bifurcation of a Delayed Equation Describing the Lag Effect of Media Impact on the Spread of Infectious Disease. Journal of Dynamics and Differential Equations, 76, 1249-1267. https://doi.org/10.1007/s00285-017-1173-y
[25]	Allen, L.J. and van den Driessche, P. (2008) The Basic Reproduction Number in Some Discrete-Time Epidemic Models. Journal of Dynamics and Differential Equations, 14, 1127-1147. https://doi.org/10.1080/10236190802332308
[26]	Mokni, K. and Elaydi, S. (2020) Discrete Evolutionary Population Models: A New Approach. Journal of Biological Dynamics, 14, 454-478. https://doi.org/10.1080/17513758.2020.1772997
[27]	van den Driessche, P. and Yakubu, A.A. (2019) Disease Extinction versus Persistence in Discrete-Time Epidemic Models. Bulletin of Mathematical Biology, 81, 4412-4446. https://doi.org/10.1007/s11538-018-0426-2
[28]	van den Driessche, P. and Yakubu, A. A. (1999) SIS Nonlinear Discrete-Time Models with Two Competing Strains. MTBI Cornell University Technical Report.
[29]	Lutscher, F. (2019) Integrodifference Equations in Spatial Ecology. Springer, Cham. https://doi.org/10.1007/978-3-030-29294-2