Dynamical analysis of a heroin–cocaine epidemic model with nonlinear incidence and spatial heterogeneity

Abstract

In this paper, we investigated a new heroin–cocaine epidemic model which incorporates spatial heterogeneity and nonlinear incidence rate. The main project of this paper is to explore the threshold dynamics in terms of the basic reproduction number ${\mathcal{R}}_0$, which was defined by applying the next-generation operator. The threshold type results shown that if ${\mathcal{R}}_0<1$, then the drug-free steady state is globally asymptotically stable. If ${\mathcal{R}}_0>1$, then heroin–cocaine spread is uniformly persistent. Furthermore, the globally asymptotic stability of the drug-free steady state has been established for the critical case of ${\mathcal{R}}_0=1$ by analysing the local asymptotic stability and global attractivity.

Keywords:

• Heroin–cocaine spread
• spatial heterogeneity
• the basic reproduction number
• uniformly persistent
• global asymptotically stable

1. Introduction

Over the past several decades, mathematical modelling as an important tool has been successfully applied to investigate infectious diseases which have posed imminent health threats to human. However, the abuse of illicit drugs is also detrimental on society, such as heroin, cocaine, which need to obtain more attention from the public health agency. When modelling heroin addiction and spread, Mackintosh and Stewart [24] have shown that classical epidemic models can also be applied to study drug addiction and informing control strategies. Based on the research, White and Comiskey [41] proposed an ODE model to consider the heroin epidemics by separating the population into three parts: susceptible individuals (S), drug users not in treatment ($U_1$) and drug users in treatment ($U_2$). The model takes the following form: (1) $\begin{array}{r}\{\begin{array}{l}{S}^{\prime }(t)=\lambda -{\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-\mu S,\\ U_1^{\prime }(t)={\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-p{U}_1+{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_1)U_1,\\ U_2^{\prime }(t)=p{U}_1-{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_2)U_2.\end{array}\end{array}$(1) The details of parameters of model (1(1) $\begin{array}{r}\{\begin{array}{l}{S}^{\prime }(t)=\lambda -{\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-\mu S,\\ U_1^{\prime }(t)={\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-p{U}_1+{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_1)U_1,\\ U_2^{\prime }(t)=p{U}_1-{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_2)U_2.\end{array}\end{array}$(1) ) are described in [41]. The authors defined the basic reproduction number ${\mathcal{R}}_0$ of the model and the sensitivity analysis are carried out. Moreover, the existence of backward bifurcation at the endemic equilibrium investigated for ${\mathcal{R}}_0=1$, and the authors concluded that effective prevention is indeed better than cure for blocking drugs abuse. Motivated by this work, Ma et al. [21–23] proposed and studied some models which incorporating relapse, media impact and psychological addicts, and the global dynamics and existence of bifurcation have been analysed.

When taking time delay for interpreting the time needed for relapse from treatment individuals ($U_2$) turn to the untreated state ($U_1$) into consideration on the base of model (1(1) $\begin{array}{r}\{\begin{array}{l}{S}^{\prime }(t)=\lambda -{\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-\mu S,\\ U_1^{\prime }(t)={\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-p{U}_1+{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_1)U_1,\\ U_2^{\prime }(t)=p{U}_1-{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_2)U_2.\end{array}\end{array}$(1) ), some DDE models have been proposed and studied (see [6,9,17,18,20,30,46,47]) and the authors mainly focused on the global dynamics and hopf bifurcation of these models. Moreover, to investigate the impact of treat age or susceptibility age on the heroin abuse, age-structured models have been proposed on the base of model (1(1) $\begin{array}{r}\{\begin{array}{l}{S}^{\prime }(t)=\lambda -{\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-\mu S,\\ U_1^{\prime }(t)={\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-p{U}_1+{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_1)U_1,\\ U_2^{\prime }(t)=p{U}_1-{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_2)U_2.\end{array}\end{array}$(1) ) by researchers (see [2,3,5,10,11,19,39]). However, the above-mentioned models are all involved in a homogeneous environment. Actually, the mobility of host population and spatial heterogeneity are also important and meaningful factors when modelling disease spreading. For example, Allen et al. [1] have proposed a diffusive SIS epidemic model and studied the effect of spatial heterogeneity on the persistence and extinction of the disease. Song et al. [33] have investigated an SEIR model in heterogeneous environment and pointed that the movement and spatial heterogeneous of the exposed and recovered individuals have a great effect on the disease dynamics. For more details about such models, one can refer literature [8,13,29,36,37,45] and references cited in. Hence, it is necessary to incorporate spatial heterogeneity for modelling. When considering the effect of spatial heterogeneity on heroin spreading, motivated by model (1(1) $\begin{array}{r}\{\begin{array}{l}{S}^{\prime }(t)=\lambda -{\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-\mu S,\\ U_1^{\prime }(t)={\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-p{U}_1+{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_1)U_1,\\ U_2^{\prime }(t)=p{U}_1-{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_2)U_2.\end{array}\end{array}$(1) ), Duan et al. proposed and studied an age-structured heroin epidemic model with reaction–diffusion and the threshold dynamics of the model has been investigated in term of the corresponding basic reproduction number ${\mathcal{R}}_0$ [8]. Wang and Sun in [37] proposed and investigated a heroin epidemic model with reaction–diffusion and homogeneous Neumann boundary conditions in a bounded domain (Ω) based on model (1(1) $\begin{array}{r}\{\begin{array}{l}{S}^{\prime }(t)=\lambda -{\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-\mu S,\\ U_1^{\prime }(t)={\displaystyle \frac{{\beta }_1{U}_{1}S}{N}}-p{U}_1+{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_1)U_1,\\ U_2^{\prime }(t)=p{U}_1-{\displaystyle \frac{{\beta }_2{U}_1{U}_2}{N}}-(\mu +{\delta }_2)U_2.\end{array}\end{array}$(1) ). By using the theory of the next-generation operator, the basic reproduction number ${\mathcal{R}}_0$ has been defined, which is regard as the threshold that determines whether the heroin extinct or not. The threshold dynamics in the light of ${\mathcal{R}}_0$ are achieved. Furthermore, the dynamics for a critical case of ${\mathcal{R}}_0$ equal to 1 has been improved by Xu and Geng in [44].

The above-mentioned studies of drug abuse dynamics models mostly focus on one kind of drugs such as heroin. In fact, heroin and cocaine users are overlapping populations, i.e. many heroin users are also use cocaine at the same time [15]. When considering the effect of interaction between the epidemic of heroin and cocaine spread, Chekroun et al. [3] have proposed and studied an age-structured model with both the two kinds of drugs, the threshold dynamics have been studied in term of the basic reproduction number. However, spatial heterogeneity and movement of the host population have not been taken into consideration in [3]. To the best of our knowledge, the heroin–cocaine co-infection models in which population density depends on both time and space variables are rarely studied. Motivated by [3,37], we propose an reaction–diffusion model which describing the dynamics of heroin and cocaine epidemics, the model takes the following form: (2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) with the homogeneous Neumann boundary conditions (3) $\begin{array}{r}\frac{\mathrm{\partial }S}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }{I}_2}{\mathrm{\partial }\nu }=\frac{\mathrm{\partial }R}{\mathrm{\partial }\nu }=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega },\end{array}$(3) and initial conditions (4) $\begin{array}{r}(S(0,x),I_1(0,x),I_2(0,x),R(0,x))=(S_0(x),I_{10}(x),I_{20}(x),R_0(x)),x\in \mathrm{\Omega }.\end{array}$(4) Here $\mathrm{\Omega }\subset {\mathbb{R}}^n$ is a bounded domain with smooth boundary denoted by $\mathrm{\partial }\mathrm{\Omega }$. $S(t,x)$, $I_1(t,x)$, $I_2(t,x)$, $R(t,x)$ denote the density of susceptible, heroin users, cocaine users and recovered individuals at time t and location x, respectively. $d_i(i=1,2,3,4)$ represent the corresponding diffusion rates of $S(t,x)$, $I_1(t,x)$, $I_2(t,x)$, $R(t,x)$. ${\beta }_1(x)$ and ${\beta }_2(x)$ are the incidence rates of heroin users and cocaine users, respectively. $\lambda (x)$ is the recruitment rate. $k_1$ is the rate that an individual recovered from the consumption of heroin becomes a cocaine user, $k_2$ is the rate that an individual recovered from the consumption of cocaine becomes a heroin user. $\mu (x)$, ${\mu }_1(x)$, ${\mu }_2(x)$ and ${\mu }_3(x)$ are the death rates of $S(t,x)$, $I_1(t,x)$, $I_2(t,x)$ and $R(t,x)$, respectively. ${\theta }_1(x)$ and ${\theta }_2(x)$ are the recovery rates of heroin users and cocaine users, respectively. ${\delta }_1(x)$ and ${\delta }_2(x)$ are the rates at which a recovery individual becomes a heroin user and cocaine user, respectively. All the location-dependent parameters are continuous and strictly positive defined on $\overline{\mathrm{\Omega }}$. The incidence rates of heroin users and cocaine users are given by general nonlinear functions $f(S,I_1)$ and $g(S,I_2)$ satisfy the following conditions:

(A1)

$f(0,I_1)=f(S,0),g(0,I_2)=g(S,0)=0,f(S,I_1)>0,g(S,I_2)>0,f(S,I_1)\le S{I}_1,g(S,I_2)\le S{I}_2$, for $S,I_1,I_2\ge 0$;

(A2)

$\frac{\mathrm{\partial }f(S,I_1)}{\mathrm{\partial }S}>0$, $\frac{\mathrm{\partial }f(S,I_1)}{\mathrm{\partial }{I}_1}>0$, $\frac{{\mathrm{\partial }}^{2}f(S,I_1)}{{\mathrm{\partial }}^2{I}_1}\le 0$; $\frac{\mathrm{\partial }g(S,I_2)}{\mathrm{\partial }S}>0$, $\frac{\mathrm{\partial }g(S,I_2)}{\mathrm{\partial }{I}_2}>0$, $\frac{{\mathrm{\partial }}^{2}g(S,I_2)}{{\mathrm{\partial }}^2{I}_2}\le 0$.

The rest of this paper is organized as follows. The well-posedness of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) is established and the basic reproduction number ${\mathcal{R}}_0$ is defined in Section 2. In Section 3, we investigate the threshold dynamics of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) in term of ${\mathcal{R}}_0=1$. A brief conclusion ends the paper.

2. Well-posedness

Let $\mathbb{X}:=C(\overline{\mathrm{\Omega }},{\mathbb{R}}^4)$ be the Banach space with the supremum norm $\Vert \cdot {\Vert }_{\mathbb{X}}$ and ${\mathbb{X}}_{+}:=C(\overline{\mathrm{\Omega }},{\mathbb{R}}_{+}^4)$. Then $(\mathbb{X},{\mathbb{X}}_{+})$ is a ordered Banach space. Denote $\overline{\phi }=\underset{x\in \overline{\mathrm{\Omega }}}{max}\{\phi (x)\}$, $\underset{\_}{\phi }=\underset{x\in \overline{\mathrm{\Omega }}}{min}\{\phi (x)\}$, $f_x(x,y)=\frac{\mathrm{\partial }f(x,y)}{\mathrm{\partial }x}$, $f_y(x,y)=\frac{\mathrm{\partial }f(x,y)}{\mathrm{\partial }y}$.

Denote ${\rho }_1(x)=\mu (x)$, ${\rho }_2(x)={\mu }_1(x)+{\theta }_1(x)$, ${\rho }_3(x)={\mu }_2(x)+{\theta }_2(x)$, ${\rho }_4(x)={\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x)$, and ${\mathrm{\Gamma }}_i(t,x,y)(i=1,2,3,4)$ are the Green function associated with $d_i\mathrm{\Delta }-{\rho }_i(x)(i=1,2,3,4)$ subject to the Neumann boundary condition. Define ${\mathbb{T}}_i(t)(i=1,2,3,4):C(\overline{\mathrm{\Omega }},\mathbb{R})\to C(\overline{\mathrm{\Omega }},\mathbb{R})$ are the $C_0$ semigroup associated with $d_i\mathrm{\Delta }-{\rho }_i(x)(i=1,2,3,4)$ subject to the Neumann boundary condition. Then we have $({\mathbb{T}}_i(t)\psi )(x)={\int }_{\mathrm{\Omega }}{\mathrm{\Gamma }}_i(t,x,y)\psi (y)\mathrm{d}y,t>0,\psi \in C(\overline{\mathrm{\Omega }},\mathbb{R}).$It follows from [31] that ${\mathbb{T}}_i(t)$ is compact and strongly positive for all t>0, then there exists a constant M>0 such that $\Vert {\mathbb{T}}_i(t)\Vert \le M{e}^{w_{i}t}$ where $w_i<0$ are the principle eigenvalue of $d_i\mathrm{\Delta }-{\rho }_i(x)(i=1,2,3,4)$ subject to the Neumann boundary condition.

Define operators $\mathbb{F}=({\mathbb{F}}_1,{\mathbb{F}}_2,{\mathbb{F}}_3,{\mathbb{F}}_4):={\mathbb{X}}^{+}\to \mathbb{X}$ by $\begin{array}{rl}{\mathbb{F}}_1(t,\psi )& =\lambda (x)-{\beta }_1(x)f({\psi }_1(x),{\psi }_2(x))-{\beta }_2(x)g({\psi }_1(x),{\psi }_3(x)),\\ {\mathbb{F}}_2(t,\psi )& ={\beta }_1(x)f({\psi }_1(x),{\psi }_2(x))+k_2{\theta }_2(x){\psi }_3(x)+{\delta }_1(x){\psi }_4(x),\\ {\mathbb{F}}_3(t,\phi )& ={\beta }_2(x)g({\psi }_1(x),{\psi }_3(x))+k_1{\theta }_1(x){\psi }_2(x)+{\delta }_2(x){\psi }_4(x),\\ {\mathbb{F}}_4(t,\phi )& =(1-k_1){\theta }_1(x){\psi }_2(x)+(1-k_2){\theta }_2(x){\psi }_3(x),\end{array}$where $\psi (x)=({\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4{)}^T\in \mathbb{X}$. Let $\mathbb{T}(t)=diag({\mathbb{T}}_i(t))$, then we can rewrite model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) as (5) $\begin{array}{r}u(t,x)=(\mathbb{T}(t)\psi )(x)+{\int }_0^t\mathbb{T}(t-s)\mathbb{F}(u(s,x))\mathrm{d}s,\end{array}$(5) where $u(t,x)=(S(x,t),I_1(t,x),I_2(t,x),R(t,x){)}^T$ and $\psi (x)=(S_0(x),I_{10}(x),I_{20}(x),R_0(x){)}^T$. It then follows from [27] that the following lemma holds.

Lemma 2.1

For any $\psi \in {\mathbb{X}}_{+}$, model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) has a unique mild solution $u(t,\cdot ,\psi )\in {\mathbb{X}}_{+}$ on the interval of $[0,{\tau }_{\mathrm{\infty }})$ where ${\tau }_{\mathrm{\infty }}\le \mathrm{\infty }$. Moreover, this solution is a classical solution.

Based on Lemma 2.1, we are ready to present the following result.

Theorem 2.2

Model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) admits a unique solution $u(t,x,\psi )=(S(t,x),I_1(t,x),I_2(t,x),R(t,x))\in {\mathbb{X}}_{+}$ on $[0,+\mathrm{\infty })$ for any $\psi \in {\mathbb{X}}_{+}$ with $u(0,x,\psi )=\psi$. Furthermore, the solution semiflow ${\mathrm{\Phi }}_t=u(t,\cdot ):{\mathbb{X}}_{+}\to {\mathbb{X}}_{+}$ of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) is defined by $\mathrm{\Phi }(t)(\psi )=u(t,\cdot ,\psi ),t\ge 0$ admits a global compact attractor.

Proof.

Based on Lemma 2.1, assume ${\tau }_{\mathrm{\infty }}<\mathrm{\infty }$, then we have $\Vert u(t,\cdot ,\psi ){\Vert }_{\mathbb{X}}\to \mathrm{\infty }$ as $t\to \mathrm{\infty }$ by Theorem 2 in [27]. It follows from the first equation of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that $\frac{\mathrm{\partial }S}{\mathrm{\partial }t}\le d_1\mathrm{\Delta }S(t,x)+\overline{\lambda }-\underset{\_}{\mu }S,t\in [0,{\tau }_{\mathrm{\infty }}),x\in \mathrm{\Omega }.$By the comparison principle and Lemma 2 in [13], there exists a constant $M_1>0$ such that $S(t,x)\le M_1$, for $t\in [0,{\tau }_{\mathrm{\infty }}),x\in \overline{\mathrm{\Omega }}$. Then we have $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}\le d_2\mathrm{\Delta }{I}_1(t,x)+\overline{{\beta }_1}{f}_{I_1}(M_1,0)I_1+k_2\overline{{\theta }_2}{I}_2-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})I_1+\overline{{\delta }_1}R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}\le d_3\mathrm{\Delta }{I}_2(t,x)+\overline{{\beta }_2}{g}_{I_2}(M_1,0)I_2+k_1\overline{{\theta }_1}{I}_1-(\underset{\_}{{\mu }_2}+\underset{\_}{{\theta }_2})I_2+\overline{{\delta }_2}R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}\le d_4\mathrm{\Delta }R(t,x))+(1-k_1)\overline{{\theta }_1}{I}_1(t,x)+(1-k_2)\overline{{\theta }_2}{I}_2\\ -(\underset{\_}{{\mu }_3}+\underset{\_}{{\delta }_1}+\underset{\_}{{\delta }_2})R,t>0,x\in \mathrm{\Omega }.\end{array}$Considering the following system: (6) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{w}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{w}_1(t,x)+\overline{{\beta }_1}{f}_{w_1}(M_1,0)w_1+k_2\overline{{\theta }_2}{w}_2\\ -(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})w_1+\overline{{\delta }_1}{w}_3,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{w}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{w}_2(t,x)+\overline{{\beta }_2}{g}_{w_2}(M_1,0)w_2+k_1\overline{{\theta }_1}{w}_1\\ -(\underset{\_}{{\mu }_2}+\underset{\_}{{\theta }_2})w_2+\overline{{\delta }_2}{w}_3,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{w}_3(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }{w}_3(t,x)+(1-k_1)\overline{{\theta }_1}{w}_1(t,x)+(1-k_2)\overline{{\theta }_2}{w}_2\\ -(\underset{\_}{{\mu }_3}+\underset{\_}{{\delta }_1}+\underset{\_}{{\delta }_2})w_3,t>0,x\in \mathrm{\Omega }.\end{array}\end{array}$(6) It follows from the standard Krein–Rutman theorem [7] that model (6(6) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{w}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{w}_1(t,x)+\overline{{\beta }_1}{f}_{w_1}(M_1,0)w_1+k_2\overline{{\theta }_2}{w}_2\\ -(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})w_1+\overline{{\delta }_1}{w}_3,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{w}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{w}_2(t,x)+\overline{{\beta }_2}{g}_{w_2}(M_1,0)w_2+k_1\overline{{\theta }_1}{w}_1\\ -(\underset{\_}{{\mu }_2}+\underset{\_}{{\theta }_2})w_2+\overline{{\delta }_2}{w}_3,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{w}_3(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }{w}_3(t,x)+(1-k_1)\overline{{\theta }_1}{w}_1(t,x)+(1-k_2)\overline{{\theta }_2}{w}_2\\ -(\underset{\_}{{\mu }_3}+\underset{\_}{{\delta }_1}+\underset{\_}{{\delta }_2})w_3,t>0,x\in \mathrm{\Omega }.\end{array}\end{array}$(6) ) admits one principle eigenvalue ${\lambda }_0$ which corresponding to a strongly positive eigenfunction $\phi =({\phi }_1,{\phi }_2,{\phi }_3)$, and hence, model (6(6) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{w}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{w}_1(t,x)+\overline{{\beta }_1}{f}_{w_1}(M_1,0)w_1+k_2\overline{{\theta }_2}{w}_2\\ -(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})w_1+\overline{{\delta }_1}{w}_3,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{w}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{w}_2(t,x)+\overline{{\beta }_2}{g}_{w_2}(M_1,0)w_2+k_1\overline{{\theta }_1}{w}_1\\ -(\underset{\_}{{\mu }_2}+\underset{\_}{{\theta }_2})w_2+\overline{{\delta }_2}{w}_3,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{w}_3(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }{w}_3(t,x)+(1-k_1)\overline{{\theta }_1}{w}_1(t,x)+(1-k_2)\overline{{\theta }_2}{w}_2\\ -(\underset{\_}{{\mu }_3}+\underset{\_}{{\delta }_1}+\underset{\_}{{\delta }_2})w_3,t>0,x\in \mathrm{\Omega }.\end{array}\end{array}$(6) ) has a solution $\zeta e^{{\lambda }_{0}t}\phi (x)$ for $t\ge 0$, where $\zeta >0$ is a constant satisfying $\zeta \phi =(w_1(0,x),w_2(0,x),w_3(0,x))\ge (I_1(0,x),I_2(0,x),R(0,x))$ for $x\in \overline{\mathrm{\Omega }}$. Then it follows from the comparison principle that $(I_1(t,x),I_2(t,x),R(t,x))\le \zeta e^{{\lambda }_{0}t}\phi (x),t\in [0,{\tau }_{\mathrm{\infty }}),x\in \overline{\mathrm{\Omega }}.$Thus there exists a constant $M_2>0$ such that $I_1(t,x)\le M_2$, $I_2(t,x)\le M_2$,$R(t,x)\le M_2$, for $t\in [0,{\tau }_{\mathrm{\infty }}),x\in \overline{\mathrm{\Omega }}$, which leads to a contradiction.

