Transmission dynamics of Zika virus incorporating harvesting

1.   Introduction

The Zika virus (ZIKV) infection, a vector-borne disease carried by Aedes africanus, is caused by the sting of Aedes aegypti. ZIKV was first discovered in Uganda in 1947 ^[1]. In 2007, it was reported ZIKV occurred on Yap Island (Federated States of Micronesia). It then spread rapidly to Asia, Africa, the United States and Brazil ^[2,3].

Aedes aegyptis (or the yellow fever mosquitos) are the main source of ZIKV transmission and the cause of dengue infection as well. The ZIKV infected by humans spreads through the bite of infected Aedes aegyptis. If one partner of a couple is infected with ZIKV, the virus can be transmitted through unprotected sexual activities. People infected with ZIKV would display mild symptoms initially because they feel uncomfortable inside and then might develop or trigger serious conditions ^[4].

As it is known to all, mathematical modelling plays an important role in epidemiology. It is an important tool to analyze and understand the spread of virus infection in a community. It can assist to strengthen health policies to control or avoid the outbreak of infection. Among the mathematical models of infectious disease, compartment model ^[5,6] is the most common. For research on ZIKV, which is considered as a mosquito-borne virus, the established compartment model classifies not only the population but also the mosquitoes. Funk et al. ^[7] considered the incubation period of human ZIKV. An SEIR model was formulated for the population and an SEI one was formulated for mosquitoes. They made full use of the commonalities of dengue and ZIKV, compared with two dengue and ZIKV outbreaks in different island environments (Yappa Main Island and Fayes) in Micronesia. They found that the proportion of reported ZIKV cases was smaller than that of dengue cases. Kucharski et al. ^[8] applied a mathematical compartment model similar to the one above to analyze the outbreaks in the six major archipelagos of French Polynesia during 2013–2014. Ndaïrou et al. ^[9] established a compartment model to simulate the mother-to-child transmission and the spread of ZIKV in Brazil. Considering seasonal effects, Suparit et al. ^[10] formulated a compartment model with a time-dependent mosquito biting rate and used a computational parameter estimation algorithm to estimate the value of unknown insect parameters.

In addition to the study of the infection dynamics of ZIKV, focus should also be given to how to control the disease. At present, in addition to the use of insecticide spray control vectors and destruction of larval breeding grounds, there is currently no established ZIKV treatment. So, if humans have targeted a series of killings on mosquitoes, does this act influence the reduction or eradication of infectious diseases? Harvesting is defined as the continuous removal (via killing or capture) of a population species ^[11]. In the harvesting model of population dynamics, the emphasis is on the killing or capture of a specific species by humans. So, how to kill mosquitoes? From a human intervention point of view, there are two main types of methods: one is the use of synthetic insecticides, such as pyrethroids and organophosphates; the other is the use of growth blocking techniques, such as the application of the larvicidal bacterium Bti (Bacillus thuringiensis var. israelensis), pyriproxyfen, larvicidal oil and so on. A study on managing Aedes aegypti populations in the first Zika transmission zones in the continental United States ^[12] shows researchers have used some practical control methods aimed at killing mosquitoes including aerial and truck sprays of adulticides and larvicides. They found Bti larvicide greatly depressed urban Aedes aegypti populations when applied weekly. Wang et al. ^[13] observed the efficacy of control agents against small larvae, large larvae, and pupae of Aedes aegypti to determine an appropriate larvicide regime to employ in emergency dengue control programs. They concluded that larvicides that kill the pupal stage (Aquatain AMF or larvicidal oil) should be included in the effectively interrupting the dengue transmission program in addition to Bti, pyriproxyfen, or temephos. Cornel et al. ^[14] believed long-term high density placements of Autocidal-Gravid-Ovitraps (AGO-B trips) could be used as an environmentally friendly trap-kill control strategy. From the studies on mosquitoes, one can find that there are methods which are used (such as BG Sentinel traps and AGO-B traps) to surveil Ae. aegypti in locations where is in high-risk. Therefore, it is possible to achieve an approximate quantitative capture of Ae. aegypti.

