Global dynamics analysis of a Zika transmission model with environment transmission route and spatial heterogeneity

1.   Introduction

Zika virus, a mosquito-borne flavivirus, was first isolated in monkeys a rhesus in Uganda in 1947. Later, it was detected in humans in Uganda and the United Republic of Tanzania in 1952 ^[1]. From the 1960s to 1980s, rare sporadic cases of human infections were found throughout Africa and Asia. The first recorded outbreak was reported from the Island of Yap in 2007 ^[2]. In March 2015, Brazil reported a large outbreak of rash illness, soon identified as Zika virus infection, and later found to be associated with Guillain-Barr $\acute{e}$ syndrome and microcephaly ^[3]. On February 1, 2016, WHO declared Zika as a "Public Health Emergency of International Concern" ^[4]. Since the outbreak of Zika in Brazil, the expansion of the Zika outbreak has seemed unstoppable. It spread rapidly from Brazil to northern Europe ^[5], Australia ^[6], through Canada ^[7], the USA ^[8], subsequently, arrived to reach Japan ^[9], China ^[10]. Zika cases have been reported in 90 countries and territories by November 4, 2019 ^[11]. Recently, November 6, 2021, 13 new cases were reported in Uttar Pradesh's Kanpur district, which took the case tally to 79 in the state during the past two weeks. Authorities in the Indian capital region said they were on alert in the wake of a spike in Zika virus cases in the neighbouring state of Uttar Pradesh. It is natural to ask how the previously unknown Zika virus spreads rapidly in the short term.

Since the outbreak of Zika in Brazil, many models have been proposed to study spread, impact, and control of Zika disease and dynamic behaviors. Zhang et al. ^[12] employed a SEIR (Susceptible-Exposed-Infected-Removed)-SEI (Susceptible-Exposed-Infected) human-vector model to estimate the time of first introduction of Zika to Brazil. Zhao et al. ^[13] considered the limited medical resources in Zika model and obtained rich bifurcation phenomena, such as, backward bifurcation, Hopf bifurcation, Bogdanov-Takens bifurcation of codimension 2 and discontinuous bifurcation. Various ordinary differential equations (ODEs) models and dynamics analyses had been applied to study Zika outbreak ^{[14,15,16,17,18,19,20]}. However, the above models ignore the effects of the spatial factors. In fact, the spread of the disease concerns not only the time, but also the spatial location. For this purpose, some researchers began to describe spatio-temporal transmission of Zika disease through partial differential equations (PDEs) ^[21,22]. Miyaoka et al. ^[21] considered spatial movement of humans and vectors and formulated a reaction diffusion model to research the effect of vaccination on the transmission and control of Zika disease. Yamazaki ^[22] added diffusive terms in Zika model in order to capture the movement of human hosts and mosquitoes, considering the unique threat of the sexual transmission of Zika disease. In the above PDEs models, all the coefficients are positive constants. That is, the dynamics of humans and vectors are described in spatially homogeneous environments. However, the diffusion dynamics of the disease is affected by the natural landscapes, the urban and rural distribution, even cultural geographical factors ^[23]. To make the model more consistent with the spread laws of the disease, the spatial heterogeneity must be considered. Hence, it is necessary to understand the transmission dynamics of the Zika disease influenced by the spatial heterogeneity ^[24]. However, in the above studies, the contaminated aquatic environments, an important transmission route for Zika virus, seem to have received little attention.

Recently, experimental result ^[25] demonstrated that $Aedes$ mosquitoes can acquire and transmit Zika virus by breeding in contaminated aquatic environments. It implies that $Aedes$ mosquitoes are infected by Zika virus not only through biting infectious hosts but also through urine excreted by Zika patients. This new transmission route makes the transmission cycle of Zika virus much shorter. It may be one of the major causes of rapid spread of Zika virus in nature. Therefore, it is more reasonable to introduce environment transmission route (That is, human-environment-mosquito-human transmission route) into Zika model ^[26]. However, few Zika models incorporate environment transmission route and spatial heterogeneity simultaneously.

The paper is organized as follows. In Section 2, we propose a novel reaction-diffusion Zika model with environment transmission route in a spatial heterogeneous environment, which is different from all Zika models mentioned earlier. In Section 3, the well-posedness and some properties of the model are also discussed. The basic offspring number $R_{0}^{m}$ and basic reproduction number $R_{0}$ for our spatial model are established in Section 4. In Section 5, by using comparison arguments and the theory of uniform persistence, the threshold dynamics for the model in terms of $R_{0}^{m}$ and $R_{0}$ are analysed. A brief conclusion is given in Section 7.

2.   Model formulation

In this section, we propose a reaction-diffusion Zika model with environment transmission route. Considering the effects of individual mobility, we assume that a host population lives in a spatial heterogeneous environment, which is represented with a bounded domain $\Omega$ with smooth boundary $\partial\Omega$. Mosquitoes are classified in aquatic and adult mosquitoes. Here we combine the egg, larval and pupal stages as one aquatic stage. Aquatic mosquitoes are divided into susceptible and infectious compartments, and their spatial densities at location $x$ and time $t$ are represented by $S_{a}(x, t)$ and $I_{a}(x, t)$, respectively. Adult mosquitoes are divided into susceptible and infectious compartments with spatial densities $S_{m}(x, t)$ and $I_{m}(x, t)$, respectively. We divide the density of total human population at location $x$ and time $t$, denoted by $N_{h}(x, t)$, into three categories: susceptible humans $S_{h}(x, t)$, infectious humans $I_{h}(x, t)$ and recovered humans $R_{h}(x, t)$. So $N_{h}(x, t) = S_{h}(x, t)+I_{h}(x, t)+R_{h}(x, t).$ $V(x, t)$ represents the density at location $x$ and time $t$ of Zika virus within the contaminated aquatic environments. Implication of $V(x, t)$ is similar to avian influenza virus concentration in water ^[27,28].