Denote ${\pi }_n$ be the eigenvalue of $d_2\mathrm{\Delta }-({\mu }_1(x)+{\theta }_1(x))$ with eigenfunction of ${\phi }_n(s)$ which satisfying ${\pi }_1>{\pi }_2\ge {\pi }_3\ge \cdots \ge {\pi }_n\ge \cdots$. It follows from [12] that ${\mathrm{\Gamma }}_2(t,x,y)=\sum _{n\ge 1}{e}^{{\pi }_{n}t}{\phi }_n(x){\phi }_n(y)$. Since ${\phi }_n(x)$ is uniformly bounded, then there exists a $\sigma >0$ such that ${\mathrm{\Gamma }}_2(t,x,y)\le \sigma \sum _{n\ge 1}{e}^{{\pi }_{n}t}$, for t>0.

In addition, assume that ${\tau }_i$ be the eigenvalue of $d_2\mathrm{\Delta }-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})$ subject to the Neumann boundary condition, which satisfies ${\tau }_1=-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})>{\tau }_2\ge {\tau }_3\ge \cdots \ge {\tau }_n\ge \cdots$. It then follows from [35] that ${\tau }_i\ge {\pi }_i$ for $i\in {\mathbb{N}}_{+}$. Furthermore, for some ${\sigma }_{\ast }>0$ we have ${\mathrm{\Gamma }}_2(t,x,y)\le \sigma \sum _{n\ge 1}{e}^{{\pi }_{n}t}\le {\sigma }_{\ast }{e}^{{\pi }_{1}t}={\sigma }_{\ast }{e}^{-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})t},\mathrm{\forall }t>0.$It follows from the comparison principle [13] that there exists a $t_0>0$ and $M_3>0$ such that $S(t,x)\le M_3$, $t\ge t_0$, $x\in \overline{\mathrm{\Omega }}$. Now, we define $P(t)={\int }_{\mathrm{\Omega }}(S(t,x)+I_1(t,x)+I_2(t,x)+R(t,x))\mathrm{d}x,$then we have $\frac{\mathrm{d}P}{\mathrm{d}t}\le {\int }_{\mathrm{\Omega }}\lambda (x)\mathrm{d}x-{\mu }_{0}P,$where ${\mu }_0=min\{\underset{\_}{\mu },{\underset{\_}{\mu }}_1,{\underset{\_}{\mu }}_2,{\underset{\_}{\mu }}_3\}$. Thus there exist $t_1>0$ and $M_4>0$ such that $P(t)\le M_4$ for $t\ge t_1$. Let $t_2=min\{t_0,t_1\}$, for all $t\ge t_2$, utilizing the comparison principle and (5(5) $\begin{array}{r}u(t,x)=(\mathbb{T}(t)\psi )(x)+{\int }_0^t\mathbb{T}(t-s)\mathbb{F}(u(s,x))\mathrm{d}s,\end{array}$(5) ), we obtain $\begin{array}{rl}{I}_1(t,x)& ={\mathbb{T}}_2(t)I_1(t_2,x)+{\int }_{t_2}^t{\mathbb{T}}_2(t-\tau )[{\beta }_1(x)f(S(\tau ,x),I_1(\tau ,x))\\ & +k_2{\theta }_2(x)I_2(\tau ,x)+{\delta }_1(x)R(\tau ,x)]\mathrm{d}\tau \\ & \le M{e}^{w_2(t-t_2)}\Vert I_1(t_2,x)\Vert +{\int }_{t_2}^t{\int }_{\mathrm{\Omega }}{\mathrm{\Gamma }}_2(t,x,y)[{\beta }_1(y)f(S(\tau ,y),I_1(\tau ,y))\\ & +k_2{\theta }_2(y)I_2(\tau ,y)+{\delta }_1(y)R(\tau ,y)]\mathrm{d}y\mathrm{d}\tau \\ & \le M{e}^{w_2(t-t_2)}\Vert I_1(t_2,x)\Vert +{\int }_{t_2}^t{\int }_{\mathrm{\Omega }}{\sigma }_{\ast }{e}^{-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})(t-\tau )}[\overline{{\beta }_1}{f}_{I_1}(M_3,0)I_1(\tau ,y)\\ & +k_2\overline{{\theta }_2}{I}_2(\tau ,y)+\overline{{\delta }_1}R(\tau ,y)]\mathrm{d}y\mathrm{d}\tau \\ & \le M{e}^{w_2(t-t_2)}\Vert I_1(t_2,x)\Vert +{\sigma }_{\ast }{M}_4[\overline{{\beta }_1}{f}_{I_1}(M_3,0)+k_2\overline{{\theta }_2}+\overline{{\delta }_1}]{\int }_{t_2}^t{e}^{-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})(t-\tau )}\mathrm{d}\tau \\ & =M{e}^{w_2(t-t_2)}\Vert I_1(t_2,x)\Vert +{\sigma }_{\ast }{M}_4[\overline{{\beta }_1}{f}_{I_1}(M_3,0)+k_2\overline{{\theta }_2}+\overline{{\delta }_1}]\frac{1-e^{-(\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1})(t-t_2)}}{\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1}}\\ & \le M{e}^{w_2(t-t_2)}\Vert I_1(t_2,x)\Vert +\frac{{\sigma }_{\ast }{M}_4[\overline{{\beta }_1}{f}_{I_1}(M_3,0)+k_2\overline{{\theta }_2}+\overline{{\delta }_1}]}{\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1}},\end{array}$which implies that $\underset{t\to \mathrm{\infty }}{lim sup}\Vert I_1(t,x)\Vert \le \frac{{\sigma }_{\ast }{M}_4[\overline{{\beta }_1}{f}_{I_1}(M_3,0)+k_2\overline{{\theta }_2}+\overline{{\delta }_1}]}{\underset{\_}{{\mu }_1}+\underset{\_}{{\theta }_1}}.$Similarly, there exists two positive constants $M_5>0,M_6>0$ such that $\underset{t\to \mathrm{\infty }}{lim sup}\Vert I_2(t,x)\Vert \le M_5$, $\underset{t\to \mathrm{\infty }}{lim sup}\Vert R(t,x)\Vert \le M_6$.

