Dynamics of an impulsive reaction-diffusion mosquitoes model with multiple control measures

1.   Introduction

Mosquito-borne diseases (MBDs) that spread by the bite of the infected mosquitoes, have become increasingly serious worldwide in recent decades. Some well-known such diseases include dengue, malaria, Zika, and so on ^[1,2]. The death of humans caused by MBDs is more than one million annually ^[3]. Prevention and control of MBDs have always been a focus of researchers. So far, there has no therapeutic drug or effective vaccine for the majority of MBDs. Controlling the mosquito population, i.e., reducing or killing the natural mosquitoes, becomes the principal means to prevent the transmission of MBDs. Spraying insecticides was the main method to effectively control mosquitoes in the early years. However, the excessive use of insecticides may lead to serious harm to human health and environment, which stimulates researchers to search for alternative control methods.

A novel and environmentally friendly way, implanting Wolbachia into mosquitoes artificially, is found by researchers to lessen the natural mosquito population and impede mosquitoes from transmitting some MBDs ^[4,5,6,7]. Wolbachia, a maternally transmitted endosymbiotic bacterium, lives in the reproductive organs of its hosts. The hosts' reproductive mechanisms can be interfered by it in diverse ways, such as cytoplasmic incompatibility (CI, the embryos do not hatch when infected males mate with uninfected females) and maternal transmission (Wolbachia still exist in the offsprings of infected females) ^{[8,9,10,11,12]}. In some insect and mosquito species, the imperfect maternal transmission was often observed, which means that the offsprings of infected females may be uninfected ^[13,14]. In particular, the same phenomenon was found in Anopheles stephensi and Aedes aegypti ^[11,15].

To understand the influence of Wolbachia on reducing the natural mosquito population and the infection of MBDs, there are plenty of different mathematical models in recent years ^{[16,17,18,19,20,21,22]}. In 1959, in order to study the impact of CI on mosquito species, Caspari and Watson ^[16] first put forward a discrete-time model and derived a condition for the infection establishment. Yu and Zheng ^[17] proposed a discrete-time model with CI effect and imperfect maternal transmission to study the population replacement. Zheng et al. ^[19,20] established ordinary differential equations to analyze the influence of imperfect maternal transmission on Wolbachia infection. The imperfect maternal transmission results in the rising of natural mosquitoes. However, in these studies, removing the negative influence seems to receive little attention. Therefore, we consider that another control measure (spraying insecticides) should be taken at the same time to get command of natural mosquitoes based on the insights of Xi and Joshi ^[23]. Two approaches are mainly used to spray insecticides: continuously and impulsively. The continuously spraying may lead to the waste of material, manpower and financial resources, which can be cut down by the pulse method. Li and Liu ^[21] constructed a state-dependent impulsive model to discuss the strategies of mosquito extinction or replacing natural mosquitoes with Wolbachia-infected mosquitoes and analyze the integrated mosquito control strategy. However, the state-dependent impulse control is harder to operate than the fixed-time impulse control. As a consequence, we adopt the method of spraying insecticides at fixed time in this paper.

It is worth noting that Wolbachia-infected mosquitoes may suffer fitness costs and fitness advantages ^[12]. Thus, the reproduction rate and mortality rate of Wolbachia-infected mosquitoes are different from those of natural mosquitoes. Mosquitoes compete with each other for food in order to survive, suggesting a density-dependent death rate. Hence, considering these factors mentioned above seems to be more practical. In addition, the laboratory experiments and field demonstrate that the reproduction rate and the mortality rate of both Aedes albopictus and Aedes aegypti are affected by temperature ^[24,25,26], which are functions of time $t$. Accordingly, it seems reasonable and unavoidable to incorporate the seasonality into the model. Moreover, what we should actually do is considering the diffusion of the population and spatial heterogeneity, which are not only important factors in epidemic modeling ^[27], but also affect the evolution of mosquitoes. Hence, it is essential to explore the dynamical behaviours of mosquito population in a diffusive spatially heterogeneous environment.

In the current paper, we employ an impulsive reaction-diffusion model with all the critical factors, including CI effect, imperfect maternal transmission and fitness effect of Wolbachia, fixed-time impulsively spraying of insecticides, and the density-dependent decay rate and spatiotemporal heterogeneity of mosquito population, to study the evolutions of Wolbachia-infected mosquitoes and natural mosquitoes. To the best of our knowledge, few studies seem to incorporate the seasonality, spatial heterogeneity, and impulsively regulating of Wolbachia-infected mosquitoes and natural mosquitoes in the research of the mosquito population, simultaneously. There are two goals in this work. One is to investigate how Wolbachia can be established in natural mosquito species in the context of spatiotemporal heterogeneity and impulsive effect of insecticides. Another is to understand the effectiveness of various control measures (Wolbachia and impulse control of insecticides), and the impacts of seasonality and environmental heterogeneity on the developments of mosquito populations. Our results will throw new insights on mosquito control and the prevention of MBDs.

The paper is organized as follows. In Section 2, we put forward an impulsive reaction-diffusion model with insecticides and seasonality in a heterogeneous environment in terms of Wolbachia-infected mosquitoes and natural mosquitoes. Then the existence and ultimate boundedness of the solution are analyzed in Section 3. Section 4 presents the dynamic behavior of this model. Section 5 takes advantage of the numerical simulations to verify the theoretical conclusions and reflect the effects of some critical factors on the evolution of mosquitoes. In Section 6, we discuss and summarize the obtained results.

2.   The model

Inspired by the remarkable experimental studies ^[28,29], it is assumed that each mosquito only mates once in its life and the sex ratio of Wolbachia-infected mosquitoes and natural mosquitoes is identical. Then, we divide the mosquitoes into two subclasses, Wolbachia-infected mosquitoes and natural mosquitoes. Motivated by the interaction of Wolbachia-infected mosquitoes, natural mosquitoes and insecticides, the reaction-diffusion impulsive hybrid model with seasonality and spatial heterogeneity takes the form as follows:

where $\Omega\in \mathbb{R}^{m}$ ($m\geq1$) and $\mathbb{N}$ are a bounded domain with smooth boundary $\partial \Omega$ and the set of all positive integers, respectively. Here the subclasses $W(t, x)$ and $U(t, x)$ denote the numbers of Wolbachia-infected mosquitoes and natural mosquitoes at time $t$ and position $x$, respectively. The $\Delta$ manifests the Laplacian operator. All parameters are strictly-positive bounded and their biological meanings are listed in Table 1. We suppose that $\rho b_{W}(t, x)\geq d_{W}(t, x)+\delta(t, x)/K(x)$ and $b_{U}(t, x)\geq d_{U}(t, x)+\delta(t, x)/K(x)$ since the natural or Wolbachia-infected mosquitoes do not die out naturally in the wild. A part of the zygote produced by mating natural mosquitoes and Wolbachia-infected mosquitoes do not survive due to CI mechanism. The term $W(t, x)/[U(t, x)+W(t, x)]$ indicates the probability of mating with Wolbachia-infected mosquitoes. Then the loss of offsprings from mating between natural mosquitoes and Wolbachia-infected mosquitoes is $qb_{U}(t, x)U(t, x)W(t, x)/[U(t, x)+W(t, x)]$. $\{t_{n}, n\in \mathbb{N}\}$ is a sequence satisfying $0 = t_{0} < t_{1} < t_{2} < \cdots < t_{n} < \cdots$ and $\lim\limits_{n\rightarrow \infty}t_{n} = +\infty$. Further, at a fixed moment $t_{n}$, the use of insecticides leads that $U(t, x)$ and $W(t, x)$ suffer transient disturbances. The initial values of systems (2.1)–(2.4) and the homogeneous Neumann boundary conditions comply with the form

and

where $\nu$ represents the outward normal unit vector on $\partial\Omega$ and the derivative along $\nu$ to the boundary $\partial\Omega$ is denoted by $\frac{\partial }{\partial \nu}$. Condition (2.6) means that there is no mosquitoes flux crosses $\partial \Omega$.

We introduce some assumptions as follows:

(P1) Functions $b_{W}(t, x), d_{W}(t, x), b_{U}(t, x), d_{U}(t, x), \delta(t, x)$ are bounded positive-valued functions on $\mathbb{R}\times\overline\Omega$, continuously differentiable in $t$ and $x$, and $T$-periodic in $t$ with a period $T > 0$.

(P2) Function $K(x)$ is continuous and bounded positive-valued function on $\overline\Omega$.

(P3) For impulsive sequence $\{t_{n}, n\in \mathbb{N}\}$, we assume that $t_{n+\omega} = t_{n}+T$ for all $n$, where $\omega\in \mathbb{N}$ corresponds to the number of spraying insecticides in a period $T$.

(P4) Sequences $\{1-\mu_{Wn}(x), n\in\mathbb{N}, x\in\Omega \}$ and $\{1-\mu_{Un}(x), n\in\mathbb{N}, x\in\Omega\}$ satisfy $1-\mu_{W(n+\omega)}(x) = 1-\mu_{Wn}(x) > 0$ and $1-\mu_{U(n+\omega)}(x) = 1-\mu_{Un}(x) > 0$ for all $n$, $\omega$ and $x$, where $\omega\in \mathbb{N}$ is the number of spraying insecticides in a period $T$.

3.   The well-posedness

The well-posedness of systems (2.1)–(2.4) is mainly derived in this section. Before stating the main results, we make some notations firstly. For a bounded function $\varphi(t, x)$, $(t, x)\in\mathbb{R}_{+}\times\overline{\Omega}$, we denote $\varphi^{I}: = \inf\limits_{(t, x)\in\mathbb{R}_{+}\times\overline{\Omega}}\varphi(t, x)$, $\varphi^{S}: = \sup\limits_{(t, x)\in\mathbb{R}_{+}\times\overline{\Omega}}\varphi(t, x)$. Define $Y: = C(\overline{\Omega}, \mathbb{R}^{2})$ and $Y^{+}: = C(\overline{\Omega}, \mathbb{R}_{+}^{2})$, here $Y$ is a Banach space of continuous functions from $\overline{\Omega}$ to $\mathbb{R}^{2}$ with the supremum norm $\left\|\cdot\right\|_Y$ and $Y^{+}$ is the positive cone of $Y$.

For $(t, s)\in [0, t_{1}]\times[0, t_{1}]$ ($t_{1}$ is the first fixed impulsive moment after 0), let $\Upsilon_{1}(t, s), \Upsilon_{2}(t, s) :C(\overline{\Omega}, \mathbb{R})\rightarrow C(\overline{\Omega}, \mathbb{R}), t\geq s$ be the evolution operators associated with

and

subject to (2.6), respectively. From [30,Corollary 7.2.3], $\Upsilon_{i}(t, s)$ ($i = 1, 2$) are compact and strongly positive for $(t, s)\in [0, t_{1}]\times[0, t_{1}]$ with $t > s$. Moreover, $\Upsilon(t, s): = {\rm{diag}}{\{\Upsilon_{1}(t, s), \Upsilon_{2}(t, s)}\}:Y \rightarrow Y$, $(t, s)\in [0, t_{1}]\times[0, t_{1}]$ with $t\geq s$, is a semigroup generated by the linear operator $\mathcal{A}(t): = {\rm{diag}}{\{\mathcal{A}_{1}(t), \mathcal{A}_{2}(t)}\}$ defined on the domain $\mathcal{D}(\mathcal{A}(t)) = \mathcal{D}(\mathcal{A}_{1}(t))\times \mathcal{D}(\mathcal{A}_{2}(t))$, where $\mathcal{A}_{i}(t)$ and $\mathcal{D} (\mathcal{A}_{i}(t))$, $i = 1, 2$, are given by

and

After that, define the nonlinear operator $\mathcal{F}:[0, t_{1}]\times Y^{+}\rightarrow Y$ by

where $\psi(x) = (\psi_{1}(x), \psi_{2}(x))^{\mathbb{T}}\in Y^{+},$ $\mathbb{T}$ represents the transpose. Therefore, for $t\in[0, t_{1}]$, systems (2.1)–(2.6) can be transformed as follows:

here $\vartheta = (\vartheta_{1}, \vartheta_{2}) ^{\mathbb{T}} = (W, U)^{\mathbb{T}}$ and $\psi = (\psi_{1}, \psi_{2})^{\mathbb{T}} = (W_{0}, U_{0})^{\mathbb{T}}$.

For any $(t, \psi(x))\in(0, t_{1}]\times Y^{+}$ and small $h > 0$, it can be obtained that

where

When $t\in(0, t_{1}]$, the solutions of systems (2.1)–(2.6) are defined as the solutions of (2.1), (2.2), (2.5) and (2.6) from ^[31]. Impulsive conditions (2.3) and (2.4) display that the functions $W(t_{1}^{+}, \cdot), U(t_{1}^{+}, \cdot)$ are continuously differentiable in $x$ and satisfy homogeneous Neumann boundary conditions. Thus, the solution on $t\in(t_{1}, t_{2}]$ can be derived by letting $(W(t_{1}^{+}, \cdot), U(t_{1}^{+}, \cdot))$ as a new initial function. In the same way, we can obtain the following theorem.

Theorem 3.1. For every initial data $\psi(x)\in Y^{+}$, systems (2.1)–(2.6) admit a unique nonnegative solution $\vartheta(t, x; \psi(x))$ on the maximal existence interval $[0, \tau)$, where $\tau\leq \infty$.

Theorem 3.2. For every $\psi(\cdot)\in Y^{+}$, the solution $\vartheta(t, \cdot; \psi(\cdot))$ of systems (2.1)–(2.6) with $\vartheta(0, \cdot; \psi(\cdot)) = \psi(\cdot)$ exists globally for $t\geq0$ and is ultimately bounded.

Proof. From Theorem 3.1, we know that systems (2.1)–(2.6) have a unique solution $\vartheta(t, \cdot; \psi(\cdot))$ on the interval $[0, \tau)$ with $\vartheta(0, \cdot; \psi(\cdot)) = \psi(\cdot)$. Firstly, according to (2.1), it follows that

Let $M_{W} = \max\limits_{x\in\overline{\Omega}}\mid W_{0}(x)\mid$. Then $W(t, x)\leq \tilde{W}(t)$, $t\in[0, t_{1}], x\in\overline\Omega$ holds with the comparison principle and the uniqueness theorem, where $\tilde{W}(t)$ is the solution of the following ordinary differential equations

So, $W(t_{1}^{+}, x) = (1-\mu_{W1}(x))W(t_{1}, x)\leq (1-\mu_{W1}^{I})\tilde{W}(t_{1})\; {\rm{for}}\; x\in\overline\Omega.$ Likewise, we can conclude that $W(t, \cdot)\leq \tilde{W}(t)$, $t\geq0, t\neq t_{n}^{+}$ and $W(t_{n}^{+}, \cdot) = (1-\mu_{Wn}(\cdot))W(t_{n}, \cdot)\leq (1-\mu_{Wn}^{I})\tilde{W}(t_{n})\; {\rm{for}}\; n\in \mathbb{N}.$ The solution of the following corresponding impulsive model (3.1) is bounded from [32,Lemma 1].

It can be obtained that $W(t, \cdot)\leq K_{W}^{*}, t\geq0$, where $K_{W}^{*} = \max{\{M_{W}, \frac{(\rho b_{W}^{S}-d_{W}^{I})K^{S}}{\delta^{I}(1-\exp(-(\rho b_{W}^{S}-d_{W}^{I})\gamma))}\}}, \gamma = \min_{i = 0, 1, 2, \cdots, \omega}(t_{i+1}-t_{i})$. Accordingly, there have $K_{1} > 0$ and $\tilde{t}_{1} > 0$, such that $W(t, \cdot)\leq K_{1}$, $t\geq \tilde{t}_{1}$, where $K_{1}$ is independent of initial data.

The number of total mosquitoes at time $t$ and position $x$ is denoted by $N(t, x)$. Let $N(t, x) = W(t, x)+U(t, x)$ and $N(t) = \int_{\Omega}N(t, x)dx$, by virtue of systems (2.1)–(2.6), we can get

