# A two-strain reaction-diffusion malaria model with seasonality and vector-
bias††thanks:  This research is supported by the NSF of China (No. 11971369)
and the Fundamental Research Funds for the Central Universities (No.
JB210711).

Huijie Chu, Zhenguo Bai
School of Mathematics and Statistics, Xidian University, Xi’an, 710126, China
E-mail<EMAIL_ADDRESS>
###### Abstract

To investigate the combined effects of drug resistance, seasonality and
vector-bias, we formulate a periodic two-strain reaction-diffusion model. It
is a competitive system for resistant and sensitive strains, but the single-
strain subsystem is cooperative. We derive the basic reproduction number
$\mathcal{R}_{i}$ and the invasion reproduction number $\mathcal{\hat{R}}_{i}$
for strain $i~{}(i=1,2)$, and establish the transmission dynamics in terms of
these four quantities. More precisely, (i) if $\mathcal{R}_{1}<1$ and
$\mathcal{R}_{2}<1$, then the disease is extinct; (ii) if
$\mathcal{R}_{1}>1>\mathcal{R}_{2}$ ($\mathcal{R}_{2}>1>\mathcal{R}_{1}$),
then the sensitive (resistant) strains are persistent, while the resistant
(sensitive) strains die out; (iii) if $\mathcal{R}_{i}>1$ and
$\mathcal{\hat{R}}_{i}>1~{}(i=1,2)$, then two strains are coexistent and
periodic oscillation phenomenon is observed. We also study the asymptotic
behavior of the basic reproduction number with respect to small and large
diffusion coefficients. Numerically, we demonstrate the phenomena of
coexistence and competitive exclusion for two strains and explore the
influences of seasonality and vector-bias on disease spreading.

Key words: Malaria model; Seasonality; Vector-bias; Two strains; Reproduction
numbers.

AMS Subject Classification:  92D30, 37N25, 34K13

##  1 Introduction

Malaria, one of the most common vector-borne diseases, is endemic in over 100
countries worldwide and causes serious public health problems and a
significant economic burden worldwide [1]. Human malaria infection is caused
by the genus Plasmodium parasite, which can be transmitted to humans by the
effective bites of adult female Anopheles mosquitoes (after taking a blood
meal from humans) [2]. According to the 2020 WHO report [3], the global tally
of malaria cases was 229 million in 2019, claiming some 409 000 lives compared
to 411 000 in 2018. Therefore, a deep understanding of malaria transmission
mechanisms will undoubtedly contribute to disease control.

Mathematical models have been proposed to study the dynamics of malaria
outbreaks in different parts of the world, the earliest model dates back to
the Ross-Macdonald model [4, 5]. Since then, various mathematical models have
been designed to describe and predict the spreading of malaria (see, e.g., [6,
7, 8, 9, 10, 11, 12, 13, 14]). However, few studies consider the following
three biological factors for malaria transmission simultaneously.

Vector-bias effect. The vector-bias describes that mosquitoes prefer biting
infectious humans to susceptible ones. Kingsolver [15] first introduced a
vector-bias model for the dynamics of malarial transmission. Following
Kingsolver’s work, Hosack et al. [16] included the incubation time in
mosquitoes to study the dynamics of the disease concerning the reproduction
number. Further, Chamchod and Britton [7] extended the model from previous
authors by defining the attractiveness in a different way. Motivated by these
works, Wang and Zhao incorporated the seasonality into a vector-bias model
with incubation period [11]. Bai et al. formulated a time-delayed periodic
reaction-diffusion model with vector-bias effect [12] and found that the
ignorance of the vector-bias effect will underestimate the infection risk. All
these results show that the vector-bias has an important impact on the
epidemiology of malaria.

Drug-resistance. Currently, due to the lack of effective and safe vaccine, the
main strategy in controlling malaria is drugs. However, the use of anti-
malarial drugs such as chloroquine, malaraquine, nivaquine, aralen and
fansidar results in the appearance and spread of resistance in the parasite
population [2, 17, 18]. This poses a significant challenge to the global
control of malaria transmission or eradication of the disease. Therefore, it
is essential to investigate the resistance in malaria transmission.

Seasonality. It is generally believed that climatic factors such as
temperature, rainfall, humidity, wind, and duration of daylight greatly
influence the transmission and distribution of vector-borne diseases [19, 20,
21]. For example, rising temperatures will reduce the number of days required
for breeding, and thereby increase mosquito development rates [22]. There have
been some mathematical models and field observations suggesting that the
strength and mechanisms of seasonality can change the pattern of infectious
diseases [22, 23]. These results are beneficial for forecasting the mosquito
abundance and further effectively controlling the disease.

Except these considerations above, human and vector populations have also
contributed to the spread of vector-borne diseases [6, 9]. Therefore, this
paper will investigate a periodic two-strain malaria model with diffusion,
which is an extension of autonomous limiting system in [24]. In view of the
intrinsic mathematical structure of the model, we choose a time-varying phase
space to carry out dynamical analysis. This idea has also been used in [25].
In particular, we prove that no subset forms a cycle on the boundary with the
aim of using uniform persistence theory. Its proof is nontrivial (see Theorem
4.3).

The rest of this paper is organized as follows. In the next section, we
formulate the model and study its well-posedness. In Section 3, we define the
basic reproduction number $\mathcal{R}_{i}$ and the invasion reproduction
number $\mathcal{\hat{R}}_{i}~{}(i=1,2)$ for the sensitive and resistant
strains, respectively. In Section 4, we investigate the uniform persistence
and extinction in terms of the reproduction numbers. In Section 5, we analyze
the asymptotic behavior of the basic reproduction number concerning small and
large diffusion coefficients. In Section 6, we conduct numerical study for our
model. And the paper ends with a brief discussion.

##  2 Model formulation

Motivated by [12, 24], we consider the model with no immunity; that is,
individuals who recovered from malaria cannot resist reinfection of the
disease and can become susceptible directly. We assume that no susceptible
individual or mosquito can be infected by two virus strains. The total human
population $N_{h}(t,x)$ is divided into three groups: susceptible
$S_{h}(t,x)$, infected individuals with drug sensitive strain $I_{1}(t,x)$ and
infected individuals with drug resistant strain $I_{2}(t,x)$. For the vector
population, only adult female mosquitoes can contract the virus due to adult
males and immature mosquitoes do not take blood. Thereby, we consider only
adult female mosquitoes in our model. The vector population $M(t,x)$ has the
epidemiological classes denoted by $S_{v}(t,x)$, $I_{v1}(t,x)$ and
$I_{v2}(t,x)$ for the susceptible, infected with sensitive and resistant
strains, respectively.

Assume that all populations remain confined to a bounded domain
$\Omega\subset\mathbb{R}^{m}(m\geq 1)$ with smooth boundary $\partial\Omega$
(when $m\geq 1$). Following the line in [12], we suppose that the density of
total human population $N_{h}(t,x)=S_{h}(t,x)+I_{1}(t,x)+I_{2}(t,x)$ satisfies
the following reaction-diffusion equation:

$\left\\{\begin{split}\frac{\partial N_{h}(t,x)}{\partial t}&=D_{h}\Delta
N_{h}(t,x)+B(x,N_{h})N_{h}(t,x)-dN_{h}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial
N_{h}(t,x)}{\partial\nu}&=0,~{}~{}~{}&&t>0,~{}x\in\partial\Omega,\end{split}\right.$
(2.1)

with

$\begin{array}[]{l}B(x,u)=\left\\{\begin{array}[]{ll}b\left[1-\frac{u}{K(x)}\right],&0\leq
u\leq K(x),~{}x\in\bar{\Omega},\\\
0,&u>K(x),~{}x\in\bar{\Omega},\end{array}\right.\end{array}$

where $\Delta$ is the usual Laplacian operator. $D_{h}>0$ is the diffusion
coefficient of humans, $b$ and $d~{}(0<d<b)$ are respectively the maximal
birth rate and the nature mortality rate of humans, and $K(x)$ denotes the
local carrying capacity, which is supposed to be a positive continuous
function of location $x$. By employing [26, Theorems 3.1.5 and 3.1.6], we
arrive at that system (2.1) admits a globally attractive positive steady state
$N(x)$ in $C(\bar{\Omega},\mathbb{R}_{+})\setminus\\{0\\}$.

We also assume that the equation of the total mosquito population
$M(t,x)=S_{v}(t,x)+I_{v1}(t,x)+I_{v2}(t,x)$ is of the form:

$\left\\{\begin{split}\frac{\partial M(t,x)}{\partial t}&=D_{v}\Delta
M(t,x)+\Lambda(t,x)-\eta(t,x)M(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial
M(t,x)}{\partial\nu}&=0,~{}~{}~{}&&t>0,~{}x\in\partial\Omega,\end{split}\right.$
(2.2)

where $D_{v}>0$ is the diffusion coefficient of mosquitoes, $\Lambda(t,x)$ is
the recruitment rate at which adult female mosquitoes emerge from larval at
time $t$ and location $x$, and $\eta(t,x)$ is the natural death rate of
mosquitoes at time $t$ and location $x$. Functions $\Lambda(t,x)$ and
$\eta(t,x)$ are Hölder continuous and nonnegative nontrivial on
$\mathbb{R}\times\bar{\Omega}$, and $\omega$-periodic in $t$ for some
$\omega>0$. It easily follows that system (2.2) admits a globally stable
positive $\omega$-periodic solution $M^{*}(t,x)$ in
$C(\bar{\Omega},\mathbb{R}_{+})$ (see, e.g., [27, Lemma 2.1]). Biologically,
we may suppose that the total human and mosquito density at time $t$ and
location $x$ respectively stabilize at $N(x)$ and $M^{*}(t,x)$, that is,
$N_{h}(t,x)\equiv N(x)$ and $M(t,x)\equiv M^{*}(t,x)$ for all $t\geq 0$ and
$x\in\Omega$.

For model parameters, since the impact of climate change on mosquitoes
activities is much more than that on humans, the parameters corresponding to
mosquitoes are assumed to be time-dependent. To incorporate a vector-bias term
into the model, we use the parameters $p$ and $l$ to describe the
probabilities that a mosquito arrives at a human at random and picks the human
if he is infectious and susceptible, respectively [7, 11]. Since infectious
humans are more attractive to mosquitoes, we assume $p\geq l>0$. Let
$\beta(t,x)$ be the biting rate of mosquitoes at time $t$ and location $x$;
$c_{1}$($\alpha_{1}$) be the transmission probability per bite from infectious
mosquitoes (humans) with sensitive strain to susceptible humans (mosquitoes),
and $c_{2}$($\alpha_{2}$) be the transmission probability per bite from
infectious mosquitoes (humans) with resistant strain to susceptible humans
(mosquitoes). According to the induction in [24], we obtain

$\displaystyle J_{1}(I_{1}(t,x),I_{2}(t,x)):=$
$\displaystyle\frac{c_{1}\beta(t,x)l(N(x)-I_{1}(t,x)-I_{2}(t,x))}{p(I_{1}(t,x)+I_{2}(t,x))+l(N(x)-I_{1}(t,x)-I_{2}(t,x))},$
$\displaystyle J_{2}(I_{1}(t,x),I_{2}(t,x)):=$
$\displaystyle\frac{\alpha_{1}\beta(t,x)pI_{1}(t,x)}{p(I_{1}(t,x)+I_{2}(t,x))+l(N(x)-I_{1}(t,x)-I_{2}(t,x))},$
$\displaystyle J_{3}(I_{1}(t,x),I_{2}(t,x)):=$
$\displaystyle\frac{c_{2}\beta(t,x)l(N(x)-I_{1}(t,x)-I_{2}(t,x))}{p(I_{1}(t,x)+I_{2}(t,x))+l(N(x)-I_{1}(t,x)-I_{2}(t,x))},$
$\displaystyle J_{4}(I_{1}(t,x),I_{2}(t,x)):=$
$\displaystyle\frac{\alpha_{2}\beta(t,x)pI_{2}(t,x)}{p(I_{1}(t,x)+I_{2}(t,x))+l(N(x)-I_{1}(t,x)-I_{2}(t,x))},$

where $J_{1}(J_{3})$ represents the number of newly infectious humans with
sensitive (resistant) strain caused by an infected mosquito with sensitive
(resistant) strain per unit time at time $t$ and location $x$; and
$J_{2}$($J_{4}$) means the force of infection on mosquitoes due to the contact
with infectious humans with sensitive (resistant) strain.

Taking into account all of these assumptions, we obtain the following periodic
reaction-diffusion model:

$\left\\{\footnotesize\begin{split}&\frac{\partial I_{1}(t,x)}{\partial
t}=D_{h}\Delta
I_{1}(t,x)-(d+\gamma_{1})I_{1}(t,x)+J_{1}(I_{1}(t,x),I_{2}(t,x))I_{v1}(t,x),~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{v1}(t,x)}{\partial t}=D_{v}\Delta
I_{v1}(t,x)-\eta(t,x)I_{v1}(t,x)\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+J_{2}(I_{1}(t,x),I_{2}(t,x))(M^{*}(t,x)-I_{v1}(t,x)-I_{v2}(t,x)),~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{2}(t,x)}{\partial t}=D_{h}\Delta
I_{2}(t,x)-(d+\gamma_{2})I_{2}(t,x)+J_{3}(I_{1}(t,x),I_{2}(t,x))I_{v2}(t,x),~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{v2}(t,x)}{\partial t}=D_{v}\Delta
I_{v2}(t,x)-\eta(t,x)I_{v2}(t,x)\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+J_{4}(I_{1}(t,x),I_{2}(t,x))(M^{*}(t,x)-I_{v1}(t,x)-I_{v2}(t,x)),~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{1}(t,x)}{\partial\nu}=\frac{\partial
I_{v1}(t,x)}{\partial\nu}=\frac{\partial
I_{2}(t,x)}{\partial\nu}=\frac{\partial
I_{v2}(t,x)}{\partial\nu}=0,~{}&&t>0,~{}x\in\partial\Omega,\\\
&I_{1}(0,x)=I_{1}^{0}(x)\geq 0,~{}~{}I_{v1}(0,x)=I_{v1}^{0}(x)\geq 0,\\\
&I_{2}(0,x)=I_{2}^{0}(x)\geq 0,~{}~{}I_{v2}(0,x)=I_{v2}^{0}(x)\geq
0,~{}&&x\in\bar{\Omega}.\end{split}\right.$ (2.3)