The above discussions imply that the system is point dissipative. In addition, it follows from [42] that the solution semiflow $\mathrm{\Phi }(t)$ is compact for t>0. Then, it follows from [14] that $\mathrm{\Phi }(t)$ admits a global compact attractor.

According to [13] that model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) admits a unique drug-free steady state $E_0(x)=(S_0^{\ast }(x),0,0,0)$, linearizing model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) yields that (7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0)I_1+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0)I_2+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(7) Let $(I_1(t,x),I_2(t,x),R(t,x))=e^{\eta t}({\psi }_2(x),{\psi }_3(x),{\psi }_4(x))$, it then follows from model (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0)I_1+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0)I_2+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(7) ) that (8) $\begin{array}{r}\{\begin{array}{l}\lambda {\psi }_2(x)=d_2\mathrm{\Delta }{\psi }_2(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0){\psi }_2+k_2{\theta }_2(x){\psi }_3\\ -({\mu }_1(x)+{\theta }_1(x)){\psi }_2+{\delta }_1(x){\psi }_4,t>0,x\in \mathrm{\Omega },\\ \lambda {\psi }_3(x)=d_3\mathrm{\Delta }{\psi }_3(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0){\psi }_3+k_1{\theta }_1(x){\psi }_2\\ -({\mu }_2(x)+{\theta }_2(x)){\psi }_3+{\delta }_2(x){\psi }_4,t>0,x\in \mathrm{\Omega },\\ \lambda {\psi }_4(x)=d_4\mathrm{\Delta }{\psi }_4(t,x)+(1-k_1){\theta }_1(x){\psi }_2(t,x)+(1-k_2){\theta }_2(x){\psi }_3\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x)){\psi }_4,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{\psi }_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\psi }_3(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\psi }_4(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(8) It follow from Theorem 7.6.1 in [31] that eigenvalue problem (8(8) $\begin{array}{r}\{\begin{array}{l}\lambda {\psi }_2(x)=d_2\mathrm{\Delta }{\psi }_2(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0){\psi }_2+k_2{\theta }_2(x){\psi }_3\\ -({\mu }_1(x)+{\theta }_1(x)){\psi }_2+{\delta }_1(x){\psi }_4,t>0,x\in \mathrm{\Omega },\\ \lambda {\psi }_3(x)=d_3\mathrm{\Delta }{\psi }_3(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0){\psi }_3+k_1{\theta }_1(x){\psi }_2\\ -({\mu }_2(x)+{\theta }_2(x)){\psi }_3+{\delta }_2(x){\psi }_4,t>0,x\in \mathrm{\Omega },\\ \lambda {\psi }_4(x)=d_4\mathrm{\Delta }{\psi }_4(t,x)+(1-k_1){\theta }_1(x){\psi }_2(t,x)+(1-k_2){\theta }_2(x){\psi }_3\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x)){\psi }_4,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{\psi }_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\psi }_3(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\psi }_4(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(8) ) admits a principle eigenvalue ${\lambda }_0(S_0^{\ast })$ with a strictly positive eigenfunction. According to the idea of [38], we will define the basic reproduction number for model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ). To this end, let $F(x)=\left(\begin{array}{ccc}{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0)& 0& 0\\ 0& {\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0)& 0\\ 0& 0& 0\end{array}\right),$and $V(x)=\left(\begin{array}{ccc}{\mu }_1(x)+{\theta }_1(x)& -k_2{\theta }_2(x)& -{\delta }_1(x)\\ -k_1{\theta }_1(x)& {\mu }_2(x)+{\theta }_2(x)& -{\delta }_2(x)\\ -(1-k_1){\theta }_1(x)& -(1-k_2){\theta }_2(x)& {\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x)\end{array}\right).$Assume $v=(I_1,I_2,R{)}^T$, $\mathrm{d}\mathrm{\Delta }v=(d_2\mathrm{\Delta }{I}_1,d_3\mathrm{\Delta }{I}_2,d_4\mathrm{\Delta }R{)}^T$ and $\mathrm{\Psi }(t):C(\overline{\mathrm{\Omega }},{\mathbb{R}}^3)\to C(\overline{\mathrm{\Omega }},{\mathbb{R}}^3)$ be the $C_0$-semigroup generated by the following system: (9) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }t}}=\mathrm{d}\mathrm{\Delta }v-V(x)v,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }\nu }}={\displaystyle \frac{\mathrm{\partial }{I}_2}{\mathrm{\partial }\nu }}={\displaystyle \frac{\mathrm{\partial }R}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(9) Assume that the heroin and cocaine drug users are introduced at time t = 0, where the distribution of initial drug users is $\psi =({\psi }_2,{\psi }_3,{\psi }_4)$, then $\mathrm{\Psi }(t)\psi$ is the distribution of those drug users as time evolves to time t under the synthetical influences of mobility, mortality and transfer of individuals in infected compartments, and hence, $F(x)\mathrm{\Psi }(t)\psi (x)$ represents the distribution of new drug users at time t. Then, by [38], the total number of new drug users is given by $\mathcal{L}(\psi )(x)={\int }_0^{\mathrm{\infty }}F(x)\mathrm{\Psi }(t)\psi (x)\mathrm{d}t,$it follows from [38] that $\mathcal{L}$ is a continuous and positive operator which maps the initial distribution of drug users to the total drug users produced during the infection period. Furthermore, the basic reproduction number is defined by ${\mathcal{R}}_0=r(\mathcal{L}),$where $r(\mathcal{L})$ represents the spectral radius of the operator $\mathcal{L}$. Consequently, the following result holds followed by [38].

Lemma 2.3

The principal eigenvalue ${\lambda }_0(S_0^{\ast }(x))$ and ${\mathcal{R}}_0-1$ have the same sign, and the drug-free steady state $E_0=(S_0^{\ast }(x),0,0,0)$ is asymptotically stable when ${\mathcal{R}}_0<1$.

3. Stability and persistence

Based on the above discussions, we first to investigate the global stability of disease-free steady state $E_0(x)$.

Theorem 3.1

If ${\mathcal{R}}_0<1$, then the drug-free steady state $E_0(x)$ of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) is globally asymptotically stable.

