A diffusive SIS epidemic model in a heterogeneous and periodically evolvingenvironment

1.   Introduction

It has been recognized that environmental heterogeneity and individual motility are significant factors that can affect the dynamics of infectious diseases. To investigate the roles of diffusion and spatial heterogeneity on disease dynamics, Allen et al. ^[1] proposed a frequency-dependent SIS reaction diffusion model

where $S(x, t)$ and $I(x, t)$ denote the numbers of susceptible and infected individuals at location $x$ and time $t$, respectively. $d_{S}$ and $d_{I}$ are the diffusion rates of the susceptible and infected individuals. $\beta(x)$ and $\gamma(x)$ are positive bounded H$\ddot{o}$lder continuous functions on $\overline{\Omega}$ that represent the spatially dependent rates of contact transmission and disease recovery at $x$, respectively. They obtained the explicit formula of the basic reproduction number ${\mathcal{R}}_0$ and characterized whether or not the domain is high (low) risk. They also showed that in high-risk domains (${\mathcal{R}}_0 > 1$) the disease-free equilibrium is always unstable, and there is a unique endemic equilibrium, while in low-risk domains (${\mathcal{R}}_0 < 1$), the disease-free equilibrium is stable if and only if the infected individuals have mobility above a threshold value (see ^[1] for more details). Recently, Li etc. ^[16] considered a SIS epidemic reaction-diffusion model governed by a mass action infection mechanism and linear birth-death growth. They studied the stability of the disease-free equilibrium, uniform persistence property in terms of the basic reproduction number and investigated the asymptotic profile of endemic equilibria in a heterogeneous environment when the movement rate of the susceptible and infected populations is small. Their results showed that factors such as infection mechanism, variation of total population, and population movement play vital but subtle roles in the transmission dynamics of diseases.

Considering an environment with the hostile boundary for the survival of population, such as extremely cold or hot temperature, the lack of resource, and so on, Huang, Han and Liu ^[12] modified the model (1.1) under the null Dirichlet boundary condition,

where $f(\bar{S}, \bar{I}) = \frac{\bar{S}}{\bar{S}+\bar{I}}$, $\beta(x), \gamma(x) \; \rm{and}\; \Lambda(x)$ is positive and continuous on $\overline{\Omega}$ and $\partial\Omega$ is $C^2$ smooth. They showed that the disease dies out when ${\mathcal{R}}_0 < 1$ and persists if ${\mathcal{R}}_0 > 1$.

When ecological phenomena are described by mathematical models, reaction diffusion systems are usually considered and the domains involved are fixed. However, the changing of domain plays a significance role in the survival of species and the transmission of disease, related problems attract much attention. One of them is the problem with free boundary, which is caused by behaviors of species themselves. In ^[8], Du and Lin proposed the logistic reaction diffusion model, and gave an spreading-vanishing dichotomy, that is, the population either successfully expends to the entire new environment, or eventually becomes extinct, see also some recent work (^{[10,14,15,23,24,25,29]}) for the spreading of species and (^[9,17]) for the transmission of disease. Another problem with regional change is that with evolving domain (^{[2,3,4,5,6,13]}), which is possibly caused by environment and climate. For example, according to monitoring meteorological satellites, Poyang lake in China covered an area of 1407 square kilometers on October 22, 2013, compared with 2022 square kilometers on August 7, its water area in summer is significantly larger than that in winter, the same is true of Dongting lake in China. In addition, in a biological context, the movement of cells is typically modelled as a diffusion-like process. For example, in the urodele amphibian axolotl ^[19] the pronephric duct extends caudally from the level of somite 7 to the cloaca. This is controlled by migratory cells at the advancing tip of the duct. Over the approximately 20 h this takes to complete, the length of the path in which duct-tip migration takes place increases from about 0.9 to 1.4 mm (lengths estimated from ^[19]). For infectious diseases, there are similar phenomenons. For example, lakes which habitat infected fishes are periodically evolving. The area infected with Japanese encephalitis (JE) is also periodically evolving. In fact, JE is an acute infectious disease caused by JE virus and transmitted by mosquitoes. Culex is the major vector of JE, and its survival, development and reproduction are influenced heavily by temperature and precipitation. JE virus in mosquitoes lost ability of infection under a temperature below $20^{o}C$. In winter temperatures are low and mosquitoes are inactive, there is less JE infection. When warm days come, the area infected with JE expands gradually.