here, $b^{\max} = \max\{b_{W}^{S}, b_{U}^{S}\}$, $d^{\min} = \min\{d_{W}^{I}, d_{U}^{I}\}$, $\mu_{n}^{\min} = \min\{\mu_{Wn}^{I}, \mu_{Un}^{I}\}$ and $M_{U} = \max\limits_{x\in\overline{\Omega}}\mid U_{0}(x)\mid$. Further, in view of [32,Lemma 1], it is easy to find that there has a $C_{N} > 0$, depending on $N(0)$, so that

To proceed, we show the ultimate boundedness of $U(t, \cdot)$. Due to $U(t_{1}^{+}, \cdot) = (1-\mu_{U1}(\cdot))U(t_{1}, \cdot) \leq U(t_{1}, \cdot)$, hence, for $t\in[t_{1}^{+}, t_{2}]$, it is straightforward to see that the solution of (2.2) and (2.4) is less than or equal to the solution of (2.2). Repeating the analysis above with $U(t_{n}^{+}, \cdot) = (1-\mu_{Un}(\cdot))U(t_{n}, \cdot) \leq U(t_{n}, \cdot)$, $n\in\mathbb{N}$, we can get that the solution of (2.2) with impulsive disturbance is less than or equal to that without impulsive disturbance for $t > 0$. In order to study the ultimate boundedness of term $U(t, \cdot)$, one only needs to verify the following Claim.

Claim. For $\varsigma\geq0$, there exist $B_{2^{\varsigma}} > 0$ which is independent of $\psi(\cdot)\in Y^{+}$, and $\tilde{t}_{2}\geq\tilde{t}_{1}$, such that

In the forthcoming, the method of induction is used to prove (3.3) step by step. For $\varsigma = 0$, (3.3) holds from (3.2). Suppose (3.3) is true for $\varsigma-1$, that is,

Multiplying both sides of (2.2) by $U^{2^{\varsigma}-1}(t, x)$ and integrating over $\Omega$, it is obtained that

By the ultimate boundedness of $W(t, x)$ and Young's inequality, one has

where $\mid\Omega\mid$ is the volume of $\Omega$, $B^{*} = (1-\rho)b_{W}^{S}(K_{1}+1)^{q^{*}}$ and $q^{*} = 2^{\varsigma}/(2^{\varsigma}-1)$. Let $G^{*} = (2^{\varsigma}-1)D_{U}/(2^{2\varsigma-2})$, $H^{*} = B^{*}+b_{U}^{S}$ and $F^{*} = (1-\rho)b_{W}^{S}\mid\Omega\mid$. Hence, (3.5) can be estimated by

Set $\varepsilon = G^{*}/(2H^{*})$. Using the interpolation inequality, there exists a $B^{\varepsilon} > 0$ such that

By virtue of (3.4), then we have

Thus,

Meaning, Claim holds. Then we know that there has a constant $B_{p} > 0$ independent of initial conditions such that $\limsup\limits_{t\rightarrow \infty}\|U(t, \cdot)\|_{p}\leq B_{p}$ for any $p\geq1$. Further, by the same analysis as those in [33,Lemma 2.4], there exist a positive constant $K_{2}$ independent of initial data, and $\tilde{t}_{2}\geq\tilde{t}_{1}$, such that $U(t, \cdot)\leq K_{2}$, $t\geq \tilde{t}_{2}$. It is obvious that $W(t, \cdot)$ and $U(t, \cdot)$ are ultimately bounded.

Therefore, the solution of systems (2.1)–(2.6) with initial data $\vartheta(0, \cdot; \psi(\cdot)) = \psi(\cdot) \in Y_{+}$ globally exists for $t > 0$ and is also ultimately bounded from the comparison argument. This finishes the proof.

Theorem 3.3. Let $\vartheta(t, \cdot; \psi(\cdot))$ be the solution of systems (2.1)–(2.6) satisfying $\psi(\cdot)\in Y^{+}$ on $[0, \infty)$. If there is $\tilde{t}_{0}\geq0$ such that $\vartheta(\tilde{t}_{0}, \cdot; \psi(\cdot))\not\equiv(0, 0)^{\mathbb{T}}$, then $\vartheta(t, \cdot; \psi(\cdot)) > (0, 0)^{\mathbb{T}}$ for $t > \tilde{t}_{0}$.

Proof. From Theorem 3.2, there exists a positive constant $C_{N}$ such that $W(t, x)+U(t, x)\leq C_{N}$, $x\in\overline\Omega$, $t\geq0$. With the aid of systems (2.1)–(2.6), then it is apparent to find that

and

where $t > 0, t\neq t_{n}^{+}, n\in \mathbb{N}, x\in \Omega$. On account of $W(\tilde{t}_{0}, \cdot)\not\equiv0$ and $U(\tilde{t}_{0}, \cdot)\not\equiv0$, by maximum principle, it can be found that $W(t, \cdot) > 0$ and $U(t, \cdot) > 0$ for $t\in(\tilde{t}_{0}, t_{\tilde{n}}]$, where $\tilde{n}$ is the first fixed impulse moment after $\tilde{t}_{0}$. In the same way, for $t\in(t_{\tilde{n}}, t_{\tilde{n}+1}]$, the positiveness of $W(t, \cdot)$ and $U(t, \cdot)$ can be proved with the positiveness of the functions $(1-\mu_{W\tilde{n}}(\cdot))$ and $(1-\mu_{U\tilde{n}}(\cdot))$. Obviously, by using the analogous argument, it can be obtained that $W(t, \cdot) > 0$ and $U(t, \cdot) > 0$ for $t\in(\tilde{t}_{0}, \infty)$. This completes the proof.

4.   Dynamic analysis

In this section, we focus on the dynamics of systems (2.1)–(2.6) by applying comparison principle of differential equations and constructing an appropriate auxiliary function.

4.1.   Extinction and permanence of the mosquitoes

Theorem 4.1. When systems (2.1)–(2.6) satisfy

then $\lim\limits_{t\rightarrow \infty}W(t, x) = 0$ and $\lim\limits_{t\rightarrow \infty}U(t, x) = 0$, $x\in\overline\Omega$.

Proof. First of all, we consider the extinction of Wolbachia-infected mosquitoes. Take into account the following inequality based on (2.1):

for $t > 0, t\neq t_{n}^{+}, n\in \mathbb{N}, x\in \Omega$. Choose a positive constant $M_{1}$ which satisfies $M_{1}\geq W_{0}(\cdot)$ and denote by $\hat{W}(t)$ the solution of the initial value problem $\frac{d \hat{W}(t)}{d t} = \hat{W}(t)(\rho b_{W}^{S}-d_{W}^{I}), \hat{W}(0) = M_{1}.$ By employing the comparison principle, it can be found that $W(t, \cdot)\leq \hat{W}(t)$ for $0\leq t\leq t_{1}$. Further, impulsive condition (2.3) implies that $W(t_{1}^{+}, \cdot) \leq(1-\mu_{W1}^{I})\hat{W}(t_{1})$. Analogously, the corresponding solutions of the following linear system (4.2) with impulse are bounded from below by solutions of (2.1) and (2.3).

It follows from condition (4.1) that all solutions of the impulsive ODE (4.2) tend to zero as $t\rightarrow \infty$. Hence, with (4.1), we can find that $W(t, \cdot)\rightarrow0$ as $t\rightarrow \infty$.

After that, the extinction of the natural mosquitoes is taken into account in event of the elimination of Wolbachia-infected mosquitoes. Combining with (4.1), we know that for any constant $\tilde{\epsilon}$, there exist $0 < \hat{\epsilon}_{1} < \tilde{\epsilon}$, $\tilde{\varrho} > 0$ and $\hat{t}_{1} > 0$ such that $W(t, \cdot) < \hat{\epsilon}_{1}$, $t\geq\hat{t}_{1}$ and