We find that some researchers have applied the idea of hunting to disease control ^[15,16,17]. Yusof et al. ^[16] applied the idea of population harvesting to the control of Hantavirus spreading. They used three harvesting methods, one for constant harvesting, one for proportional harvesting, and another for seasonal harvesting. From numerical simulation results, it was found that harvesting did not eliminate the disease but could reduce the spread of the disease and the proportional-based harvesting was the most significant. One can find that seasonal harvesting is actually a periodic concussion in a neighborhood of constant harvesting. In a large time range, this is in fact similar to that of the constant. In this paper, we focus on the first two harvesting models.

In this paper, we analyze an ZIKV-compartment model and investigate the transmission dynamics of ZIKV with two types of harvesting: proportional harvesting and constant harvesting. Firstly, the existence and stability of the equilibrium point of the model are calculated. Next, the relevant numerical simulations are conducted to verify the conclusions obtained before. Finally, the impact of the harvesting on the spread of disease is discussed.

2.   Model derivations

Thanks to the insightful work of Bonyah ^[18]. Divide the human population into four sub-classes, susceptible humans ${S_H}(t)$, exposed humans ${E_H}(t)$, infected humans ${I_H}(t)$ and recovered humans ${R_H}(t)$. The total human population is represented as

Similarly, ${N_M}(t)$ is the total number of mosquitoes which is partitioned into susceptible mosquitoes ${S_M}(t)$, exposed mosquitoes ${E_M}(t)$ and infected mosquitoes ${I_M}(t)$. Hence

The compartmental mathematical model is given by the following system:

where all the parameters are positive, which are described in Table 1 ^[18].

Incorporating the factor of harvesting to control the ZIKV results in a model given by the following system:

where ${H_i}(t), {\text{ }}i = 1, 2, 3$ represents the population harvesting or examine of mosquitoes. In this study, we propose two strategies for harvesting: proportional harvesting and constant harvesting, and the functions are given as follow:

(a) Proportional harvesting: ${H_1}(t) = e{S_M}, {H_2}(t) = e{E_M}, {H_3}(t) = e{I_M}$.