As in ^[2,22], let $\Omega (t)\subset {\bf R}^n$ be a simply connected bounded shifting domain at time $t\ge 0$ with its changing boundary $\partial\Omega(t)$. What calls for special attention is that we considered $x\in\Omega(t)\subset {\bf R}^n$ with $n\leq2$ in this paper. For $n = 1$, the evolving interval $[0, x(t)]$ can be regarded as a simplified form of a lake, $0$ represents the top of the water column and $x(t)$ is the average depth of the lake, see ^[11]. Certainly, a lake is actually 3-dimensional and its water area $\Omega(t)\subset {\bf R}^3$. For $n = 2$, the evolving domain $\Omega(t)$ can be used to describe the area infected with Japanese encephalitis and the temperature of the area is above $20^{o}C$. For any point

we assume that $\bar{S}(x(t), t)$ and $\bar{I}(x(t), t)$ are the density of susceptible and infected species at position $x(t)$ and time $t\ge 0$. By Reynolds transport theorem (^[21]), we have

where $f_1(\bar{S}, \bar{I}, t) = -\beta(x)f(\bar{S}, \bar{I})\bar{I}+\gamma(x)\bar{I}+\Lambda(x), \; f_2(\bar{S}, \bar{I}, t) = \beta(x)f(\bar{S}, \bar{I})\bar{I}-\gamma(x)\bar{I}$ and ${\bf a} = \dot{x}(t)$, $\nabla \bar{S} \cdot {\bf a}$ and $\nabla \bar{I} \cdot {\bf a}$ are called advection terms while $(\nabla \cdot{\bf a})\bar{S}$ and $(\nabla \cdot{\bf a})\bar{I}$ are called dilution terms. In order to circumvent the difficulty induced by the evolving domain, we have to modify equations in (1.3). Let $y_1, y_2, \ldots, y_n$ be fixed cartesian coordinates in a fixed domain $\Omega(0)$ such that

Then $(\bar{S}, \bar{I})$ is mapped into the new vector $(S, I)$ defined as

Thus equations (1.3) can be translated to another form which are defined on the fixed domain $\Omega(0)$ with respect to $y = (y_1, y_2, \ldots, y_n)$. However, the new equations are still very complicated. To further simplify the model equations (1.3), we assume that domain evolution is uniform and isotropic. That is, the evolution of the domain takes place at the same proportion in all directions as time elapses. Mathematically, $x(t) = (x_1(t), x_2(t), \ldots, x_n(t))$ can be described as follows:

where the positive continuous function $\rho(t)$ is called evolving rate subject to $\rho(0) = 1$. Furthermore, if $\rho(t) = \rho(t+T)$ for some $T > 0$, the domain is periodically evolving, which has been discussed in ^[13,20]. If $\dot \rho(t) \ge0$, the domain is then called growing one (^[21,22]), and if $\dot \rho(t) \le0$, the domain is shrinking, see (^[27]) and references therein. Using (1.5) yields

Then (1.3) becomes

Now we transform the SIS epidemic model on the periodically evolving domain $\Omega(t)$ into the following problem in a fixed domain $\Omega(0)$:

with the initial condition

for later application, we also consider problem (1.7) with the periodic condition

where $f(S, I)$ is monotonically decreasing with respect to $I$ and increasing with respect to $S$ and $\lim\limits_{I\rightarrow0}f(S, I) = 1$, see for example, $f(S, I) = \frac {S}{S+I}$ for the standard incidence rate $\frac{\beta (x) SI}{S+I}$ in ^[12].

The remaining work is organized as follows. In Section 2, we focus on the existence and uniqueness of disease-free equilibrium (DFE). We define the basic reproduction number and analyze the stability of DFE in Section 3. The paper ends with some simulations and epidemiological explanations for our analytical findings.

2.   The existence and uniqueness of DFE

We first present the existence and uniqueness of the disease-free equilibrium $(S^{*}(y, t), 0)$. When $I = 0$, (1.7), (1.9) becomes the following problem

Let $u(y, t) = e^{n\ln\rho(t)}S(y, t)$, then problem (2.1) turns to

In order to find the positive solution of problem (2.2), we define the upper solution $\overline{u} = MW(y)$ and the lower solution $\underline{u} = \varepsilon W(y)$, where $W(y)$ satisfies the following equations

Since $e^{n\ln\rho(t)}\Lambda(\rho(t)y)$ is bounded function, then we can choose sufficiently large $M$ and small $\varepsilon$ such that

which implies

we easily see that $\overline{u}$ and $\underline{u}$ are the ordered upper and lower solution of problem (2.2). As a result, we conclude that there exists an $u^*(y, t)\in[\underline{u}, \overline{u}]$ satisfying problem (2.2). So problem (2.1) admits a positive solution $S^*(y, t)$.