In this paper, we extend our previous model ^[26] to consider mosquitoes and humans in spatially heterogeneous environments. So, the transmission path of Zika virus is similar to literature ^[26]. In order to incorporate the multiple factors of diffusion and spatial heterogeneity in the spatial domain $\Omega$, we assume that the parameters $K(x)$, $\vartheta(x)$, $\omega(x)$, $\mu_{a}(x)$, $\mu_{m}(x)$, $\Lambda(x)$, $\mu_{h}(x)$, $r(x)$, $a(x)$, $\beta_{1}(x)$, $\beta_{2}(x)$, $\beta_{v}(x)$, $\theta(x)$, $\delta(x)$ are functions of the spatial location $x$ where the contact occurs, and these space dependent parameters are continuous and strictly positive. Mathematically, we assume that all aquatic mosquitoes do not diffuse, and all adult mosquitoes have the same diffusion rate, denoted by $d_m > 0$, while all humans have the same diffusion rate, denoted by $d_h > 0$. The biological meanings of all parameters are shown in Table 1.

On the basis of above assumptions, following the flow diagram in Figure 1, we will focus on the spatiotemporal reaction-diffusion Zika model with environment transmission route as follows:

where $\triangle$ represents the Laplacian operator. The density of total human population $N_{h}(x, t)$ can be determined by the following equation

From Lemma 1 in ^[29], system (2.2) admits a globally attractive positive steady state $H(x)$ in $C(\bar{\Omega}, \; \mathbb{R}_{+})$. For simplicity, we assume that the density of total human population at location $x$ and time $t$ stabilizes at $H(x)$. That is, $N_{h}(x, t) \equiv H(x)$, $\forall\; t \geq0, x\in \Omega.$ Therefore, it suffices to consider the following reduced system:

where $\rho(x) = ls\vartheta(x).$

3.   Well-posedness

In this section, we will study the well-posedness of system (2.3). Let $\mathbb{X}: = C(\bar{\Omega}, \; \mathbb{R}^{7})$ be the Banach space with the supremum norm $\parallel \cdot\parallel$. Define $\mathbb{X_{+}}: = C(\bar{\Omega}, \; \mathbb{R}_{+}^{7})$, then $(\mathbb{X}, \; \mathbb{X_{+}})$ is a strongly ordered Banach space. Let $\mathbb{X}_{K}$ be the subset in $\mathbb{X}$ defined by

In order to simplify notations, we set $u = (u_{1}, u_{2}, u_{3}, u_{4}, u_{5}, u_{6}, u_{7})^{T} = (S_{a}, I_{a}, S_{m}, I_{m}, S_{h}, I_{h}, V)^{T}$, and the initial data satisfies $u^{0} = (u_{1}^{0}, u_{2}^{0}, u_{3}^{0}, u_{4}^{0}, u_{5}^{0}, u_{6}^{0}, u_{7}^{0})^{T} = (S_{a0}, I_{a0}, S_{m0}, I_{m0}, S_{h0}, I_{h0}, V_{0})^{T}$. Throughout, for any $w\in C(\bar{\Omega}, \; \mathbb{R})$, we denote $\overline{w} : = \max\limits_{x\in\bar{\Omega}}w(x)$, $\underline{w} : = \min\limits_{x\in\bar{\Omega}}w.$

We define $T_{i}(t),$ $T_{j}(t),$ $T_{5}(t),$ $T_{6}(t),$ $T_{7}(t)$ : $C(\bar{\Omega}, \; \mathbb{R})\rightarrow C(\bar{\Omega}, \; \mathbb{R})$ as the $C_{0}$ semigroups associated with

subject to the Neumann boundary condition, respectively, $i = 1, \; 2,$ $j = 3, \; 4.$ Then it follows that for any $\psi\in C(\bar{\Omega}, \; \mathbb{R}), \; t\geq 0$,

in which $\Gamma_{i}$, $\Gamma_{j}$, $\Gamma_{6}$, $\Gamma_{7}$ and $\Gamma_{7}$ are the Green functions associated with

subject to the Neumann boundary condition, respectively, $i = 1, \; 2,$ $j = 3, \; 4.$

It is well-known that for all $t > 0$ and $k = 1, \; 2, \; 3, \; 4, \; 5, \; 6, 7,$ $T_{k}$ is compact and strongly positive (see ^[30], Section 7.1 and Corollary 7.2.3). Moreover, $T(t) = (T_{1}(t), \; T_{2}(t), \; T_{3}(t), \; T_{4}(t), \; T_{5}(t), \; T_{6}(t), \; T_{7}(t))^{T}$, $t\geq 0,$ is a $C_{0}$ semigroup. For $\forall\; x\in\partial\Omega$ and $u = (u_{1}, \; u_{2}, \; u_{3}, \; u_{4}, \; u_{5}, \; u_{6}, \; u_{7})^{T}\in\mathbb{X}_{K},$ the nonlinear operator $F = (F_{1}, \; F_{2}, \; F_{3}, \; F_{4}, \; F_{5}, \; F_{6}, \; F_{7})^{T}$ : $\mathbb{X}_{K}\rightarrow \mathbb{X}$ is defined by

Then system (2.3) can be rewritten as the following integral equation

For any $\phi\in\mathbb{X}_{K}$ and $h\geq 0$, then we have

and

This means that

It then follows from Corollary 8.1.3 in ^[31] (see also Corollary 4 in ^[32]), we have the following result.

Lemma 3.1. For every initial value function $\phi\in \mathbb{X}_{K},$ system (2.3) admits a unique mild solution, denoted by

on its maximal existence interval $[0, \; b_{\phi})$ with $u^{0} = \phi$, where $b_{\phi}\leq \infty$.Moreover, $u(\cdot, t, \phi)\in \mathbb{X}_{K}$ for $\forall\; t\in (0, \; b_{\phi})$ and $u(\cdot, t, \phi)$ is a classical solution of system (2.3).

Next, we will show that solutions of system (2.3) exist globally on $[0, \; \infty)$, and admit a global compact attractor on $\mathbb{X}_{K}.$

Lemma 3.2. For every initial value function $\phi\in\mathbb{X}_{K}$, system (2.3) has a unique solution, denoted by