Here, the positive constants $\gamma_{1}$ and $\gamma_{2}$ denote the recovery
rate of the sensitive and resistant strains for humans, respectively. The
function $\beta(t,x)$ is Hölder continuous and nonnegative but not zero
identically on $\mathbb{R}\times\bar{\Omega}$, and $\omega$-periodic in $t$.
Other parameters are the same as above.

Let $\mathbb{X}:=C(\bar{\Omega},\mathbb{R}^{4})$ be the Banach space with
supremum norm $\|\cdot\|$ and
$\mathbb{X}^{+}:=C(\bar{\Omega},\mathbb{R}_{+}^{4})$. For each $t\geq 0$, we
define

$\begin{split}X(t):=\\{\varphi=(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})\in\mathbb{X}^{+}:~{}&0\leq\varphi_{1}(x)+\varphi_{3}(x)\leq
N(x),\\\ &0\leq\varphi_{2}(x)+\varphi_{4}(x)\leq M^{*}(t,x),~{}\forall
x\in\bar{\Omega}\\}.\end{split}$

Let $\mathbb{Y}:=C(\bar{\Omega},\mathbb{R})$ and
$\mathbb{Y}^{+}:=C(\bar{\Omega},\mathbb{R}_{+})$. Let
$T_{1}(t,s),~{}T_{2}(t,s),~{}T_{3}(t,s):\mathbb{Y}\rightarrow\mathbb{Y},t\geq
s$, be the linear evolution operators associated with

$\begin{split}\frac{\partial v_{1}(t,x)}{\partial t}&=D_{h}\Delta
v_{1}(t,x)-(d+\gamma_{1})v_{1}(t,x):=A_{1}v_{1}(t,x),\\\ \frac{\partial
v_{2}(t,x)}{\partial t}&=D_{h}\Delta
v_{2}(t,x)-(d+\gamma_{2})v_{2}(t,x):=A_{2}v_{2}(t,x),\\\ \frac{\partial
v_{3}(t,x)}{\partial t}&=D_{v}\Delta
v_{3}(t,x)-\eta(t,x)v_{3}(t,x):=A_{3}v_{3}(t,x),\end{split}$

subject to the Neumann boundary condition, respectively. Noting that
$T_{j}(t,s)=T_{j}(t-s),j=1,2$, we have $T_{j}(t+\omega,s+\omega)=T_{j}(t,s)$
for $(t,s)\in\mathbb{R}^{2}$ with $t\geq s,j=1,2$. Since $\eta(t,x)$ is
$\omega$-periodic in $t$, [28, Lemma 6.1] implies that
$T_{3}(t+\omega,s+\omega)=T_{3}(t,s)$ for $(t,s)\in\mathbb{R}^{2}$ with $t\geq
s$. Moreover, for $(t,s)\in\mathbb{R}^{2}$ with $t>s$,
$T_{j}(t,s),~{}j=1,2,3$, are compact and strongly positive. Set $T(t,s)={\rm
diag}\\{T_{1}(t,s),T_{3}(t,s),T_{2}(t,s),T_{3}(t,s)\\}$ and $A(t)={\rm
diag}\\{A_{1},A_{3}(t),A_{2},A_{3}(t)\\}$. Define
$F=(F_{1},F_{2},F_{3},F_{4}):[0,\infty)\times\mathbb{X}^{+}\rightarrow\mathbb{X}$
by

$\displaystyle F_{1}(t,\varphi)$
$\displaystyle=\frac{c_{1}\beta(t,\cdot)l(N(\cdot)-\varphi_{1}(\cdot)-\varphi_{3}(\cdot))}{p(\varphi_{1}(\cdot)+\varphi_{3}(\cdot))+l(N(\cdot)-\varphi_{1}(\cdot)-\varphi_{3}(\cdot))}\varphi_{2}(\cdot),$
$\displaystyle F_{2}(t,\varphi)$
$\displaystyle=\frac{\alpha_{1}\beta(t,\cdot)p\varphi_{1}(\cdot)}{p(\varphi_{1}(\cdot)+\varphi_{3}(\cdot))+l(N(\cdot)-\varphi_{1}(\cdot)-\varphi_{3}(\cdot))}(M^{*}(t,\cdot)-\varphi_{2}(\cdot)-\varphi_{4}(\cdot)),$
$\displaystyle F_{3}(t,\varphi)$
$\displaystyle=\frac{c_{2}\beta(t,\cdot)l(N(\cdot)-\varphi_{1}(\cdot)-\varphi_{3}(\cdot))}{p(\varphi_{1}(\cdot)+\varphi_{3}(\cdot))+l(N(\cdot)-\varphi_{1}(\cdot)-\varphi_{3}(\cdot))}\varphi_{4}(\cdot),$
$\displaystyle F_{4}(t,\varphi)$
$\displaystyle=\frac{\alpha_{2}\beta(t,\cdot)p\varphi_{3}(\cdot)}{p(\varphi_{1}(\cdot)+\varphi_{3}(\cdot))+l(N(\cdot)-\varphi_{1}(\cdot)-\varphi_{3}(\cdot))}(M^{*}(t,\cdot)-\varphi_{2}(\cdot)-\varphi_{4}(\cdot))$

for all $t\geq 0$ and
$\varphi=(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})\in\mathbb{X}$. Then
system $\eqref{model}$ becomes

$\left\\{\begin{split}&\frac{du}{dt}=A(t)u+F(t,u),~{}~{}~{}&t>0,\\\
&u(0)=\varphi\in\mathbb{X}^{+},\end{split}\right.$

which can be written as an integral equation

$\begin{split}u(t,\varphi)=T(t,0)\varphi+\int_{0}^{t}T(t,s)F(s,u)ds,~{}~{}~{}\forall
t\geq 0,~{}\varphi\in\mathbb{X}^{+},\end{split}$ (2.4)

where

$u(t,x)=(u_{1}(t,x),u_{2}(t,x),u_{3}(t,x),u_{4}(t,x))=(I_{1}(t,x),I_{v1}(t,x),I_{2}(t,x),I_{v2}(t,x)).$

As usual, solutions of (2.4) are called mild solutions to system (2.3).

######  Lemma 2.1.

For any $\varphi\in X(0)$, system (2.3) has a unique solution
$u(t,\cdot,\varphi)$ with $u(0,\cdot,\varphi)=\varphi$ such that
$u(t,\cdot,\varphi)\in X(t)$ for all $t\in[0,\infty)$. Moreover, system (2.3)
generates an $\omega$-periodic semiflow $Q(t):X(0)\rightarrow X(t)$ defined by
$Q(t)\varphi=u(t,\cdot,\varphi),t\geq 0$. In addition, $Q:=Q(\omega)$ admits a
strong global attractor in $X(0)$.

###### Proof.

From the expression of $F$, we see that $F$ is locally Lipschitz continuous.
For any $(t,\psi)\in\mathbb{R}_{+}\times\mathbb{X}^{+}$ and $h>0$, in view of
$p\geq l>0$, we have

$\begin{array}[]{c}\begin{aligned} &\psi(x)+hF(t,\psi)(x)\\\
&=\left(\begin{aligned}
&\psi_{1}(x)+h\frac{c_{1}\beta(t,x)l(N(x)-\psi_{1}(x)-\psi_{3}(x))}{p(\psi_{1}(x)+\psi_{3}(x))+l(N(x)-\psi_{1}(x)-\psi_{3}(x))}\psi_{2}(x)\\\
&\psi_{2}(x)+h\frac{\alpha_{1}\beta(t,x)p\psi_{1}(x)(M^{*}(t,x)-\psi_{2}(x)-\psi_{4}(x))}{p(\psi_{1}(x)+\psi_{3}(x))+l(N(x)-\psi_{1}(x)-\psi_{3}(x))}\\\
&\psi_{3}(x)+h\frac{c_{2}\beta(t,x)l(N(x)-\psi_{1}(x)-\psi_{3}(x))}{p(\psi_{1}(x)+\psi_{3}(x))+l(N(x)-\psi_{1}(x)-\psi_{3}(x))}\psi_{4}(x)\\\
&\psi_{4}(x)+h\frac{\alpha_{2}\beta(t,x)p\psi_{3}(x)(M^{*}(t,x)-\psi_{2}(x)-\psi_{4}(x))}{p(\psi_{1}(x)+\psi_{3}(x))+l(N(x)-\psi_{1}(x)-\psi_{3}(x))}\end{aligned}\right)\\\
&\geq\left(\begin{aligned}
&\psi_{1}(x)\left(1-h\frac{c_{1}\beta(t,x)}{N(x)}\psi_{2}(x)\right)-h\frac{c_{1}\beta(t,x)}{N(x)}\psi_{3}(x)\psi_{2}(x)\\\
&\psi_{2}(x)\left(1-h\frac{\alpha_{1}\beta(t,x)p}{lN(x)}\psi_{1}(x)\right)-h\frac{\alpha_{1}\beta(t,x)p}{lN(x)}\psi_{1}(x)\psi_{4}(x)\\\
&\psi_{3}(x)\left(1-h\frac{c_{2}\beta(t,x)}{N(x)}\psi_{4}(x)\right)-h\frac{c_{2}\beta(t,x)}{N(x)}\psi_{1}(x)\psi_{4}(x)\\\
&\psi_{4}(x)\left(1-h\frac{\alpha_{2}\beta(t,x)p}{lN(x)}\psi_{3}(x)\right)-h\frac{\alpha_{2}\beta(t,x)p}{lN(x)}\psi_{3}(x)\psi_{2}(x)\end{aligned}\right).\end{aligned}\end{array}$

This implies that

$\lim_{h\rightarrow 0^{+}}\frac{1}{h}{\rm
dist}(\psi+hF(t,\psi),~{}\mathbb{X}^{+})=0,~{}~{}~{}\forall(t,\psi)\in\mathbb{R}_{+}\times\mathbb{X}^{+}.$

In addition, $T(t,s)\mathbb{X}^{+}\subseteq\mathbb{X}^{+},\forall t\geq s\geq
0$. Therefore, by [29, Corollary 4] with $K=\mathbb{X}^{+}$ and
$S(t,s)=T(t,s)$, system (2.3) admits a unique non-continuable mild solution
$u(t,\cdot,\psi)$ on its maximal existence interval $[0,t_{\psi})$ with
$u(0,\cdot,\psi)=\psi$, and $u(t,\cdot,\psi)\geq 0$ for all
$t\in[0,t_{\psi})$, where $t_{\psi}\leq\infty$.

Based on the above analysis, we obtain that for any $\varphi\in
X(0)\subset\mathbb{X}^{+}$, system (2.3) has a unique solution
$u(t,\cdot,\varphi)\in\mathbb{X}^{+}$ on $[0,t_{\varphi})$ with
$u(0,\cdot,\varphi)=\varphi$, where $t_{\varphi}\leq\infty$. Next we want to
show that $u(t,x,\varphi)$ is bounded for all $t\in[0,t_{\varphi})$, which
then implies $t_{\varphi}=\infty$. To this end, we set

$I_{h}(t,x)=I_{1}(t,x)+I_{2}(t,x),~{}~{}~{}I_{v}(t,x)=I_{v1}(t,x)+I_{v2}(t,x).$

It turns out that $N(x)$ and $M^{*}(t,x)$ are respectively the upper solutions
of the following two equations

$\displaystyle\frac{\partial I_{h}(t,x)}{\partial t}=$ $\displaystyle
D_{h}\Delta I_{h}(t,x)-dI_{h}(t,x)-\gamma_{1}I_{1}(t,x)-\gamma_{2}I_{2}(t,x)$
$\displaystyle+J_{1}(I_{1}(t,x),I_{2}(t,x))I_{v1}(t,x)+J_{3}(I_{1}(t,x),I_{2}(t,x))I_{v2}(t,x),$

and

$\displaystyle\frac{\partial I_{v}(t,x)}{\partial t}=$ $\displaystyle
D_{h}\Delta I_{v}(t,x)-\eta(t,x)I_{v}(t,x)$
$\displaystyle+J_{2}(I_{1}(t,x),I_{2}(t,x))(M^{*}(t,x)-I_{v}(t,x))$
$\displaystyle+J_{4}(I_{1}(t,x),I_{2}(t,x))(M^{*}(t,x)-I_{v}(t,x))$

for $t\in(0,t_{\varphi})$ and $x\in\bar{\Omega}$. Thus, the comparison
principle implies that solutions of (2.3) are bounded on $[0,t_{\varphi})$,
and thus, $t_{\varphi}=\infty.$ In addition, we also have that
$u(t,\cdot,\varphi)\in X(t)$ for all $t\geq 0$, and it is classic for $t>0$ in
light of the analyticity of $T(t,s)$ when $t>s.$

Define a family of operators $\\{Q(t)\\}_{t\geq 0}$ from $X(0)$ to $X(t)$ by

$[Q(t)\varphi](x)=u(t,x,\varphi),~{}~{}~{}\forall\varphi\in
X(0),~{}x\in\bar{\Omega}.$

By the proof of [27, Lemma 2.1], we can show that $Q(t)$ is an
$\omega$-periodic semiflow, and thus $Q:=Q(\omega):X(0)\rightarrow
X(\omega)=X(0)$ is the Poincaré map associated with system (2.3). The fact
that $u(t,\cdot,\varphi)\in X(t)$ for all $t\geq 0$ when $\varphi\in X(0)$
also implies that solutions of (2.3) are ultimately bounded. Hence, by [30,
Theorem 2.9], $Q$ has a strong global attractor in $X(0)$. ∎

######  Lemma 2.2.

For any $\varphi\in X(0)$, let $u(t,x,\varphi)$ be the solution of system
(2.3). If there exists some $t_{0}\geq 0$ such that
$u_{i}(t_{0},x,\varphi)\not\equiv 0,~{}i=1,2,3,4$, then

$u_{i}(t,x,\varphi)>0,~{}~{}~{}i=1,2,3,4,~{}\forall t\geq
t_{0},~{}x\in\bar{\Omega}.$

###### Proof.

For any given $\varphi\in X(0)$, one easily sees

$\left\\{\begin{split}\frac{\partial u_{1}(t,x)}{\partial t}&\geq D_{h}\Delta
u_{1}(t,x)-(d+\gamma_{1})u_{1}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial u_{2}(t,x)}{\partial t}&\geq D_{v}\Delta
u_{2}(t,x)-\bar{\eta}u_{2}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial u_{3}(t,x)}{\partial t}&\geq D_{h}\Delta
u_{3}(t,x)-(d+\gamma_{2})u_{3}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial u_{4}(t,x)}{\partial t}&\geq D_{v}\Delta
u_{4}(t,x)-\bar{\eta}u_{4}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial u_{1}(t,x)}{\partial\nu}&=\frac{\partial
u_{2}(t,x)}{\partial\nu}=\frac{\partial
u_{3}(t,x)}{\partial\nu}=\frac{\partial
u_{4}(t,x)}{\partial\nu}=0,~{}~{}~{}&&t>0,~{}x\in\partial\Omega,\end{split}\right.$