Assume that $M_{2} > 0$, satisfying $M_{2}\geq U_{0}(\cdot)$, is a fixed constant, and $\tilde{U}(t)$ is the solution of the problem $\frac{d \tilde{U}(t)}{d t} = \tilde{U}(t)(b_{U}^{S}-d_{U}^{I})+(1-\rho) b_{W}^{S}\hat{\epsilon}_{1}$ with initial value $\tilde{U}(0) = M_{2}.$ Further, for $t\geq\hat{t}_{1}, t\neq t_{n}^{+}, n\in \mathbb{N}$ and $x\in\overline\Omega$, from (2.3), we yield the following inequality:

By the comparison principle, one has that $U(t, x)\leq \tilde{U}(t)$ for $t\in[\hat{t}_{1}, t_{\hat{n}}], \hat{n}$ is the first fixed impulse moment after $\hat{t}_{1}$. Moreover, we acquire that $U(t_{\hat{n}}^{+}, \cdot) \leq (1-\mu_{U\hat{n}}^{I})\tilde{U}(t_{\hat{n}})$ by using impulsive condition (2.4). Proceeding in this fashion, for $t\geq\hat{t}_{1}$, it is concluded that solutions of (2.2) and (2.4) are bounded from above by the homologous solutions of the following impulsive system

In fact, in order to prove $U(t, \cdot)\rightarrow0$ as $t$ tends to $\infty$, we only need to consider the following three cases for any positive solution $\tilde{U}(t)$ of system (4.4).

Case Ⅰ: There is a $\tilde{t}_{21} \geq \hat{t}_{1}$ such that $\tilde{U}(t)\geq \tilde{\epsilon}$ for all $t\geq \tilde{t}_{21}$.

Let $t = \tilde{t}_{21}+lT$, $l\geq0$ is any positive integer. By integrating system (4.4) from $\tilde{t}_{21}$ to $t$ with (4.3), we can get

Therefore, $\tilde{U}(t)\rightarrow0$ as $l$ tends to $\infty$, which leads to a contradiction.

Case Ⅱ: $\tilde{U}(t)$ is oscillatory about $\tilde{\epsilon}$ for all $t\geq \hat{t}_{1}$.

We can select two sequences $\{\zeta_{j}, j\in\mathbb{N}\}$ and $\{\zeta_{j}^{*}, j\in\mathbb{N}\}$, satisfying $\lim\limits_{j\rightarrow \infty}\zeta_{j} = \lim\limits_{j\rightarrow \infty}\zeta_{j}^{*} = \infty$ and $\hat{t}_{1} < \zeta_{1} < \zeta_{1}^{*} < \cdots < \zeta_{j} < \zeta_{j}^{*} < \cdots$, such that

For any $t\geq\zeta_{1}$, if $t\in(\zeta_{j}, \zeta_{j}^{*}]$ for some integer $j$, then we can choose integer $l$ and constant $0\leq\upsilon < T$ such that $t = \zeta_{j}+lT+\upsilon$. Because $\frac{d \tilde{U}(t)}{d t}\leq \tilde{U}(t)(b_{U}^{S}-d_{U}^{I}+\frac{(1-\rho) b_{W}^{S}\hat{\epsilon}_{1}}{\tilde{\epsilon}}) \; {\rm for\; all}\; t\in(\zeta_{j}, \zeta_{j}^{*}), t\neq t_{n}^{+},$ then integrating this inequality from $\zeta_{j}$ to $t$, we have

where $e = \max\{|b_{U}^{S}-d_{U}^{I}|+(1-\rho) b_{W}^{S}\}, G = \sum\limits_{i = 1}^{p}|\ln(1-\mu_{Ui}^{I})|$. If there exists an integer $j$ such that $t\in(\zeta_{j}^{*}, \zeta_{j+1}]$, then we obviously have $\tilde{U}(t)\leq \tilde{\epsilon} < \tilde{\epsilon}\exp(eT+G)$. Thus, for all $t\geq \zeta_{1}$, we obtain $\tilde{U}(t)\leq\tilde{\epsilon}\exp(eT+G)$ for Case Ⅱ.

Case Ⅲ: there is a $\tilde{t}_{21} \geq \hat{t}_{1}$ such that $\tilde{U}(t)\leq \tilde{\epsilon}$ for all $t\geq \tilde{t}_{21}$.

When Case Ⅲ holds, one yields that $\tilde{U}(t)\leq\tilde{\epsilon}\exp(eT+G)$. In brief, due to the arbitrary of $\tilde{\epsilon}$, it follows that all solutions of (4.4) tend to 0 as $t$ tends to $\infty$. Namely, it is directly gotten under (4.1) that $U(t, \cdot)\rightarrow0$ as $t\rightarrow \infty$. Therefore, the mosquitoes will go to extinction with (4.1). This ends the proof.

Remark 4.1. When $\mu_{W1}^{I} = \mu_{W2}^{I} = \cdots = \mu_{W\omega}^{I}$ and $\mu_{U1}^{I} = \mu_{U2}^{I} = \cdots = \mu_{U\omega}^{I}$, Theorem 4.1 indicates that the number of mosquitoes may eventually reduce to zero if $\omega > \omega_{*}$, $\omega_{*} = \max\{\frac{-T(\rho b_{W}^{S}-d_{W}^{I})}{\ln(1-\mu_{Wi}^{I})}, \frac{-T(b_{U}^{S}-d_{U}^{I})}{\ln(1-\mu_{Ui}^{I})}\}$. It means that mosquitoes will die out when the spraying number of insecticides $\omega$ exceeds critical value $\omega_{1}^{*}$ and the other parameters remain unchanged.

Theorem 4.2. For systems (2.1)–(2.6) with nonnegative nontrivial initial value, if

then $\lim\limits_{t\rightarrow \infty}W(t, \cdot) = 0$ and $\liminf\limits_{t\rightarrow \infty}U(t, \cdot)\geq\sigma$.

Proof. Firstly, we think about that the Wolbachia-infected mosquitoes will become extinct with some conditions. Taking the first inequality of (4.5), it is distinctly testified that $W(t, \cdot)\rightarrow0$ as $t\rightarrow \infty$ via same argument as in the proof of Theorem 4.1.

Afterwards, we consider the evolution of natural mosquitoes under the extinction of Wolbachia-infected mosquitoes. By Theorem 3.3, we know that $W(t, x) > 0$ and $U(t, x) > 0$ for all $t > 0$ and $x\in\overline{\Omega}$ while $W_{0}(x)\geq0$ and $U_{0}(x)\geq0$ that are not identically zero. It is result that the solution separates from zero on $[\varepsilon, \infty)$ with some sufficiently small $\varepsilon > 0$. Therefore, it can be assumed, in general, that $\min\limits_{x\in \overline{\Omega}}W_{0}(x): = m_{1} > 0$ and $\min\limits_{x\in \overline{\Omega}}U_{0}(x): = m_{2} > 0$. From Theorem 3.2 and (4.5), we have that there exist $\tilde{\epsilon}_{11} > 0$ small enough, $\tilde{t}_{12} > 0$ large enough and $\check{\gamma} > 0$ such that $W(t, \cdot) < \tilde{\epsilon}_{11}$, $t\geq\tilde{t}_{12}$ and

Further, for $t\geq\tilde{t}_{12}, t\neq t_{n}^{+}, n\in \mathbb{N}$ and $x\in\overline\Omega$, considering the inequality as follows:

and applying the comparison principle, $U(t, x)\geq \check{U}(t)$ for $t\in[\tilde{t}_{12}, t_{\check{n}}], \check{n}$ is the first fixed impulse moment after $\tilde{t}_{12}$, where $\check{U}(t)$ is the solution of the initial value problem $\frac{d \check{U}(t)}{d t} = \check{U}(t)(b_{U}^{I}-d_{U}^{S}+\frac{\delta^{S}(\check{U}(t)+\tilde{\epsilon}_{11})}{K^{I}}) +b_{U}^{S}\check{U}(t)q, \check{U}(0) = m_{2}.$ Moreover, impulsive condition (2.4) results in $U(t_{\check{n}}^{+}, \cdot) \geq (1-\mu_{U\check{n}}^{S})\check{U}(t_{\check{n}})$. Proceeding in this fashion, $U(t, \cdot)$, the solution of (2.2) and (2.4), is bounded from below by the corresponding solution of the following logistic model with impulse for $t\geq\tilde{t}_{12}$