on $[0, \; \infty)$ with $u^{0} = \phi$. Moreover, define the semiflow $\Phi(t):\mathbb{X}_{K}\rightarrow \mathbb{X}_{K}$ associated withsystem (2.3) by

Then the semiflow $\Phi(t):\mathbb{X}_{K}\rightarrow \mathbb{X}_{K}$ admits a global compact attractor on $\mathbb{X}_{K}, \; \forall t\geq 0.$

Proof. Clearly, for $\forall\; \phi\in\mathbb{X}_{K}$, we have $0\leq S_{h}(\cdot, t, \phi), \; I_{h}(\cdot, t, \phi)\leq H(\cdot)$ for all $t \geq 0$. The comparison principle (^[30], Theorem 7.3.4) implies that $S_{h}(x, t, \phi)$ and $I_{h}(x, t, \phi)$ are uniformly bounded and ultimately bounded. It then follows from the seventh equation of (2.3) that

Thus, the comparison principle shows that $V(\cdot, t, \phi)$ is uniformly bounded on $[0, b_{\phi})$, and

More precisely, there exists a $t_{1} > 0$ such that

Letting $N_{1} = \max\left\{\max\limits_{t\in[0, t_{1}], \; x\in\bar{\Omega}}V(x, t), \; 2\frac{\overline{\theta} \overline{H} }{\underline {\delta} }\right\} < \infty$, we deduce

Let $A(x, t) = S_{a}(x, t)+I_{a}(x, t),$ $M(x, t) = S_{m}(x, t)+I_{m}(x, t).$ Then it follows from the first four equations of (2.3) that $(A(x, t), \; M(x, t))$ satisfies

where $A _{0}(x) = \phi_{1}(x)+\phi_{2}(x),$ $M_{0}(x) = \phi_{3}(x)+\phi_{4}(x)$. It is easy to see that there exists a positive vector $\nu = (\nu_{1}, \; \nu_{2}): = \left(\overline{K}, \; \frac{f\overline{\omega}\overline{K}}{\underline {\mu_{m}}}\right)$ such that

Thus, $\nu$ is an upper solution of (3.8). The comparison principle implies that solutions of (3.8) are uniformly bounded on $[0, b_{\phi})$. Hence, so are $S_{a}(x, t), \; I_{a}(x, t)$, $S_{m}(x, t)$ and $I_{m}(x, t)$. Then, we can extend the local unique solution from Lemma 3.1 to global in time via a standard a priori estimates and continuation of local theory. That is, solutions of (2.3) exist on $[0, \; \infty)$. Next, we show that $S_{a}(x, t), \; I_{a}(x, t)$, $S_{m}(x, t)$ and $I_{m}(x, t)$ are ultimately bounded.

From Lemma 3.1, we have $S_{a}(x, t)+I_{a}(x, t)\leq \overline{K}$ for all $x\in \bar{\Omega}, \; t\geq0$. This implies that $S_{a}(x, t)$ and $I_{a}(x, t)$ are ultimately bounded and

It then follows from the third equation of (2.3) that

Consider

From Lemma 1 in ^[29], $\frac{\overline{\omega} \overline{K} }{\underline {\mu_m} }$ is a unique positive steady state that is globally attractive in $C(\bar{\Omega}, \; R_{+})$. Hence there exists $t_{2}\geq 0$ such that $W(x, t)\leq 2 \frac{\overline{\omega} \overline{K} }{\underline {\mu_m} }$. By comparison principle, $S_{m}(\cdot, t)\leq 2\frac{\overline{\omega} \overline{K} }{\underline {\mu_m} }$ when $t > t_{2}$.

Then,

Similarly, $I_{m}(x, t)$ is ultimately bounded. More precisely, there exists a $t_{3} > 0$ such that $I_{m}(\cdot, t)\leq \frac{2\overline{\omega}\overline{K}}{\underline {\mu_m}}\left(1+2 \frac{\overline{a} \overline{\beta_1} }{\underline {\mu_m} } \right),$ for $\forall t > t_{3}$. Thus, we can obtain

where $N_{2} = \max\left\{\max\limits_{t\in[0, t_{3}], \; x\in\bar{\Omega}}I_{m}(x, t), \; \frac{2\overline{\omega}\overline{K}}{\underline {\mu_m}}\left(1+2 \frac{\overline{a} \overline{\beta_1} }{\underline {\mu_m} } \right)\right\} < \infty.$

In addition, since the first two equations and the last equation of (2.3) have no diffusion term, the solution semiflow $\Phi(t)$ is not compact. However, due to $-\omega(x) -\mu_{a}(x) < 0, -\delta(x) < 0, \forall x\in\bar{\Omega},$ using similar arguments from Theorem 4.1 in ^[33] (also Lemma 4.1 in ^[34] and Theorem 2.6 in ^[35]) that the semiflow $\Phi(t):\mathbb{X}_{K}\rightarrow \mathbb{X}_{K}$ has a global compact attractor on $\mathbb{X}_{K}, \; \forall t\geq 0.$ This completes the proof of Lemma 3.2.

The following result shows that the solution of system (2.3) is strictly positive.

Lemma 3.3. Let $(S_{a}(\cdot, t, \phi), I_{a}(\cdot, t, \phi), S_{m}(\cdot, t, \phi), I_{m}(\cdot, t, \phi), S_{h}(\cdot, t, \phi), I_{h}(\cdot, t, \phi), V(\cdot, t, \phi))^{T}$ be the solution of system (2.3) with the initial value $\phi\in\mathbb{X}_{K}$. If there exists some $t_{0}\geq 0$ such that $I_{a}(\cdot, t_{0}, \phi)\not\equiv 0$, $I_{m}(\cdot, t_{0}, \phi)\not\equiv 0$, $I_{h}(\cdot, t_{0}, \phi)\not\equiv 0$, $V(\cdot, t_{0}, \phi)\not\equiv 0$, then the solution of system (2.3)satisfies

Moreover, for any initial value $\phi\in\mathbb{X}_{K}$, there exists some positive constant $\zeta_0$ such that