Proof.

Since ${\mathcal{R}}_0<1$, according to Lemma 2.3 that ${\lambda }_0<0$. By continuity, there exists a small $ϵ>0$ such that ${\lambda }_0^{ϵ}<0$, where ${\lambda }_0^{ϵ}$ is the principle eigenvalue (10) $\begin{array}{r}\{\begin{array}{l}{\lambda }_0^{ϵ}{\psi }_2(x)=d_2\mathrm{\Delta }{\psi }_2(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x)+ϵ,0){\psi }_2+k_2{\theta }_2(x){\psi }_3\\ -({\mu }_1(x)+{\theta }_1(x)){\psi }_2+{\delta }_1(x){\psi }_4,t>0,x\in \mathrm{\Omega },\\ {\lambda }_0^{ϵ}{\psi }_3(x)=d_3\mathrm{\Delta }{\psi }_3(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x)+ϵ,0){\psi }_3+k_1{\theta }_1(x){\psi }_2\\ -({\mu }_2(x)+{\theta }_2(x)){\psi }_3+{\delta }_2(x){\psi }_4,t>0,x\in \mathrm{\Omega },\\ {\lambda }_0^{ϵ}{\psi }_4(x)=d_4\mathrm{\Delta }{\psi }_4(t,x)+(1-k_1){\theta }_1(x){\psi }_2(t,x)+(1-k_2){\theta }_2(x){\psi }_3\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x)){\psi }_4,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{\psi }_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\psi }_3(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\psi }_4(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(10) It follows from the first equation (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that ${\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }t}}\le d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x),t>0,x\in \mathrm{\Omega }.$Then, by the comparison principle, there exists a $t_3>0$ such that $S(t,x)\le S_0(x)+ϵ$ for all $t\ge t_3$ and $x\in \overline{\mathrm{\Omega }}$. Furthermore, it follows from model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that (11) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}\le d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x)+ϵ,0)I_1+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}\le d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x)+ϵ,0)I_2+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}\le d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(11) In addition, for any $\psi \in {\mathbb{X}}_{+}$, assume that there exists a ${\xi }_{\ast }>0$ such that $(I_1(t_3,x,\psi ),I_2(t_3,x,\psi ),R_1(t_3,x,\psi ))\le {\xi }_{\ast }({\psi }_{I_1}^{ϵ}(x),{\psi }_{I_2}^{ϵ}(x),{\psi }_R^{ϵ}(x)),$where $({\psi }_{I_1}^{ϵ}(x),{\psi }_{I_2}^{ϵ}(x),{\psi }_R^{ϵ})(x)$ is a strongly positive eigenfunction corresponding to ${\lambda }_0^{ϵ}$. By the comparison principle, we can obtain that $(I_1(t,x,\psi ),I_2(t,x,\psi ),R_1(t,x,\psi ))\le {\xi }_{\ast }{e}^{{\lambda }_0^{ϵ}(t-t_3)}({\psi }_{I_1}^{ϵ}(x),{\psi }_{I_2}^{ϵ}(x),{\psi }_R^{ϵ}(x)),t\ge t_3.$Since ${\lambda }_0^{ϵ}<0$, it then follows that $\underset{t\to \mathrm{\infty }}{lim}(I_1(t,x,\psi ),I_2(t,x,\psi ),R_1(t,x,\psi ))=(0,0,0).$Thus the equation of S in model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) is asymptotic to $\frac{\mathrm{\partial }S}{\mathrm{\partial }t}=d_1\mathrm{\Delta }S+\lambda (x)-\mu (x)S.$It follows from Lemma 2 in [13] and Corollary 4.3 in [34] that $\underset{t\to \mathrm{\infty }}{lim}S(t,x,\psi )=S_0^{\ast }(x)$, and hence, it completes the proof.

Lemma 3.2

Let $u(t,x,\psi )$ be the solution of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) with $u_0=\psi \in {\mathbb{X}}_{+}$. If there exists some $t_0>0$ such that $I_1(t_0,x,\psi )\not\equiv 0$, $I_2(t_0,x,\psi )\not\equiv 0$ and $R(t_0,x,\psi )\not\equiv 0$, then $I_1(t,x,\psi )>0$, $I_2(t,x,\psi )>0$ and $R(t,x,\psi )>0$ for $t>t_0$ and $x\in \overline{\mathrm{\Omega }}$. Furthermore, there exists some ${\eta }_{\ast }>0$ such that $\underset{t\to \mathrm{\infty }}{lim inf}S(t,x,\psi )\ge {\eta }_{\ast }$ for any $\psi \in {\mathbb{X}}_{+}$.

