Bifurcation analysis of a wild and sterile mosquito model

1.   Introduction

Mosquitoes not only buzz around us, but also transmit various diseases by biting susceptible persons. Diseases, such as Dengue fever, malaria, yellow fever, Zika viruses, West Nile virus and Japanese encephalitis are carried and transmitted by mosquitoes ^{[1,2,3,4,5,6]}. Mosquito-borne diseases have been a considerable public health concern worldwide ^[7]. According to the estimate of WHO report, 2.5 billion of people are living in dengue risk area and 50 to 100 million cases of dengue occur every year ^[8].

One of the ways to prevent the transmissions of mosquito-borne diseases is to control the number of the wild mosquito population. An important effort to reduce or eradicate wild mosquitoes is insect sterilization, also called sterile insect technique (SIT). The male mosquitoes are modified to be sterile ones by the application of radical or other chemical/physical methods and then released into the environment to mate with wild mosquitoes. A wild female mosquito that mates with a sterile male will either not reproduce, or produce eggs which will not hatch. Recently, a mosquito eradication project named "Debug Fresno" was carried out by the life sciences arm of Google's parent company ^[9]. The scientists released up to 20 million sterile male mosquitoes infected with Wolbachia and let them mate with wild female mosquitoes to reduce the impact that disease-carrying mosquitoes have on human health. The theoretical basis of project "Debug Fresno" is sterile insect technique since Wolbachia bacteria can cause a form of conditional sterility.

Before SIT is carried out, it is critical to choose the optimal release strategies and estimate the release number of sterile mosquitoes. Mathematical models have been proven to be useful in overcoming such challenging problems in population dynamics and epidemiology. A number of mathematical models have been formulated to investigate the interactive dynamics and the control of the wild and sterile mosquito populations ^{[10,11,12,7,13,14,15,16,17,18,19,20]}. Recently, other models are also proposed to investigate the spreading dynamics by incorporating stage-structure ^[21,22] or delay in differential equations ^[23,24,25], spatial diffusion in reaction-diffusion equations ^[26,27], and the environmental heterogeneity in stochastic equations ^[28]. In ^[24], the number of released mosquitoes is assumed to be a variable satisfying an independent dynamical equation. This new idea seems to be the first in such a direction in the mathematical modeling of mosquito-borne diseases. It is found that the model displays bistability with two stable steady-states and one unstable steady-state when the abundance of released males is smaller than an exact value of the threshold releasing intensity. Simulation also suggests the existence of one or more stable periodic solutions.

In ^[7], the authors proposed the following continuous-time ordinary differential equation model to study the interactive dynamics of the wild and sterile mosquitoes:

Here $w(t)$ and $g(t)$ are the numbers of wild and sterile mosquitoes at time $t$, respectively; $C(N)$ is the number of mating per individual per unit of time with $N = w+g$; $a_0$ is the number of wild offspring produced per mate; $R(\cdot)$ is the release rate of the sterile mosquitoes; $\mu_i$ and $\xi_i$, $i = 1, 2$, are the density independent and dependent death rates of wild and sterile mosquitoes, respectively. In ^[7], the authors focused on the impact of the different release strategies of the sterile insect technique on disease transmission and consider the following three cases:

(ⅰ) $C(N) = c$, $R(w) = b$;

(ⅱ) $C(N) = \frac{c_0 N}{1+N}$, $R(w) = b w$;

(ⅲ) $C(N) = \frac{c_0 N}{1+N}$, $R(w) = \frac{b w}{1+w}$.

Here, $c, c_0, b$ are all constants. In case (ⅰ), both the mating rate and the release rate are constants. In case (ⅱ), the release rate is a linear function, and in case (ⅲ), the release rate form takes into account an Allee effect as in ^[29,7,30]. The mating rate is a Holling-Ⅱ type function both in case (ⅱ) and (ⅲ). The authors explored the existence of all possible equilibria and the stability of these equilibria. They found that the positive equilibrium is a stable node and the system (1.1) has no closed orbits in case (ⅰ). An explicit threshold release value for $b$ was also established for this case. For case (ⅱ), they gave some results about the stability of the positive equilibrium when $\mu_1\leq \mu_2$ and $\xi_1\leq 2\xi_2$. A threshold release value for $b$ was also obtained implicitly because of the complexity of the system. They also gave some results about the existence and stability of the positive equilibrium when $\mu_1 = \mu_2$ and $\xi_1 = \xi_2$ for case (ⅲ).

In the present paper, we will give a detailed follow-up work focused on the case (ⅱ). Let $a = a_0c_0$, the dynamics of the interacting mosquitoes in case (ⅱ) is governed by the following system:

We aim to investigate the impact of the mating rate and the release strategy on the dynamics of (1.2) without the restriction $\xi_1\leq 2\xi_2$. Using $a$ and $b$ as the bifurcation parameters, we find that the system (1.2) possesses rich dynamics in different parameter ranges.