Proof. For a give $\phi\in\mathbb{X}_{K}$, it is easy to see that $I_{a}(x, t, \phi)$, $I_{m}(x, t, \phi)$, $I_{h}(x, t, \phi)$ and $V(x, t, \phi)$ satisfy

If there exists some $t_{0}\geq 0$ such that $I_{a}(\cdot, t_{0}, \phi)\not\equiv 0$, $I_{m}(\cdot, t_{0}, \phi)\not\equiv 0$, $I_{h}(\cdot, t_{0}, \phi)\not\equiv 0$, $V(\cdot, t_{0}, \phi)\not\equiv 0$, it then follows from the strong maximum principle (see Proposition 13.1 in ^[33]) that $I_{a}(x, t, \phi) > 0, I_{m}(x, t, \phi) > 0, I_{h}(x, t, \phi) > 0, V(x, t, \phi) > 0, \; \; \forall \; t > t_{0}, \; x\in\bar{\Omega}.$

Next, we will prove $S_{a}(x, t, \phi) > 0, \; \forall t > t_{0}, \; x\in\bar{\Omega}.$ To this end, we first show $A(x, t) = S_{a}(x, t)+I_{a}(x, t) < K(x)$ for all $x\in \bar{\Omega}, \; t\geq0$. If not, then there exists $x_1\in \bar{\Omega}, \; t_1\geq0$, such that $A(x_1, t_1) = K(x_1)$. Since $A(x, t)$ satisfies

It implies that $A(x_1, t_1) = 0$. Then $K(x_1) = 0$ which contradicts that $K(x)$ is strictly positive. So, $S_{a}(x, t), I_{a}(x, t) < K(x)$ for all $x\in \bar{\Omega}, \; t\geq0$. Next, we show $S_{a}(x, t, \phi) > 0, \; \forall t > t_{0}, \; x\in\bar{\Omega}.$ Suppose not, there exists $x_2\in \bar{\Omega}, \; t_2 > t_{0}$, such that $S_{a}(x_2, t_2) = 0$. From the first equation of system (2.3), we have

It is contradictory. Thus, $S_{a}(x, t, \phi) > 0, \; \forall t > t_{0}, \; x\in\bar{\Omega}.$ Similarly, we can show $S_{m}(x, t, \phi) > 0, S_{h}(x, t, \phi) > 0, \; \forall t > t_{0}, \; x\in\bar{\Omega}.$

Moreover, for a give $\phi\in\mathbb{X}_{K}$, it follows from the fifth equation of system (2.3) and (3.13) that $S_{h}(x, t, \phi)$ satisfies

According to the comparison principle, we have

This completes the proof of Lemma 3.3.

4.   Threshold index

One of the most important concepts in epidemiology is the basic reproduction number $R_{0}$ which is a threshold index to determine the disease invasion. It is defined to be the average number of secondary cases produced in a completely susceptible population, by a typical infective individual, during its lifetime as infectious. By using the concept of next generation operators, Diekmann, Heesterbeek and Metz ^[36] presented a general approach to $R_{0}$ for autonomous epidemic models. Van den Driessche and Watmough ^[37] gave a computation formula of $R_{0}$ for compartmental models of ordinary differential equations with constant coefficients. For the reaction-diffusion system with spatially dependent coefficients, Wang and Zhao ^[38] defined the basic reproduction ratio as the spectral radius of the next infection operator. According to the above methods, in this section, in order to obtain the basic reproduction number, we should first find the disease free equilibrium (DFE) (infection-free steady state). System (2.3) admits two possible DFEs: $E_{01}(x) = (0, 0, 0, 0, H(x), 0, 0)$ and $E_{02}(x) = (A(x)^{*}, 0, M^{*}(x), 0, H(x), 0, 0).$ $E_{01}(x)$ is characterized by the absence of mosquitoes. $E_{02}(x)$ represents an eradication of Zika in the presence of mosquito population. Firstly, we give the basic offspring number which determines whether mosquito population persists, corresponds to the stability of $E_{01}(x)$.

4.1.   Basic offspring number $R_{0}^{m}$

Since $(A(x, t), \; M(x, t))$ satisfies system (3.8), it suffices to consider system (3.8). Obviously, system (3.8) always admits a DFE $(0, 0).$ Linearizing system (3.8) at $(0, 0),$ we have

The eigenvalue problem associated with (4.1) is as follows

where $\psi_{1}$ and $\psi_{2}$ are both positive for $x\in\Omega$. Let $B_{m}$ be defined as follows

Denote the basic offspring number as $R_{0}^{m}$. According to Lemma 4.2 in ^[38], $\frac{1}{ R_{0}^{m}}$ is the unique positive eigenvalue of the eigenvalue problem

By using the variational characterization of principal eigenvalue ^[39], we can obtain

By Lemma 2.2 and Lemma 2.3 in ^[40], we have the following observation.

Lemma 4.1. Let $\lambda_{m}^{*} = s(B_{m})$ be the spectral bound of $B_{m}$.

(A1) If $\lambda_{m}^{*}\geq 0,$ then $\lambda_{m}^{*}$ is the principal eigenvalue of (4.2) with a strongly positive eigenfunction.

(A2) $R_{0}^{m}-1$ and $\lambda_{m}^{*}$ have the same sign.

Let $\mathbb{Y}_{K}: = \{(A_{0}, M_{0})^{T}\in C(\bar{\Omega}, \mathbb{R}_{+}^{2})\; :\; 0\leq A_{0}(x)\leq K(x)), \; \forall\; x\in \bar{\Omega} \}$. The following result is concerned with the global dynamics of system (3.8).

Lemma 4.2. (^[40], Lemma 2.5) Suppose that $R_{0}^{m} > 1$.System (3.8) admits a unique steady state $(A^{*}(x), M^{*}(x))$ which is globally asymptotically stable in $\mathbb{Y}_{K}\backslash \{(0, 0)\}$. Moreover, $0 < A^{*}(x) < K(x), \; \forall\; x\in \bar{\Omega}.$

Below we use the method proposed in Wang and Zhao ^[38] to introduce the basic reproduction number.

4.2.   Basic reproduction number $R_{0}$