To illustrate the uniqueness of the solution, let $S_1$ and $S_2$ be two solutions. Set

Clearly $\Lambda$ contains a neighbourhood of $0$. We claim that $1\in\Lambda$. Suppose not, then

Therefore

Recalling that $f(S, \ t)+K^\ast S = \Lambda(\rho(t)y)+(K^\ast-\frac{n\dot \rho(t)}{\rho(t)})S$ is increasing on $[0, \ \max S_2]$ for $K^\ast = n \max\limits_{[0, T]}\frac{\dot \rho(t)}{\rho(t)}$. Then

for $y\in\Omega(0), t > 0$. On the other hand, for $y\in\partial \Omega(0), t > 0$, $S_2(y, t)-h_0S_1(y, t) = 0$. Using the strong maximum principle we have assertions as follows.

(ⅰ) $\ S_2-h_0 S_1\geq0(\not\equiv0)$, using the strong maximum principle gives $\ S_2-h_0 S_1 > 0 \ \rm{in}\ \Omega(0)\times[0, T]$ with $\frac{\partial}{\partial S }(S_2-h_0 S_1) < 0 \ \rm{on}\ \partial\Omega(0)\times[0, T]$. Then, clearly there is some $\varepsilon > 0$ such that $S_2-h_0 S_1 \geq \varepsilon S_1$. Thus $h_0 + \varepsilon \in \Lambda$, which contradicts the maximality of $h_0$.

(ⅱ) $S_2-h_0 S_1\equiv0$ in $\overline \Omega(0)\times[0, T]$. This case is also impossible since we would have the equation $f(S_2, t) = h_0 f(S_1, t)$, but $f(S_2, t) = f(h_0S_1, t) > h_0 f(S_1, t)$.

Therefore, problem (2.1) admits only a positive periodic solution.

3.   The basic reproduction number

First, we define the basic reproduction number (${\mathcal{R}}_0$), and investigate its properties and implications for the reaction-diffusion system (1.7). Usually, the basic reproduction number is used as threshold for the transmission mechanism of the disease. Biologically, ${\mathcal{R}}_0$ is the expected number of secondary infections due to an infected individual over its infection period ^[7]. We know that for epidemic models described by spatially-independent systems, ${\mathcal{R}}_0$ can be obtained by the second generation matrix method ^[26].

First, a routine computation gives rise to the corresponding linearized system of problem (1.7) about the disease free equilibrium $(S^*(y, t), 0)$,

with the same periodic condition (1.9), where

Let $V(t, s)$ be the evolution operator of the problem

By the standard semigroup theory, it is easily seen that there exist positive constants $K$ and $c_0$ such that

Let $C_T$ be the ordered Banach space consisting of all $T-$ periodic and continuous function from $R$ to $C(\overline{\Omega}(0), R)$ with the maximum norm $\|\cdot\|$ and the positive cone $C_T^+: = \{\xi\in C_T:\xi(t)y\geq0, \forall \ t\in R, y\in\overline{\Omega}(0)\}$. The notation $\xi(y, t): = \xi(t)y$ will be adopted for any given $\xi\in C_T.$ After supposing that $\eta = (\xi, \zeta)\in C_T\times C_T$ is the density distribution of $w$ at the spatial locaton $y\in \Omega(0)$ and time $s$, we introduce the linear operator as in ^[28], which may be called as the next infection operator

It is easily seen that $L$ is positive, continuous and compact on $C_T\times C_T$. We define the spectral radius of $L$

as the basic reproduction number for periodic system (1.7), (1.9). Besides, we have the following results.

Lemma 3.1. $(i)$ ${\mathcal{R}}_0 = \mu_0$, where $\mu_0$ is the principle eigenvalue of the following periodic-parabolic eigenvalue problem

$(ii)\ \rm{sign} (1-{\mathcal{R}}_0) = \rm{sign}\lambda_0$, where $\lambda_0$ is the principal eigenvalue of the following reaction-diffusion problem

Particularly, assume that the coefficients in problem (1.7) are all positive constant, that is $\beta(\rho(t)y)\equiv\beta^*$ and $\gamma(\rho(t)y)\equiv\gamma^*$, then

where $\lambda^*$ is the principal eigenvalue of the following problem

It is easy to see that ${\mathcal{R}}_0$ is decreasing with respect to $\overline {\rho^{-2}}$ (: = $\frac 1T\int^T_0 \frac 1{\rho^2(t)}dt$).

Theorem 3.2. The following statements are valid:

$(i)$ If ${\mathcal{R}}_0 < 1$, then the disease-free equilibrium $(S^{*}(y, t), 0)$ is globally asymptotically stable for system (1.7), (1.8), that is to say, for any nonnegative solutions to problem (1.7), (1.8), we can deduce that $\lim\limits_{t\to\infty}{I(y, t)} = 0$ for $y\in\overline \Omega(0)$ and $\lim\limits_{m\rightarrow+\infty}S(y, t+mT) = S^{*}(y, t)$ for $(y, t)\in \overline \Omega(0)\times[0, T].$

$(ii)$ If ${\mathcal{R}}_0 > 1$, then there exists $\varepsilon_0 > 0$ such that any positive solution of system (1.7), (1.8) satisfies $\limsup\limits_{t\to\infty}\|(S(y, t), I(y, t))-(S^*(y, t), 0)\|\geq\varepsilon_0$.

Proof: $(i)$ It follows from Lemma 3.1 that problem (3.3) admits an eigen-pair $({\mathcal{R}}_0; \varphi, \phi)$ such that $\psi(y, t), \phi(y, t) > 0$ for $(y, t)\in\overline \Omega(0)\times[0, T]$. Letting $\bar{I}(y, t) = Me^{-\lambda t}\phi(y, t),$ where $0 < \lambda\leq \beta(\rho(t)y)(\frac{1}{{\mathcal{R}}_0}-1)$ for $(y, t)\in\overline \Omega(0)\times[0, T]$. Recalling that $0\leq f(S, I)\leq 1$, we then have

therefore $\bar{I}$ is the upper solution of the following problem

if $M$ is large enough. Since $\lim\limits_{t\rightarrow+\infty}\bar{I}(y, t) = 0,$ then $\lim\limits_{t\rightarrow+\infty}{I}(y, t) = 0$ uniformly for $y\in\overline \Omega(0)$.

The above limit implies that for any $\varepsilon > 0$, there exists $T_\varepsilon > 0$ such that $0\leq I(y, t)\leq\varepsilon$ for $y\in \overline \Omega(0)$ and $t > T_\varepsilon$, we then have

where $M^* = \max\limits_{(y, t)\in\overline \Omega(0)\times[0, T]}\{\beta(\rho(t)y)+\gamma(\rho(t)y)\}$. Assume that $\overline{S}_\varepsilon$ and $\underline{S}_\varepsilon$ are solutions of the following problems

and

We can deduce that $\overline{S}_\varepsilon$ and $\underline{S}_\varepsilon$ are the upper and lower solution of problem (1.7), (1.8), respectively. So the solution $(S(y, t), I(y, t))$ of problem (1.7), (1.8) satisfies $\underline{S}_\varepsilon\leq S(y, t)\leq \overline{S}_\varepsilon$ in $\Omega(0)\times[0, +\infty)$. Let $\overline{S}_\varepsilon^{(m)}$ and $\underline{S}_\varepsilon^{(m)}$ be the maximal and minimal sequences obtained from the following problem with initial iterations $\overline{S}_\varepsilon^{(0)} = \overline{S}_\varepsilon$ and $\underline{S}_\varepsilon^{(0)} = \underline{S}_\varepsilon$,

where $m = 1, 2, \cdots$ and

According to Lemma 3.1 in ^[18], it follows that the sequences $\overline{S}_\varepsilon^{(m)}$ and $\underline{S}_\varepsilon^{(m)}$ admit the monotone property

and the limits exist,

which means that

Recalling that

and letting $S_m(y, t) = S(y, t+mT)$ yields

Considering the system (1.7) with the initial condition $S_0(y) = S_1(y, 0)$, since by the initial condition in (3.10) for $m = 1$,

and

we see that

and using comparison principle gives that

Assume, by induction, that

we can deduce by the comparison principle that

and therefore, for $(y, t)\in \overline \Omega(0) \times [0, +\infty)$,

On the other hand, for $(y, t)\in \overline \Omega(0) \times [0, +\infty)$,

where $\underline{S}_\varepsilon^*(y, t)$ satisfies

and $\overline{S}_\varepsilon^*(y, t)$ satisfies

Due to the uniqueness of the solution to problem (2.1), we have

and then

$(ii)$ Since $\lim\limits_{I\rightarrow0}f(S, I) = 1$, take $\delta_0 = \frac{1}{2}(1-\frac{1}{{\mathcal{R}}_0}) > 0$, there exists $\varepsilon_0 > 0$ such that

if $0\leq I(y, t)\leq\varepsilon_0$.