Proof.

By applying the similar method of [13,45], it follows from Lemma 2.1 and equation of $I_1,I_2$ and R that $\{\begin{array}{l}{d}_2\mathrm{\Delta }{I}_1-{\displaystyle \frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }t}}-({\mu }_1(x)+{\theta }_1(x))I_1\le 0,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega },\end{array}$and $\{\begin{array}{l}{d}_3\mathrm{\Delta }{I}_1-{\displaystyle \frac{\mathrm{\partial }{I}_2}{\mathrm{\partial }t}}-({\mu }_2(x)+{\theta }_2(x))I_2\le 0,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega },\end{array}$and $\{\begin{array}{l}{d}_4\mathrm{\Delta }R-{\displaystyle \frac{\mathrm{\partial }R}{\mathrm{\partial }t}}-({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R\le 0,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$It follows from the strong maximum principle and Hopf boundary lemma in [28] that $I_1(t,x,\psi )>0$, $I_2(t,x,\psi )>0$, $R(t,x,\psi )>0$. Furthermore, according to the proof of Theorem 2.2, there exists a positive constant $M_7>0$ and $t_3>0$ such that $I_1(t,x)$, $I_2(t,x)$, $R(t,x)\le M_7$ for $t\ge t_3$. Then, from model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) we have $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }t}}\ge d_1\mathrm{\Delta }S+\lambda (x)-(\mu (x)+{\beta }_1(x)M_7+{\beta }_2(x)M_7)S,t\ge t_3,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }\nu }}=0,t\ge t_3,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$Considering the following comparison system: (12) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }v+\lambda (x)-(\mu (x)+{\beta }_1(x)M_7+{\beta }_2(x)M_7)v,t\ge t_3,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }v(t,x)}{\mathrm{\partial }\nu }}=0,t\ge t_3,x\in \mathrm{\partial }\mathrm{\Omega },\end{array}\end{array}$(12) it follows from the Lemma 2 in [13] that system (12(12) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }v+\lambda (x)-(\mu (x)+{\beta }_1(x)M_7+{\beta }_2(x)M_7)v,t\ge t_3,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }v(t,x)}{\mathrm{\partial }\nu }}=0,t\ge t_3,x\in \mathrm{\partial }\mathrm{\Omega },\end{array}\end{array}$(12) ) admits a unique positive steady state $v^{\ast }(x)$ which is globally asymptotically stable. Thus, the standard parabolic comparison implies that $\underset{t\to \mathrm{\infty }}{lim inf}S(t,x,\psi )\ge v^{\ast }(x)$, and hence, the proof is complete.

Theorem 3.3

For every initial value function $\psi \in \mathbb{X}$, assume model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) admits a solution $u(t,x,\psi )=(S(t,x),I_1(t,x),I_2(t,x),R(t,x))$ on $[0,\mathrm{\infty })$. If ${\mathcal{R}}_0>1$, then there exists a constant $\rho >0$ such that for any $\psi \in {\mathbb{X}}_{+}$ with ${\psi }_2\ne 0$, ${\psi }_2\ne 0$ and ${\psi }_4\ne 0$, then we have $\underset{t\to \mathrm{\infty }}{lim inf}(S(t,x),I_1(t,x),I_2(t,x),R(t,x))\ge (\rho ,\rho ,\rho ,\rho ).$Moreover, model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) admits at least one positive steady state.

Proof.

Let ${\mathbb{X}}_0=\{\psi =({\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4)\in {\mathbb{X}}_{+}:{\psi }_2(\cdot )\not\equiv 0,{\psi }_3(\cdot )\not\equiv 0\mathrm{and}{\psi }_4(\cdot )\not\equiv 0\}$and $\mathrm{\partial }{\mathbb{X}}_0:={\mathbb{X}}_{+}\mathrm{\setminus }{\mathbb{X}}_{+}=\{{\psi }_2(\cdot )=0\mathrm{or}{\psi }_3(\cdot )=0\mathrm{or}{\psi }_4(\cdot )=0\}.$For $\psi \in {\mathbb{X}}_0$, it follows from Lemma 3.2 that $I_1(t,x,\psi )>0$, $I_2(t,x,\psi )>0$ and $R(t,x,\psi )>0$. That is, $\mathrm{\Phi }(t){\mathbb{X}}_0\subseteq {\mathbb{X}}_0$. Moreover, assume that $M_{\mathrm{\partial }}=\{\psi \in \mathrm{\partial }{X}_0:\mathrm{\Phi }(t)(\psi )\in \mathrm{\partial }{X}_0,\mathrm{\forall }t\ge 0\}$ and $\omega (\psi )$ is the omega limit set of orbit ${\gamma }^{+}(\psi ):=\{\mathrm{\Phi }(t)\psi :\mathrm{\forall }t\ge 0\}$.

Claim 1. $\omega (\psi )=\{(S_0^{\ast }(x),0,0,0)\}$ for $\psi \in M_{\mathrm{\partial }}$.

Since $\psi \in M_{\mathrm{\partial }}$, then it follows that $\mathrm{\Phi }(t)\psi \in \mathrm{\partial }{\mathbb{X}}_0$, and hence, $I_1(t,\cdot ,\psi )\equiv 0$ or $I_2(t,\cdot ,\psi )\equiv 0$ or $R(t,\cdot ,\psi )\equiv 0$. For $I_1(t,\cdot ,\psi )\equiv 0$, it then follows from model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that $k_2{\theta }_2(x)I_2(t,\cdot ,\psi )+{\delta }_1(x)R(t,\cdot ,\psi )\equiv 0$, then we have $I_1(t,\cdot ,\psi )\equiv R(t,\cdot ,\psi )\equiv 0$. Furthermore, it follows from model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that the equation of S is asymptotic to $\frac{\mathrm{\partial }S}{\mathrm{\partial }t}=d_1\mathrm{\Delta }S+\lambda (x)-\mu (x)S,$which implies that $\underset{t\to \mathrm{\infty }}{lim}S(t,x,\psi )=S_0^{\ast }(x)$. If $I_1(t,\cdot ,\psi )\ne 0$ for some $\hat{t}>0$ and $x\in \overline{\mathrm{\Omega }}$, then Lemma 3.2 implies that $I_1(t,\cdot ,\psi )>0$ for $t\ge \hat{t}$ and $x\in \overline{\mathrm{\Omega }}$, and hence, $I_2(t,\cdot ,\psi )\equiv 0$ or $R(t,\cdot ,\psi )\equiv 0$. For the case of $I_2(t,\cdot ,\psi )\equiv 0$, then model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) implies that $k_1{\theta }_1(x)I_1(t,\cdot ,\psi )+{\delta }_2(x)R(t,\cdot ,\psi )\equiv 0$, and thus $I_1(t,\cdot ,\psi )\equiv R(t,\cdot ,\psi )\equiv 0$, which leads a contradiction. Furthermore, considering the case where $I_2(t,\cdot ,\psi )\ne 0$ for some $\breve{t}>0$ then $I_2(t,\cdot ,\psi )>0$ for $t\ge \breve{t}$, and thus $R(t,\cdot ,\psi )\equiv 0$. It follows from (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that $I_1(t,\cdot ,\psi )\equiv I_2(t,\cdot ,\psi )\equiv 0$, which is also a contradiction. Thus the claim is valid.

Recall that ${\mathcal{R}}_0>1$, then Lemma 2.3 implies that ${\lambda }_0>0$. By continuity, there exists a sufficient small ${\eta }_{\ast }>0$ such that ${\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}>0$, where ${\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}$ is the principle eigenvalue of (13) $\begin{array}{r}\{\begin{array}{l}{\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}{\tilde{\psi }}_2(x)=d_2\mathrm{\Delta }{\tilde{\psi }}_2(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x)-{\eta }_{\ast },0){\tilde{\psi }}_2+k_2{\theta }_2(x){\tilde{\psi }}_3\\ -({\mu }_1(x)+{\theta }_1(x)){\tilde{\psi }}_2+{\delta }_1(x){\tilde{\psi }}_4,t>0,x\in \mathrm{\Omega },\\ {\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}{\tilde{\psi }}_3(x)=d_3\mathrm{\Delta }{\tilde{\psi }}_3(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x)-{\eta }_{\ast },0){\tilde{\psi }}_3+k_1{\theta }_1(x){\tilde{\psi }}_2\\ -({\mu }_2(x)+{\theta }_2(x)){\tilde{\psi }}_3+{\delta }_2(x){\tilde{\psi }}_4,t>0,x\in \mathrm{\Omega },\\ {\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}{\tilde{\psi }}_4(x)=d_4\mathrm{\Delta }{\tilde{\psi }}_4(t,x)+(1-k_1){\theta }_1(x){\tilde{\psi }}_2(t,x)+(1-k_2){\theta }_2(x){\tilde{\psi }}_3\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x)){\tilde{\psi }}_4,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{\tilde{\psi }}_2(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\tilde{\psi }}_3(t,x)}{\mathrm{\partial }t}}={\displaystyle \frac{\mathrm{\partial }{\tilde{\psi }}_4(t,x)}{\mathrm{\partial }t}}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}\end{array}$(13) Claim 2. $\{(S_0^{\ast },0,0,0)\}$ is uniformly weak repeller in the sense that $\underset{t\to \mathrm{\infty }}{lim sup}\Vert \mathrm{\Phi }(t)\psi -(S_0^{\ast },0,0,0)\Vert \ge {\eta }^{\ast }$, $\mathrm{\forall }\psi \in {\mathbb{X}}_0$.