This paper is organized as follows. In Section 2, we give some result about the existence, and local stability of equilibria especially when $\xi_1 > 2\xi_2$. In Section 3, we show the existence of Hopf bifurcations and homoclinic bifurcations. The Bogdanov-Takens bifurcation analysis is given in Section 4. We end with some discussions in Section 5.

2.   Equilibria and stability

First, for the convenience of readers, we show some known results about the existence of a positive equilibrium for model (1.2) which have been obtained in ^[7] and give some more comprehensive results about its stability.

Define

Then the system (1.2) admits a positive equilibrium if and only if $F(N) = 0$ admits a positive root.

Define

Then $F(N) = 0$ admits a positive root if and only if there exists a positive solution $N$ to the equation $b = B(N)$. It is useful to note that for any $N > 0$, there is at most one $b > 0$ which is $b = B(N)$ such that $(w, g)$ is a positive equilibrium of (1.2), where

For the function $b = B(N)$ defined in (2.2), if $a > \mu_1+\xi_1+2\sqrt{\mu_1\xi_1}$, there is a unique $\bar{N} = \bar{N}(a)$ such that $B'(\bar{N}) = 0$. We define

Then $\bar{b} > 0$ and we have the following results.

Proposition 2.1. Suppose all parameters $a, b, \mu_1, \mu_2, \xi_1, \xi_2$ are positive and $a > \mu_1+\xi_1+2\sqrt{\mu_1\xi_1}$. Let $F(N)$ and $\bar{b}$ be defined as in (2.1) and (2.3), respectively. Then the equation $F(N) = 0$ has no positive roots if $b > \bar{b}$, one positive root $\bar{N}$ if $b = \bar{b}$ and two positive roots $N^{\pm}$ with $N^- < N^+$ if $0 < b < \bar{b}$. Furthermore, $N^{\pm}$ depend on $b$ continually, and $N^{-'}(b) > 0$, $N^{+'}(b) < 0$.

Proof. The first part of the proposition is obvious, and we only need to prove the last part. It is easy to know that $N^- < \bar{N} < N^+$ whenever $N^{\pm}$ and $\bar{N}$ exist. Define a two-variable function

By the implicit function theorem on the intervals $[N_1, \bar{N}]$ and $[\bar{N}, N_2]$, two functions $N^- = N^-(b)$ and $N^+ = N^+(b)$ satisfying $N_1 < N^- < \bar{N}$ and $\bar{N} < N^+ < N_2$ can be obtained from the equation $\tilde{F}(N, b) = 0$, where

Then we have $N'(b) = -\frac{\partial \tilde{F}}{\partial b}/\frac{\partial \tilde{F}}{\partial N}$. Simple calculations show that $\frac{\partial \tilde{F}}{\partial b} = (1+N)(\mu_1+\xi_1N) > 0$ for any positive values of $N$ and $\frac{\partial \tilde{F}}{\partial N} = F'(N)$. Note that $F'(N^-) < 0$ and $F'(N^+) > 0$. Then $N^{-'}(b) > 0$ and $N^{+'}(b) < 0$.

Then the existence and multiplicity of positive equilibria for system (1.2) follow immediately from Proposition 2.1.

Lemma 2.2. Suppose all parameters $a, b, \mu_1, \mu_2, \xi_1, \xi_2$ are positive and $a > \mu_1+\xi_1+2\sqrt{\mu_1\xi_1}$. Define $\bar{b}$ as in (2.3). Then

(1) If $b > \bar{b}$, then the system (1.2) has no positive equilibria.

(2) If $b = \bar{b}$, then the system (1.2) has a unique positive equilibrium $(\bar{w}, \bar{g})$, where $\bar{w} = (1+\bar{N})(\mu_1+\xi_1\bar{N})/a, \bar{g} = \bar{b}(1+\bar{N})(\mu_1+\xi_1\bar{N})/(a(\mu_2+\xi_2 \bar{N}))$.

(3) If $0 < b < \bar{b}$, then the system (1.2) has two positive equilibria $(w^{\pm}, g^{\pm})$, where $w^{\pm} = (1+N^{\pm})(\mu_1+\xi_1N^{\pm})/a, g^{\pm} = b(1+N^{\pm})(\mu_1+\xi_1N^{\pm})/(a(\mu_2+\xi_2 N^{\pm}))$.

Note that if $b$ is regarded as a bifurcation parameter, then $b = \bar{b}$ is a saddle-node bifurcation point marked by $b_{LP}$. From Lemma 2.2, the system (1.2) may have three nonnegative equilibria: $(0, 0), (w^-, g^-)$ and $(w^+, g^+)$. In the following, $(w^+, g^+)$ is referred as the higher wild mosquitoes state, while $(w^-, g^-)$ is the lower wild mosquitoes state.

Next we investigate the local stability of these equilibria. From ^[7], $(0, 0)$ is locally asymptotically stable and $(w^-, g^-)$ is an unstable saddle whenever it exists. Then we only need to consider the stability of $(w^+, g^+)$. The associated Jacobian matrix at the equilibrium $(w^+, g^+)$ is given by

Then we have

and the trace of $J(w^+, g^+)$ is given by

where $T(N)$ is defined as