Actually, system (4.7) possesses a unique strictly positive and piece-wise continuous periodic solution $\check{U}^{*}(t)$ by using [34,Theorem 2.1] and condition (4.6). Every solution $\check{U}(t)$ of system (4.7) with $m_{2} > 0$ possesses the property $\lim\limits_{t\rightarrow \infty}\check{U}(t) = \check{U}^{*}(t)$. Together with $U(t, x)\geq\check{U}(t)$, hence, there exists a positive constant $\sigma$ such that $\liminf\limits_{t\rightarrow \infty}U(t, x)\geq\sigma$. Thus, Wolbachia-infected mosquitoes will be wiped out and natural mosquitoes will permanent with (4.5) when time $t$ is sufficiently large. The proof is finished.

Remark 4.2. Assume $\mu_{W1}^{I} = \mu_{W2}^{I} = \cdots = \mu_{W\omega}^{I}$ and $\mu_{U1}^{S} = \mu_{U2}^{S} = \cdots = \mu_{U\omega}^{S}$. It follows from Theorem 4.2 that Wolbachia-infected mosquitoes will die out and natural mosquitoes are permanent if $\underline{\omega} < \omega < \overline{\omega}$, $\underline{\omega} = \frac{-T(\rho b_{W}^{S}-d_{W}^{I})}{\ln(1-\mu_{Wi}^{I})}, \overline{\omega} = \frac{-T(b_{U}^{I}-d_{U}^{S}-b_{U}^{S}q)}{\ln(1-\mu_{Ui}^{S})}$.

For the sake of simplicity, let

Theorem 4.3. Let $(W, U)^{\mathbb{T}}$ be the solution of systems (2.1)–(2.6) with non-negative nontrivial initial value. If

then $\liminf\limits_{t\rightarrow \infty}W(t, x)\geq\sigma_{1}$ and $\liminf\limits_{t\rightarrow \infty}U(t, x)\geq\sigma_{2}$, $x\in\overline\Omega$, where $\sigma_{1}$ and $\sigma_{2}$ are positive constants.

Proof. From Theorem 3.3, we have that if $W_{0}(\cdot)\geq0$, $U_{0}(\cdot)\geq0$, and $W_{0}(\cdot)\not\equiv0$, $U_{0}(\cdot)\not\equiv0$, then $W(t, \cdot) > 0$ and $U(t, \cdot) > 0$ for all $t > 0$. That implies $W(\varepsilon, \cdot), U(\varepsilon, \cdot)$ separate from zero for some small $\varepsilon > 0$. Therefore, it is relatively reasonable to take $\min_{x\in \overline{\Omega}}W_{0}(x): = m_{1} > 0$, and $\min_{x\in \overline{\Omega}}U_{0}(x): = m_{2} > 0$. Based on the proof of Theorem 3.2 and Young's inequality, it can be obtained that

Then, it is deduced that $U(t, x)\leq K_{U}^{*}$ for $t\geq0, x\in\overline\Omega$ by using the same analysis process as the proof of Theorem 3.2. Subsequently, for $t > 0, t\neq t_{n}^{+}, n\in \mathbb{N}$ and $x\in\overline\Omega$, the following inequality is given by (2.1):

Combining with the condition (4.8) and using the same analysis as $U(t, \cdot)$ in Theorem 4.2, it is easy to yield that $\liminf\limits_{t\rightarrow \infty}W(t, \cdot)\geq\sigma_{1}$. Then, there are $\sigma_{1} > 0$ and $\bar{t}_{1} > 0$, such that $W(t, \cdot)\geq \sigma_{1}$ for $t\geq \bar{t}_{1}$.

Later, with regard to the following inequality:

for $t\geq\bar{t}_{1}, t\neq t_{n}^{+}, n\in \mathbb{N}$, $x\in\overline\Omega$, by the comparison principle, we can conclude that $U(t, x)\geq \bar{U}(t)$ for $t\in[\bar{t}_{1}, t_{\bar{n}}]$, where $\bar{n}$ is the first fixed impulse moment after $\bar{t}_{1}$ and $\bar{U}(t)$ represents the solution of the initial value problem $\frac{d \bar{U}(t)}{d t} = \bar{U}(t)(-d_{U}^{S}-\frac{\delta^{S}(K_{W}^{*}+K_{U}^{*}) }{K^{I}}-b_{U}^{S}q)+(1-\rho) b_{W}^{I}\sigma_{1}, \bar{U}(0) = m_{2}.$ Furthermore, combining with (2.4), an impulsive condition, one gets that $U(t_{\bar{n}}^{+}, \cdot) \geq (1-\mu_{U\bar{n}}^{S})\bar{U}(t_{\bar{n}})$. In a similar way, according to the impulsive equations as follows:

it is achieved that the corresponding solutions of system (4.9) are bounded from above by solutions of (2.2) and (2.4) for $t\geq\bar{t}_{1}$. As a matter of fact, system (4.9) has a unique and strictly positive solution $\bar{U}^{*}(t)$, which is periodic and piece-wise continuous. Then $U(t, \cdot)\geq\bar{U}(t)$ and $\bar{U}(t)\rightarrow \bar{U}^{*}(t)$ as $t\rightarrow \infty$ for any solution $\bar{U}(t)$ of system (4.9). In a words, there has $\sigma_{2} > 0$ such that $\liminf\limits_{t\rightarrow \infty}U(t, \cdot)\geq\sigma_{2}$. As a result, Wolbachia-infected mosquitoes and natural mosquitoes are permanent under conditions. This accomplishes the proof.

Remark 4.3. Theorem 4.3 shows that increasing maternal transmission rate $\rho$ to pass $\rho_{*}$ will ensure the persistence of the Wolbachia-infected mosquitoes in the habitat, where $\rho_{*} = -\frac{1}{T}\sum\limits_{i = 1}^{\omega}\ln(1-\mu_{Wi}^{S}) +d_{W}^{S}+\frac{\delta^{S}K_{U}^{*}}{K^{I}}$. This means that large $\rho$ is better chance for the establishment of Wolbachia in natural mosquito population.

Remark 4.4. When $\mu_{W1}^{S} = \mu_{W2}^{S} = \cdots = \mu_{W\omega}^{S}$, if $\omega < \omega^{*}$, Wolbachia-infected mosquitoes and natural mosquitoes coexist, where $\omega^{*} = \frac{-T(\rho b_{W}^{I}-d_{W}^{S}-\frac{\delta^{S}K_{U}^{*}}{K^{I}})}{\ln(1-\mu_{Wi}^{S})}$.

Remark 4.5. Although the eradication of mosquitoes can be came true by the use of large quantities of highly effective insecticides (see Theorem 4.1 and Figure 2), it is unrealistic because the overuse of insecticides will cause a lot of pollution to the environment on which we live, the mosquito resistance to insecticides, and may cause cancer, nerve paralysis and other human diseases. Therefore, we should use insecticides appropriately to reduce the number of natural mosquitoes under the premise of ensuring the establishment of Wolbachia in natural mosquito population (see Theorem 4.3, Figures 1 and 6). This is exactly what we pay attention to and is one of the main research objectives of this work.

4.2.   Periodic solution

The existence, uniqueness and stability of periodic solution for systems (2.1)–(2.6) are basic and vital problem. Thus, we study this problem by constructing an appropriate auxiliary function in this subsection.

Based on Theorems 3.2 and 4.3, it is obtained that there have $\sigma^{*} > 0$ and $K^{*} > 0$ such that the solutions of systems (2.1)–(2.6) with non-negative nontrivial initial value satisfy $(W(t, \cdot), U(t, \cdot))\in \Pi = \{(W(t, \cdot), U(t, \cdot)):\sigma^{*}\leq W(t, \cdot)\leq K^{*}, \sigma^{*}\leq U(t, \cdot)\leq K^{*} \}$ for sufficiently large $t$.

Theorem 4.4. Assume that the condition of Theorem 4.3 holds. If