where $\bar{\eta}=\max_{(t,x)\in[0,\omega]\times\bar{\Omega}}\eta(t,x)$. If
there exists $t_{0}\geq 0$ such that $u_{i}(t_{0},x,\varphi)\not\equiv 0$ for
some $i\in\\{1,2,3,4\\}$, it then follows from the parabolic maximum principle
[31, Proposition 13.1] that $u_{i}(t,x,\varphi)>0$ for all $t>t_{0}$ and
$x\in\bar{\Omega}$. ∎

##  3 Reproduction numbers

In this section, we first define the basic reproduction number
$\mathcal{R}_{0}$ of (2.3), and then introduce the invasion reproduction
number $\hat{\mathcal{R}}_{i}$ for strain $i~{}(i=1,2)$.

###  3.1 Basic reproduction number

In order to derive the basic reproduction number of (2.3), we first consider
subsystems: one involves sensitive strains alone and the other involves
resistant strains alone. We fix $i\in\\{1,2\\}$ and let $I_{j}(t,x)\equiv
0,I_{vj}(t,x)\equiv 0,~{}\forall t\geq 0,~{}x\in\bar{\Omega},~{}j=1,2$ and
$j\neq i$. Then system (2.3) reduces to the following single-strain model:

$\left\\{\begin{split}\frac{\partial I_{i}(t,x)}{\partial t}=&D_{h}\Delta
I_{i}(t,x)-(d+\gamma_{i})I_{i}(t,x)\\\
&+\frac{c_{i}\beta(t,x)l(N(x)-I_{i}(t,x))}{pI_{i}(t,x)+l(N(x)-I_{i}(t,x))}I_{vi}(t,x),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{vi}(t,x)}{\partial t}=&D_{v}\Delta
I_{vi}(t,x)-\eta(t,x)I_{vi}(t,x)\\\
&+\frac{\alpha_{i}\beta(t,x)pI_{i}(t,x)}{pI_{i}(t,x)+l(N(x)-I_{i}(t,x))}(M^{*}(t,x)-I_{vi}(t,x)),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{i}(t,x)}{\partial\nu}=&\frac{\partial
I_{vi}(t,x)}{\partial\nu}=0,~{}&&t>0,~{}x\in\partial\Omega,\end{split}\right.$
(3.1)

Let $\mathbb{E}:=C(\bar{\Omega},\mathbb{R}^{2})$ and
$\mathbb{E}^{+}:=C(\bar{\Omega},\mathbb{R}_{+}^{2})$. Linearizing (3.1) at
$(0,0)$ yields

$\left\\{\begin{split}\frac{\partial I_{i}(t,x)}{\partial t}=&D_{h}\Delta
I_{i}(t,x)-(d+\gamma_{i})I_{i}(t,x)\\\
&+c_{i}\beta(t,x)I_{vi}(t,x),~{}&&t>0,~{}x\in\Omega,\\\ \frac{\partial
I_{vi}(t,x)}{\partial t}=&D_{v}\Delta I_{vi}(t,x)-\eta(t,x)I_{vi}(t,x)\\\
&+\frac{\alpha_{i}\beta(t,x)pM^{*}(t,x)}{lN(x)}I_{i}(t,x),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{i}(t,x)}{\partial\nu}=&\frac{\partial
I_{vi}(t,x)}{\partial\nu}=0,~{}&&t>0,~{}x\in\partial\Omega.\\\
\end{split}\right.$ (3.2)

Define the operator $\mathcal{F}_{i}(t):\mathbb{E}\rightarrow\mathbb{E}$ by

$\mathcal{F}_{i}(t)\left(\begin{array}[]{c}\psi_{1}\\\
\psi_{2}\end{array}\right)=\left(\begin{array}[]{c}c_{i}\beta(t,\cdot)\psi_{2}(\cdot)\\\
\frac{\alpha_{i}\beta(t,\cdot)pM^{*}(t,\cdot)}{lN(\cdot)}\psi_{1}(\cdot)\end{array}\right),~{}~{}~{}\forall
t\in\mathbb{R},~{}\psi=(\psi_{1},\psi_{2})\in\mathbb{E}.$

Let $-V_{i}(t)v=D\Delta v-W_{i}(t)v$, where $D={\rm diag}(D_{h},D_{v})$ and

$-[W_{i}(t)](x)=\left(\begin{array}[]{cc}-(d+\gamma_{i})&0\\\
0&-\eta(t,x)\end{array}\right),~{}~{}~{}\forall
t\in\mathbb{R},~{}x\in\bar{\Omega}.$

Then $\Psi_{i}(t,s)={\rm diag}(T_{i}(t,s),T_{3}(t,s)),~{}t\geq s$, is the
evolution operator on $\mathbb{E}$ associated with the following system

$\frac{dv}{dt}=-V_{i}(t)v$

subject to the Neumann boundary condition. The exponential growth bound of
$\Psi_{i}(t,s)$ is defined as

$\bar{\omega}(\Psi_{i})=\inf\\{\tilde{\omega}_{i}:~{}\exists M\geq
1~{}\mbox{such that}~{}\|\Psi_{i}(t+s,s)\|_{\mathbb{E}}\leq
Me^{\tilde{\omega}_{i}t},~{}\forall s\in\mathbb{R},~{}t\geq 0\\}.$

By the Krein-Rutman Theorem and [31, Lemma 14.2], we have

$0<r(\Psi_{i}(\omega,0))=\max\\{r(T_{i}(\omega,0)),r(T_{3}(\omega,0))\\}<1,$

where $r(\Psi_{i}(\omega,0))$ is the spectral radius of $\Psi_{i}(\omega,0)$.
Then, it follows from [32, Proposition 5.5] with $s=0$ that
$\bar{\omega}(\Psi_{i})<0$. Note that $\Psi_{i}(t,s)$ is a positive operator
in the sense that $\Psi_{i}(t,s)\mathbb{E}^{+}\subseteq\mathbb{E}^{+}$ for all
$t\geq s$. Therefore, $\mathcal{F}_{i}(t)$ and $\Psi_{i}(t,s)$ satisfy

  * (H1)

For each $t\geq 0$, $\mathcal{F}_{i}(t)$ is a positive operator on
$\mathbb{E}$.

  * (H2)

For any $t\geq s$, $\Psi_{i}(t,s)$ is a positive operator on $\mathbb{E}$, and
$\bar{\omega}(\Psi_{i})<0$.

Let $C_{\omega}(\mathbb{R},\mathbb{E})$ be the Banach space of all
$\omega$-periodic and continuous functions from $\mathbb{R}$ to $\mathbb{E}$
equipped with the maximum norm. Following the theory developed in [33, 34], we
define two linear operators on $C_{\omega}(\mathbb{R},\mathbb{E})$ by

$[\mathcal{L}_{i}v](t):=\int^{\infty}_{0}\Psi_{i}(t,t-s)\mathcal{F}_{i}(t-s)v(t-s)ds,~{}~{}~{}\forall
t\in\mathbb{R},~{}v\in C_{\omega}(\mathbb{R},\mathbb{E}),~{}i=1,2.$

Motivated by the concept of next generation operators [32, 35], we define the
basic reproduction number as $\mathcal{R}_{i}:=r(\mathcal{L}_{i})$, where
$r(\mathcal{L}_{i})$ is the spectral radius of $\mathcal{L}_{i}$.

The disease-free state of (2.3) is $(0,0,0,0)$ and the corresponding
linearized system is

$\left\\{\footnotesize\begin{split}&\frac{\partial I_{1}(t,x)}{\partial
t}=D_{h}\Delta
I_{1}(t,x)-(d+\gamma_{1})I_{1}(t,x)+c_{1}\beta(t,x)I_{v1}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{v1}(t,x)}{\partial t}=D_{v}\Delta
I_{v1}(t,x)-\eta(t,x)I_{v1}(t,x)+\frac{\alpha_{1}\beta(t,x)pM^{*}(t,x)}{lN(x)}I_{1}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{2}(t,x)}{\partial t}=D_{h}\Delta
I_{2}(t,x)-(d+\gamma_{2})I_{2}(t,x)+c_{2}\beta(t,x)I_{v2}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{v2}(t,x)}{\partial t}=D_{v}\Delta
I_{v2}(t,x)-\eta(t,x)I_{v2}(t,x)+\frac{\alpha_{2}\beta(t,x)pM^{*}(t,x)}{lN(x)}I_{2}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
&\frac{\partial I_{1}(t,x)}{\partial\nu}=\frac{\partial
I_{v1}(t,x)}{\partial\nu}=\frac{\partial
I_{2}(t,x)}{\partial\nu}=\frac{\partial
I_{v2}(t,x)}{\partial\nu}=0,~{}~{}~{}&&t>0,~{}x\in\partial\Omega.\\\
\end{split}\right.$ (3.3)

Similarly, we can derive the basic reproduction number of (2.3), which is
given by

$\mathcal{R}_{0}=\max\\{\mathcal{R}_{1},\mathcal{R}_{2}\\}.$