Then $(w^+, g^+)$ is locally asymptotically stable if $T(N^+) < 0$ and unstable if $T(N^+) > 0$ ^[7].

In the present paper, we always assume that the sterile mosquitoes are prone to death for their adaptability and competitive abilities are diminished as compared to the wild populations, i.e. $\mu_1\leq \mu_2$. In ^[7], the authors have shown that the equilibrium $(w^+, g^+)$ is locally asymptotically stable if $\xi_1\leq 2\xi_2$. Now we will discuss the stability of $(w^+, g^+)$ for the case $\xi_1 > 2\xi_2$.

Define

Note that $\xi_2-2\xi_1 < 0$. Then if $a < a_0$ or $a > a_0$ and $\triangle_1 > 0$, we have $T(N) < 0$ which implies that the equilibrium $(w^+, g^+)$ is locally asymptotically stable. For the case that $a > a_0$ and $\triangle_1 > 0$, the equation $T(N) = 0$ has two positive roots, denoted by $N_*$ and $N^*$ with $N_* < N^*$, where

To discuss the stability of $(w^+, g^+)$ for this case, we first give the following proposition.

Proposition 2.3. Suppose all parameters $a, b, \mu_1, \mu_2, \xi_1, \xi_2$ are positive and $a > \mu_1+\xi_1+2\sqrt{\mu_1\xi_1}$, $\mu_1\leq \mu_2$. Then the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $b\approx 0$.

Proof. It is useful to note that $a N^+ = (1+N^+)(\mu_1+\xi_1N^+)$ if $b = 0$. Then if $b = 0$, we have

which implies that the equilibrium $(w^+, g^+)$ is locally asymptotically stable if $b = 0$. By the continuity of $T(N^+)$ with respect to $b$, we have $T(N^+) < 0$ for $b\approx 0$ which implies the local stability of $(w^+, g^+)$ when $b\approx 0$.

When $0 \ll b < \bar{b}$, the discussion is divided into the following cases:

(ⅰ) $N_2 < N_*$ or $N^* < \bar{N} < N_2$. For these cases, we have $T(N^+) < 0$ which implies that the equilibrium $(w^+, g^+)$ is locally asymptotically stable.

(ⅱ) $N_* < \bar{N} < N^* < N_2$. For this case, there exists $b = b_1$ such that $N^+(b_1) = N^*$. Then $T(N^+) < 0$ when $0 < b < b_1$ and $T(N^+) > 0$ when $b_1 < b < \bar{b}$ which implies that the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < b_1$ and unstable when $b_1 < b < \bar{b}$.

(ⅲ) $\bar{N} < N_* < N^* < N_2$. For this case, there exist $b = b_1$ and $b = b_2$ with $b_1 < b_2$ such that $N^+(b_1) = N^*$ and $N^+(b_2) = N_*$. Then $T(N^+) < 0$ when $0 < b < b_1$ or $b_2 < b < \bar{b}$, and $T(N^+) > 0$ when $b_1 < b < b_2$, which implies that the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < b_1$ or $b_2 < b < \bar{b}$, and unstable when $b_1 < b < b_2$.

Then we have the following results about the stability of the equilibrium $(w^+, g^+)$.

Theorem 2.4. Assume all parameters $a, b, \mu_1, \mu_2, \xi_1, \xi_2$ are positive and $a > \mu_1+\xi_1+2\sqrt{\mu_1\xi_1}$, $\mu_1\leq \mu_2$. Then the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $\xi_1\leq2\xi_2$. When $\xi_1 > 2\xi_2$, define $\bar{b}$ as in (2.3), $a_0, \triangle_1$ as in (2.5), and $N_*, N^*$ as in (2.6). Then the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $a < a_0$ or $a > a_0$ and $\triangle_1 > 0$. When $a > a_0$ and $\triangle_1 > 0$, the following cases may occur:

(1) If $N_2 < N_*$ or $N^* < \bar{N} < N_2$, then the equilibrium $(w^+, g^+)$ is locally asymptotically stable.

(2) If $N_* < \bar{N} < N^* < N_2$, then there exists $b = b_1$ such that $N^+(b_1) = N^*$ and the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < b_1$ and unstable when $b_1 < b < \bar{b}$.

(3) If $\bar{N} < N_* < N^* < N_2$, then there exist $b = b_1$ and $b = b_2$ with $b_1 < b_2$ such that $N^+(b_1) = N^*$ and $N^+(b_2) = N_*$. The equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < b_1$ or $b_2 < b < \bar{b}$, and unstable when $b_1 < b < b_2$.

Now we give some examples to illustrate our results and show how the existence and stability of the positive equilibrium $(w^+, g^+)$ vary with $b$.

Example 2.5. Choosing the values of the parameters as $\mu_1 = 1, \mu_2 = 1, \xi_1 = 4, \xi_2 = 0.51, a = 20$, we can numerically obtain $N_2 = 3.6821, N_* = 0.1517, N^* = 0.9882, \bar{N} = 0.6863$ which satisfy $N_* < \bar{N} < N^* < N_2$. Furthermore, we have $\bar{b} = 1.5841$. According to Theorem 2.4, there exists $b = b_1$ such that $N^+(b_1) = N^*$. Actually, we obtain $b_1 = 1.5146$. Then the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < 1.5146$ and unstable when $1.5146 < b < 1.5841$.