Assume, for the sake of contradiction, that there exists a positive solution $(S, I)$ of problem (1.7), (1.8) such that

For the above given $\varepsilon_0$, there exists $T_{\varepsilon_0}$ such that

Then we have

for $y\in\Omega(0), \; t\geq T_{\varepsilon_0}$. We now choose a sufficiently small number $\eta > 0$ such that

where $\phi(y, t) > 0$ for $(y, t)\in\overline \Omega(0)\times[0, T]$ satisfies (3.3) with ${\mathcal{R}}_0 > 1$. Set $0 < \lambda_0\leq\frac{1}{2}(1-\frac{1}{{\mathcal{R}}_0})\beta(\rho(t)y)$, and direct calculations show that, $\underline{I}(y, t) = \eta e^{\lambda_0t}\phi(y, t)$ satisfies

It follows from (3.15) and the comparison principle that

therefore, $I(y, t) \rightarrow\infty$ as $t\rightarrow\infty$, which contradicts (3.14). This proves statement $(ii)$.

4.   Simulation and discussion

In this section, we first carry out numerical simulations for problem (1.7), (1.8) to illustrate the theoretical results by using Matlab. Let us fix some coefficients. Assume that

and subsequently $\lambda^* = \pi^2$, then the asymptotic behaviors of the solution to problem are shown by choosing different $\rho(t)$ and $\beta^*$.

Example 1. Fix $\beta^*_1 = 0.27$. We first choose $\rho_1(t)\equiv 1$, which means that the habitat is fixed. Direct calculations show that

It is easily seen from Figure 1 that the infected individual $I$ decays to zero.

Now we choose $\rho_2(t) = e^{0.1(1-\cos(4t))}$, it follows from (3.5) that

and

It is easily seen from Figure 2 that $I$ stabilizes to a positive periodic steady state.

One can see from the example that the infected individual vanishes in a fixed domain, but persist in a periodically evolving domain.

Example 2. Fix $\beta^*_2 = 0.3$. We first choose $\rho_3(t)\equiv 1$ and consider the corresponding problem in the fixed domain. Calculations show that

It is easily seen from Figure 3 that $I$ stabilizes to a positive periodic steady state.

If we choose $\rho_4(t) = e^{-0.2(1-\cos(4t))}$, it follows from (3.5) that

and

It is easily seen from Figure 4 that $I$ decays to zero and the infected individual vanishes eventually.

Results in the example imply that the infected individual spreads in a fixed domain, but vanishes in a periodically evolving domain.

Shifting of habitat for species or expending of infected domain for disease plays considerable biological significance, related problems have been attracting much attention. To explore the impact of the periodic evolution in habitats on the prevention and control of the infectious disease, we study a SIS reaction-diffusion model with periodical and isotropic domain evolution.

We first transform the SIS epidemic model with periodical evolving domain into a reaction-diffusion system on a fixed domain with time-dependent diffusion term, and then introduce the spatial-temporal risk index ${\mathcal{R}}_0(\rho)$ by using the next infection operator. ${\mathcal{R}}_0(\rho)$ depends on the domain evolution rate $\rho(t)$ and its average value $\overline{\rho^{-2}}: = \frac{1}{T}\int_0^T\frac{1}{\rho^2(t)}dt$ plays an important role, see the explicit formula (3.5). It is proved in Theorem 3.2 that If ${\mathcal{R}}_0 < 1$, the disease-free equilibrium $(S^{*}(y, t), 0)$ is globally asymptotically stable for system (1.7), (1.8), while for ${\mathcal{R}}_0 > 1$, there exists $\varepsilon_0 > 0$ such that any positive solution of system (1.7), (1.8) satisfies $\limsup\limits_{t\to\infty}\|(S(y, t), I(y, t))-(S^*(y, t), 0)\|\geq\varepsilon_0$, which means the disease-free equilibrium $(S^{*}(y, t), 0)$ is unsatble. Moreover, our numerical simulations show that the periodical domain evolution with large evolution rate has a negative effect on the control of the disease (see Figures 1 and 2), and that with small evolution rate has a positive effect on the control of the disease (see Figures 3 and 4). However, mathematically, we can not derive the property of the endemic equilibrium at present, which deserves further study.

Acknowledgments

We are very grateful to the anonymous referee for careful reading and helpful comments which led to improvements of our original manuscript. The first author is supported by Research Foundation of Young Teachers of Hexi University(QN2018013) and the second author is supported by the NNSF of China (Grant No. 11771381).

Conflict of interest

The authors declare there is no conflict of interest.