then, systems (2.1)–(2.6) satisfying $W_{0}\geq(\not\equiv)0, U_{0}\geq (\not\equiv)0$ has a unique, strictly positive, and global attractive, piece-wise continuous $T$-periodic solution, here $S_{i} = \max\limits_{x\in\overline\Omega}\{(1-\mu_{Wi}(x))^{2}, (1-\mu_{Ui}(x))^{2}\}$, and $\lambda_{M}$ delegates for the maximal eigenvalue of the matrix $E$:

Proof. We will prove this result in three steps. It should illustrate that $b_{W}, d_{W}, b_{U}, d_{U}$, $\delta, K$, $\mu_{Wn}$, and $\mu_{Un}$ are the abbreviations for $b_{W}(t, \cdot), d_{W}(t, \cdot), b_{U}(t, \cdot), d_{U}(t, \cdot)$, $\delta(t, \cdot), K(\cdot)$, $\mu_{Wn}(\cdot)$ and $\mu_{Un}(\cdot)$ for the convenience of marking below, respectively.

Step 1: At first, the periodic solutions for systems (2.1)–(2.6) exist. From Theorems 3.2 and 4.3, it is easy to know that the permanent of systems (2.1)–(2.6) is ensured. Defined operator $\Psi:\Pi\rightarrow \Pi$ by $\Psi(W(t_{0}^{+}), U(t_{0}^{+})) = (W(t_{n}^{+}), U(t_{n}^{+}))$. Then the operator $\Psi$ has at least one fixed point $(W^{*}, U^{*})\in\Pi$ on account of the Brouwer's fixed point theorem. Hence, systems (2.1)–(2.6) have at least one strictly positive and piecewise continuous $T$-periodic solution.

Step 2: Assume that the $(\bar{W}, \bar{U})^{\mathbb{T}}$ is a periodic solution of systems (2.1)–(2.6) and $(W, U)^{\mathbb{T}}$ is another solution in $\Pi$. Next, we need to prove the global attractively of the solution, that is, it is verified that

Constructing the auxiliary function $\mathcal{V}(t) = \int_{\Omega}[(W-\bar{W})^{2}+(U-\bar{U})^{2}]dx,$ the derivative of $\mathcal{V}(t)$ at $t\neq t_{n}^{+}$ is as follows:

After some simple calculations, we can obtain

and

Therefore, it follows from Theorems 3.2 and 4.3 that

Thus, we get $\mathcal{V}(t_{i+1})\leq \mathcal{V}(t_{i}^{+})\exp(\lambda_{M}(t_{i+1}-t_{i}))$ and

Let $S^{*} = \prod\limits_{i = 1}^{\omega} S_{i}\exp(\lambda_{M}T)$. It is obvious to find that $\mathcal{V}(t+T)\leq S^{*}\mathcal{V}(t) = \prod\limits_{i = 1}^{\omega} S_{i}\mathcal{V}(t)\exp(\lambda_{M}T)$. The condition (4.10) implies $S^{*} < 1$. Accordingly, $\mathcal{V}(\varpi T+s)\leq S^{*^{\varpi}}\mathcal{V}(s)\rightarrow0$ as $\varpi\rightarrow \infty$, which means that $\parallel W-\bar{W}\parallel_{L_{2}(\Omega)}\rightarrow0$ and $\parallel U-\bar{U}\parallel_{L_{2}(\Omega)}\rightarrow0$ as $t$ tends to infinity. Due to the boundary of the solution of systems (2.1)–(2.6) in the space $Y$ from [32,Theorem 9], we can get

Therefore, the solutions are globally attractive.

Step 3: It is vital to demonstrate the uniqueness of the periodic solutions. Let us consider the sequence $\{v(nT, v_{0}) = (W(nT, x), U(nT, x)), n\in\mathbb{N}, x\in\Omega\}$. We know that the sequence is compact in the space $Y$ from [32,Theorem 9]. Let $\bar{v}$ be a limit point of this sequence, that is, $\bar{v} = \lim\limits_{\bar{k}\rightarrow \infty}v(n_{\bar{k}}T, v_{0})$. Since $v(T, v(n_{\bar{k}}T, v_{0})) = v(n_{\bar{k}}T, v(T, v_{0}))$ and $\lim\limits_{\bar{k}\rightarrow \infty} v(n_{\bar{k}}T, v(T, v_{0}))-v(n_{\bar{k}}T, v_{0}) = 0$, then

as $\bar{k}\rightarrow \infty$. This implies $\bar{v} = v(T, \bar{v})$. For $\{ v(nT, v_{0}) \}$, the limit point is unique. On the contrary, for the sequence, suppose that there exist two limit points $\bar{v} = \lim\limits_{\bar{k}\rightarrow \infty}v(n_{\bar{k}}T, v_{0})$ and $\tilde{v} = \lim\limits_{\bar{k}\rightarrow \infty}v(n_{\bar{k}}T, v_{0})$. Thus, according to (4.11) and $\tilde{v} = v(n_{\bar{k}}T, \tilde{v})$, then $\parallel \tilde{v}-\bar{v}\parallel_{Y} \leq \parallel \tilde{v}-v(n_{\bar{k}}T, v_{0}) \parallel_{Y}+\parallel v(n_{\bar{k}}T, v_{0})-\bar{v} \parallel_{Y}\rightarrow0$ when $\bar{k}$ tends to $\infty$. Therefore, $\tilde{v} = \bar{v}$, which shows that the solution $(\bar{W}, \bar{U})$ is the unique periodic solution of systems (2.1)–(2.6). This ends the proof of Theorem 4.4.

Corollary 4.1. For systems (2.1)–(2.6), if there are no impulsive controls, the condition (4.8) can be reduced to: $\rho b_{W}^{I}-d_{W}^{S}-\frac{\delta^{S}K_{U}^{*}}{K^{I}} > 0$, which is sufficient to have the permanence for the system.

Corollary 4.2. Suppose that there are no impulsive controls in systems (2.1)–(2.6), if $\rho b_{W}^{I}-d_{W}^{S}-\frac{\delta^{S}K_{U}^{*}}{K^{I}} > 0$ and $\lambda_{M} < 0$, where $\lambda_{M}$ is the same definition as Theorem 4.4, then systems (2.1)–(2.6) without impulsive controls have a unique and strictly positive $T$-periodic solution which is globally attractive.

5.   Numerical simulations

In this section, some numerical simulations will be utilized to derive some important insights about the developments of mosquito populations from our analytic conclusions. More importantly, we will present the effects of some relevant critical factors on the evolution of mosquito species and the influences of the two control policies for reducing the number of natural mosquitoes infected by some viruses.

5.1.   Long term behavior

In this part, some numerical simulations are presented to substantiate the analytic results of Section 4. For simplicity, suppose $\Omega = (0, 6)$. With the reference to ^[35,36], fix $T = 12\; {\rm months}$,

Motivated by recent work ^[37], we assume that the insecticides are sprayed on average $5$ times a month, then $\omega = 60$ and $t_{n} = 0.2n\; {\rm month}$. All the parameters and their values are listed in Table 2. For illustration, we take $b_{W}(t, x) = 0.9b_{U}(t, x)$, $d_{W}(t, x) = 0.85d_{U}(t, x)$ and select the initial data as

where $W_{00} = 1,500,000, U_{00} = 2,500,000$. Applying this set of parameters, we numerically calculate condition (4.8) to obtain

Clearly, as shown in Figure 1, the Wolbachia-infected mosquitoes and natural mosquitoes will coexist. This is coincident with the consequence of Theorem 4.3. That is to say, the partial substitution with Wolbachia is a feasible strategy.

Further, in order to simulate the result of Theorem 4.1, we use $\omega = 432$ and keep other parameter values unchanged. By simple calculations, we acquire that parameters satisfy the conditions of Theorem 4.1, i.e.,

From Figure 2, it is not difficult to find that $W(t, x)$ and $U(t, x)$ go to zero eventually, which implies that the strategy of mosquito extinction can be successful when the effective rate of insecticides unchanged. In the numerical simulation of Theorem 4.2, all parameter values in systems (2.1)–(2.6) are line with the parameters in Figure 1, except for $\rho$ and $q$. Changing $\rho = 0.2$ and $q = 0.35$, then direct computations give that