If it is not true, then there exists $\hat{\psi }\in {\mathbb{X}}_0$ such that $\underset{t\to \mathrm{\infty }}{lim sup}\Vert \mathrm{\Phi }(t)\hat{\psi }-(S_0^{\ast },0,0,0)\Vert <{\eta }^{\ast }$, and hence, there exists a $t^{\ast }>0$ such that for $t\ge t^{\ast }$, $S(t,\cdot ,\hat{\psi })>S_0^{\ast }-{\eta }_{\ast }$, $0<I_1(t,\cdot ,\hat{\psi })<{\eta }_{\ast }$, $0<I_2(t,\cdot ,\hat{\psi })<{\eta }_{\ast }$, $0<R(t,\cdot ,\hat{\psi })<{\eta }_{\ast }$. It follows from model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }t}}\ge d_2\mathrm{\Delta }{I}_1+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x)-{\eta }_{\ast },{\eta }_{\ast })I_1+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>t^{\ast },x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2}{\mathrm{\partial }t}}\ge d_3\mathrm{\Delta }{I}_2+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x)-{\eta }_{\ast },{\eta }_{\ast })I_2+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>t^{\ast },x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R}{\mathrm{\partial }t}}\ge d_4\mathrm{\Delta }R+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>t^{\ast },x\in \mathrm{\Omega }.\end{array}$Assume $({\tilde{\psi }}_2,{\tilde{\psi }}_3,{\tilde{\psi }}_4)$ is the strongly positive eigenfunction associated with the principle eigenvalue ${\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}>0$. Suppose $(I_1(t^{\ast },x,\hat{\psi }),I_2(t^{\ast },x,\hat{\psi }),R(t^{\ast },x,\hat{\psi }))\ge {\zeta }_{\ast }({\tilde{\psi }}_2,{\tilde{\psi }}_3,{\tilde{\psi }}_4)$ for some ${\zeta }_{\ast }>0$, and hence, $(I_1(t,x,\hat{\psi }),I_2(t,x,\hat{\psi }),R(t,x,\hat{\psi }))\ge {\zeta }_{\ast }({\tilde{\psi }}_2,{\tilde{\psi }}_3,{\tilde{\psi }}_4)e^{{\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}(t-t^{\ast })},\mathrm{\forall }t\ge t^{\ast },$where ${\zeta }_{\ast }({\tilde{\psi }}_2,{\tilde{\psi }}_3,{\tilde{\psi }}_4)e^{{\lambda }_0^{S_0^{\ast }(x)-{\eta }_{\ast }}(t-t^{\ast })}$ is a solution of the following linear system: $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x)-{\eta }_{\ast },{\eta }_{\ast })I_1+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>t^{\ast },x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x)-{\eta }_{\ast },{\eta }_{\ast })I_2+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>t^{\ast },x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>t^{\ast },x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }\nu }}={\displaystyle \frac{\mathrm{\partial }{I}_2}{\mathrm{\partial }\nu }}={\displaystyle \frac{\mathrm{\partial }R}{\mathrm{\partial }\nu }}=0,t>t^{\ast },x\in \mathrm{\partial }\mathrm{\Omega }.\end{array}$Then we have $\underset{t\to \mathrm{\infty }}{lim}{I}_1(t,x)=\mathrm{\infty }$, $\underset{t\to \mathrm{\infty }}{lim}{I}_2(t,x)=\mathrm{\infty }$, $\underset{t\to \mathrm{\infty }}{lim}R(t,x)=\mathrm{\infty }$, this arrives a contradiction.

Define a function $p:{\mathbb{X}}_{+}\to [0,\mathrm{\infty })$ by $p(\psi )=min\{\underset{x\in \overline{\mathrm{\Omega }}}{min}{\psi }_2(0,x),\underset{x\in \overline{\mathrm{\Omega }}}{min}{\psi }_3(0,x),\underset{x\in \overline{\mathrm{\Omega }}}{min}{\psi }_4(0,x)\},\psi \in {\mathbb{X}}_{+}.$It follows that $p^{-1}(0,+\mathrm{\infty })\subseteq {\mathbb{X}}_0$ and have the property that if either $p(\psi )>0$ or $p(\psi )=0$ and $\psi \in {\mathbb{X}}_0$, then $p(\mathrm{\Phi }(t)\psi )>0$. This implies p is a generalized distance function for the semiflow $\mathrm{\Phi }(t):{\mathbb{X}}_{+}\to {\mathbb{X}}_{+}$ (see [32]). Moreover, it follows from the above discussion that any forward orbit of $\mathrm{\Phi }(t)$ in ${\mathcal{M}}_{\mathrm{\partial }}$ converges to $\{(S_0^{\ast }(x),0,0,0)\}$, thus $\{(S_0^{\ast }(x),0,0,0)\}$ is isolated in ${\mathbb{X}}_{+}$ and $W^s(S_0^{\ast }(x),0,0,0)\cap {\mathbb{X}}_0=\mathrm{\varnothing }$, where $W^s(S_0^{\ast }(x),0,0,0)$ is the stable set of $\{(S_0^{\ast }(x),0,0,0)\}$. Thus there is no cycle in ${\mathcal{M}}_{\mathrm{\partial }}$ from $\{(S_0^{\ast }(x),0,0,0)\}$ to $\{(S_0^{\ast }(x),0,0,0)\}$. It follows from Theorem 2.1 that the semiflow $\mathrm{\Phi }(t):{\mathbb{X}}_{+}\to {\mathbb{X}}_{+}$ has a global compact attractor in ${\mathbb{X}}_{+}$. By [32], there exists a $\hat{\eta }>0$ such that $\underset{\psi \in \omega (\psi )}{min}p(\psi )>\hat{\eta }$ for $\psi \in {\mathbb{X}}_0$, which implies that $\underset{t\to \mathrm{\infty }}{lim inf}{I}_1(t,\cdot ,\psi )\ge \hat{\eta }$, $\underset{t\to \mathrm{\infty }}{lim inf}{I}_2(t,\cdot ,\psi )\ge \hat{\eta }$, $\underset{t\to \mathrm{\infty }}{lim inf}R(t,\cdot ,\psi )\ge \hat{\eta }$. Accordingly, choosing $0<\rho <min\{{\eta }_{\ast },\hat{\eta }\}$, and by Lemma 3.2, it then follows that $\underset{t\to \mathrm{\infty }}{lim inf}{I}_1(t,\cdot ,\psi )\ge \rho$, $\underset{t\to \mathrm{\infty }}{lim inf}{I}_2(t,\cdot ,\psi )\ge \rho$, $\underset{t\to \mathrm{\infty }}{lim inf}R(t,\cdot ,\psi )\ge \rho$. Furthermore, it follows from Lemma 3.2 and [25] that $\mathrm{\Phi }(t)$ admits at least a steady state in ${\mathbb{X}}_0$, which is a positive steady state of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ). This completes the proof.

Furthermore, motivated by the idea of [4,26,43], we investigate the global stability of $E_0(x)$ in the critical case of ${\mathcal{R}}_0=1$.

Theorem 3.4

If ${\mathcal{R}}_0=1$, and $f_{I_1}(S,0)$ and $g_{I_2}(S,0)$ are Lipschitz continuous, then $E_0(x)$ is globally asymptotically stable.

Proof.

First, we prove the local stability of $E_0(x)$ of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ). Assume $\hat{\epsilon }>0$ and set $u_0=(S_0,I_{10},I_{20},R_0)$ with $\Vert u_0-(S_0^{\ast }(x),0,0,0)\Vert \le \sigma$ for a small $\sigma >0$. Define $U(t,x)=\frac{S(t,x)}{S_0^{\ast }(x)}-1\mathrm{and}\chi (t)=\underset{x\overline{\mathrm{\Omega }}}{max}\{U(t,x),0\}.$Then recall that $D_1\mathrm{\Delta }{S}_0^{\ast }(x)+\lambda (x)-\mu (x)S_0^{\ast }(x)=0$, and hence, we have $U_t-d_1\mathrm{\Delta }U-2d_1\frac{\mathrm{\nabla }{S}_0^{\ast }(x)\mathrm{\nabla }U}{S_0^{\ast }(x)}+\frac{\lambda (x)}{S_0^{\ast }(x)}U=-\frac{{\beta }_1(x)f(S,I_1)+{\beta }_2(x)g(S,I_2)}{S_0^{\ast }(x)}.$Assume the operator $D_1\mathrm{\Delta }+2D_1\frac{\mathrm{\nabla }{S}_0^{\ast }(x)\mathrm{\nabla }}{S_0^{\ast }(x)}-\frac{\lambda (x)}{S_0^{\ast }(x)}$ admits a positive semigroup ${\tilde{\mathrm{\Gamma }}}_1(t)$, then it follows from [16] that $\Vert {\tilde{\mathrm{\Gamma }}}_1\Vert \le L_1{e}^{-qt}$ for some $L_1>0,q>0$. Therefore, $U(t,x)$ can be rewritten by (14) $\begin{array}{rl}U(t,x)& ={\tilde{\mathrm{\Gamma }}}_1(t)U(0,x)\\ & -{\int }_0^t{\tilde{\mathrm{\Gamma }}}_1(t-s)\frac{{\beta }_1(x)f(S(s,x),I_1(s,x))+{\beta }_2(x)g(S(s,x),I_2(s,x))}{S_0^{\ast }(x)}\mathrm{d}s,\end{array}$(14) where $U(0,x)=\frac{S(0,x)}{S_0^{\ast }(x)}-1$.