This sub-section is devoted to formulate of the reproduction number for system (2.3) that determines invasion of Zika disease. So, it is the essential condition that guarantees the persistence of mosquito population. From Lemma 4.2, we assume that $R_{0}^{m} > 1$, and then $E_{02}(x)$ exists. Linearizing system (2.3) at $E_{02}(x)$, and then considering only the equations of infective compartments, we have

Let $\bar{T}(t)$ : $C(\bar{\Omega}, \; \mathbb{R}^{4})\rightarrow C(\bar{\Omega}, \; \mathbb{R}^{4})$ be the solution semigroup generated by system (4.6). It is easy to see that $\bar{T}(t)$ is a positive $C_{0}$ semigroup, and its generator $B$ can be written as

Further, $B$ is a closed and resolvent positive operator (see Theorem 3.12 in ^[41]). The eigenvalue problem associated with (4.6) is as follows

By a similar argument as Theorem 7.6.1 in ^[30], we have the following observation.

Lemma 4.3. Let $\lambda^{*} = s(B)$ be the spectral bound of $B$. If $\lambda^{*}\geq 0$, then $\lambda^{*}$ is the principal eigenvalue of the eigenvalue problem (4.7) with a strongly positive eigenfunction.

In the following, we will use the ideas in ^[38] to define the basic reproduction number. Let $\mathbb{T}(t)\; :$ $C(\bar{\Omega}, \; \mathbb{R}^{4})\rightarrow C(\bar{\Omega}, \; \mathbb{R}^{4})$ be the solution semigroup generated by the following linear system

It is easy to see that $\mathbb{T}(t)$ is a $C_{0}$ semigroup on $C(\bar{\Omega}, \; \mathbb{R}^{4})$.

We define

In order to define the basic reproduction number for system (2.3), we assume that the state variables are near DFE $E_{02}$, and introduce the distribution of initial infective individuals described by $\phi(x)\in C(\bar{\Omega}, \; \mathbb{R}^{4}).$ Thus, it is easy to see that $\mathbb{T}(t)\phi(x)$ represents the distribution of those infective individuals at time $t$. Thus, $\mathbb{F}(x)\mathbb{T}(t)\phi(x)$ represents the distribution of new infective individuals at time $t$.

Define $\mathbb{L}\; :\; C(\bar{\Omega}, \; \mathbb{R}^{4})\rightarrow C(\bar{\Omega}, \; \mathbb{R}^{4})$ as follows

It then follows that $\mathbb{L}(\phi)(x)$ represents the distribution of the total new population generated by initial infective individuals $\phi(x)$ during their infection period. So, $\mathbb{L}$ is the next generation operator. We define the spectral radius of $\mathbb{L}$ as the basic reproduction number of system (2.3). That is,

From ^[38], we have the following observation.

Lemma 4.4. $R_{0}-1$ and $\lambda ^{*}$ have the same sign.

The following result indicates that basic offspring number $R_{0}^m$ is a threshold index for eradication or persistence of the Zika disease.

5.   Global dynamic behavior

We firstly focus on the global dynamic behaviors of the DFEs $E_{01}$ and $E_{02}$ of system (2.3).

Theorem 5.1. If $R_{0}^{m} < 1$, then the DFE $E_{01}(x)$ is globally attractive in $\mathbb{X}_{K}$ for system (2.3).

Proof. Assume $R_{0}^{m} < 1$. It follows from Lemma 4.1 that $\lambda_{m}^{*} > 0$. $\lambda_{m }^{*}$ is the principal eigenvalue of eigenvalue problem (4.2) with a strongly positive eigenfunction $(\psi_{1}, \; \psi_{2})$. Since $(A(x, t), \; M(x, t))$ satisfies system (3.8), it follows that

For any given $\phi = (\phi_{1}, \phi_{2})\in \mathbb{Y}_{K}$, there exists some $q > 0$ such that

Note that the linear system (4.1) admits a solution $q e^{\lambda_{m}^{*} t }(\psi_{1}, \; \psi_{2}), \; \forall\; t\geq0$. Then the comparison principle implies that

Hence, $\lim\limits_{t\rightarrow \infty }(A(x, t, \phi), M(x, t, \phi)) = (0, 0),$ uniformly for all $x\in \bar{\Omega}.$ Then, from $A(x, t) = S_{a}(x, t)+I_{a}(x, t)$ and $M(x, t) = S_{m}(x, t)+I_{m}(x, t),$ together with positivity of solutions, for system (2.3), we have for every initial value function $\phi = (\phi_{1}, \; \phi_{2}, \; \phi_{3}, \; \phi_{4}, \; \phi_{5}, \; \phi_{6}, \; \phi_{7})^{T} \in\mathbb{X}_{K}$,

Then, $S_{h}(\cdot, t)$ in system (2.3) is asymptotic to the following system

By the theory for asymptotically autonomous semiflows (see Corollary 4.3 in ^[42]), together with Lemma 1 in ^[29], it follows that

Similarly, $I_{h}(\cdot, t)$ in system (2.3) is asymptotic to the following system

Therefore,

Then $V(\cdot, t)$ in system (2.3) is asymptotic to the following system

Thus,

This completes the proof of Theorem 5.1.

Remark 1. Biologically, Theorem 5.1 shows that the basic offspring number $R_{0}^m$ can be used as a control parameter which determines whether mosquito population is absent or not. It means that mosquito population can be vanished, and the Zika virus will eradicate in human population and contaminated aquatic environment by reducing $R_{0}^m$ below $1$.

Theorem 5.2. Let $u(x, t, \phi)$ be the solution of system (2.3) with $u(\cdot, 0, \phi) = \phi\in\mathbb{X}_{K}$.If $R_{0}^{m} > 1$, $R_{0} < 1$, then the DFE $E_{02}(x)$ is globally attractive for system (2.3). That is, for any $\phi\in\mathbb{X}_{K}$, if $(\phi_{1}, \phi_{3})\not\equiv (0, 0)$, then