Example 2.6. Choosing the values of the parameters as $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 48$, we can numerically obtain $N_2 = 3.2448, N_* = 0.2498, N^* = 1.0207, \bar{N} = 0.3947$ which satisfy $N_* < \bar{N} < N^* < N_2$. Furthermore, we have $\bar{b} = 3.6646$. For this case, $b_1 = 2.9187$ and the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < 2.9187$ and unstable when $2.9187 < b < 3.6646$.

Example 2.7. Choosing the values of the parameters as $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 40$, we can numerically obtain $N_2 = 2.5092, N_* = 0.3881, N^* = 0.5680, \bar{N} = 0.3783$ which satisfy $\bar{N} < N_* < N^* < N_2$. Furthermore, we have $\bar{b} = 2.6593$. According to Theorem 2.4, there exist $b = b_1$ and $b = b_2$ with $b_1 < b_2$ such that $N^+(b_1) = N^*$ and $N^+(b_2) = N_*$. Actually, we obtain $b_1 = 2.5374$ and $b_2 = 2.6588$. Then the equilibrium $(w^+, g^+)$ is locally asymptotically stable when $0 < b < 2.5374$ or $2.6588 < b < 2.6593$ and unstable when $2.5374 < b < 2.6588$.

3.   Hopf bifurcations and numerical simulations

In this section, we mainly give a simulation illustration for Hopf bifurcations of system (1.2). From the discussion in Section 2 there exist values of $b$ such that the corresponding characteristic matrix of $(w^+, g^+)$ has a pair of complex roots, denoted by

where

and $F(N)$ and $T(N)$ are defined as in (2.1) and (2.4), respectively. When $b = b_k, k = 1, 2$, we have $T(N^+) = 0$. Then there are a pair of purely imaginary eigenvalues for $J(w^+, g^+)$. This means that a Hopf bifurcation may occur. Choosing $b$ as a bifurcation parameter, the existence of Hopf bifurcations could be obtained by algebraic and logical methods ^[31]. In fact, we have

and