Then we have $\chi (t)\le \underset{x\in \overline{\mathrm{\Omega }}}{max}\{{\tilde{\mathrm{\Gamma }}}_1(t)U(0,x),0\}\le \Vert {\tilde{\mathrm{\Gamma }}}_1(t)U(0,x)\Vert \le \frac{L_1\sigma e^{-qt}}{S_{min}},$where $S_{min}=\underset{x\in \overline{\mathrm{\Omega }}}{min}\{S_0^{\ast }(x)\}$. Moreover, let ${\tilde{\mathrm{\Gamma }}}_2(t)$ be the semigroup for the following system: $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0)I_1(t,x)+{\delta }_1(x)R(t,x)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1(t,x),\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0)I_2(t,x)+k_1{\theta }_1(x)I_1(x)+{\delta }_2(x)R(t,x)\\ -({\mu }_2(x)+{\theta }_2(x))I_2(t,x),\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R(t,x).\end{array}$It then follows from [40] that there exists a positive constant ${\tilde{L}}_1$ such that $\Vert {\tilde{\mathrm{\Gamma }}}_2\Vert \le {\tilde{L}}_1$ for $t\ge 0$, and hence, we have $\begin{array}{rl}\left(\begin{array}{c}{I}_1(t,x)\\ I_2(t,x)\\ R(t,x)\end{array}\right)& ={\tilde{\mathrm{\Gamma }}}_2(t)\left(\begin{array}{c}{I}_1(0,x)\\ I_2(0,x)\\ R(0,x)\end{array}\right)\\ & +{\int }_0^t{\tilde{\mathrm{\Gamma }}}_2(t-s)\left(\begin{array}{c}{\beta }_1(x)f(S(s,x),I_1(s,x))-{\beta }_1(x)f_{I_1}(S_0^{\ast }(x),0)I_1(s,x)\\ {\beta }_2(x)g(S(s,x),I_2(s,x))-{\beta }_2(x)g_{I_2}(S_0^{\ast }(x),0)I_2(s,x)\\ 0\end{array}\right)\mathrm{d}s\\ & \le {\tilde{\mathrm{\Gamma }}}_2(t)\left(\begin{array}{c}{I}_1(0,x)\\ I_2(0,x)\\ R(0,x)\end{array}\right)\\ & +{\int }_0^t{\tilde{\mathrm{\Gamma }}}_2(t-s)\left(\begin{array}{c}{\beta }_1(x)[f_{I_1}(S(s,x),0)-f_{I_1}(S_0^{\ast }(x),0)]I_1(s,x)\\ {\beta }_2(x)[g_{I_2}(S(s,x),0)-g_{I_2}(S_0^{\ast }(x),0)]I_2(s,x)\\ 0\end{array}\right)\mathrm{d}s\\ & \le {\tilde{\mathrm{\Gamma }}}_2(t)\left(\begin{array}{c}{I}_1(0,x)\\ I_2(0,x)\\ R(0,x)\end{array}\right)+{\int }_0^t{\tilde{\mathrm{\Gamma }}}_2(t-s)\left(\begin{array}{c}{\beta }_1(x)Q_1|S(s,x)-S_0^{\ast }(x)|I_1(s,x)\\ {\beta }_2(x)Q_2|S(s,x)-S_0^{\ast }(x)|I_2(s,x)\\ 0\end{array}\right)\mathrm{d}s,\end{array}$where $Q_1$ and $Q_2$ are Lipschitz constants for $f_{I_1}$ and $g_{I_2}$, respectively. Denote $\beta =\underset{x\in \overline{\mathrm{\Omega }}}{max}\{{\beta }_1(x),{\beta }_2(x)\},Q=max\{Q_1,Q_2\}.$Then, we have $\begin{array}{rl}& max\{\Vert I_1(t,x)\Vert ,\Vert I_2(t,x)\Vert ,\Vert R(t,x)\Vert \}\\ & \le {\tilde{L}}_1\sigma +{\tilde{L}}_1\beta Q\Vert S_0^{\ast }(x)\Vert \frac{L_1\sigma }{S_{min}}{\int }_0^t{e}^{-qs}(\Vert I_1(s,x)\Vert +\Vert I_2(s,x)+R(s,x)\Vert )\mathrm{d}s\\ & ={\tilde{L}}_1\sigma +{\tilde{L}}_2\sigma {\int }_0^t{e}^{-qs}(\Vert I_1(s,x)\Vert +\Vert I_2(s,x)+R(s,x)\Vert )\mathrm{d}s,\end{array}$where ${\tilde{L}}_2=\frac{{\tilde{L}}_1{L}_1\beta Q\Vert S_0^{\ast }(x)\Vert }{S_{min}}$, thus we have $\begin{array}{rl}& \Vert I_1(t,x)\Vert +\Vert I_2(t,x)\Vert +\Vert R(t,x)\Vert \\ & \le 3{\tilde{L}}_1\sigma +3{\tilde{L}}_2\sigma {\int }_0^t{e}^{-qs}(\Vert I_1(s,x)\Vert +\Vert I_2(s,x)+R(s,x)\Vert )\mathrm{d}s,\end{array}$according to Gronwall's inequality, it follows that $\Vert I_1(t,x)\Vert +\Vert I_2(t,x)\Vert +\Vert R(t,x)\Vert \le 3{\tilde{L}}_1\sigma e^{\frac{3{\tilde{L}}_2\sigma }{q}}.$Thus, from the first equation of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ), we can obtain $\frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}\ge d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-3{\tilde{L}}_1\sigma e^{\frac{3\sigma {\tilde{L}}_2}{q}}({\beta }_1(x)+{\beta }_2(x))S(t,x).$Assume $\hat{u}(t,x)$ is the solution of the following system: (15) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }\hat{u}(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }\hat{u}(t,x)+\lambda (x)-\mu (x)\hat{u}(t,x)-3{\tilde{L}}_1\sigma e^{{\displaystyle \frac{3\sigma {\tilde{L}}_2}{q}}}({\beta }_1(x)\\ +{\beta }_2(x))\hat{u}(t,x),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\hat{u}}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega },\\ \hat{u}(0,x)=S_0(x),x\in \mathrm{\Omega }.\end{array}\end{array}$(15) Then, the comparison principle implies that $S(t,x)>\hat{u}(t,x)$ for $t\ge 0$, $x\in \overline{\mathrm{\Omega }}$. Let $S_{\sigma }$ be the positive steady state of (15(15) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }\hat{u}(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }\hat{u}(t,x)+\lambda (x)-\mu (x)\hat{u}(t,x)-3{\tilde{L}}_1\sigma e^{{\displaystyle \frac{3\sigma {\tilde{L}}_2}{q}}}({\beta }_1(x)\\ +{\beta }_2(x))\hat{u}(t,x),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\hat{u}}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega },\\ \hat{u}(0,x)=S_0(x),x\in \mathrm{\Omega }.\end{array}\end{array}$(15) ) and $\hat{S}(t,x)=\hat{u}(t,x)-S_{\sigma }(x)$, then $\hat{S}(t,x)$ satisfies $\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }\hat{S}(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }\hat{S}(t,x)+\lambda (x)-(\mu (x)+3{\tilde{L}}_1\sigma e^{{\displaystyle \frac{3\sigma {\tilde{L}}_2}{q}}}({\beta }_1(x)\\ +{\beta }_2(x)))\hat{S}(t,x),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\hat{S}}{\mathrm{\partial }\nu }}=0,t>0,x\in \mathrm{\partial }\mathrm{\Omega },\\ \hat{S}(0,x)=\hat{u}(0,x)-S_{\sigma }(x)=S_0(x)-S_{\sigma }(x),x\in \mathrm{\Omega }.\end{array}$It then follows that $\hat{S}(t,x)={\mathbb{T}}_1(t)\hat{S}(0,x)-{\int }_0^t{\mathbb{T}}_1(t-s)3{\tilde{L}}_1\sigma e^{\frac{3\sigma {\tilde{L}}_2}{q}}({\beta }_1(x)+{\beta }_2(x))\hat{S}(s,x)\mathrm{d}s,$which implies that $\hat{S}(t,x)\le M{e}^{{\omega }_{1}t}\Vert S_0(x)-S_{\sigma }(x)\Vert +{\tilde{L}}_3{\int }_0^t{e}^{{\omega }_1(t-s)}\Vert \hat{S}(s,x)\Vert \mathrm{d}s,$where ${\tilde{L}}_3=3M{\tilde{L}}_1\sigma e^{\frac{3\sigma {\tilde{L}}_2}{q}}(\Vert {\beta }_1(x)\Vert +\Vert {\beta }_2(x)\Vert )$. By Gronwall's inequality, we have $\hat{S}(t,x)=\Vert \hat{u}(t,x)-S_{\sigma }(x)\Vert \le M\Vert S_0(x)-S_{\sigma }(x)\Vert e^{{\tilde{L}}_{3}t+{\omega }_{1}t}.$By selecting $\sigma >0$ sufficiently small such that ${\tilde{L}}_3<-\frac{{\omega }_1}{2}$, then we have (16) $\begin{array}{r}\Vert \hat{u}(t,x)-S_{\sigma }(x)\Vert \le M\Vert S_0(x)-S_{\sigma }(x)\Vert e^{\frac{{\omega }_{1}t}{2}}.\end{array}$(16) It follows that $\begin{array}{rl}S(t,x)-S_0^{\ast }(x)& \ge \hat{u}(t,x)-S_0^{\ast }(x)=\hat{u}(t,x)-S_{\sigma }(x)+S_{\sigma }(x)-S_0^{\ast }(x)\\ & \ge -M\Vert S_0^{\ast }(x)-S_{\sigma }(x)\Vert e^{\frac{{\omega }_{1}t}{2}}+S_{\sigma }(x)-S_0^{\ast }(x)\\ & \ge -M(\Vert S_0(x)-S_0^{\ast }(x)\Vert +\Vert S_0^{\ast }(x)-S_{\sigma }(x)\Vert )-\Vert S_{\sigma }(x)-S_0^{\ast }(x)\Vert \\ & \ge -M\sigma -(M+1)\Vert S_{\sigma }(x)-S_0^{\ast }(x)\Vert .\end{array}$Recall that $\chi (t)\le \frac{L_1\sigma e^{-qt}}{S_{min}}$, then we have $S(t,x)-S_0^{\ast }(x)\le \Vert S_0^{\ast }(x)\Vert \frac{L_1\sigma }{S_{min}}$, and hence, $\Vert S(t,x)-S_0^{\ast }(x)\Vert \le max\{\Vert S_0^{\ast }(x)\Vert \frac{L_1\sigma }{S_{min}},M\sigma +(M+1)\Vert S_{\sigma }(x)-S_0^{\ast }(x)\Vert \}.$Since $\underset{\sigma \to 0}{lim}{S}_{\sigma }(x)=S_0^{\ast }(x)$, then we can selecting a sufficiently small σ such that $\Vert S(t,x)-S_0^{\ast }(x)\Vert ,\Vert I_1(t,x)\Vert ,\Vert I_2(t,x)\Vert ,\Vert R(t,x)\Vert \le \hat{\epsilon },\mathrm{\forall }t>0,$which implies the local asymptotic stability of $(S_0^{\ast }(x),0,0,0)$.