Proof. Suppose $R_{0} < 1$. By Lemma 4.4, we have $\lambda^{*} < 0$. So, there exists a sufficiently small positive number $\epsilon_{0}$ such that $\lambda^{*}_{\epsilon_{0}} < 0$, where $\lambda^{*}_{\epsilon_{0}} < 0$ is the principal eigenvalue of the following eigenvalue problem

with a strongly positive eigenfunction $(\varphi_{2}^{\epsilon_{0}}, \varphi_{4}^{\epsilon_{0}}, \varphi_{6}^{\epsilon_{0}}, \varphi_{7}^{\epsilon_{0}})$. It follows from the condition $R_{0}^{m} > 1$, Lemma 4.2 and the positivity of solutions that there exists a $t_{0} > 0$ such that

for all $x\in\bar{\Omega}$, $t\geq t_{0}.$ Hence, by the $I_{a}$, $I_{m}$, $I_{h}$ and $V$ equations of system (2.3), it follows that

For any given $\phi \in \mathbb{X}_{K}$, there exists some $q_{1} > 0$ such that

Note that the following linear system

admits a solution $q_{1} e^{\lambda_{\epsilon_{0}}^{*}(t-t_{0})}(\varphi_{2}^{\epsilon_{0}}, \varphi_{4}^{\epsilon_{0}}, \varphi_{6}^{\epsilon_{0}}, \varphi_{7}^{\epsilon_{0}}), \; \forall\; t\geq t_{0}.$ Then the comparison principle implies that

Hence, $\lim\limits_{t\rightarrow \infty }(I_{a}(x, t, \phi), \; I_{m}(x, t, \phi), \; I_{h}(x, t, \phi), \; V(x, t, \phi)) = (0, 0, 0, 0),$ uniformly for all $x\in \bar{\Omega}.$ Then, $(S_{a}(\cdot, t), S_{m}(\cdot, t))$ in system (2.3) is asymptotic to system (3.8). By the theory for asymptotically autonomous semiflows, together with Lemma 1 in ^[29], it follows that

Similarly, $S_{h}(\cdot, t)$ in system (2.3) is asymptotic to system (5.2). That is,

This completes the proof of Theorem 5.2.

Remark 2. Biologically, Theorem 5.2 shows that mosquito population is present when the basic offspring number $R_{0}^m > 1$. Under this premise, the basic reproduction number $R_{0}$ can be used as a control parameter which determines whether the disease will eventually die out or not. It means that the disease can be eradicated by reducing $R_{0}$ below $1$.

Before giving the disease persistence, we first give the following lemma.

Lemma 5.1. Suppose that $R_m > 1,$ and $\phi_i\equiv 0, i = 2, 4.$ If there exists some $\zeta_1 > 0$ such that

then there exists some $\zeta_2 > 0$ such that

Proof. From $\liminf\limits_{t\rightarrow +\infty}I_{h}(x, t, \phi)\geq \zeta_1, \; uniformly \; for \; all \; x\in \bar{\Omega},$ we have that there exists $t_{11} > 0$ such that

It follows from the last equation of system (2.3) that $V(x, t)$ satisfies

By comparison principle, we have

Thus, there is a $t_{12} > t_{11}$ such that

Due to $R_m > 1$ and $\phi_i\equiv 0, i = 2, 4,$ and from Lemma 4.2, we can obtain that there exists $t_{13} > t_{12}$ such that

From the first equation of system (2.3) and (3.7), one has

Then, we can obtain

Thus, there is a $t_{14} > t_{13}$ such that

From the second equation of system (2.3), we have

Then

which implies that there is a $t_{15} > t_{14}$ such that

Similarly, it follows from the third and fourth equations of system (2.3) that

So, there is a $t_{16} > t_{15}$ such that

By Lemma 3.3, we know

Hence, there is a $t_{17} > t_{16}$ such that

Therefore, letting $\zeta_2 = \frac{1}{3}\max\{\zeta_0, \zeta_1, e_1, e_2, e_3, e_4, e_5\}$, we can obtain that (5.8) holds. This completes the proof of Lemma 5.1.

Theorem 5.3. Let

If $R_{0}^{m} > 1$ and $R_{0} > 1$, then system (2.3) is uniformly persistent, i.e., and there is a constant $\varsigma > 0$ such that, for any initial value $\phi\in\mathbb{M}_{0}$, we can obtain

Proof. The following four steps are taken to prove this result.

Step I $M_0$ is invariant under $\Psi(t).$

For any initial value $\phi\in M_0$, from Lemma 3.3, we can obtain

Then $\Psi(t)\phi\in M_0.$ So $M_0$ is invariant under $\Psi(t).$

Step II For any $\phi \in \partial M_0$, one obtains the $\omega-$ limit set $\omega(\phi)$ = $\{E_{01}\}\cup\{E_{02}\},$ where $\omega(\phi)$ is the omega limit set of the forward orbit $\gamma^+: = \{\Psi(t)\phi:t\geq 0\}.$

Define

For any given $\phi\in\Gamma_{\partial}$, we have $\Psi(t)\phi\in\partial M_0, \forall \; t\geq0.$ That is, for every $t\geq0$, we have

We first consider the case $I_{m}(\cdot, t, \phi) \equiv 0$ for all $t\geq0$. From the sixth equation of system (2.3), we can get that $I_{h}(x, t, \phi)$ satisfies system (5.3). Thus,

From the seventh equation of system (2.3), and according to Corollary 4.3 in ^[42], one has

Similarly, we can obtain

In addition, $S_{h}(\cdot, t, \phi)$ in system (2.3) is asymptotic to system (5.2). Thus,

In the case $I_{m}(\cdot, t, \phi)\not \equiv 0$ for all $t\geq0$, $S_{a}(\cdot, t, \phi)$ and $S_{m}(\cdot, t, \phi)$ are as follows.

(i) $S_{a}(\cdot, t, \phi) \equiv 0$ and $S_{m}(\cdot, t, \phi) \equiv 0$ for all $t\geq0$.

(ii) $S_{a}(\cdot, t, \phi) \equiv 0$ for all $t\geq0$, and $S_{m}(\cdot, t_{21}, \phi)\not\equiv 0$ for some $t_{21} > 0.$