For any given $t\geq 0$, let $P_{i}(t)$ be the solution map of (3.2) on
$\mathbb{E}$. Then $P_{i}:=P_{i}(\omega)$ is the associated Poincaré map. Let
$r(P_{i})$ be the spectral radius of $P_{i}$. By [34, Theorem 3.7] with
$\tau=0$, we have the following nice property.

######  Lemma 3.1.

$\mathcal{R}_{i}-1$ has the same sign as $r(P_{i})-1,i=1,2$, and thus
$\mathcal{R}_{0}-1$ has the same sign as $r(P)-1$, where
$r(P)=\max\\{r(P_{1}),r(P_{2})\\}$ is the spectral radius of the Poincaré map
$P$ associated with (3.3).

###  3.2 Invasion reproduction number

In this subsection, we define the invasion reproduction number for each
strain. The invasion reproduction number gives the ability of strain
$i~{}(i=1,2)$ to invade strain $j~{}(j=1,2,~{}j\neq i)$ measured as the number
of secondary infections strain $i$ one-infected individual can produce in a
population where strain $j$ is at an endemic state [36]. We express it by
$\hat{\mathcal{R}_{i}}~{}(i=1,2)$ and give their definition by analyzing the
boundary $\omega$-periodic solution of (2.3), that is, the sensitive strain
$\omega$-periodic solution or resistant strain $\omega$-periodic solution.

For each $t\geq 0$, let $E(t)$ be subset in $\mathbb{E}$ defined by

$E(t):=\\{\psi=(\psi_{1},\psi_{2})\in\mathbb{E}^{+}:0\leq\psi_{1}(x)\leq
N(x),~{}0\leq\psi_{2}(x)\leq M^{*}(t,x),~{}\forall x\in\bar{\Omega}\\}.$

After a similar process in [25, Lemma 3], we obtain that for any $\psi\in
E(0)$, system (3.1) has a unique solution
$v_{i}(t,\cdot,\psi)=(I_{i}(t,x),I_{vi}(t,x))$ with $v_{i}(0,\cdot,\psi)=\psi$
such that $v_{i}(t,\cdot,\psi)\in E(t)$ for all $t\geq 0$. Moreover, by
employing the arguments in [25, Theorem 1], one immediately obtains the
following result.

######  Theorem 3.1.

The following statements are valid:

  * (i)

If $\mathcal{R}_{i}\leq 1$, then $(0,0)$ is globally asymptotically stable for
(3.1) in $E(0)$;

  * (ii)

If $\mathcal{R}_{i}>1$, then system (3.1) admits a unique positive
$\omega$-periodic solution $(I_{i}^{*}(t,x),I_{vi}^{*}(t,x))$, and it is
globally asymptotically stable for (3.1) in $E(0)\setminus\\{(0,0)\\}$.

For ease of presentation, we introduce the following notations:
$\bullet$ $E_{0}=(0,0,0,0)$: The disease-free state of (2.3).
$\bullet$ $E_{1}(t,x)=(I_{1}^{*}(t,x),I_{v1}^{*}(t,x),0,0)$: The sensitive
strain $\omega$-periodic solution of (2.3).
$\bullet$ $E_{2}(t,x)=(0,0,I_{2}^{*}(t,x),I_{v2}^{*}(t,x))$: The resistant
strain $\omega$-periodic solution of (2.3).

By Theorem 3.1, we see that when $\mathcal{R}_{i}>1~{}(i=1,2)$, system (2.3)
admits a unique semitrivial boundary $\omega$-periodic solution $E_{i}(t,x)$.
Linearizing (2.3) at the $E_{j}(t,x),~{}j\neq i,~{}i,j=1,2$, and considering
only the equations for $I_{i}(t,x)$ and $I_{vi}(t,x)$, we get

$\left\\{\begin{split}\frac{\partial I_{i}(t,x)}{\partial t}=&D_{h}\Delta
I_{i}(t,x)-(d+\gamma_{i})I_{i}(t,x)\\\
&+\frac{c_{i}\beta(t,x)l(N(x)-I_{j}^{*}(t,x))}{pI_{j}^{*}(t,x)+l(N(x)-I_{j}^{*}(t,x))}I_{vi}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{vi}(t,x)}{\partial t}=&D_{v}\Delta
I_{vi}(t,x)-\eta(t,x)I_{vi}(t,x)\\\
&+\frac{\alpha_{i}\beta(t,x)p(M^{*}(t,x)-I_{vj}^{*}(t,x))}{pI_{j}^{*}(t,x)+l(N(x)-I_{j}^{*}(t,x))}I_{i}(t,x),~{}~{}~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{i}(t,x)}{\partial\nu}=&\frac{\partial
I_{vi}(t,x)}{\partial\nu}=0,~{}~{}~{}&&t>0,~{}x\in\partial\Omega.\\\
\end{split}\right.$ (3.4)

Similar to Section 3.1, we can define the invasion reproduction numbers
$\hat{\mathcal{R}}_{i}~{}(i=1,2)$. Further, we have the following
characterization of $\hat{\mathcal{R}}_{i}$.

######  Lemma 3.2.

$\hat{\mathcal{R}}_{i}-1$ has the same sign as $r(\hat{P}_{i})-1$, where
$\hat{P}_{i}$ is the Poincaré map associated with (3.4), and $r(\hat{P}_{i})$
is the spectral radius of $\hat{P}_{i}$.

##  4 Disease extinction and uniform persistence

In this section, we establish the dynamics of (2.3) in terms of
$\mathcal{R}_{i}$ and $\hat{\mathcal{R}}_{i},i=1,2$.

###  4.1 Global extinction

######  Theorem 4.1.

If $\mathcal{R}_{1}<1$ and $\mathcal{R}_{2}<1$, then $E_{0}$ is globally
attractive for (2.3) in $X(0)$.

###### Proof.

Let $(I_{1}(t,x),I_{v1}(t,x),I_{2}(t,x),I_{v2}(t,x))$ be the solution of (2.3)
with initial data $\varphi\in X(0)$. It is easily seen that

$\left\\{\begin{split}\frac{\partial I_{1}(t,x)}{\partial t}\leq&D_{h}\Delta
I_{1}(t,x)-(d+\gamma_{1})I_{1}(t,x)\\\
&+\frac{c_{1}\beta(t,x)l(N(x)-I_{1}(t,x))}{pI_{1}(t,x)+l(N(x)-I_{1}(t,x))}I_{v1}(t,x),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{v1}(t,x)}{\partial t}\leq&D_{v}\Delta
I_{v1}(t,x)-\eta(t,x)I_{v1}(t,x)\\\
&+\frac{\alpha_{1}\beta(t,x)pI_{1}(t,x)}{pI_{1}(t,x)+l(N(x)-I_{1}(t,x))}(M^{*}(t,x)-I_{v1}(t,x)),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{1}(t,x)}{\partial\nu}=&\frac{\partial
I_{v1}(t,x)}{\partial\nu}=0,~{}&&t>0,~{}x\in\partial\Omega.\end{split}\right.$
(4.1)

and

$\left\\{\begin{split}\frac{\partial I_{2}(t,x)}{\partial t}\leq&D_{h}\Delta
I_{2}(t,x)-(d+\gamma_{2})I_{2}(t,x)\\\
&+\frac{c_{2}\beta(t,x)l(N(x)-I_{2}(t,x))}{pI_{2}(t,x)+l(N(x)-I_{2}(t,x))}I_{v2}(t,x),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{v2}(t,x)}{\partial t}\leq&D_{v}\Delta
I_{v2}(t,x)-\eta(t,x)I_{v2}(t,x)\\\
&+\frac{\alpha_{2}\beta(t,x)pI_{2}(t,x)}{pI_{2}(t,x)+l(N(x)-I_{2}(t,x))}(M^{*}(t,x)-I_{v2}(t,x)),~{}&&t>0,~{}x\in\Omega,\\\
\frac{\partial I_{2}(t,x)}{\partial\nu}=&\frac{\partial
I_{v2}(t,x)}{\partial\nu}=0,~{}&&t>0,~{}x\in\partial\Omega.\end{split}\right.$
(4.2)

When $\mathcal{R}_{1}<1$, Theorem 3.1 implies that $(0,0)$ is globally stable
for (3.1). Hence the comparison principle applies to (4.1) and ensures that

$\lim_{t\rightarrow\infty}(I_{1}(t,x),I_{v1}(t,x))=(0,0)~{}\mbox{uniformly
for}~{}x\in\bar{\Omega}.$

In the case where $\mathcal{R}_{2}<1$, by using the similar procedure as above
to (4.2), one attains

$\lim_{t\rightarrow\infty}(I_{2}(t,x),I_{v2}(t,x))=(0,0)~{}\mbox{uniformly
for}~{}x\in\bar{\Omega}.$

Therefore, the desired result is established. ∎

###  4.2 Competitive exclusion and coexistence

######  Theorem 4.2.

Let
$(I_{1}(t,\cdot,\varphi),I_{v1}(t,\cdot,\varphi),I_{2}(t,\cdot,\varphi),I_{v2}(t,\cdot,\varphi))$
be the solution of (2.3) through $\varphi\in X(0)$. Then the following
assertions hold.

  * (1)

If $\mathcal{R}_{1}>1>\mathcal{R}_{2}$ and $\varphi_{1},\varphi_{2}\not\equiv
0$, then

$\lim_{t\rightarrow\infty}((I_{1}(t,x,\varphi),I_{v1}(t,x,\varphi),I_{2}(t,x,\varphi),I_{v2}(t,x,\varphi))-E_{1}(t,x))=(0,0,0,0)$

uniformly for $x\in\bar{\Omega}$.

  * (2)

If $\mathcal{R}_{2}>1>\mathcal{R}_{1}$ and $\varphi_{3},\varphi_{4}\not\equiv
0$, then

$\lim_{t\rightarrow\infty}((I_{1}(t,x,\varphi),I_{v1}(t,x,\varphi),I_{2}(t,x,\varphi),I_{v2}(t,x,\varphi))-E_{2}(t,x))=(0,0,0,0)$

uniformly for $x\in\bar{\Omega}$.

###### Proof.

We only prove statement (1), since statement (2) can be treated similarly. In
the case where $\mathcal{R}_{2}<1$, one immediately has that
$\lim_{t\rightarrow\infty}(I_{2}(t,\cdot,\varphi),I_{v2}(t,\cdot,\varphi))=(0,0)$.
Then the limiting system of (2.3) is the system (3.1) with $i=1$. Moreover, by
employing the theory of internally chain transitive sets (see, e.g., [26]), we
conclude that

$\lim_{t\rightarrow\infty}(I_{1}(t,x,\varphi),I_{12}(t,x,\varphi)-(I_{1}^{*}(t,x),I_{v1}^{*}(t,x))=(0,0)$

uniformly for $x\in\bar{\Omega}$. Hence, statement (1) is established. ∎

For each $t\geq 0$, define

$X_{0}(t):=\\{(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})\in\
X(t):~{}\varphi_{i}(\cdot)\not\equiv 0,~{}i=1,2,3,4\\},$

and

$\displaystyle\partial X_{0}(t):=X(t)\backslash
X_{0}(t)=\\{(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4})\in\ X(t):~{}$
$\displaystyle\varphi_{i}(\cdot)\equiv 0~{}\mbox{at least for one}~{}i\\}.$

In order to study the coexistence of strains, we first give the following
lemma for our subsequent coexistence result.

######  Lemma 4.1.

Let
$(I_{1}(t,\cdot,\varphi),I_{v1}(t,\cdot,\varphi),I_{2}(t,\cdot,\varphi),I_{v2}(t,\cdot,\varphi))$
be the solution of (2.3) with the initial value $\varphi\in X_{0}(0)$. If
$\mathcal{R}_{1}>1$ and $\mathcal{R}_{2}>1$, then there exists $\delta>0$ such
that

$\limsup_{t\rightarrow\infty}\|(I_{1}(t,\cdot,\varphi),I_{v1}(t,\cdot,\varphi),I_{2}(t,\cdot,\varphi),I_{v2}(t,\cdot,\varphi))-E_{0}\|\geq\delta.$

###### Proof.

Suppose, by contradiction, that there exists some $\psi\in X_{0}(0)$ such that

$\limsup_{t\rightarrow\infty}\|(I_{1}(t,\cdot,\psi),I_{v1}(t,\cdot,\psi),I_{2}(t,\cdot,\psi),I_{v2}(t,\cdot,\psi))-E_{0}\|<\delta.$

Then there exists a $t_{1}>0$ such that

$0<I_{i}(t,x,\psi)<\delta,~{}~{}~{}0<I_{vi}(t,x,\psi)<\delta,~{}~{}~{}i=1,2$

for all $t\geq t_{1}$ and $x\in\bar{\Omega}$. Then $I_{1}(t,\cdot,\psi)$ and
$I_{v1}(t,\cdot,\psi)$ satisfy