Next, we prove the $(S_0^{\ast }(x),0,0,0)$ is globally attractive. According to Theorem 2.2, ${\mathrm{\Phi }}_t$ has a global attractor ${\mathcal{M}}_0$. Define $\mathrm{\partial }{\mathbb{X}}_1=\{(S,I_1,I_2,R)\in {\mathbb{X}}_{+}:I_1=I_2=R=0\}.$Claim 1. $u_0=(S_0,I_{10},I_{20},R_0)\in {\mathcal{M}}_0$, the omega limit set $\omega (u_0)\subset \mathrm{\partial }{\mathbb{X}}_1$.

Similar to [4,26,36,43,45], we define $c(t,u_0):=inf\{\tilde{c}\in \mathbb{R}:I_1(t,\cdot )\le \tilde{c}{\psi }_2,I_2(t,\cdot )\le \tilde{c}{\psi }_3,R(t,\cdot )\le \tilde{c}{\psi }_4\},$and hence, $c(t,u_0)>0$ for t>0. It claims that $c(t,u_0)$ is strictly decreasing. To this end, fixed a $t_0>0$ and let ${\tilde{I}}_1(t,\cdot )=c(t_0,u_0){\psi }_2$, ${\tilde{I}}_2(t,\cdot )=c(t_0,u_0){\psi }_3$ and $\tilde{R}(t,\cdot )=c(t_0,u_0){\psi }_4$ for $t\ge t_0$, such that $({\tilde{I}}_1(t,\cdot ),{\tilde{I}}_2(t,\cdot ),\tilde{R}(t,\cdot ))\ge (I_1(t,\cdot ),I_2(t,\cdot ),R(t,\cdot ))$ for $(x,t)\in \overline{\mathrm{\Omega }}\times [t_0,\mathrm{\infty })$, where $({\tilde{I}}_1(t,\cdot ),{\tilde{I}}_2(t,\cdot ),\tilde{R}(t,\cdot ))$ satisfies (17) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{\tilde{I}}_1}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{\tilde{I}}_1+{\beta }_1(x)S{\tilde{I}}_1+k_2{\theta }_2(x){\tilde{I}}_2+{\delta }_1(x)\tilde{R}-({\mu }_1(x)+{\theta }_1(x)){\tilde{I}}_1,t>t_0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{\tilde{I}}_2}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{\tilde{I}}_2+{\beta }_2(x)S{\tilde{I}}_2+k_1{\theta }_1(x){\tilde{I}}_1+{\delta }_2(x)\tilde{R}-({\mu }_2(x)+{\theta }_2(x)){\tilde{I}}_2,t>t_0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }\tilde{R}}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }\tilde{R}+(1-k_1){\theta }_1(x){\tilde{I}}_1+(1-k_2){\theta }_2(x){\tilde{I}}_2-({\mu }_3(x)\\ +{\delta }_1(x)+{\delta }_2(x))\tilde{R},t>t_0,x\in \mathrm{\Omega },\\ {\tilde{I}}_1(t_0,x)\ge I_1(t_0,x),{\tilde{I}}_2(t_0,x)\ge I_2(t_0,x),\tilde{R}(t_0,x)\ge R(t_0,x).\end{array}\end{array}$(17) It follows from the comparison principle that $c(t_0,u_0){\psi }_2={\tilde{I}}_1(t,x)>I_1(t,x)$, $c(t_0,u_0){\psi }_3=\tilde{I}(t,x)>I_2(t,x)$ and $c(t_0,u_0){\psi }_4=\tilde{R}(t,x)>R(t,x)$ for $t\ge t_0$, $x\in \overline{\mathrm{\Omega }}$. Due to the arbitraryness of $t_0>0$, thus $c(t,u_0)$ is strictly decreasing.

Denote $c_{\ast }=\underset{t\to \mathrm{\infty }}{lim}c(t,u_0)$ and $v=(v_1,v_2,v_3,v_4)\in \omega (u_0)$, then there exists a sequence $\{t_k\}$ with $t_k\to \mathrm{\infty }$ such that $\mathrm{\Phi }(t_k)u_0\to v$. Since, $\underset{t_k\to \mathrm{\infty }}{lim}\mathrm{\Phi }(t+t_k)u_0=\mathrm{\Phi }(t)\underset{t_k\to \mathrm{\infty }}{lim}\mathrm{\Phi }(t_k)u_0=\mathrm{\Phi }(t)v,$which implies $c(t,v)=c_{\ast }$ for $t\ge 0$. If $v_2\ne 0$ or $v_3\ne 0$ or $v_4\ne 0$, a similar argument procedures to the above leads to that $c(t,v)$ is strictly decreasing. This is a contradiction and hence $v_2=v_3=v_4=0$, then we have $(I_1(t,x),I_2(t,x),R(t,x))\to (0,0,0)$. Consequently, it follows from the first equation of model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ) that $S(t,x)\to S_0^{\ast }(x)$ as $t\to \mathrm{\infty }$.

Claim 2. ${\mathcal{M}}_0=\{(S_0^{\ast }(x),0,0,0)\}$.

It follows from the discussions above that $\{(S_0^{\ast }(x),0,0,0)\}$ is globally attractive in $\mathrm{\partial }{\mathbb{X}}_1$, and $\{(S_0^{\ast }(x),0,0,0)\}$ forms the only compact invariant subset in $\mathrm{\partial }{\mathbb{X}}_1$. Since the omega limits set $\omega (u_0)$ is compact invariant and $\omega (u_0)\subset \mathrm{\partial }{\mathbb{X}}_1$ for $u_0\in {\mathcal{M}}_0$, thus $\omega (u_0)=\{(S_0^{\ast }(x),0,0,0)\}$. Since the global attractor ${\mathcal{M}}_0$ is compact invariant, it follows from Lemma 3.11 in [43] that ${\mathcal{M}}_0=\{(S_0^{\ast }(x),0,0,0)\}$. Global attractivity and local asymptotic stability immediately lead to the global asymptotically stability of $\{(S_0^{\ast }(x),0,0,0)\}$. This completes the proof.

4. Conclusion

In this paper, a diffusion epidemic model for heroin and cocaine was formulated. The basic reproduction number was defined by using the method of next generation operators. The threshold type dynamics of the model in terms of ${\mathcal{R}}_0$ have performed by applying the comparison arguments and persistence theory: if ${\mathcal{R}}_0<1$, then the drug-free steady state is globally asymptotically stable, while the drug spread is persistent if ${\mathcal{R}}_0>1$. Moreover, by means of Gronwall's inequality, comparison theorem and the properties of semigroup, the global stability of the drug-free steady state in the critical case ${\mathcal{R}}_0=1$ was investigated.

Though the spatial heterogeneity was introduced in model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }S(t,x)}{\mathrm{\partial }t}}=d_1\mathrm{\Delta }S(t,x)+\lambda (x)-\mu (x)S(t,x)-{\beta }_1(x)f(S,I_1)\\ -{\beta }_2(x)g(S,I_2),t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_1(t,x)}{\mathrm{\partial }t}}=d_2\mathrm{\Delta }{I}_1(t,x)+{\beta }_1(x)f(S,I_1)+k_2{\theta }_2(x)I_2\\ -({\mu }_1(x)+{\theta }_1(x))I_1+{\delta }_1(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }{I}_2(t,x)}{\mathrm{\partial }t}}=d_3\mathrm{\Delta }{I}_2(t,x)+{\beta }_2(x)g(S,I_2)+k_1{\theta }_1(x)I_1\\ -({\mu }_2(x)+{\theta }_2(x))I_2+{\delta }_2(x)R,t>0,x\in \mathrm{\Omega },\\ {\displaystyle \frac{\mathrm{\partial }R(t,x)}{\mathrm{\partial }t}}=d_4\mathrm{\Delta }R(t,x)+(1-k_1){\theta }_1(x)I_1(t,x)+(1-k_2){\theta }_2(x)I_2\\ -({\mu }_3(x)+{\delta }_1(x)+{\delta }_2(x))R,t>0,x\in \mathrm{\Omega },\end{array}\end{array}$(2) ), the age structure was not been considered in the model. As mentioned in [3,8,10,11,19], treat age or susceptibility age on drug abuse is also should be taken into consideration. Thus a model with age structure and spatial heterogeneity is an interesting project, we leave this in future study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China (#11701445, #11971379, #11801439), by Natural Science Basic Research Program in Shaanxi Province of China (2022JM-042, 2020JQ-831, 2021JM-320).