Here, we use the facts that $T'(N^+)|_{b = b_1} = T'(N^*) < 0$, $T'(N^+)|_{b = b_2} = T'(N_*) > 0$ and $N^{+'}(b_k) < 0$ which result in $\alpha'(b_1) > 0$ and $\alpha'(b_2) < 0$. Then we obtain the following result.

Theorem 3.1. Assume all parameters $a, b, \mu_1, \mu_2, \xi_1, \xi_2$ are positive and $b = b_k, k = 1, 2$ are defined as in Theorems 2.4. Then system (1.2) undergoes a Hopf bifurcation at $b_1$ and $b_2$ if $b_1$ or $b_2$ exists.

In order to determine the type of the Hopf bifurcation, the direction of the Hopf bifurcations at $b = b_1$ and $b = b_2$ is calculated (see Appendix A for details). We focus on the exhibition of the existence and the direction of the Hopf bifurcations by numerical simulation. The Hopf bifurcations are classified into the following four different cases to show that the system (1.2) may admit one or two Hopf bifurcation points and both supercritical and subcritical Hopf bifurcations may occur for different parameter ranges. In all of the following bifurcation diagrams, the points labelled "H" denote the points at which the trace of the corresponding characteristic equation equal to zero, and the points labelled "LP" denote the limit points. In all of our phase portraits, the red points denote the equilibria and the green curves are the stable and unstable manifolds of the positive equilibrium $(w^-, g^-)$ which is a saddle.

Case 1. Existence of a subcritical Hopf bifurcation and one unstable limit cycle

We choose the values of the parameters as $\mu_1 = 1, \mu_2 = 1, \xi_1 = 4, \xi_2 = 0.51, a = 20$ (see Example 2.5 in Section 2). In Figure 1 (a)-(c), the bifurcation diagram, the maximum and the minimum of the cycles and the periods of limit cycles for different values of $b$ are shown. Here the bifurcation diagram (produced with MatCont) shows there exists a saddle-node bifurcation point at $b = b_{LP}$, which is labelled "LP" and one Hopf bifurcation points at $b = b_1$ (labelled a black "H") in Figure 1 (a). The point (labelled a red "H") represents a neutral saddle which is not a Hopf bifurcation point. Numerically calculated first Lyapunov coefficients at the Hopf bifurcation point is positive which indicates that the Hopf bifurcation point is subcritical. The branch of limit cycles from $b_1$ is not a closed loop, but is an "open-ended" one. Then from the global Hopf bifurcation theorem ^[32,33], the period of the periodic orbits on the bifurcating branch must be unbounded. In Figure 1 (c), one can observe that the period of cycles is decreasing in $b$ for the bifurcation diagram corresponding to (a). Furthermore, the period approach to $\infty$ when $b$ tends to some $b = b_{HL}$, which is a homoclinic bifurcation point where the limit cycles converge to a homoclinic orbit based on the saddle equilibrium $(w^-, g^-)$. This can be observed in Figure 1 (b) as the limit of minimum of the cycles is the lower positive equilibrium. We have called this bifurcation diagram a "lotus" in ^[34] with the Hopf bifurcation point as the base and the homoclinic bifurcation point the top of the lotus.

Figure 2 shows the evolution of phase portraits of system (1.2) with $\mu_1 = 1, \mu_2 = 1, \xi_1 = 4, \xi_2 = 0.51, a = 20$ when $b$ increases. One can numerically calculate that the bifurcation values $b_1 = 1.51462$ (Hopf bifurcation point), and $b_{LP} = 1.584068$ (saddle-node bifurcation point). Here the Hopf bifurcation point is subcritical so the bifurcating periodic orbits are unstable. Moreover the homoclinic bifurcation point is $b_{HL} = 1.4964786$. Apparently when $b > b_{LP}$, the equilibrium $(0, 0)$ is the attractor and attracts all positive orbits. When $0 < b < b_{HL}$, the positive equilibria $(w^+, g^+)$ and $(0, 0)$ are both stable, and the basins of attraction for the two locally stable equilibria are separated by the stable manifold of the saddle equilibrium $(w^-, g^-)$ (see Figure 2 (a)); when $b = b_{HL}$ a homoclinic orbit exists (see Figure 2 (b)); when $b$ increases to $b_{HL} < b < b_1$, the equilibrium $(w^+, g^+)$ is still stable and an unstable limit cycle emerges from the subcritical Hopf bifurcation (see Figure 2 (c)); when $b$ increases to $b_1 < b < b_{LP}$, the equilibrium $(w^+, g^+)$ becomes unstable and almost all the positive orbits except for positive equilibria $(w^{\pm}, g^{\pm})$ and stable manifold (the green orbits) of $(w^-, g^-)$ tend to the trivial equilibrium $(0, 0)$ in Figure 2 (d).

Case 2. Existence of a supercritical Hopf bifurcation and one stable limit cycle

In Figure 3 (a)-(c), the bifurcation diagram, the maximum and the minimum of the cycles and periods of limit cycles for different values of $b$ are shown when we choose $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 48$ (see Example 2.6 in Section 2). Here the bifurcation diagram (produced with MatCont) shows there exists a saddle-node bifurcation point at $b = b_{LP}$, which is labeled "LP" and one Hopf bifurcation points at $b = b_1$ (labelled a black "H") in Figure 3 (a). The point (labelled a red "H") represents a neutral saddle which is not a bifurcation point. Unlike the case in Figure 1, numerically calculated first Lyapunov coefficients at the bifurcation point is negative which indicates that this Hopf bifurcation point is supercritical. The branch of limit cycles from $b_1$ is also called a "lotus" in ^[34]. In Figure 3 (c), one can also observe that the period of cycles is increasing in $b$ for the bifurcation diagram corresponding to (a). Furthermore, the period approach to $\infty$ when $b$ tends to some $b = b_{HL}$, which is a homoclinic bifurcation point where the limit cycles converge to a homoclinic orbit based on the saddle equilibrium $(w^-, g^-)$. This can be observed in Figure 3 (b) as the limit of minimum of the cycles is the lower positive equilibrium.

Figure 4 shows the evolution of phase portraits of system (1.2) with $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 48$ when $b$ increases. One can numerically calculate that the bifurcation values $b_1 = 2.918737$ (Hopf bifurcation point), and $b_{LP} = 3.664595$ (saddle-node bifurcation point). Here the Hopf bifurcation point is supercritical so the bifurcating periodic orbits are stable. Moreover the homoclinic bifurcation point is $b_{HL} = 2.9568$. Apparently when $b > b_{LP}$, the equilibrium $(0, 0)$ is the attractor and attracts all positive orbits. When $0 < b < b_1$, the positive equilibria $(w^+, g^+)$ and $(0, 0)$ are both stable, and the basins of attraction for the two locally stable equilibria are separated by the stable manifold of the saddle equilibrium $(w^-, g^-)$ (see Figure 4 (a)); when $b$ increases to $b_1 < b < b_{HL}$, the equilibrium $(w^+, g^+)$ becomes unstable and a stable limit cycle emerges from the supercritical Hopf bifurcation (see Figure 4 (b)). At $b = b_{HL}$, a homoclinic orbit exists (see Figure 4 (c)). When $b$ increases to $b_{HL} < b < b_{LP}$, the homoclinic orbit is broken and the equilibrium $(w^+, g^+)$ is still unstable. All of the positive orbits except for positive equilibria $(w^{\pm}, g^{\pm})$ and stable manifold (the green orbits) of $(w^-, g^-)$ tend to the trivial equilibrium $(0, 0)$ in Figure 4 (d).

Case 3. Existence of a subcritical Hopf bifurcation and two limit cycles

An interesting phenomenon occurs when we choose the values of the parameters as $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 52$. In Figure 5 (a)-(c), the bifurcation diagram, the maximum and the minimum of the cycles and periods of limit cycles for different values of $b$ are shown. Here the bifurcation diagram (produced with MatCont) shows there exist a saddle-node bifurcation point at $b = b_{LP}$, which is labelled "LP" and one Hopf bifurcation points at $b = b_1$ (labelled a black "H") in Figure 5 (a). The point (labelled a red "H") is a neutral saddle which is not a bifurcation point. Numerically calculated first Lyapunov coefficients at the Hopf bifurcation point is positive which indicates this Hopf bifurcation point is subcritical. In Figure 5 (c), one can observe that the period approach to $\infty$ when $b$ tends to some $b = b_{HL}$, which is a homoclinic bifurcation point where the limit cycles converge to a homoclinic orbit based on the saddle equilibrium $(w^-, g^-)$. This can be observed in Figure 5 (b) as the limit of minimum of the cycles is the lower positive equilibrium. Note that, in Figure 5 (c), one can also observe that the period of cycles is not increasing in $b$ and there exists some range of $b$, $b_{LPC} < b < b_{HL}$, that the period of cycles is not a single-valued function of $b$. This implies that the system (1.2) admits multiple periodic orbits for the same $b$ value for the bifurcation diagram corresponding to (a). Unlike the case in 1, the branch of limit cycles from $b_1$ is an "open-ended" one which is called a "pepper" for there is another saddle-node bifurcation of cycles on the branch ^[34].

Figure 6 shows the evolution of phase portraits of system (1.2) with $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 52$ when $b$ increases. One can numerically calculate that the bifurcation values $b_1 = 3.130339$ (Hopf bifurcation point), and $b_{LP} = 4.168163$ (saddle-node bifurcation point). Here the Hopf bifurcation point is subcritical so the bifurcating periodic orbits are unstable. Moreover the homoclinic bifurcation point is $b_{HL} = 3.1562665$ and a limit cycle saddle-node bifurcation point is $b_{LPC} = 3.129271$.

Apparently when $b > b_{LP}$, the equilibrium $(0, 0)$ is the attractor and attracts all positive orbits. When $0 < b < b_{LPC}$, the positive equilibria $(w^+, g^+)$ and $(0, 0)$ are both stable, and the basins of attraction for the two locally stable equilibria are separated by the stable manifold of the saddle equilibrium $(w^-, g^-)$ (see Figure 6 (a)); when $b$ increases to $b_{LPC} < b < b_1$, the equilibrium $(w^+, g^+)$ is also stable and an unstable limit cycle emerges from the subcritical Hopf bifurcation which implies the existence of two limit cycles (see Figure 6 (b)). when $b$ increases to $b_1 < b < b_{LP}$, the equilibrium $(w^+, g^+)$ becomes unstable and there exists one stable limit cycle (see Figure 6 (c)). At $b = b_{HL}$, a homoclinic orbit exists (see Figure 6 (d)). When $b$ increases to $b_{HL} < b < b_{LP}$, the homoclinic orbit is broken and the equilibrium $(w^+, g^+)$ is still unstable. All of the positive orbits except for positive equilibria $(w^{\pm}, g^{\pm})$ and stable manifold (the green orbits) of $(w^-, g^-)$ tend to the trivial equilibrium $(0, 0)$ in Figure 6 (e).

Case 4. Existence of two supercritical Hopf bifurcations and one stable limit cycle

We choose the values of the parameters as $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 40$ (see Example 2.7 in Section 2). In Figure 7 (a)-(c), the bifurcation diagram, the maximum and the minimum of the cycles and periods of limit cycles for different values of $b$ are shown. Here the bifurcation diagram (produced with MatCont) shows there exists a saddle-node bifurcation point at $b = b_{LP}$, which is labeled "LP" and two Hopf bifurcation points at $b = b_1$ and $b = b_2$ with $b_1 < b_2$ (labelled a black "H") in Figure 7 (a). Numerically calculated first Lyapunov coefficients at the bifurcation points are both negative which indicates the two Hopf bifurcation points are supercritical. The branch of limit cycles from $b_1$ is also a "lotus". In Figure 7 (c), one can also observe that the period of cycles is increasing in $b$ for the bifurcation diagram corresponding to (a). Furthermore, the period approach to $\infty$ when $b$ tends to some $b = b_{HL}$, which is a homoclinic bifurcation point where the limit cycles converge to a homoclinic orbit based on the saddle equilibrium $(w^-, g^-)$. This can be observed in Figure 7 (b) as the limit of minimum of the cycles is the lower positive equilibrium. The branch of limit cycles from $b_1$ is also a "lotus" which is unlike the cases $(b), (d)$ and $(e)$ which we called a "bubble" or a "heart" in ^[34] when there exist two supercritical Hopf bifurcation points.

Figure 8 shows the evolution of phase portraits of system (1.2) with $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95, a = 40$ when $b$ increases. One can numerically calculate that the bifurcation values $b_1 = 2.537385$ (Hopf bifurcation point), $b_2 = 2.658842$ (Hopf bifurcation point) and $b_{LP} = 2.659327$ (saddle-node bifurcation point). Here the Hopf bifurcation points are supercritical so the bifurcating periodic orbits are stable. Moreover the homoclinic bifurcation point is $b_{HL} = 2.59096$. Apparently when $b > b_{LP}$, the equilibrium $(0, 0)$ is the attractor and attracts all positive orbits. When $0 < b < b_1$, the positive equilibria $(w^+, g^+)$ and $(0, 0)$ are both stable, and the basins of attraction for the two locally stable equilibria are separated by the stable manifold of the saddle equilibrium $(w^-, g^-)$ (see Figure 8 (a)); when $b$ increases to $b_1 < b < b_{HL}$, the equilibrium $(w^+, g^+)$ becomes unstable and a stable limit cycle emerges from the supercritical Hopf bifurcation (see Figure 8 (b)). At $b = b_{HL}$, a homoclinic orbit exists (see Figure 8 (c)). When $b$ increases to $b_{HL} < b < b_2$, the homoclinic orbit is broken and the equilibrium $(w^+, g^+)$ is still unstable and all of the positive orbits except for positive equilibria $(w^{\pm}, g^{\pm})$ and stable manifold (the green orbits) of $(w^-, g^-)$ tend to the trivial equilibrium $(0, 0)$ in Figure 8 (d). When $b$ increases to $b_2 < b < b_{LP}$, the positive equilibrium $(w^+, g^+)$ becomes stable and and the basins of attraction for the two locally stable equilibria $(w^+, g^+)$ and $(0, 0)$ are separated by the stable manifold of the saddle equilibrium $(w^-, g^-)$ (see Figure 8 (e)).

From the simulations above, we show that the system (1.2) may admit one or two Hopf bifurcation points and the Hopf bifurcation points may be supercritical and subcritical. At a critical value $b = b_{HL}$, a homoclinic orbit exists. Thus, we conclude that the system may undergo a sequence of bifurcations including saddle-node bifurcation, supercritical Hopf bifurcation, subcritical Hopf bifurcation, homoclinic bifurcation. Furthermore, the system may admit monostable, bistable or tristable dynamics. There exist bistable regions in which the non-mosquito equilibrium coexists with a positive equilibrium, i.e. the extinction of both mosquitoes coexists with the "steady state persistence" of both mosquitoes, or the non-mosquito equilibrium coexists with a stable limit cycle, i.e. the extinction of both mosquitoes coexists with the "oscillatory persistence" of both mosquitoes. There also exist tristable regions in which two equilibria coexist with one limit cycle, i.e. the extinction of both mosquitoes coexists with the "steady state persistence" and the "oscillatory persistence" of both mosquitoes.

4.   Bogdanov-Takens bifurcation

In Section 3, we show the existence of the Hopf bifurcations through a single bifurcation parameter $b$. By using two parameters $a$ and $b$ as bifurcation parameters, we now show that a Bogdanov-Takens bifurcation may occur in the system (1.2).

To that end, we fix $\mu_1, \mu_2, \xi_1, \xi_2 (\xi_1 > 2\xi_2)$, and choose $a$ and $b$ satisfying the following conditions:

(ⅰ) $a = a^*$ such that $T(\bar{N}(a^*), a^*) = 0$ where $\bar{N}$ is the unique critical point of $B(N)$ which is defined as (2.2);

(ⅱ) $b = \bar{b}(\bar{N}(a^*))$.

Then system (1.2) admits a unique positive equilibrium $(\bar{w}, \bar{g})$ (at the saddle-node bifurcation point) defined by

and

Therefore, if conditions (ⅰ) and (ⅱ) are satisfied then the Jacobian matrix at $(\bar{w}, \bar{g})$ has a zero eigenvalue with multiplicity $2$. This suggests that (1.2) may admit a Bogdanov-Takens bifurcation. With similar procedures to those in ^[35,36,34], we confirm that $(w^*, b^*)$ is a cusp singularity of codimension $2$. Now we show that the system (1.2) may admit a Bogdanov-Takens bifurcation by simulations.

Figure 9 shows the two-parameter bifurcation diagrams on the $b-a$ plane for $\mu_1 = 1, \mu_2 = 1, \xi_1 = 4, \xi_2 = 0.51$ in Figure 9 (a) and $\mu_1 = 1, \mu_2 = 2, \xi_1 = 11, \xi_2 = 0.95$ in Figure 9 (b). In each panel of Figure 9, the blue curve represents the saddle-node bifurcation curve and the cyan curve represents the Hopf bifurcation curve. The saddle-node curve $\Gamma_1$ and the Hopf curve $\Gamma_2$ divide the $b-a$ plane into three parts, marked by "nonexistence" (no positive equilibrium), "stable" ($(w^+, g^+)$ is stable) and "unstable" ($(w^+, g^+)$ is unstable), respectively. The intersection point of the saddle-node curve $\Gamma_1$ and the Hopf curve $\Gamma_2$ is the Bogdanov-Takens bifurcation point marked by "BT" (again by using package MatCont). A degenerate bifurcation point marked by "GH" is the generalized Hopf bifurcation point where the first Lyapunov coefficient vanishes while the second Lyapunov coefficient does not vanish, and that is where the Hopf bifurcation changes from subcritical to supercritical. In Figure 9 (a), the Hopf bifurcation curve is a monotone curve, which implies that there is only one Hopf bifurcation point $b = b_1$ when using $b$ as bifurcation parameter (for example fixed $a = 20$). On the other hand, in Figure 9 (b) the Hopf bifurcation curve is not monotone which means that there may exist two Hopf bifurcation points $b = b_1$ and $b = b_2$ (for example fixed $a = 40$).

In Figure 9, the three different parts "nonexistence", "stable" and "unstable" on the $b-a$ plane, which are divided by the saddle-node curve $\Gamma_1$ and the Hopf curve $\Gamma_2$, may give us some insights into the release strategy of the sterile mosquitoes. Of course, it must be the best result to eradicate wild mosquitoes. From the point of economic benefits and the control of wild mosquito population, it is practical to choose an appropriate release rate $b$ for fixed values of mating rate $a$ and other parameters to assure that the positive equilibrium lies in the region of "nonexistence". Another advantage is that the release rate can be adjust or the mating rate can be disturbed according to the observation of the current states of the wild mosquitoes and the sterile mosquitoes. For example, the release sterile mosquitoes can be increased appropriately if a cyclical phenomenon is observed.

5.   Conclusions

Sterile insect technique (SIT) is an important method to control and eliminate the number of wild mosquito population. It is significant to use mathematical model to investigate the releasing strategy and estimate the number of sterile mosquitoes in SIT. In this paper, an ordinary differential equation model (1.2) for the dynamics of the wild mosquitoes and the sterile mosquitoes is investigated. By the bifurcation analysis, we find rich bifurcation structure in model (1.2), including saddle-node bifurcations, (subcritical/supercritical) Hopf bifurcations, homoclinic bifurcation, and Bogdanov-Takens bifurcations.

For different release rates of sterile mosquito and circumstance parameters, we find the system (1.2) admits monostable, bistable and tristable dynamics. When the release rate of sterile mosquitoes is large, only the non-mosquito equilibrium exists and it is globally stable, which implies the extinction of both mosquitoes. With the decrease of the release rate of sterile mosquitoes, there exist bistable regions in which the non-mosquito equilibrium coexists with a positive equilibrium, i.e. the extinction of both mosquitoes coexists with the "steady state persistence" of both mosquitoes, or the non-mosquito equilibrium coexists with a stable limit cycle, i.e. the extinction of both mosquitoes coexists with the "oscillatory persistence" of both mosquitoes. There also exist tristable regions in which two equilibria coexist with one limit cycle. In the bistable and tristable parameter regions, the asymptotic states of both mosquitoes are determined by the initial densities of the wild and sterile mosquitoes.

From the viewpoint of biology, the wild mosquitoes will be eliminated if the release rate of sterile mosquitoes is large, which is unrealistic sometime due to the cost. When the release rate of sterile mosquitoes lies in a intermediate level, whether the wild mosquito population is eliminated or not depends on the initial wild mosquitoes and sterile mosquitoes densities. Then, first and foremost, the densities of wild mosquitoes should be investigated before sterile mosquitoes are released to insure the numbers of both populations lie asymptotically in the region "nonexistence".

Acknowledgments

The authors wish to express their grateful thanks to the anonymous referee for their careful reading and valuable comments and suggestions which greatly improved this work. This work is partially supported by grants from National Science Foundation of China (11671327, 11701472, 11871060, 11871403), China Scholarship Council and US-NSF grants DMS-1715651.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendices

A.   Direction of Hopf bifurcation

When $b = b_k, k = 1, 2,$ we have $T(N^+) = 0$. Then the eigenvalues of the Jacobian matrix at $(w^+, g^+)$ are $\lambda_1 = \omega_0 i$ and $\lambda_2 = -\omega_0i$, where $\omega_0: = \omega(b_k) = \sqrt{\mbox{Det}(J(w^+, g^+))} = \sqrt{\frac{w^+}{1+N^+}F'(N^+)}.$

By using the transformation of $x = w-w^+, y = g-g^+$ to (1.2), we obtain

where

and

By applying the translation $u = x, v = -\frac{1}{\omega_0}(a_{11}x+a_{12}y)$ to (9.1), we obtain

where

Set

By the results in ^[37], the direction of the Hopf bifurcation is determined by the sign of $\sigma$. Therefore, we have the following result.

Theorem A.1. Assume all parameters $a, b, \mu_1, \mu_2, \xi_1, \xi_2$ are positive and $b = b_k, k = 1, 2$ are defined as in Theorems 2.4. Then a family of periodic solutions bifurcate from the equilibrium $(w^+, g^+)$. Furthermore, system (1.2) undergoes a supercritical Hopf bifurcation if $\sigma < 0$ and a subcritical Hopf bifurcation if $\sigma > 0$, where $\sigma$ is defined as in (A.3).