$\left\\{\begin{split}\frac{\partial I_{1}(t,x)}{\partial t}\geq&D_{h}\Delta
I_{1}(t,x)-(d+\gamma_{1})I_{1}(t,x)\\\
&+\frac{c_{1}\beta(t,x)l(N(x)-2\delta)}{2p\delta+lN(x)}I_{v1}(t,x),~{}~{}~{}&&t\geq
t_{1},~{}x\in\Omega,\\\ \frac{\partial I_{v1}(t,x)}{\partial
t}\geq&D_{v}\Delta I_{v1}(t,x)-\eta(t,x)I_{v1}(t,x)\\\
&+\frac{\alpha_{1}\beta(t,x)p(M^{*}(t,x)-2\delta)}{2p\delta+lN(x)}I_{1}(t,x),~{}~{}~{}&&t\geq
t_{1},~{}x\in\Omega,\\\ \frac{\partial
I_{1}(t,x)}{\partial\nu}=&\frac{\partial
I_{v1}(t,x)}{\partial\nu}=0,~{}~{}~{}&&t\geq t_{1},~{}x\in\partial\Omega.\\\
\end{split}\right.$

Let $P_{1}^{\delta}:X(0)\rightarrow X(0)$ be the Poincaré map associated with
the following system:

$\left\\{\begin{split}\frac{\partial I_{1}(t,x)}{\partial t}=&D_{h}\Delta
I_{1}(t,x)-(d+\gamma_{1})I_{1}(t,x)\\\
&+\frac{c_{1}\beta(t,x)l(N(x)-2\delta)}{2p\delta+lN(x)}I_{v1}(t,x),~{}~{}~{}&&t\geq
0,~{}x\in\Omega,\\\ \frac{\partial I_{v1}(t,x)}{\partial t}=&D_{v}\Delta
I_{v1}(t,x)-\eta(t,x)I_{v1}(t,x)\\\
&+\frac{\alpha_{1}\beta(t,x)p(M^{*}(t,x)-2\delta)}{2p\delta+lN(x)}I_{1}(t,x),~{}~{}~{}&&t\geq
0,~{}x\in\Omega,\\\ \frac{\partial I_{1}(t,x)}{\partial\nu}=&\frac{\partial
I_{v1}(t,x)}{\partial\nu}=0,~{}~{}~{}&&t\geq
0,~{}x\in\partial\Omega.\end{split}\right.$ (4.3)

In view of Lemma 3.1, we have that $\mathcal{R}_{1}>1$ is equivalent to
$r(P_{1})>1$. By continuity, we see that $\lim_{\delta\rightarrow
0}r(P_{1}^{\delta})=r(P_{1})>1$. Thus, we can fix a sufficiently small number
$\delta>0$ such that

$0<\delta<\min\\{\min_{(t,x)\in[0,\omega]\times\bar{\Omega}}\frac{M^{*}(t,x)}{2},~{}\min_{x\in\bar{\Omega}}\frac{N(x)}{2}\\}~{}~{}\mbox{and}~{}~{}r(P_{1}^{\delta})>1.$

Since $P_{1}^{\delta}$ is compact and strongly positive on $E(0)$, then Krein-
Rutman Theorem implies that $r(P_{1}^{\delta})$ is a simple eigenvalue of
$P_{1}^{\delta}$ having a strongly positive eigenvector. It then follows from
[37, Theorem 2.16 and Remark 2.20] that there is a positive $\omega$-periodic
function $\vartheta^{\delta}(t,x)$ such that
$e^{\mu^{\delta}t}\vartheta^{\delta}(t,x)$ is a positive solution of (4.3),
where $\mu^{\delta}=\frac{\ln r(P_{1}^{\delta})}{\omega}>0$. From Lemma 2.2,
we know that

$I_{1}(t,x,\psi)>0~{}~{}\mbox{and}~{}~{}I_{v1}(t,x,\psi)>0,~{}~{}~{}\forall
t>0,x\in\bar{\Omega}.$

Thus, we may choose a $c>0$ such that

$(I_{1}(t_{1},x,\psi),I_{v1}(t_{1},x,\psi))\geq
ce^{\mu^{\delta}t_{1}}\vartheta^{\delta}(t_{1},x),~{}~{}~{}\forall
x\in\bar{\Omega}.$

A simple comparison leads to

$(I_{1}(t,\cdot,\psi),I_{v1}(t,\cdot,\psi))\geq
ce^{\mu^{\delta}t}\vartheta^{\delta}(t,x),~{}~{}~{}\forall t\geq
t_{1},~{}x\in\bar{\Omega},$

Since $\mu^{\delta}>0$, it follows that

$\lim_{t\rightarrow\infty}I_{1}(t,x,\psi)=\infty,~{}\lim_{t\rightarrow\infty}I_{v1}(t,x,\psi)=\infty~{}~{}{\rm
uniformly~{}for}~{}x\in\bar{\Omega}.$

By performing a similar analysis on $(I_{2}(t,x,\psi),I_{v2}(t,x,\psi))$, when
$\mathcal{R}_{2}>1$,

$\lim_{t\rightarrow\infty}I_{2}(t,x,\psi)=\infty,~{}\lim_{t\rightarrow\infty}I_{v2}(t,x,\psi)=\infty~{}~{}{\rm
uniformly~{}for}~{}x\in\bar{\Omega}.$

This contradicts the boundedness of $I_{i}$ and $I_{vi}$, $i=1,2$. ∎

######  Theorem 4.3.

Suppose that
$\mathcal{R}_{1}>1,~{}\mathcal{R}_{2}>1,~{}\hat{\mathcal{R}}_{1}>1,$ and
$\hat{\mathcal{R}}_{2}>1$, then system (2.3) admits at least one positive
$\omega$-periodic solution, and there exists a constant $\delta^{*}>0$ such
that for any $\varphi\in X_{0}(0)$, we have

$\liminf_{t\rightarrow\infty}\min_{x\in\bar{\Omega}}I_{i}(t,x)\geq\delta^{*},~{}~{}~{}\liminf_{t\rightarrow\infty}\min_{x\in\bar{\Omega}}I_{vi}(t,x)\geq\delta^{*},~{}~{}~{}i=1,2.$

###### Proof.

For any $\varphi\in X_{0}(0)$, by Lemma 2.2, we have

$I_{i}(t,x,\varphi)>0,~{}I_{vi}(t,x,\varphi)>0,~{}~{}~{}i=1,2,~{}\forall
t>0,~{}x\in\bar{\Omega}.$

Thus, $Q^{n}(X_{0}(0))\subset X_{0}(0),~{}\forall n\in\mathbb{N}$.
Furthermore, $Q$ admits a global attractor on $X(0)$.

Next we prove that $Q$ is uniformly persistent with respect to
$(X_{0}(0),\partial X_{0}(0))$. Recalling the definitions of
$E_{0},E_{1}(t,x),E_{2}(t,x)$ in Section 3.2, we let

$M_{0}=E_{0},~{}~{}M_{1}=E_{1}(0,\cdot),~{}~{}M_{2}=E_{2}(0,\cdot).$

Then we have the following claims.

Claim 1. There is a $\delta_{1}>0$ such that

$\limsup_{n\rightarrow\infty}\|Q^{n}(\varphi)-M_{0}\|\geq\delta_{1},~{}~{}~{}\forall\varphi\in
X_{0}(0).$

This claim directly follows from Lemma 4.1.

Consider an auxiliary system with parameter $\varepsilon$:

$\left\\{\begin{split}\frac{\partial I_{2}(t,x)}{\partial t}=&D_{h}\Delta
I_{2}(t,x)-(d+\gamma_{2})I_{2}(t,x)\\\
&+\frac{c_{2}\beta(t,x)l(N(x)-I_{1}^{*}(t,x)-2\varepsilon)}{p(I_{1}^{*}(t,x)+2\varepsilon)+l(N(x)-I_{1}^{*}(t,x)-2\varepsilon)}I_{v2}(t,x),\\\
\frac{\partial I_{v2}(t,x)}{\partial t}=&D_{v}\Delta
I_{v2}(t,x)-\eta(t,x)I_{v2}(t,x)\\\
&+\frac{\alpha_{2}\beta(t,x)p(M^{*}(t,x)-I_{v1}^{*}(t,x)-2\varepsilon)}{p(I_{1}^{*}(t,x)+2\varepsilon)+l(N(x)-I_{1}^{*}(t,x)-2\varepsilon)}I_{2}(t,x),\\\
\frac{\partial I_{2}(t,x)}{\partial\nu}=&\frac{\partial
I_{v2}(t,x)}{\partial\nu}=0,\\\ \end{split}\right.$ (4.4)

for all $t>0$. Let
$\hat{P}_{2}^{\varepsilon}:=\hat{P}_{2}^{\varepsilon}(\omega)$ be the Poincaré
map of (4.4). Since $\lim_{\varepsilon\rightarrow
0}r(\hat{P}_{2}^{\varepsilon})=r(\hat{P}_{2})>1$, we can fix a small number
$\varepsilon>0$ such that $r(\hat{P}_{2}^{\varepsilon})>1$. As discussed in
Lemma 4.1, there is a positive $\omega$-periodic function
$\vartheta^{\varepsilon}(t,x)$ such that
$e^{\mu^{\varepsilon}t}\vartheta^{\varepsilon}(t,x)$ is a positive solution of
(4.4), where $\mu^{\varepsilon}=\frac{\ln
r(\hat{P}_{2}^{\varepsilon})}{\omega}>0$. For $\varepsilon>0$ above, by the
continuous dependence of solutions on the initial value, there exists
$\delta_{2}=\delta_{2}(\varepsilon)>0$ such that for all $\varphi\in X_{0}(0)$
with $\|\varphi-M_{1}\|\leq\delta_{2}$, we have
$\|Q(t)\varphi-Q(t)M_{1}\|<\varepsilon,~{}\forall t\in[0,\omega]$.

Claim 2.
$\limsup_{n\rightarrow\infty}\|Q^{n}(\varphi)-M_{1}\|\geq\delta_{2},~{}\forall\varphi\in
X_{0}(0).$

Suppose the claim is false, then
$\limsup_{n\rightarrow\infty}\|Q^{n}(\psi)-M_{1}\|<\delta_{2}$ for some
$\psi\in X_{0}(0)$. Then there exists an integer $N_{1}\geq 1$ such that
$\|Q^{n}(\psi)-M_{1}\|<\delta_{2}$ for $n\geq N_{1}$. For any $t\geq
N_{1}\omega$, letting $t=n\omega+\tilde{t}$ with $n=[t/\omega]$ and
$\tilde{t}\in[0,\omega)$, we have

$\|Q(t)\psi-Q(t)M_{1}\|=\|Q(\tilde{t})(Q^{n}(\psi))-Q(\tilde{t})M_{1}\|<\varepsilon.$

According to the above inequality and Lemma 2.2, we infer that

$\displaystyle
0<I_{1}(t,x,\psi)<I_{1}^{*}(t,x,\psi)+\varepsilon,~{}0<I_{v1}(t,x,\psi)<I_{v1}^{*}(t,x,\psi)+\varepsilon,$
$\displaystyle
0<I_{2}(t,x,\psi)<\varepsilon,~{}0<I_{v2}(t,x,\psi)<\varepsilon,~{}~{}~{}\forall
t\geq N_{1}\omega,~{}x\in\bar{\Omega}.$

As a result, $I_{2}(t,x,\psi)$ and $I_{v2}(t,x,\psi)$ satisfy

$\left\\{\begin{split}\frac{\partial I_{2}(t,x)}{\partial t}\geq&D_{h}\Delta
I_{2}(t,x)-(d+\gamma_{2})I_{2}(t,x)\\\
&+\frac{c_{2}\beta(t,x)l(N(x)-I_{1}^{*}(t,x)-2\varepsilon)}{p(I_{1}^{*}(t,x)+2\varepsilon)+l(N(x)-I_{1}^{*}(t,x)-2\varepsilon)}I_{v2}(t,x),\\\
\frac{\partial I_{v2}(t,x)}{\partial t}\geq&D_{v}\Delta
I_{v2}(t,x)-\eta(t,x)I_{v2}(t,x)\\\
&+\frac{\alpha_{2}\beta(t,x)p(M^{*}(t,x)-I_{v1}^{*}(t,x)-2\varepsilon)}{p(I_{1}^{*}(t,x)+2\varepsilon)+l(N(x)-I_{1}^{*}(t,x)-2\varepsilon)}I_{2}(t,x),\\\
\frac{\partial I_{2}(t,x)}{\partial\nu}=&\frac{\partial
I_{v2}(t,x)}{\partial\nu}=0,\\\ \end{split}\right.$ (4.5)

for all $t\geq N_{1}\omega$. Since $\psi\in X_{0}(0)$, $I_{2}(t,x,\psi)>0$ and
$I_{v2}(t,x,\psi)>0$ for all $t\geq 0$ and $x\in\bar{\Omega}$, there exists a
$k>0$ such that

$(I_{2}(N_{1}\omega,x,\psi),I_{v2}(N_{1}\omega,x,\psi))\geq
ke^{\mu^{\varepsilon}N_{1}\omega}\vartheta^{\varepsilon}(N_{1}\omega,x),~{}~{}~{}x\in\bar{\Omega}.$

An application of the comparison theorem to (4.5) yields

$(I_{2}(t,x,\psi),I_{v2}(t,x,\psi))\geq
ke^{\mu^{\varepsilon}t}\vartheta^{\varepsilon}(t,x),~{}~{}~{}\forall t\geq
N_{1}\omega,~{}x\in\bar{\Omega}.$

Since $\mu^{\varepsilon}>0$, one sees that
$I_{2}(t,\cdot,\psi)\rightarrow\infty,~{}I_{v2}(t,\cdot,\psi)\rightarrow\infty$
as $t\rightarrow\infty$. This gives rise to a contradiction, and thereby, the
above claim is true.

In a similar way, we can prove the following claim.

Claim 3. There exists a $\delta_{3}>0$ such that

$\limsup_{n\rightarrow\infty}\|Q^{n}(\varphi)-M_{2}\|\geq\delta_{3},~{}\forall\varphi\in
X_{0}(0).$

With the above three claims, we see that $M_{0},M_{1}$ and $M_{2}$ are
isolated invariant sets for $Q$ in $X(0)$ and $W^{s}(M_{i})\bigcap
X_{0}(0)=\emptyset,i=0,1,2$, where $W^{s}(M_{i})$ is the stable set of $M_{i}$
for $Q$. Set

$M_{\partial}:=\\{\varphi\in\partial X_{0}(0):~{}Q^{n}(\varphi)\in\partial
X_{0}(0),~{}\forall n\in\mathbb{N}\\}.$

We now show that $M_{\partial}=\mathcal{M}_{0}$, where

$\mathcal{M}_{0}:=\\{\varphi\in\partial
X_{0}(0):~{}\varphi_{1}(\cdot)=\varphi_{2}(\cdot)\equiv
0~{}\mbox{or}~{}\varphi_{3}(\cdot)=\varphi_{4}(\cdot)\equiv 0\\}.$

Obviously, it suffices to prove $M_{\partial}\subset\mathcal{M}_{0}$. For any
given $\varphi\in M_{\partial}$, we have $I_{1}(n\omega,\cdot,\varphi)\equiv
0$ or $I_{v1}(n\omega,\cdot,\varphi)\equiv 0$ or
$I_{2}(n\omega,\cdot,\varphi)\equiv 0$ or $I_{v2}(n\omega,\cdot,\varphi)\equiv
0,\forall n\in\mathbb{N}$. Assume that $\varphi\notin\mathcal{M}_{0}$, then
there are eight possibilities as below:

  * (i)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)>0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)>0$.

  * (ii)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)>0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)>0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)\equiv