In this case, from the third equation of system (2.3), we can get

(iii) $S_{m}(\cdot, t, \phi) \equiv 0$ for all $t\geq0$, and $S_{a}(\cdot, t_{22}, \phi)\not\equiv 0$ for some $t_{22} > 0.$

In this case, from the first equation of system (2.3), we can obatin

(iv) $S_{a}(\cdot, t_{23}, \phi)\not\equiv 0$ and $S_{m}(\cdot, t_{23}, \phi)\not\equiv 0$ for some $t_{23} > 0.$

In this case, $(S_{a}(x, t, \phi), S_{m}(x, t, \phi))$ in system (2.3) is asymptotic to system (3.8). From Lemma 4.2, we have

Thus, we obtain $\omega(\phi)$ = $\{E_{01}\}\cup\{E_{02}\}.$

Next, we assume $I_{m}(\cdot, t_{24}, \phi)\not \equiv 0$ for some $t_{24} > 0.$ Form Lemma 4.2, one has $I_{m}(\cdot, t, \phi) > 0$ for all $t > t_{24}.$ Then, we get

Here, we only show the case $I_{a}(\cdot, t, \phi) \equiv 0$ for all $t > t_{24}$. It follows from the second equation of system (2.3) that

If $S_{a}(\cdot, t, \phi) \equiv 0, \; for\; all \; t > t_{24},$ then, from the first equation of system (2.3), we have $I_{m}(\cdot, t, \phi) \equiv 0, \; for\; all \; t > t_{24},$ which contradicts our assumption.

If $V(\cdot, t, \phi) \equiv 0, \; for\; all \; t > t_{24},$ then, from the sixth equation of system (2.3), one has $I_{m}(\cdot, t, \phi) \equiv 0, \; for\; all \; t > t_{24},$ which contradicts our assumption. Thus, the step II is proved.

Step III $W^s(\{E_{01}(x)\})\cap M_0 = \emptyset$ and $W^s(\{E_{02}(x)\})\cap M_0 = \emptyset$. In this step, we will show the following two claims.

$Claim 1$ $E_{01}(x)$ is a uniform weak repeller for $M_0.$ That is, there exists $\varepsilon_1 > 0$ such that

$Claim 2$ $E_{02}(x)$ is a uniform weak repeller for $M_0.$ That is, there exists $\varepsilon_2 > 0$ such that

Here, we just prove $Claim 1$. $Claim 2$ can be similarly proven.

If (5.10) does not hold, then

That is, there exists $t_{25} > 0$ such that

Thus, $S_a(x, t, \tilde{\phi})$ and $S_m(x, t, \tilde{\phi})$ satisfy

Consider the following auxiliary linear system

where $\nu = (\nu_1, \nu_2)^T$, and

We know $\tilde{B}_{m}(0) = B_m$, where $B_m$ is defined by (4.3). From Lemma 4.1, if $R_{0}^{m} > 1,$ then $\lambda_m^* = s(B_m) > 0.$ $\tilde{B}_{m}(\varepsilon_1)$ is continuous for small $\varepsilon_1$. So, when $\varepsilon_1$ is small enough, we have $s(\tilde{B}_{m}(\varepsilon_1)) > 0.$ Denote $\lambda_{m\varepsilon_1}^*: = s(\tilde{B}_{m}(\varepsilon_1))$. Obviously, $\lambda_{m\varepsilon_1}^* > 0$.

Let $(\tilde{\varphi}_1, \tilde{\varphi}_2)$ be the positive eigenfunction corresponding to $\lambda_{m\varepsilon_1}^*$. Then auxiliary system (5.14) admits a solution $(\nu_1(x, t), \nu_2(x, t)) = e^{\lambda_{m\varepsilon_1}^*}(\tilde{\varphi}_1, \tilde{\varphi}_2)$. Due to $S_a(x, t, \tilde{\phi}) > 0$, $S_m(x, t, \tilde{\phi}) > 0$ for $\forall\; t\geq t_{25},$ there exists $\varrho_1 > 0$ such that

According to the comparison principle, we can obtain

Since $\lambda_{m\varepsilon_1}^* > 0$, we get

which contradicts with (5.13).

The above discussion implies that $\{E_{01}(x)\}$ and $\{E_{02}(x)\}$ are isolated invariant sets in $M_0$, and $W^s(\{E_{01}(x)\})\cap M_0 = \emptyset$, $W^s(\{E_{02}(x)\})\cap M_0 = \emptyset$. Clearly, every orbit in $\Gamma_{\partial}$ converges to either $E_{01}(x)$ or $E_{02}(x)$, and there are no subsets of $\{ E_{01}(x), E_{02}(x)\}$ forms a cycle in $\Gamma_{\partial}$. From Theorem 4.1 in ^[43], system (2.3) is uniformly persistent if $R^m_{0} > 1$ and $R_{0} > 1$.

Step IV Define a continues function $\bf{g}:\; \mathbb{X}_{K}\rightarrow [0, +\infty)$ with

It follows from Lemma 3.3 that $\bf{g}^{-1}(0, +\infty)\subseteq M_0.$ If $\bf{g}(\phi) > 0$, or $\phi\in M_0$ with $\bf{g}(\phi) = 0$, then $\bf{g}(\Psi(t)\phi) > 0$, $\forall\; t > 0$. That is, $\bf{g}$ is a generalized distance fuction for the semiflow $\Psi(t)$ (see ^[44]). From Lemma 3.2 and according to Theorem 3 in ^[44], there exists a $\varrho_2 > 0$ such that

Thus, $\liminf\limits_{t\rightarrow +\infty}I_{h}(\cdot, t, \phi)\geq \varrho_2, \forall\; \phi\in M_0.$ From Lemma 5.1, there exists some $\varsigma > 0$ such that (5.9) holds. This completes the proof.

Remark 3. Biologically, Theorem 5.3 shows that mosquito population and the disease will persist when the basic offspring number $R_{0}^m > 1$ and basic reproduction number $R_{0} > 1$.

6.   Numerical simulations