Figure 3 displays that the Wolbachia-infected mosquitoes become extinct and the natural mosquitoes become permanent, which presents the conclusion of Theorem 4.2. Hence, Wolbachia fails to invade in natural mosquito species.

5.2.   The effects of the parameters related to control measures on the mosquito populations dynamics

In this subsection, aiming to explore control measures, it is necessary to understand the importance of some critical parameters for the development of natural mosquitoes and Wolbachia-infected mosquitoes. Since we take into account the two control measures-Wolbachia and insecticides, then we numerically present the effects of $\omega$, $\mu_{n}(x)$, $q$, $\rho$, $d(t, x)$ and $W_{0}(x)$ at location $x = 3$ in our model, which are summarized in Figure 4. $\mu_{n}(x) = \mu_{Wn}(x)/\mu_{Un}(x)$ expounds the ratio of the effectiveness of insecticides to Wolbachia-infected mosquitoes to the effectiveness of insecticides to natural mosquitoes at location $x$, $d(t, x) = d_{W}(t, x)-d_{U}(t, x)$ expresses the fitness effect of Wolbachia-infected mosquitoes at location $x$ and time $t$. $d(t, x) < 0$ or $d(t, x) > 0$ represents the fitness advantageous or fitness cost, which depends on the Wolbachia strains.

First of all, Figure 4(a)–(d) reflect that the dynamic evolution trend of population size can be affected by the impulsive parameters $\omega$ and $\mu_{n}(x)$. Figure 4(a), (b) illustrate the influence of $\mu_{n}(x)$ for $W(t, x)$ and $U(t, x)$. The yellow, cyan and blue curves indicate that the case of $\mu_{n}(x)\leq1$ (i.e., $\mu_{Wn}(x)\leq\mu_{Un}(x)$), is more in favor of reducing natural mosquito population than $\mu_{n}(x) > 1$ (i.e., $\mu_{Wn}(x) > \mu_{Un}(x)$). That is to say, it is better to achieve that the insecticides response to Wolbachia-infected mosquitoes is less than or equal to the insecticides response to natural mosquitoes. In addition, it is straightforward to find that $W(t, x)$ increases and $U(t, x)$ decreases when $\mu_{n}(x)$ ($\mu_{n}(x)\leq1$) decreases (see the black, purple, yellow and cyan in Figure 4(a), (b)). For example, when Wolbachia-infected mosquitoes are resistant to insecticides (fixing $\mu_{Wn}(x) = 0.6$), if $\mu_{Un}(x)$ heightens from 0.6 to 0.8 by enhancing the dosage of insecticides, i.e., $\mu_{n}(x)$ changes from 1 to 0.75, we see that the highest and lowest values of $U(t, x)$ decline to 94.01% (1172410-1102188) and 47.01% (468964-220438), those of $W(t, x)$ enhance to 102.52% (6660516-6828118) and 102.52% (2664207-2731247), respectively. So, it seems imperative to improve the resistance of Wolbachia-infected mosquitoes and the effectiveness of insecticides for the natural mosquito population with the purpose of establishing a high level of population substitution and reducing the natural mosquitoes. As can be seen from Figure 4(c), (d), it is obvious to survey that the bigger impulsive frequency $\omega$ is, the smaller the numbers of Wolbachia-infected and natural mosquitoes are. While from a realistic point, it is necessary to think about and choose the appropriate insecticide and its spraying times in combination with economic and environmental factors.

Further, Figure 4(e), (f) depict the influences of $q$ and $\rho$ on the solution of systems (2.1)–(2.6). We observe that $W(t, x)$ descends and $U(t, x)$ rises as the decrease of $q$ or $\rho$. What is more, we find that $\rho$ has high sensitivity for $W(t, x)$ and $U(t, x)$. Thus, the vertical propagation and CI effect caused by Wolbachia play an important role in mosquito population. As given in Figure 4(g), (h), it is set to different constants for $d(t, x)$ to evaluate the impact of the fitness effect on $W(t, x)$ and $U(t, x)$. From the blue and red lines, we see that if $d(t, x) < 0$, that is, there has the fitness advantageous for Wolbachia-infected mosquitoes, the smaller $d(t, x)$ is, the larger $W(t, x)$ is and the smaller $U(t, x)$ is. And if Wolbachia leads to fitness cost, i.e., $d(t, x) > 0$, similar conclusions can be obtained (see cyan, purple and green curves). These simulations show that the larger beneficial fitness is or smaller fitness cost is, the higher the level of population substitution is. With the above description, decreasing $d_{W}(t, x)$ or improving $q$, $\rho$ can lead to the increase of $W(t, x)$ and the reduction of $U(t, x)$. Thus, choosing a more effective Wolbachia strain, i.e., satisfying higher of maternal transmission rate, probability of CI effect and resistance for insecticides and lower death rate, is beneficial to obtaining a higher level of population replacement earlier and easier, helping control the mosquitoes and the propagation of MBDs as much as possible.

Besides, we investigate the effects of the initial value of Wolbachia-infected mosquitoes on the numbers of Wolbachia-infected mosquitoes and natural mosquitoes. One can see from Figure 4(i), (j) that $W_{0}(x)$ also plays a significant impact on the evolution of systems (2.1)–(2.6). Obviously, from the green curve, if there is a very small part of Wolbachia-infected mosquitoes initially in the natural mosquitoes, then Wolbachia-infected mosquitoes will be extinct, i.e., there cannot achieve population replacement. If $W_{0}(x)$ increases, the population replacement strategy could be set up (see the comparison of brow and green curves). In addition, all lines except green line in Figure 4(i), (j) show that when $W_{0}(x)$ is larger, the time it takes for the mosquito to first reach a stable periodic oscillation is shorter. Such as, by the purple line and the red line, if $W_{00}$ decreases by 60% (2500000-1000000), the population substitution strategy could remain. Nevertheless, the time for mosquitoes to first reach a stable periodic oscillation over time delays 30 days. The fact analyzed above gives that whether the Wolbachia can invade successfully or not is determined by the initial proportion of Wolbachia-infected mosquitoes. And this shows that the more the increase of the initial Wolbachia-infected mosquitoes, the easier and quicker it is to realize the replacement. We suggest that adequate Wolbachia-infected mosquitoes should be satisfied in nature to make the substitution strategy come true as soon as possible.

5.3.   The effects of control measures on the mosquito populations

In systems (2.1)–(2.6), biological control and chemical control are considered. But is the integrated approach better than the single control? In order to investigate how the integrated control affects the dynamics of the population, we employ numerical simulations to compare situations when the integrated control does occur or not.

In the first place, the impact of Wolbachia on the development of the mosquito population will be shown in the installment. If the Wolbachia-infected mosquitoes are absent ($W_{00} = 0$), the natural mosquitoes are permanent, as shown in Figure 5. From Figure 1, the maximum level and the minimum level of natural mosquitoes are $1.84*10^{6}$ and $4.73*10^{5}$, respectively. It is not different to see from Figure 5 that the peak value of natural mosquitoes reaches about $15.57*10^{6}$ and the lowest increases to $4.16*10^{6}$ approximately when $W_{00}$ declines from 1,500,000 to 0. According to Figures 1 and 5, it follows that the size of natural mosquitoes depends on whether Wolbachia-infected mosquitoes are already present in the population or not. These simulations indicate that, for the purpose of reducing the number of natural mosquitoes, it is a good choice to use comprehensive control measures, rather than only to spray insecticides taking the form of pulse in nature.