0$.

  * (iii)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)>0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)>0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)>0$.

  * (iv)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)>0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)>0$.

  * (v)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)>0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)>0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)\equiv
0$.

  * (vi)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)>0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)>0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)>0$.

  * (vii)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)>0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)>0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)\equiv
0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)>0$.

  * (viii)

$\varphi_{1}(\cdot)=I_{1}(0,\cdot,\varphi)>0,~{}\varphi_{2}(\cdot)=I_{v1}(0,\cdot,\varphi)>0,~{}\varphi_{3}(\cdot)=I_{2}(0,\cdot,\varphi)>0,~{}\varphi_{4}(\cdot)=I_{v2}(0,\cdot,\varphi)\equiv
0$.

By Lemma 2.2, in case (i), we obtain that
$I_{v1}(t,x,\varphi)>0,I_{v2}(t,x,\varphi)>0$ for all $t>0$ and
$x\in\bar{\Omega}$. Further, using the first and third equation of (2.3), one
obtains that $I_{1}(t,x,\varphi)>0,~{}I_{2}(t,x,\varphi)>0,\forall
t>0,x\in\bar{\Omega}$, which contradicts with the fact $\varphi\in
M_{\partial}$. By performing a similar analysis, we can show that (ii)-(viii)
are impossible. Hence, $\varphi\in\mathcal{M}_{0}$, and hence,
$M_{\partial}\subset\mathcal{M}_{0}$. This proves
$M_{\partial}=\mathcal{M}_{0}$.

Let $\omega(\varphi)$ be the omega limit set of the forward orbit
$\gamma^{+}(\varphi):=\\{Q^{n}(\varphi):~{}\forall n\in\mathbb{N}\\}$. We
further have the following claims.

Claim 4. $\cup_{\varphi\in
M_{\partial}}\omega(\varphi)\subset\\{M_{0},M_{1},M_{2}\\}$.

Obviously, there are three possibilities for $\varphi\in
M_{\partial}=\mathcal{M}_{0}:$

$\displaystyle\mbox{\bf Case
1}:~{}~{}\varphi_{1}(\cdot)=\varphi_{2}(\cdot)\equiv
0,~{}~{}\varphi_{3}(\cdot)>0~{}\mbox{or}~{}\varphi_{4}(\cdot)>0;$
$\displaystyle\mbox{\bf Case
2}:~{}~{}\varphi_{3}(\cdot)=\varphi_{4}(\cdot)\equiv
0,~{}~{}\varphi_{1}(\cdot)>0~{}\mbox{or}~{}\varphi_{2}(\cdot)>0;$
$\displaystyle\mbox{\bf Case
3}:~{}~{}\varphi_{1}(\cdot)=\varphi_{2}(\cdot)=\varphi_{3}(\cdot)=\varphi_{4}(\cdot)\equiv
0.$

In what follows, we aim to show that claim 4 holds for each of the above three
cases.

If Case 1 happens, then $I_{1}(t,x,\varphi)=I_{v1}(t,x,\varphi)\equiv 0$ for
all $t\geq 0$ and $x\in\bar{\Omega}$. In view of system (2.3),
$I_{2}(t,x,\varphi),I_{v2}(t,x,\varphi)$ satisfy system (3.1) with $i=2,j=1$.
Since $\mathcal{R}_{2}>1$, it follows from Theorem 3.1 that

$\lim_{t\rightarrow\infty}\big{(}I_{2}(t,x,\varphi)-I_{2}^{*}(t,x)\big{)}=0,~{}~{}\lim_{t\rightarrow\infty}\big{(}I_{v2}(t,x,\varphi)-I_{v2}^{*}(t,x)\big{)}=0$

uniformly for $x\in\bar{\Omega}$. Hence, $\omega(\varphi)=M_{2}$ for any
$\varphi\in M_{\partial}$.

For Case 2, by repeating arguments similar to Case 1, we can show that
$\omega(\varphi)=M_{1}$ for any $\varphi\in M_{\partial}$. For Case 3, one
immediately finds that

$(I_{1}(t,x,\varphi),I_{v1}(t,x,\varphi),I_{2}(t,x,\varphi),I_{v2}(t,x,\varphi))=(0,0,0,0),~{}~{}~{}\forall
t\geq 0,x\in\bar{\Omega}.$

This implies that $\omega(\varphi)=M_{0}$ for any $\varphi\in M_{\partial}$.
Thus claim 4 is obtained.

Claim 5. $M_{1}$ and $M_{2}$ are locally stable, and $M_{0}$ is unstable for
$Q$ in $M_{\partial}$.

Suppose that $\varphi\in M_{\partial}$, we have that
$M_{\partial}=M_{\partial}^{1}\cup M_{\partial}^{2}$, where

$M_{\partial}^{1}:=\\{\varphi\in\partial
X_{0}(0):~{}\varphi_{1}(\cdot)=\varphi_{2}(\cdot)\equiv
0\\},~{}~{}M_{\partial}^{2}:=\\{\varphi\in\partial
X_{0}(0):~{}\varphi_{3}(\cdot)=\varphi_{4}(\cdot)\equiv 0\\}.$

If $\varphi\in M_{\partial}^{1}$, then system (2.3) restricted on
$M_{\partial}^{1}$ is a monotone system. Thus, $M_{2}$ is locally Lyapunov
stable for $Q$ in $M_{\partial}^{1}$ due to [26, Lemma 2.2.1], and $M_{0}$ is
unstable in $M_{\partial}^{1}$. In a similar manner, if $\varphi\in
M_{\partial}^{2}$, we can prove that $M_{1}$ is locally Lyapunov stable for
$Q$ in $M_{\partial}^{2}$, and $M_{0}$ is unstable in $M_{\partial}^{2}$.

The claim 5 implies that no subset of $\\{M_{0},M_{1},M_{2}\\}$ forms a cycle
in $\partial X_{0}(0)$. Based on the above analysis, it follows from the
acyclicity theorem on uniform persistence for maps [26, Theorem 1.3.1 and
Remark 1.3.1] that $Q:X(0)\rightarrow X(0)$ is uniformly persistent with
respect to $(X_{0}(0),\partial X_{0}(0))$ in the sense that there exists
$\tilde{\delta}>0$ such that

$\liminf_{n\rightarrow\infty}d(Q^{n}(\varphi),\partial
X_{0}(0))\geq\tilde{\delta},~{}~{}~{}\forall\varphi\in X_{0}(0).$

By [30, Theorem 4.5] with $\rho(\phi)=d(\phi,\partial X_{0}(0))$, $Q$ admits a
global attractor $A_{0}$ in $X_{0}(0)$, and $Q$ has a fixed point in
$\varphi^{*}\in A_{0}$. Clearly, $u(t,\cdot,\varphi^{*})$ is an
$\omega$-periodic solution of (2.3) and it is strictly positive due to Lemma
2.2.

Finally, we use the arguments in [26, Section 11.2] to obtain the practical
uniform persistence. Since $A_{0}=Q(\omega)A_{0}=Q(A_{0})$, we have
$\varphi_{i}>0,i=1,2,3,4,\forall\varphi\in A_{0}$. Let
$B_{0}:=\cup_{t\in[0,\omega]}Q(t)A_{0}$. Then [26, Theorem 3.1.1] implies that
$B_{0}\subset X_{0}(0)$, and $\lim_{t\rightarrow\infty}d(Q(t)\varphi,B_{0})=0$
for all $\varphi\in X_{0}(0).$ Define a continuous function
$p:X(0)\rightarrow\mathbb{R}_{+}$ by

$p(\varphi):=\min\\{\min_{x\in\bar{\Omega}}\varphi_{1}(x),~{}\min_{x\in\bar{\Omega}}\varphi_{2}(x),\min_{x\in\bar{\Omega}}\varphi_{3}(x),\min_{x\in\bar{\Omega}}\varphi_{4}(x)\\},~{}\forall\varphi\in
X(0).$

Since $B_{0}$ is compact, it follows that $\inf_{\varphi\in
B_{0}}p(\varphi)=\min_{\varphi\in B_{0}}p(\varphi)>0$. Therefore, there exists
a $\delta^{*}>0$ such that for any $\varphi\in X_{0}(0)$,

$\displaystyle\liminf_{t\rightarrow\infty}{\rm
min}_{x\in\bar{\Omega}}I_{i}(t,x,\varphi)\geq\delta^{*},~{}~{}~{}\liminf_{t\rightarrow\infty}{\rm
min}_{x\in\bar{\Omega}}I_{vi}(t,x,\varphi)\geq\delta^{*},~{}~{}~{}i=1,2.$

The proof is complete. ∎

##  5 Asymptotic behavior of $\mathcal{R}_{0}$

In this section, we use the recent theory developed in [38] to study the
asymptotic behavior of the basic reproduction number as the diffusion
coefficients go to zero and infinity. To do this, we write

$\mathcal{R}_{0}(D_{h},D_{v})=\max\\{\mathcal{R}_{1}(D_{h},D_{v}),\mathcal{R}_{2}(D_{h},D_{v})\\}.$

Observe that for each $x\in\bar{\Omega}$, the equation

$\frac{\partial M(t,x)}{\partial t}=\Lambda(t,x)-\eta(t,x)M(t,x)$

admits a globally stable positive $\omega$-periodic solution $M_{0}(t,x)$, and
it is continuous on $\mathbb{R}\times\bar{\Omega}$. Define
$\widetilde{g}(t):=|\Omega|^{-1}\int_{\Omega}g(t,x)dx$. One immediately sees
that the following scalar periodic equation

$\frac{\partial M(t,x)}{\partial
t}=\widetilde{\Lambda}(t)-\widetilde{\eta}(t)M(t,x)$

has a unique positive $\omega$-periodic solution