In this section, we implement numerical simulations to confirm the analytic results. For the sake of convenience, we concentrate on one dimensional domain $\Omega$, which can be taken, without loss of generality, to be $(0, \pi)$. For the sake of simplicity, we focus on model (2.3) and fix some coefficients and functions as follows: $H(x) = 100$, $K(x) = 500$, $\vartheta(x) = 5$, $l = 1/2$, $\omega(x) = 0.05$, $\mu_{a}(x) = 0.15$, $\mu_{m}(x) = 0.05$, $\mu_{h}(x) = \frac{1}{75\times365}$, $\Lambda(x) = H(x)\times\mu_{h}(x)$, $r(x) = \frac{1}{7}$, $a(x) = 0.3$, $\theta(x) = 0.1$, $\delta(x) = 0.3$, $d_m = 0.01,$ $d_h = 0.1,$

We set initial data

Next, we will change the parameters $s$ and $\beta_v(x)$ and then observe the longtime behavior of the solution to model (2.3).

Example 6.1. Choose $s = 0.04$ and $\beta_{v}(x) = 0.00035(1+\cos(2x))$. It follows from (4.5) that $R_{0}^{m} = 0.52 < 1.$ Theorem 5.1 shows that the disease free equilibrium $E_{01}(x)$ is globally attractive when $R_{0}^{m} < 1$. We can confirm this in Figure 2. Figure 2 shows that the evolutions of $S_{a}(x, t), S_{m}(x, t), I_{a}(x, t), I_{m}(x, t), I_{h}(x, t), V(x, t)$ decay to zero. Biologically, mosquito population will vanish, and Zika virus will eradicate in human population and contaminated aquatic environment.

Example 6.2. Choose $s = 0.4$ and $\beta_{v}(x) = 0.00035(1+\cos(2x))$. It follows from (4.5) that $R_{0}^{m} = 5 > 1$. According to the method of calculating the regeneration numbers in ^[38], we can get $R_{0} = 0.98 < 1$. Theorem 5.2 shows that the disease free equilibrium $E_{02}(x)$ is globally attractive when $R_{0}^{m} > 1$ and $R_{0} < 1$. We can confirm this in Figure 3. Figure 3 shows that the evolutions of $S_{a}(x, t), S_{m}(x, t)$ tend to steady state $A^*(x), M^*(x)$, and the evolutions of diseased compartments $I_{a}(x, t), I_{m}(x, t), I_{h}(x, t), V(x, t)$ decay to zero. Biologically, mosquito population is present, and Zika virus will eradicate in mosquito population, human population and contaminated aquatic environment.

Example 6.3. Choose $s = 0.4$ and $\beta_{v}(x) = 0.035(1+\cos(2x))$. It follows from (4.5) that $R_{0}^{m} = 5 > 1$. According to the method of calculating the regeneration numbers in ^[38], we can get $R_{0} = 2.8 > 1$. Theorem 5.3 shows that system (2.3) is uniformly persistent when $R_{0}^{m} > 1$ and $R_{0} > 1$. We can confirm this in Figure 4. Figure 4 shows that the evolutions of $S_{a}(x, t), S_{m}(x, t)$, $I_{a}(x, t), I_{m}(x, t), I_{h}(x, t), V(x, t)$ keep positive. Biologically, mosquito population is present, and Zika virus will persist in mosquito population, human population and contaminated aquatic environment.

7.   Conclusions

The main contribution of this study, based on experimental proving evidence ^[25], is that we propose a new Zika model by introducing the environment transmission route in a spatial heterogeneous environment. In contrast to ^[12,14,24], we consider environment transmission (human-environment-mosquito-human) route in Zika model. From Figures 3 and 4, we can get that increasing environment transmission rate $\beta_{v}(x)$ can induce the outbreak of Zika. Therefore, the environment transmission route is indispensable. Our work is an extension of previous mathematical Zika models that the transmission of Zika virus involves both mosquito-borne transmission (human-mosquito-human) routes ^[12,14,24]. In fact, environmental transmission route in our model is similar to other waterborne disease models, such as cholera ^[45]. So, our model analysis can also be applied to other waterborne diseases.

We derive a biologically meaningful threshold indexes, the basic offspring number $R_{0}^{m}$ and basic reproduction number $R_{0}$ by the theory developed by Wang and Zhao ^[38], which is characterized as the spectral radius of the next generation operator. Then, we prove that if $R_{0}^{m} < 1$, then both mosquitoes and Zika virus will become vanish. If $R_{0}^{m} > 1$ and $R_{0} < 1$, then mosquitoes will persist and Zika virus will die out. If $R_{0}^{m} > 1$ and $R_{0} > 1$, then mosquitoes and Zika virus are persistent presences. Finally, numerical simulations conform these results.

Our current study has some limitations. In our model, we do not consider sexual transmission route, but it is indeed an important route of spreading of Zika virus ^[46]. Future, we will study a Zika model which incorporates mosquito-borne transmission, sexual transmission and environment transmission threes routes. In addition, in our model, we assume same diffusive coefficients $d_m$ for both $S_m$ and $I_m$, same $d_h$ for all of $S_h$, $I_h$ and $R_h$. However, due to mobility for $S_m$ and $I_m$, $S_h$, $I_h$ and $R_h$ is different, studies of different diffusion coefficients have more realistic implications. Yin ^[47] studied a mathematical model for an infectious disease such as Cholera with different diffusion coefficients. By the delicate theory of elliptic and parabolic equations, global asymptotic behavior of the solution was obtained. This work makes progress toward the case of different diffusive coefficients. Therefore, it motivates us to consider different diffusive coefficients into our model. We leave it for future investigation.

Acknowledgements

The work is partially supported by the Startup Foundation for Introducing Talent of Anhui Polytechnic University (No. 2020YQQ069), Pre-Research National Natural Science Foundation of China of Anhui Polytechnic University (No. Xjky2020134), and Joint Special Funds for Basic Research in Local Undergraduate Universities (Part) of Yunnan Province of China (Nos 2018FH001-113).

Conflict of interest

The authors declare that they have no competing interests.