On the other hand, we will discuss the impact of impulsively spraying of insecticides. Figure 6 illustrates the evolutions of the two types of mosquitoes when $\mu_{Wn}(x) = \mu_{Un}(x) = 0$. It embodies that the persistence of natural mosquitoes and Wolbachia-infected mosquitoes may still remain if the impulse effects do not exist and this is coincident with Corollary 4.1. By Figures 2 and 6, it is easy to observe that the existence of impulsive control can significantly affect the temporal and spatial dynamics of the system. According to Figure 4(a), (b), we know that $W(t, x)$ will become zero and $U(t, x)$ is permanent if $\mu_{Wn}(x)$ is much greater than $\mu_{Un}(x)$. This is a bad influence for the establishment of Wolbachia in the natural mosquito population and the reduction of natural mosquitoes, which affirm that using the separate control of Wolbachia is better than the combined control strategy sometimes for this case. Next, we consider the case where $\mu_{Wn}(x)$ is less than or equal to $\mu_{Un}(x)$. Comparing Figures 1 and 6, if $\mu_{Wn}(x) = \mu_{Un}(x)$ decreases from 0.5 to 0, the peak values and the least values of Wolbachia-infected mosquitoes and natural mosquitoes increase. It is worth pointing out that the results of the two control measures are much better than that of only releasing Wolbachia-infected mosquitoes in terms of reducing the number of mosquito population if $\mu_{Wn}(x)\leq\mu_{Un}(x)$.

Hence, the aforementioned analysis and the comparison of Figures 1–2 and Figures 4–6 show that compared with single control strategy, combined control strategy is a better choice under $\mu_{Wn}(x)\leq\mu_{Un}(x)$.

5.4.   The effects of the spatiotemporal factors on the mosquito populations dynamics

In this part, in order to understand the role of spatial heterogeneity for the mosquito population, the evolution of mosquitoes distributed in a homogeneous environment is simulated, as shown in Figure 7. In the case with homogeneous environment, the space-averaged parameter defined by $[g(t, x)]: = \frac{1}{L}\int^{L}_{0}g(t, x)dx$, $L = 6$ is adopted, where $g$ represents the coefficients of systems (2.1)–(2.6). The selection value of $g$ is the same with Figure 1. From Figure 7, we find that the systems (2.1)–(2.6) exist spatially homogeneous periodic solutions, i.e., $W(t, x)$ and $U(t, x)$ are distributed uniformly over space at the same time. In contrast, it is easy to see that systems (2.1)–(2.6) generate a solution satisfying heterogeneous in space direction and periodic in time direction from Figure 1. We can find that the dotted lines all lie below the solid lines in boundary position and the situation is opposite in middle position (see Figure 8). This shows that the rural areas and urban areas should take different levels of control strategies to reduce mosquitoes to the same level. With the help of Figure 1 and Figures 7–8, it is clearly observed that spatial heterogeneity makes a noticeable effect on mosquitoes control.

On the other hand, we analyze the impact of seasonality on the dynamic of mosquitoes. Let the time-averaged parameter $[f(t, x)]: = \frac{1}{T}\int^{T}_{0}f(t, x)dt$, $f$ represents the coefficients of systems (2.1)–(2.6), and its selection value is consistent with Figure 1. Figure 9 depicts the numerical results of systems (2.1)–(2.6) with the time-averaged parameters. It is easily seen from Figure 9 that $W(t, x)$ and $U(t, x)$ are in a stationary state, that is, a solution without large periodic in time direction and with heterogeneous in space direction comes into being. Of which, large periodic refers to a period of 12 months. Nevertheless, as shown in Figure 1, $W(t, x)$ and $U(t, x)$ exhibit the large periodicity in time direction and the heterogeneity in space direction. Particularly, from Figure 10, the cyan curves show that solution is in a stationary situation when systems (2.1)–(2.6) is under the time-averaged parameters and the light pink curves show that the solution with the time-periodic parameters decreases first and then increases monotonically with time in a large periodic, which embodies fluctuations over time and is more realistic. Redouble, these indicate that it is hard to find periods of time when mosquitoes reach higher levels if the seasonality is ignored in the study of distribution of the mosquito population. This is not conducive to the rational allocation of resources for mosquito control. Thus, seasonality, a vital factor, cannot be neglected in the study of distribution of the mosquito population.

6.   Discussion

Controlling the vector, the mosquito population, is a very useful and efficient method for keeping humans from infecting some MBDs from area to area. In this work, we consider the integrated control strategy: biological control-Wolbachia and chemical control-insecticides for mosquito population in the spatiotemporal heterogeneous environment. Then, there are some biological questions: (a) Which is a better choice, two controls or a single control? (b) Is it possible to establish Wolbachia in natural mosquitoes under the impulsively spraying of insecticides in a spatially non-homogeneous environment? And What strain of Wolbachia is suitable to choose? (c) Whether or not the spatial heterogeneity and seasonality have an impact on the development of mosquitoes? Thus, it is worth making an inquiry into the comprehensive effects of integrated control, spatial heterogeneity and seasonal variation in the study of vector control.

To address the above questions, we propose and analyze a reaction-diffusion mosquitoes model with impulse effects and seasonality in this paper, which expounds the interactions among Wolbachia-infected mosquitoes, natural mosquitoes and insecticides in the spatially heterogeneous environment. The global existence and ultimately bounded of solutions (see Theorem 3.2) and the complex dynamic results of the systems (2.1)–(2.6) are established. Specifically, the mosquitoes will be died out when condition (4.1) holds. That is, the strategy of mosquito eradication could be set up (see Theorem 4.1 and Figure 2). And it is proved that if the sufficient condition (4.5) holds, Wolbachia-infected mosquitoes become extinction and natural mosquitoes are permanent. Meaning, the invasion of natural mosquito population by Wolbachia may not be successful (see Theorem 4.2 and Figure 3). Furthermore, we also certify the two types of mosquitoes are persistent by applying the comparison principle and condition (4.8), which implies that the strategy of mosquito partial substitution can come true (see Theorem 4.3 and Figure 1). In addition, the model has a unique, strictly positive, piecewise continuous and globally attractive $T$-periodic solution under certain conditions proved by an appropriate auxiliary function (see Theorem 4.4). Intuitively, it seems that Wolbachia-infected mosquitoes and natural could simultaneously persist when insecticides, spatial structure and seasonality are taken into account.

By Figure 4, firstly, we know that the smaller $q$, $\rho$ are, or the greater $d_{W}(t, x)$ is, the lower the number of Wolbachia-infected mosquitoes is. As a result, we suggest that Wolbachia-infected mosquitoes, which carry bacteria with high maternal transmission rate, probability of CI effect and resistance for insecticides, as well as low fitness cost, should be cultivated. And the larger the extra releases are at the initial moment, the easier and earlier the Wolbachia persistence is. Secondly, the number of impulse $\omega$ in a period and impulsive perturbations $\mu_{Wn}(x)$, $\mu_{Un}(x)$ have the energetic effect for the control of mosquitoes. In order to insure the coexistence of Wolbachia-infected mosquitoes and natural mosquitoes and reduce the number of natural mosquitoes at the same time, insecticides and the spraying times should be chosen based on the actual conditions and requirements. Furthermore, in view of the contrasting of Figures 1–2 and Figures 4–6, we suggest that the comprehensive control strategy should be considered to put into effect when the resistance of Wolbachia-infected mosquito population for insecticides is higher than those of natural mosquito population. In addition, for mosquito population, spatial heterogeneity and seasonality have a noticeable impact on the development and control of mosquitoes, which are illustrated in Figures 7–10. In summery, though our work is not perfect, our theoretical and numerical results can provide helpful information for controlling mosquitoes and MBDs.

There are other factors affecting mosquitoes control. For instance, it takes time for mosquitoes to evolve from egg state to adult, then considering the maturation delay in mathematical model is more realistic. How to study the impulsive reaction-diffusion system with a time delay? This is a great challenge and interesting work. On the other hand, MBDs are mainly transmitted between humans and mosquitoes by the bite of adult females, therefore it is necessary to extend the model by taking into account humans. We leave these issues for future investigation.

Acknowledgments

We thank the editors and reviewers for their comments and suggestions, which help improve the presentation of the paper. The research is supported by Natural Science Foundation of China (No. 11971013), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (No. KYCX20_0169) and Nanjing University of Aeronautics and Astronautics PhD short-term visiting scholar project (No. ZDGB2021026) at the University of Alberta.

Conflict of interest

All authors declare no conflicts of interest in this paper.