$\widetilde{M}_{\infty}(t)=\left[\int_{0}^{t}\widetilde{\Lambda}(s)e^{\int_{0}^{s}\widetilde{\eta}(\xi)d\xi}ds+\frac{\int_{0}^{\omega}\widetilde{\Lambda}(s)e^{\int_{0}^{s}\widetilde{\eta}(\xi)d\xi}ds}{e^{\int_{0}^{\omega}\widetilde{\eta}(s)ds}-1}\right]e^{-\int_{0}^{t}\widetilde{\eta}(s)ds},$

which is globally asymptotically stable.

It is easy to verify that assumptions (H1)-(H5) in [38] are valid. An direct
application of [38, Theorems 5.2 and 5.5] leads to

$\lim_{D_{v}\rightarrow
0}\|M^{*}(t,\cdot)-M_{0}(t,\cdot)\|_{\mathbb{Y}}=0,~{}~{}\lim_{D_{v}\rightarrow\infty}\|M^{*}(t,\cdot)-\widetilde{M}_{\infty}(t)\|_{\mathbb{Y}}=0$

hold uniformly on $t\in\mathbb{R}$. For each $x\in\bar{\Omega}$, let
$\\{\Gamma^{i}_{x,0}(t,s):t\geq s\\}~{}(i=1,2)$ be the evolution family on
$\mathbb{R}^{2}$ associated with the following system:

$\left\\{\begin{split}\frac{\partial I_{i}(t,x)}{\partial
t}&=-(d+\gamma_{i})I_{i}(t,x),~{}~{}~{}&&t>s,~{}x\in\bar{\Omega},\\\
\frac{\partial I_{vi}(t,x)}{\partial
t}&=-\eta(t,x)I_{vi}(t,x),~{}~{}~{}&&t>s,~{}x\in\bar{\Omega},\end{split}\right.$

and define

$F^{i}_{0}(t,x)\left(\begin{array}[]{c}\psi_{1}\\\
\psi_{2}\end{array}\right)=\left(\begin{array}[]{c}c_{i}\beta(t,x)\psi_{2}\\\
\frac{\alpha_{i}\beta(t,x)pM_{0}(t,x)}{lN(x)}\psi_{1}\end{array}\right),~{}\forall
t\in\mathbb{R},~{}x\in\bar{\Omega},~{}\psi=(\psi_{1},~{}\psi_{2})\in\mathbb{R}^{2}.$

Let $\\{\widetilde{\Gamma}^{i}_{\infty}(t,s):t\geq s\\}~{}(i=1,2)$ be the
evolution family on $\mathbb{R}^{2}$ of the following system:

$\left\\{\begin{split}\frac{\partial I_{i}(t,x)}{\partial
t}&=-(d+\gamma_{i})I_{i}(t,x),~{}~{}~{}&&t>s,~{}x\in\bar{\Omega},\\\
\frac{\partial I_{v2}(t,x)}{\partial
t}&=-\widetilde{\eta}(t)I_{vi}(t,x),~{}~{}~{}&&t>s,~{}x\in\bar{\Omega},\end{split}\right.$

and define

$\widetilde{F}^{i}_{\infty}(t)\left(\begin{array}[]{c}\psi_{1}\\\
\psi_{2}\end{array}\right)=\left(\begin{array}[]{c}c_{i}\widetilde{\beta}(t)\psi_{2}\\\
\widetilde{f}^{i}_{21}(t)\widetilde{M}_{\infty}(t)\psi_{1}\end{array}\right),~{}~{}~{}\forall
t\in\mathbb{R},~{}\psi=(\psi_{1},~{}\psi_{2})\in\mathbb{R}^{2},$

where

$\widetilde{f}^{i}_{21}(t):=|\Omega|^{-1}\int_{\Omega}\frac{\alpha_{i}\beta(t,x)p}{lN(x)}dx,~{}~{}~{}\forall
t\in\mathbb{R}.$

Let $C_{\omega}(\mathbb{R},\mathbb{R}^{2})$ be the Banach space of all
continuous and $\omega$-periodic functions from $\mathbb{R}$ to
$\mathbb{R}^{2}$, which is endowed with the maximum norm. For each
$x\in\bar{\Omega}$, we respectively define bounded linear positive operators
$L^{i}_{x,0}$ and $\widetilde{L}^{i}_{\infty},i=1,2$, on
$C_{\omega}(\mathbb{R},\mathbb{R}^{2})$ by

$[L^{i}_{x,0}v](t):=\int^{\infty}_{0}\Gamma^{i}_{x,0}(t,t-s)F^{i}_{0}(t-s,x)v(t-s)ds,~{}~{}~{}\forall
t\in\mathbb{R},~{}v\in C_{\omega}(\mathbb{R},\mathbb{R}^{2}),$

and

$[\widetilde{L}^{i}_{\infty}v](t):=\int^{\infty}_{0}\widetilde{\Gamma}^{i}_{\infty}(t,t-s)\widetilde{F}^{i}_{\infty}(t-s)v(t-s)ds,~{}~{}~{}\forall
t\in\mathbb{R},~{}v\in C_{\omega}(\mathbb{R},\mathbb{R}^{2}).$

Then we define $R_{i}(x,0):=r(L^{i}_{x,0}),~{}\forall
x\in\bar{\Omega},~{}i=1,2$ and
$\widetilde{R}_{i}(\infty):=r(\widetilde{L}^{i}_{\infty}),~{}i=1,2$. By [38,
Theorem 4.1] with $\kappa={\rm diag}(D_{h},D_{v}),~{}\chi=D_{v}$ and
$\chi_{0}=0$, and $\kappa={\rm diag}(D_{h},D_{v}),~{}\chi=\frac{1}{D_{v}}$ and
$\chi_{0}=0$, respectively, it follows that

$\lim_{\max(D_{h},D_{v})\rightarrow
0}\mathcal{R}_{i}(D_{h},D_{v})=\max_{x\in\bar{\Omega}}R_{i}(x,0),$
$\lim_{\min(D_{h},D_{v})\rightarrow\infty}\mathcal{R}_{i}(D_{h},D_{v})=\widetilde{R}_{i}(\infty),~{}~{}i=1,2.$

Therefore,

$\lim_{\max(D_{h},D_{v})\rightarrow
0}\mathcal{R}_{0}(D_{h},D_{v})=\max_{x\in\bar{\Omega}}\\{R_{1}(x,0),R_{2}(x,0)\\},$
$\lim_{\min(D_{h},D_{v})\rightarrow\infty}\mathcal{R}_{0}(D_{h},D_{v})=\max\\{\widetilde{R}_{1}(\infty),\widetilde{R}_{2}(\infty)\\}.$

##  6 Numerical simulations

To verify these analytic results and examine the effects of seasonality and
vector-bias on the malaria transmission, we perform illustrative numerical
investigations.

###  6.1 Competitive exclusion and coexistence

We choose the period of our model to be $T=12~{}\mbox{months}$ and concentrate
on one dimensional domain $\Omega=[0,\pi]$. For illustrative purpose, we only
let $\beta(t,x)$ be the time-dependent parameters, given by

$\displaystyle\beta(t,x)=$ $\displaystyle
4\times(5.1492-1.83692\cos(0.523599t)-0.175817\cos(1.0472t)$
$\displaystyle-0.166233\cos(1.5708t)-0.16485\cos(2.0944t)-0.17681\cos(2.61799t)$
$\displaystyle-1.37079\sin(0.523599t)+0.296267\sin(1.0472t)+0.2134\sin(1.5708t)$
$\displaystyle-0.295228\sin(2.0944t)-0.201712\sin(2.61799t))~{}\mbox{month}^{-1},$

which is adapted from [8]. Unless stated otherwise, the baseline parameters
are seen in Table 1. We use the numerical scheme proposed in [34, Lemma 2.5
and Remark 3.2] to compute the reproduction number of each strain. In order to
demonstrate the outcomes of competitive exclusion and coexistence, we consider
the following three cases.

Table 1: Parameter values. Parameter | Value (range) | Dimension | Reference
---|---|---|---
$N(x)$ | 110 | dimensionless | [12]
$M^{*}(t,x)$ | 220 | dimensionless | [24]
$d$ | $1/(72\times 12)$ | month-1 | [12]
$\eta(t,x)$ | 0.8 | month-1 | [24]
$D_{h}$ | 0.4 | km${}^{2}\cdot$ month-1 | [12]
$D_{v}$ | 0.02 | km${}^{2}\cdot$ month-1 | [12]
$p$ | 0.8 (0,1) | dimensionless | [12]
$l$ | 0.2 (0,1) | dimensionless | [12]

Case 1. $\mathcal{R}_{1}>1,\mathcal{R}_{2}>1,\hat{\mathcal{R}}_{1}>1$ and
$\hat{\mathcal{R}}_{2}>1$. We choose
$\gamma_{1}=0.096~{}\mbox{month}^{-1},\gamma_{2}=0.082~{}\mbox{month}^{-1},\alpha_{1}=0.56,\alpha_{2}=0.6,c_{1}=0.25,c_{2}=0.2$.
Then we obtain $\mathcal{R}_{1}=11.1267$, $\mathcal{R}_{2}=10.3022$,
$\hat{\mathcal{R}}_{1}=2.2919$, and $\hat{\mathcal{R}}_{2}=2.2605$. Fig. 1
shows that the disease is uniformly persistent, and periodic oscillation
phenomenon occurs, which is consistent with Theorem 4.3.

Figure 1: Two strains coexist: (a) the sensitive strains $I_{1},I_{v1}$, (b)
The resistant strains $I_{2},I_{v2}$. The initial data are chosen as
$I_{1}(0,x)=6(1+\cos(2x)),~{}I_{v1}(0,x)=10(1+\cos(2x)),~{}I_{2}(0,x)=5(1+\cos(2x)),~{}I_{v2}(0,x)=8(1+\cos(2x)),\forall~{}x\in[0,\pi]$.

Case 2. $\mathcal{R}_{1}>1,\mathcal{R}_{2}>1,\hat{\mathcal{R}}_{1}>1$ and
$\hat{\mathcal{R}}_{2}<1$. We choose
$\gamma_{1}=0.083~{}\mbox{month}^{-1},\gamma_{2}=0.082~{}\mbox{month}^{-1},\alpha_{1}=0.35,\alpha_{2}=0.55,c_{1}=0.2,c_{2}=0.1.$
Then we have $\mathcal{R}_{1}=7.8683$, $\mathcal{R}_{2}=6.9746$,
$\hat{\mathcal{R}}_{1}=1.0934$, $\hat{\mathcal{R}}_{2}=0.4729$. Fig. 2 shows
that the sensitive strains are persistent, but the resistant strains die out.

Figure 2: (a) $I_{1}$ and $I_{v1}$ persist; (b) $I_{2}$ and $I_{v2}$ die out.
The initial data are the same as Fig. 1.

Case 3. $\mathcal{R}_{1}>1,\mathcal{R}_{2}>1,\hat{\mathcal{R}}_{1}<1$ and
$\hat{\mathcal{R}}_{2}>1$. We choose
$\gamma_{1}=0.096~{}\mbox{month}^{-1},\gamma_{2}=0.082~{}\mbox{month}^{-1},\alpha_{1}=0.55,\alpha_{2}=0.45,c_{1}=0.15,c_{2}=0.2.$
Then we get $\mathcal{R}_{1}=8.5413$, $\mathcal{R}_{2}=8.9219$,
$\hat{\mathcal{R}}_{1}=0.7052$, $\hat{\mathcal{R}}_{2}=1.9601$. Fig. 3 depicts
that the resistant strains persist, but the sensitive strains go extinct.

It should be pointed out that in Figs. 1-3, we only plot the graph of
$x$-intersection with $x=0.7448$. In addition, for the second and third case,
the competitive exclusion phenomena are also observed even though
$\mathcal{R}_{0}>1$. It is a pity that we now can not prove it, which is left
for future consideration.

Figure 3: (a) $I_{1}$ and $I_{v1}$ die out; (b) $I_{2}$ and $I_{v2}$ persist.
The initial data are the same as Fig. 1.

###  6.2 Effects of parameters on $\mathcal{R}_{0}$

In order to explore the effect of seasonality, we set the biting rate
$\beta(t)\approx a_{0}(1-b_{0}\cos(0.523599t))$, where $a_{0}$ is the average
biting rate, and $b_{0}\in[0,1]$ is the strength of seasonal forcing. We use
the same parameter values as in Case 1 in Section 6.1. Fig. 4 describes the
dependence of $\mathcal{R}_{0}$ on $a_{0}$ and $b_{0}$. The More precisely,
Fig. 4(a) shows that $\mathcal{R}_{0}$ is an increasing function of $a_{0}$
for fixed $b_{0}$. Fig. 4(b) compares the influences of the time-dependent
biting rate and the time-averaged biting rate on $\mathcal{R}_{0}$. As can be
seen in Fig. 4(b), $\mathcal{R}_{0}$ increases as $b_{0}$ increases. This
implies that the use of the time-averaged biting rate may underestimate the
risk of disease transmission. It should be emphasized that this phenomenon is
not observed in all malaria models, which is dependent on model parameters.

Figure 4: The effect of seasonality on $\mathcal{R}_{0}$. (a)
$\mathcal{R}_{0}$ as a function of $a_{0}$ when $b_{0}=0.35674$; (b)
$\mathcal{R}_{0}$ as a function of $b_{0}$ when $a_{0}=5.1492$.

Next, we investigate the vector-bias effect. We use $q:=l/p$ to measure the
relative attractivity of susceptible host versus infection one. Our numerical
result in Fig. 5 shows that $\mathcal{R}_{0}$ decreases as $q$ increases,
which indicates that the ignorance of the vector-bias effect will
underestimate the value of $\mathcal{R}_{0}$.

Figure 5: The effect of vector-bias on $\mathcal{R}_{0}$.

In fact, we can analytically prove the monotonicity of $\mathcal{R}_{0}$ with
respect to $q$. Let $A^{i}$ and $B^{i}~{}(i=1,2)$ be two bounded linear
operators on $C_{\omega}(\mathbb{R},\mathbb{E})$ given by

$[A^{i}v](t):=\int^{\infty}_{0}\Psi_{i}(t,t-s)v(t-s)ds,~{}[B^{i}v](t):=\mathcal{F}_{i}(t)v,~{}\forall
t\in\mathbb{R},v\in C_{\omega}(\mathbb{R},\mathbb{E}),$

where $\Psi_{i}$ and $\mathcal{F}_{i}(t)$ are defined as in Section 3.
Inspired by Section 4.2 in [37], we write

$A^{i}v=(A^{i}_{1}v_{1},A^{i}_{2}v_{2}),~{}{\rm
and}~{}B^{i}v=(B^{i}_{1}v_{2},B^{i}_{2}v_{1}),~{}~{}~{}\forall
v=(v_{1},v_{2})\in C_{\omega}(\mathbb{R},\mathbb{E}),$

where

$\displaystyle[A^{i}_{1}v_{1}](t)=\int^{\infty}_{0}T_{i}(t,t-s)v_{1}(t-s)ds,~{}~{}[A^{i}_{2}v_{2}](t)=\int^{\infty}_{0}T_{3}(t,t-s)v_{2}(t-s)ds,$
$\displaystyle[B^{i}_{1}v_{2}](t)=c_{i}\beta(t,\cdot)v_{2}(\cdot),~{}~{}~{}[B^{i}_{2}v_{1}](t)=\frac{\alpha_{i}\beta(t,\cdot)pM^{*}(\cdot)}{lN(\cdot)}v_{1}(\cdot),~{}~{}~{}i=1,2.$

According to Section 3.1,
$\mathcal{L}_{i}(q)=A^{i}B^{i}v=(A^{i}_{1}B^{i}_{1}v_{2},A^{i}_{2}B^{i}_{2}v_{1})$,
it then follows that

$\mathcal{L}^{2}_{i}(q)v=(A^{i}_{1}B^{i}_{1}A^{i}_{2}B^{i}_{2}v_{1},A^{i}_{2}B^{i}_{2}A^{i}_{1}B^{i}_{1}v_{2})=\frac{1}{q}\mathcal{L}^{2}_{i}(1)v,~{}~{}~{}i=1,2,$

and hence, $\mathcal{L}^{2}_{i}(q)=\frac{1}{q}\mathcal{L}^{2}_{i}(1)$. In view
of $r^{2}(\mathcal{L}_{i}(q))=r(\mathcal{L}^{2}_{i}(q))$, we obtain

$\mathcal{R}_{i}(q):=r(\mathcal{L}_{i}(q))=\frac{1}{\sqrt{q}}r(\mathcal{L}_{i}(1))=\frac{1}{\sqrt{q}}\mathcal{R}_{i}(1),~{}~{}~{}i=1,2.$

Therefore,
$\mathcal{R}_{0}(q)=\max\\{\mathcal{R}_{1}(q),\mathcal{R}_{2}(q)\\}=\frac{1}{\sqrt{q}}\max\\{\mathcal{R}_{1}(1),\mathcal{R}_{2}(1)\\}$.
This supports our numerical finding.

##  7 Discussion

In this paper, we have proposed a two-strain malaria model with seasonality
and vector-bias. It is of interest to note that our model is a competitive
system for sensitive and resistent strains, but the corresponding subsystem of
each strain is cooperative. To characterize this mathematical structure, we
define a time-dependent region $X(t)$. Although the introduction of time-
varying region brings out some mathematical difficulties, the solution map
$Q(t):X(0)\rightarrow X(t)$ is an $\omega$-periodic semiflow. This nice
property makes us use uniform persistence theory for model dynamics. Our
results show that the zero solution is global attractiveness if
$\mathcal{R}_{0}=\max\\{\mathcal{R}_{1},\mathcal{R}_{2}\\}<1$ (see Theorem
4.1); sensitive (resistent) strains are uniformly persistent if
$\mathcal{R}_{1}>1>\mathcal{R}_{2}~{}(\mathcal{R}_{2}>1>\mathcal{R}_{1})$ (see
Theorem 4.2); and the model is uniformly persistent and admits a positive
periodic solution if $\mathcal{R}_{1}>1,\mathcal{R}_{2}>1$,
$\hat{\mathcal{R}}_{1}>1$ and $\hat{\mathcal{R}}_{2}>1$ (see Theorem 4.3). We
also have analyzed the asymptotic behavior of the basic reproduction number
with small and large diffusion coefficients. Numerically, we have demonstrated
the long-time behaviors of solutions: competitive exclusion and coexistence,
and revealed the influences of some key parameters on the basic reproduction
number. It is found that $\mathcal{R}_{0}$ increases as the strength of
seasonal forcing increases, but it is a decreasing function of the relative
attractivity of susceptible host versus infection one.

Finally, we mention that under certain condition, system (2.3) is a monotone
system with respect to the partial order $\leq_{K}$, which is induced by the
cone $K=\mathbb{E}^{+}\times(-\mathbb{E}^{+})$. Hence, if we can prove the
uniqueness of positive periodic solution in Theorem 4.3, then the positive
periodic solution is globally attractive in $X(0)\setminus\\{0\\}$ by the
virtue of the theory of monotone systems. This is a challenging problem and
left for future study.

## References

  * [1] J. B. Gutierrez, M. R. Galinski, S. Cantrell, et al., From within host dynamics to the epidemiology of infectious disease: Scientific overview and challenges, Math. Biosci. 270 (2015) 143–155.
  * [2] F. Forouzannia, A. B. Gumel, Mathematical analysis of an age-structured model for malaria transmission dynamics, Math. Biosci. 247 (2014) 80–94.
  * [3] The World Health Report 2020. Website: https://www.who.int/news/item/30-11-2020-who-calls-for-reinvigorated-action-to-fight-malaria.
  * [4] R. Ross, The prevention of malaria, 2nd edition, Murray, London, 1911.
  * [5] G. Macdonald, The epidemiology and control of malaria, Oxford University Press, London, 1957.
  * [6] C. Cosner, J. C. Beier, R. S. Cantrell, et al., The effects of human movement on the persistence of vector-borne diseases, J. Theor. Biol. 258 (2009) 550–560.
  * [7] F. Chamchod, N. F. Britton, Analysis of a vector-bias model on malaria transmission, Bull. Math. Biol. 73 (2011) 639–657.
  * [8] Y. Lou, X.-Q. Zhao, A climate-based malaria transmission model with structured vector population, SIAM J. Appl. Math. 70 (2010) 2023–2044.
  * [9] Y. Lou, X.-Q. Zhao, A reaction-diffusion malaria model with incubation period in the vector population, J. Math. Biol. 62 (2011) 543–568.
  * [10] Y. Xiao, X. Zou, Transmission dynamics for vector-borne diseases in a patchy environment. J. Math. Biol. 69 (2014) 113–146.
  * [11] X. Wang, X.-Q. Zhao, A periodic vector-bias malaria model with incubation period, SIAM J. Appl. Math. 77 (2017) 181–201.
  * [12] Z. Bai, R. Peng, X.-Q. Zhao, A reaction-diffusion malaria model with seasonality and incubation period, J. Math. Biol. 77 (2018) 201–228.
  * [13] R. Wu, X.-Q. Zhao, A reaction-diffusion model of vector-borne disease with periodic delays, J. Nonlinear Sci. 29 (2019) 29–64.
  * [14] B.-G. Wang, L. Qiang, Z.-C. Wang, An almost periodic Ross-Macdonald model with structured vector population in a patchy environment, J. Math. Biol. 80 (2020) 835–863.
  * [15] J. G. Kingsolver, Mosquito host choice and the epidemiology of malaria, Am. Nat. 130 (1987) 811–827.
  * [16] G. R. Hosack, P. A. Rossignol, P. van den Driessche, The control of vector-borne disease epidemics, J. Theor. Biol. 255 (2008) 16–25.
  * [17] S. J. Aneke, Mathematical modelling of drug resistant malaria parasites and vector populations, Math. Methods Appl. Sci. 25 (2002) 335–346.
  * [18] E. Y. Klein, Antimalarial drug resistance: a review of the biology and strategies to delay emergence and spread, Int. J. Antimicrob. Ag. 41 (2013) 311–317.
  * [19] F. B. Agusto, A. B. Gumel, P. E. Parham, Qualitative assessment of the role of temperature variations on malaria transmission dynamics, J. Biol. Syst. 23 (2015) 1–34.
  * [20] P. Cailly, A. Tran, T. Balenghiene, et al., A climate-driven abundance model to assess mosquito control strategies, Ecol. Model. 227 (2012) 7–17.
  * [21] M. B. Hoshen, A. P. Morse, A weather-driven model of malaria transmission, Malaria J. 3 (2004), article 32.
  * [22] D. A. Ewing, C. A. Cobbold, B. V. Purse, et al., Modelling the effect of temperature on the seasonal population dynamics of temperate mosquitoes, J. Theor. Biol. 400 (2016) 65–79.
  * [23] S. Altizer, A. Dobson, P. Hosseini, et al., Seasonality and the dynamics of infectious diseases, Ecol. Lett. 9 (2006) 467–484.
  * [24] Y. Shi, H. Zhao, Analysis of a two-strain malaria transmission model with spatial heterogeneity and vector-bias, J. Math. Biol. 82 (2021) 24.
  * [25] F. Li, X.-Q. Zhao, Global dynamics of a reaction-diffusion model of Zika virus transmission with seasonality, Bull. Math. Biol. 83 (2021) 43.
  * [26] X.-Q. Zhao, Dynamical systems in population biology, 2nd edition, Springer, New York, 2017.
  * [27] L. Zhang, Z.-C. Wang, X.-Q. Zhao, Threshold dynamics of a time periodic reaction-diffusion epidemic model with latent period, J. Differential Equations 258 (2015) 3011–3036.
  * [28] D. Daners, P. K. Medina, Abstract Evolution Equations, Periodic Problems and Applications, Pitman Research Notes in Mathematics Series, vol. 279, Longman Scientific and Technical, Harlow, UK, 1992.
  * [29] R. H. Martin, H. L. Smith, Abstract functional differential equations and reaction-diffusion systems, Trans. Amer. Math. Soc. 321 (1990) 1–44.
  * [30] P. Magal, X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM J. Math. Anal. 37 (2005) 251–275.
  * [31] P. Hess, Periodic-Parabolic Boundary Value Problems and Positivity, Pitman Research Notes in Mathematics Series, vol. 247, Longman Scientific and Technical, Harlow, UK, 1991.
  * [32] H. R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM J. Appl. Math. 70 (2009) 188–211.
  * [33] X.-Q. Zhao, Basic reproduction ratios for periodic compartmental models with time delay, J. Dynam. Differential Equations 29 (2017) 67–82.
  * [34] X. Liang, L. Zhang, X.-Q. Zhao, Basic reproduction ratios for periodic abstract functional differential equations (with application to a spatial model for Lyme disease), J. Dynam. Differential Equations 31 (2019) 1247–1278.
  * [35] N. Bacaër, S. Guernaoui, The epidemic threshold of vector-borne diseases with seasonality, J. Math. Biol. 53 (2006) 421–436.
  * [36] N. Tuncer, M. Martcheva, Analytical and numerical approaches to coexistence of strains in a two-strain SIS model with diffusion, J. Biol. Dyn. 6 (2012) 406–439.
  * [37] X. Liang, L. Zhang, X.-Q. Zhao, The principal eigenvalue for degenerate periodic reaction-diffusion systems, SIAM J. Math. Anal. 49 (2017) 3603–3636.
  * [38] L. Zhang, X.-Q. Zhao, Asymptotic behavior of the basic reproduction ratio for periodic reaction-diffusion systems, SIAM J. Math. Anal. 53 (2021) 6873–6909.