(b) Constant harvesting: ${H_1}(t) = h, {H_2}(t) = \left\{ {\begin{array}{*{20}{c}} h & {{E_M} > 0} \\ 0 & {{E_M} \leqslant 0} \end{array}} \right., {H_3}(t) = \left\{ {\begin{array}{*{20}{c}} h & {{I_M} > 0} \\ 0 & {{I_M} \leqslant 0} \end{array}} \right.$.

where $e, h$ are positive.

Lemma 2.1 Let the initial value $F(0) \geqslant {\text{0}}$, where

Then the solutions $F(t)$ of the model (2.2) are non-negative for all time $t > 0$. Furthermore, $\mathop {\lim }\limits_{t \to \infty } \sup$ ${N_H}(t) \leqslant \frac{{{\Lambda _H}}}{{{\mu _H}}}, \mathop {\lim }\limits_{t \to \infty } \sup {N_M}(t) \leqslant \frac{{{\Lambda _M}}}{{{\mu _M}}}.$

Proof Let the total dynamics of the human population is given by

then we have $0 \leqslant {N_H} \leqslant \frac{{{\Lambda _H}}}{{{\mu _H}}}$.

The total dynamics of mosquito population is given as

i.e. ${N'_M} \leqslant {\Lambda _M} - {\mu _M}{N_M}$. So same as $N_H$, when $t \to \infty$, we have $0 \leqslant {N_M} \leqslant \frac{{{\Lambda _M}}}{{{\mu _M}}}$.

It is obvious that

Hence

is positively invariant for the model (2.2) with non-negative initial conditions in $R_ + ^7$.

3.   Main results of the model with proportional harvesting

Considering the proportional harvesting, then model (2.1) can be rewritten as follows:

where $e$ is the fraction of the population removed for each time period. It is worthy to note that, if we view $e + {\mu _M}$ as ${\mu _M}$, model (3.1) becomes the basic model (2.1), which is the main research subject in reference ^[18]. Following the reference ^[18], we can get the main result of the model (3.1).

The disease-free equilibrium of the model (3.1) is ${E_0}(\frac{{{\Lambda _H}}}{{{\mu _H}}}, 0, 0, 0, \frac{{{\Lambda _M}}}{{{\mu _M} + e}}, 0, 0)$. Using the next generation operator method ^[19] on the model (3.1), the matrix ${F}$ and ${V}$ for model (3.1) are respectively given by

where ${k_1} = {\mu _H} + {\alpha _H}, {k_2} = {\mu _H} + r + \eta, {k_3} = {\mu _M} + {\delta _M}$. The basic reproduction number is given by $\rho ({F}{V}^{ - 1})$, which is ${R_0} = R{}_1 + \sqrt {{R_1}^2 + {R_2}}$, where

Before analyzing the stability, we might as well look at the biological significance of the basic reproduction number. $2{R_1} = \frac{{\rho {\beta _H}{\Lambda _H}{\alpha _H}}}{{{\mu _H}{k_1}k{}_2}}$ mainly describes the number of people who are directly infected by a human carrying Zika virus in unit space. $\frac{{{\beta _H}{\alpha _H}{\Lambda _M}}}{{{\mu _H}{k_1}k{}_2}}$ reflects the number of mosquitoes which are infected by a Zika virus-infected human in unit space. $\frac{{{\beta _M}{\delta _M}{\Lambda _H}}}{{{{(\mu _M^{} + e)}^2}({k_3} + e)}}$ gives the number of people infected by a infected mosquito biting in unit space. Therefore, during an infection period, the average number of infected human by a virus-infected person is $R_0$.

In this article, we mainly consider controlling the virus by reducing the number of mosquitoes. It doesn't affect the process of direct infection among humans. If ${\text{2}}{R_1} > 1,$ that means the virus breaks out primarily through direct human transmission, controlling mosquitoes has little effect on the spread of the virus, let alone eliminating it. Hence, we will focus on the situation ${\text{2}}{R_1} < 1$. And there are similar assumptions in section 4.

Lemma 3.1 For ${\text{2}}{R_1} < 1, 2{R_1} + {R_2} < 1$ holds if and only if ${R_0} < 1$. Furthermore, $2{R_1} + {R_2} = 1$ holds if and only if $R_0 = 1$.

Proof ${R_0} = R{}_1 + \sqrt {{R_1}^2 + {R_2}}$. Thus, if $R_0 < 1$, then $R{}_1 + \sqrt {{R_1}^2 + {R_2}} < 1$ or $\sqrt {{R_1}^2 + {R_2}} < 1 - {R_1}.$ Squaring both sides of the inequality gives $2{R_1} + {R_2} < 1$. On the other hand, $2{R_1} + {R_2} < 1$ means ${R_2} < 1 - 2{R_1}$. So ${R_2} + R_1^2 < 1 - 2{R_1} + R_1^2$, and then $\sqrt {{R_1}^2 + {R_2}} < 1 - {R_1}.$ And if $2{R_1} + {R_2} < 1$, then ${R_0} < 1$. Similarly, if ${R_0} = 1$, then $2{R_1} + {R_2} = 1$.

Based on the main results of the stability of disease-free equilibrium in ^[18], we can directly draw the following conclusions:

Theorem 3.1 ^[18] The disease-free equilibrium ${E_0}(\frac{{{\Lambda _H}}}{{{\mu _H}}}, 0, 0, 0, \frac{{{\Lambda _M}}}{{{\mu _M} + e}}, 0, 0)$ of model (3.1) is globally asymptotically stable if $R_0 < 1$, otherwise unstable.

The endemic equilibrium ${E^*} = (S_H^*, E_H^*, I_H^*, R_H^*, S_M^*, E_M^*, I_M^*)$ is the positive solution of the following system:

Then according to Eq (3.2), we can get:

And $I_H^*$ is the positive root of the following equation:

where

Case 1: If ${R_0} < 1$, which means $2{R_1} + {R_2} < 1$, then both $b > 0$ and $c > 0$ hold.

Case 2: If ${R_0} > 1$, which means $2{R_1} + {R_2} > 1$, then $c < 0$.

Case 3: If ${R_0} = 1$, then $b > 0, c = 0$.

The roots of Eq (3.3) are given as

Hence, we can get the following results:

Theorem 3.2 For model (3.1) :

(1). if ${R_0} > 1$, there exists a unique endemic equilibrium ${E^*} = (S_H^*, E_H^*, I_H^*, R_H^*, S_M^*, E_M^*, I_M^*)$, where $I_H^* = I_{H2}^*$;

(2). If ${R_0} \leqslant 1$, there is no endemic equilibrium.

Remark When $e = 0$, model (3.1) becomes model (2.1), which is model of ^[18]. From Theorem 3.2, we know that model (3.1) has only one endemic equilibrium if and only if $R_0 > 1$. That is quite different from the result of reference^[18], in which they claimed that there are two positive equilibrium points in the model.

Directly using the results of Theorem 4.2 and Theorem 5.2 in ^[18], the following results can be obtained:

Theorem 3.3 ^[18] If ${R_0} > 1$, then the unique endemic equilibrium of model (3.1) is not only locally asymptotically stable, but also globally asymptotically stable.

Therefore, the basic reproduction number is a crucial factor to determine whether the mosquito harvest can control the virus. From the expression ${R_0} = R{}_1 + \sqrt {{R_1}^2 + {R_2}}$, where

it is obvious that

Then there are two cases: one is $R_1\geq 1$ and the other is $R_1 < 1$. For the former, increasing $e$ can reduce the number of final infections but not eliminate the virus; but for the latter, it can be done.

4.   Main results of the model with constant harvesting

Considering the constant harvesting, then model (2.1) can be rewritten as follows:

where ${H_1}(t) = h, {H_2}(t) = \left\{ {\begin{array}{*{20}{c}} h & {{E_M} > 0} \\ 0 & {{E_M} \leqslant 0} \end{array}} \right., {H_3}(t) = \left\{ {\begin{array}{*{20}{c}} h & {{I_M} > 0} \\ 0 & {{I_M} \leqslant 0} \end{array}}. \right.$

4.1.   Basic reproduction number and stability of the disease-free equilibrium

It is clear that, if $0 \leqslant h \leqslant {\Lambda _M}$, model (4.1) has a disease-free equilibrium ${\overline E _0}{\text{ = }}(\frac{{{\Lambda _H}}}{{{\mu _H}}}, 0, 0, 0, \frac{{{\Lambda _M} - h}}{{{\mu _M}}}, 0,$ $0)$. Using the next generation operator method ^[19], the matrix ${F^*}$ and ${V^*}$ for model (4.1) are respectively given by

where ${k_1} = {\mu _H} + {\alpha _H}, {k_2} = {\mu _H} + r + \eta, {k_3} = {\mu _M} + {\delta _M}$. The basic reproduction number $R_0^*$ is equal to $\rho ({F^*}{V^*}^{ - 1})$, which can be written as $R_0^* = R_1^* + \sqrt {R{{_1^*}^2} + R_2^*}$, where

Next we mainly discuss in the case of $2R_1^* < 1$.

Similar to Lemma 3.1, we can get:

Lemma 4.1 $2R_1^* < 1$ and $2R_1^* + R_2^* < 1$ hold if and only if $R_0^* < 1$. Furthermore, $2R_1^* + R_2^* = 1$ if $R_0^* = 1$.

Theorem 4.1 If $0 \leqslant h \leqslant {\Lambda _M}$, the disease-free equilibrium ${\overline E _0}{\text{ = }}(\frac{{{\Lambda _H}}}{{{\mu _H}}}, 0, 0, 0, \frac{{{\Lambda _M} - h}}{{{\mu _M}}}, 0, 0)$ of model (4.1) is globally asymptotically stable if $R_0^* < 1$, otherwise unstable.

Proof The associated Jacobian matrix of the model (4.1) at ${\overline E _0}$ is given as

Clearly, $- {\mu _H}$ and $-{ \mu _M}$ are the eigenvalues of the Jacobian matrix above, and $- {\mu _H}$ is a double eigenvalue. The remaining four eigenvalues can be determined by the following equation:

where

Both $2R_1^* < 1$ and $2R_1^* + R_2^* < 1$ hold, since $R_0^* < 1$. Thus, the coefficients ${\overline G _i}$ for $i = 1, 2, 3, 4$ and all the order principal minor determinants are positive. Using the Routh Hurtwiz criteria ^[20], model (4.1) at the disease free equilibrium ${\overline E _0}$ is locally asymptotically stable if $R_0^* < 1$, otherwise unstable.

To show the global stability, we define the following Lyapunov function:

The time derivative of $\overline V$ is:

By model (4.1), we have

At $E_0$, we have

Therefore,

Choose the constants:

Notice that

we get

Thus, $\frac{{d\overline V (t)}}{{dt}} < 0$ if $R_0^* \leqslant 1$. And $\frac{{d\overline V (t)}}{{dt}}$ is zero if and only if

Therefore the largest compact invariant set in $\Omega$ is the singleton set $\{ {\overline E _0}\}$, and the model (4.1) is globally asymptotically stable in the interior of $\Omega$.

4.2.   Existence and stability of the endemic equilibrium of model (4.1)

4.2.1.   Existence of the endemic equilibrium

The endemic equilibrium ${\overline E ^*}{\text{ = }}(\overline S _H^*, \overline E _H^*, \overline I _H^*, \overline R _H^*, \overline S _M^*, \overline E _M^*, \overline I _M^*)$ is a positive solution of following system:

then according to the first four equations of (4.2), we can get

From the last three equations of (4.2), we can have

Rewrite Eq (4.4) as follows:

According to Eqs (4.3) and (4.5), let

Obviously, the positive roots of $f({I_H}) = 0$ are the key that could determine the existence of the positive equilibrium points of model (4.1). If the two curves of functions Eqs (4.6) and (4.7) (denoted by ${C_1}, {C_2}$) meet in the first quadrant, or the curve of function Eq (4.8) (denoted by $C$) and the positive half axis of the transverse axis meet, model (4.4) has positive equilibrium points.

From Eq (4.6) and Eq (4.7), the positive roots of $f({I_H}) = 0$ should fall into the interval $[{I_2}, {I_1}]$, and $h$ must meet $h < h_0$ to ensure the establishment of ${I_2} < {I_1}$, where

For curve $C$,

Then we have $f''({I_H}) > 0$ and $f'({I_H})$ monotone increasing on interval $[{I_2}, {I_1}]$, where the curve $C$ is concave. Let

When comparing the size of $y_1$ and $y_2$, the following corollaries can be obtained:

Corollary 1 If $2R_1^* + R_2^* \leqslant 1$ holds, we have ${y_1} \geqslant {y_2}$, then $f({I_H}) > 0$ for all ${I_H} > 0$.

Corollary 2 If $2R_1^* + R_2^* > 1$ and $2R_1^* < 1$ hold, which implies ${y_1} < {y_2}$, then there is a unique positive ${\hat I_H} \in ({I_2}, {I_1})$ so that $f'({\hat I_H}) = {f'_1}({\hat I_H}) - {f'_2}({\hat I_H}) = 0$. Furthermore, ${f'_{}}({I_H}) > 0$ if ${I_H} > {\hat I_H}$, and ${f'_{}}({I_H}) < 0$ if ${I_H} < {\hat I_H}$.

Curve $C_1$ has a vertical asymptote ${I_H} = {I_1}$. When ${I_H} = 0$ or ${I_H} = \frac{{{\mu _H}}}{{\rho {\beta _H}}}(2R_1^*-1)$, the curve ${C_1}$ intersects the axis ${I_H}$. Curve $C_2$ has a horizontal asymptote ${I_M} = \frac{{{\delta _M}{\Lambda _M}-(2{\delta _M} + {k_3})h}}{{{k_3}{\mu _M}}}$ and a vertical asymptote ${I_H} = -\frac{{{\mu _M}}}{{{\beta _M}}}$. The intersection of the curve and the axis $I_M$ is $(0, -\frac{{{k_3} + {\delta _M}}}{{{k_3}{\mu _M}}}h)$.

According to the analysis above, if $2R_1^* < 1$ holds, the possible location relationship of the curves within the first quadrant is shown in the following Figures:

According to the Figure 1(a) and 1(b) and corollary 2, when $f'({\hat I_H}) = f({\hat I_H}) = 0$, ${\hat I_H} \in ({I_2}, {I_1})$ holds, and a corresponding $h_1$ exists, which is

If ${h_1} < h < {h_0}$, then $f'({I_H}) > 0$ holds. And if $0 < h < {h_1}$, then $f'({I_H}) < 0$ holds.

Next, we discuss the number of intersection points. From corollary 1, it is concluded that there must be no chance for two curves $C_1$ and $C_2$ to meet in first quadrant, or curve $C$ and axis $I_H$ meet as $R_0^* \leqslant 1$. Since both $2R_1^* < 1$ and $R_0^* > 1$ hold, the results are as follows:

Case 1:The two curves (or curve $C$ and axis $I_H$) have two intersections as $0 < h < {h_1}$;

Case 2: The two curves (or curve $C$ and axis $I_H$) meet at only one point as $h = {h_1}$;

Case 3:If ${h_1} < h < {h_0}$ holds, there is no intersection of the two curves (or curve $C$ and axis $I_H$) due to $f({\hat I_H}) > 0$;

Case 4:If ${h_0} \leqslant h < \frac{{{\delta _M}{\Lambda _M}}}{{2{\delta _M} + {k_3}}}$, there is no intersection of curves $C_1$ and $C_2$ because ${I_2} \geqslant {I_1}$;

Case 5: If $\frac{{{\delta _M}{\Lambda _M}}}{{2{\delta _M} + {k_3}}} \leqslant h$, there is no intersection of curves $C_1$ and $C_2$, for the curve $C_2$ does not exist in the first quadrant according to the position of the horizontal asymptote.

Hence, we establish the following:

Theorem 4.2 Model (4.1) has no endemic equilibrium if $R_0^* \leqslant 1$. For $R_0^* > 1$ and $2R_1^* < 1$, there is no endemic equilibrium if ${h_1} < h$, only one endemic equilibrium if $h = {h_1}$ and two endemic equilibria if $0 < h < {h_1}$, where

4.2.2.   Stability of the endemic equilibrium

Here, we still mainly give the results under the condition $R_0^* > 1$ and $2R_1^* < 1$.

For $0 < h < {h_1}$, according to theorem 4.2, we can get two endemic equilibria, denoted by

where $\overline I _{H1}^* > \overline I _{H2}^*$.

Consider the Jacobian matrix of model (4.1) evaluated as

where ${\lambda _{HM}} = {\beta _H}(\rho {I_H} + {I_M}), {\lambda _{HH}} = {\beta _H}{S_H}, {\lambda _{MM}} = {\beta _M}{S_M}$.

The associate characteristic equation of $J$ is

Obviously, $- {\mu _H}$ and $- {\mu _M}$ are negative eigenvalues. Therefore, we care about the eigenvalues of the following equation:

where $E$ represents any positive equilibrium point of model (4.1), and

Set ${\lambda _1}, {\lambda _2}, {\lambda _3}, {\lambda _4}, {\lambda _5}$ are the five roots of the Eq (4.14), then according to the relationship between the roots of the algebraic equation and its coefficient, we have

From the analysis on $y = f({I_H})$ in section 4.2.1 and Figure 1, for endemic equilibrium $\overline E_1^*$, all the coefficient ${K_i}(\overline E_1^*), \ i = 1, 2, \cdots 5$ are positive. On the other hand, for another endemic equilibrium $\overline E_2^*$, ${K_5}(\overline E_2^*) < 0$ holds.

Theorem 4.3 For $R_0^* > 1$, $2R_1^* < 1$, and $0 < h < {h_1}$, the endemic equilibrium $\overline E _1^*$ of model (4.1) will be locally asymptotically stable, or unstable dimension and the number of central manifolds are both even. And the endemic equilibrium $\overline E _2^*$ of model (4.1) must be unstable.

In fact, through a lot of numerical simulations, we find that the positive equilibrium point $\overline E _1^*$ of model (4.1) should be locally asymptotically stable, rather than globally asymptotically stable.

5.   Numerical simulations and transmission dynamics comparisons

In this section, we continue to study transmission dynmamics of ZIKV model (3.1) and model (4.1) through numerical approach. We aim to investigate how harvesting affect disease spreading on the ZIKV dynamics.

5.1.   Proportional Harvesting

Firstly, we use a set of data from^[18], which are

When $\eta = 0.11, e = 0.001$, Reference ^[18] shows there is a backward bifurcation. In fact, the moment witnesses that ${R_0} = 1.011632168 > 1$ and a stable endemic equilibrium exists (see Figure 2). When $\eta = 0.2$, ${R_0} = 0.8271019214 < 1$ is hold, disease-free equilibrium for the model (3.1) is stable(See Figure 3).

5.2.   Constant harvesting

In order to better explain the local stability of the positive equilibrium, we change the value of parameter ${\Lambda _M}$ to 3, and the rest remain the same as in (5.1). We also give two sets of initial values: the initial value 1 is $(2, 4, 4, 3, 6, 4, 2)$ and the initial value 2 is $(0.01, 0.1, 0.5, 3, 0.1, 0.01, 0.1)$. It is shown in Figure 4 that the results are related not only to the value of $h$, but also to the initial value.

Figure 4 shows that the values of $I_H$ and $I_M$ departure under different initial conditions finally converge to 0 with $h = 0.004$; Figure 5 shows that when $h = 0.0035,$ the final $I_H$ values of departure from different initial conditions converge to different values: 0.2109 and 0 as well as the final $I_M$ converge to different values: 56.52 and 0.

5.3.   The effect of harvesting mosquitoes

To control the ZIKV infection, we can increase the value of $e$ (the fraction of the mosquitoes population harvesting during each time period) or $h$ (fixed harvesting constant). Theoretically speaking, $R_0$ (or $R_0^*$) decreases monotonically as $e$ (or $h$) increases, which means it is possible to control Zika virus by increasing the value of $e$ (or $h$). Here we focus on a more detailed result of this strategy on the number of infected people. we use another set of data shown in Table 2:

Assume $T$ is the peak time of disease transmission. From the second, third, sixth and seventh equations of system (3.1), we can get

Therefore we can build the relationship between $S_H(T)$ and $S_M(T)$ as follows:

Because $S_M(T)$ is decreasing with the mortality of mosquito increasing, supposing $e$ increases to $e_1$ ($e_1 > e$) and the corresponding peak time becomes $T_1$, then $S_H(T) < S_H(T_1)$. Noting $\frac{dS_H}{dt}\geq 0$, then it will take longer time to get $S_H(T_1)$. This means the peak time is postponed.

Similarly, we can get the following expression from system (4.1)

where $T_h$ is the peak time and $L(h) = \dfrac{\beta_M\delta_M}{\delta_M+k_3}S_M(T_h)-\dfrac{h}{I_H(T_h)}$. The first order approximate linear expression of $L(h)$ is

where $T_0$ is the peak time as $h = 0$. Because $S_M(T_h)$ is decreasing with the increase of $h$, $\frac{dS_M(T_h)}{dh} < 0$ holds. Then $L(h)$ decrease with $h$ increase. Take this result into Eq (5.3) and let $h$ grow to $h+\triangle h$ $(\triangle h > 0)$, then $S_H(T_h) < S_H(T_{h+\triangle h})$, which suggests that the peak time is delayed.

We increase the value of $e$ from 0.05 to 0.08, to 0.1 and then to 0.8, draw the time series diagram of the number of infected people respectively, and also give the time series diagram of the number of infected people under the four values of $h$ (0.0001, 0.00015, 0.0002, and 0.00025). Figure 6 shows that harvesting can postpone the peak of disease transmission with the mortality of mosquito increasing.

Now let's look at what happens to the number of people infected at peak time. If $\frac{dI_H}{dt}$ keeps being positive for two different values of $e$ or $h$, then we can get the longer peak time corresponding the greater the value of $e$ or $h$, and the more people will be infected at the peak time. Take Figure 6 as an example, the number of infected human reaches to about 0.23 units at peak time $T = 1847$ days when $e = 0.1$. However it is no more than 0.14 units at peak time $T = 1240$ when $e = 0.08$. If $\frac{dI_H}{dt}$ is negative for two different values of $e$ or $h$, then the less people will be infected at the peak time with $e$ or $h$ increase.

6.   Conclusions

In the present paper, we mainly study the control of ZIKV by taking advantage of continuous harvesting mosquitoes. Through the analysis of the models, the steady states of disease-free and endemic equilibrium are obtained for two types of harvesting models: proportional harvesting and constant harvesting. We find that it is possible to eliminate the virus by harvesting mosquitoes under certain conditions, no matter which of the two harvesting strategies is adopted.

Compare the two harvesting strategies, we find there is something in common:

(1) If ZIKV is primarily transmitted among humans, which means $2{R_1} > 1$ (or $2R_1^* > 1$), harvesting may only reduce the number of infections, but not eliminate the disease.

(2) Increasing harvesting (enhance the value of $e$ or $h$) may stimulate an increase in the number of virus infections at some point.

(3) Harvesting can postpone the peak of disease transmission with the mortality of mosquito increasing.

We also should pay attention to the difference between the proportional harvesting and constant harvesting. The proportional harvesting is easier to control the virus through the basic reproduction number $R_0$. If $2{R_1} < 1$, the purpose of permanently eliminating the virus will be achieved by adjusting the value of $e$ so that the value of $R_0$ is less than 1. And get the minimum value of $e$ that can kill the virus. For the constant harvesting, it is hard to confirm the minimum value of $h$ because different initial values lead to different results. In other words, for the proportional harvesting, if and only if ${R_0} < 1$, the Zika virus can be wiped out. But for the constant harvesting, the Zika virus may be controlled even though ${R_0^*} > 1$ holds.

Figure 7 shows the relationship between $R_0$ and mosquitoes recruitment rate ${\Lambda _M}$ under different capture coefficients $e$.You will find that $R_0$ gradually decreases with the increase of $e$. If the recruitment rate ${\Lambda _M}$ is 0.3, then $R_0$ will be less than one as $e = 0.6$. When ${\Lambda _M}$ grows to 0.5, taking $e$ equal to 0.8 will force $R_0$ to be less than one and ZIKV will be wiped out.

But for the constant harvesting strategy, it is quite different. As can be seen from Figure 8, under the same parameter conditions, different initial value conditions lead to different final infection numbers. Virus may eventually disappear, although $R_0^*$ is greater than 1.

In Figure 9, the red solid curve shows $R_0^*$ changes with $h$ while the black dotted line respects $R_0^* = 1$.The vertical green dotted line is $h = {h_1} \approx 0.002196696388$, which is a split line. According to theoretical reasoning that, the system has no positive equilibrium point when the selected $h$ is larger than $h_1$ (See Figure 10), which means that the disease has been eliminated. Here we get $h_1$ from the following:

where ${\hat I_H}$ is the root of $f'({\hat I_H}) = f({\hat I_H}) = 0$, ${\hat I_H} \in ({I_2}, {I_1})$. Of course, we can choose $h = {h_0}$ (See Figure 7, the vertical blue dotted line is $h = {h_0} \approx 0.002475094697$), which is larger than $h_1$ and can be figured out from

This looks easier to calculate. To sum up, both $h_1$ and $h_0$ enable us to choose an appropriate value of $h$ to control the virus. Clearly this approach has nothing to do with $R_0^*$.

Therefore, if we control the spread of the disease by means of constant harvesting of mosquitoes, we must consider the initial values and some other known parameters to select the appropriate harvesting constants, otherwise it may not control the disease, but will cause more people to be infected.

Conflict of interest

The authors declare that they have no competing interests.