Bifurcation analysis of a wild and sterile mosquito model

Abstract: The bifurcation of an ordinary differential equation model describing interaction of the wild and the released sterile mosquitoes is analyzed. It is shown that the model undergoes a sequence of bifurcations including saddle-node bifurcation, supercritical Hopf bifurcation, subcritical Hopf bifurcation, homoclinic bifurcation and Bogdanov-Takens bifurcation. We also find that the model displays monostable, bistable or tristable dynamics. This analysis suggests that the densities of the initial wild mosquitoes and the released sterile ones determine the asymptotic states of both populations. This study may give an insight into the estimation number of the released sterile mosquitoes.


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: = C(N) a 0 w w + g − (µ 1 + ξ 1 (w + g)) w, dg dt = R(w) − (µ 2 + ξ 2 (w + g))g. (1.1) 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(·) is the release rate of the sterile mosquitoes; µ i and ξ 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: Here, c, c 0 , b are all constants. In case (i), both the mating rate and the release rate are constants. In case (ii), the release rate is a linear function, and in case (iii), the release rate form takes into account an Allee effect as in [29,7,30]. The mating rate is a Holling-II type function both in case (ii) and (iii). 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 (i). An explicit threshold release value for b was also established for this case. For case (ii), they gave some results about the stability of the positive equilibrium when µ 1 ≤ µ 2 and ξ 1 ≤ 2ξ 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 µ 1 = µ 2 and ξ 1 = ξ 2 for case (iii).
In the present paper, we will give a detailed follow-up work focused on the case (ii). Let a = a 0 c 0 , the dynamics of the interacting mosquitoes in case (ii) 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 ξ 1 ≤ 2ξ 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 ξ 1 > 2ξ 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.

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 (2.1) 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 Thenb > 0 and we have the following results.
Note that if b is regarded as a bifurcation parameter, then b =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 Det(J(w + , g + )) = w + 1 + N + F (N + ) > 0. and the trace of J(w + , g + ) is given by where T (N) is defined as 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. µ 1 ≤ µ 2 . In [7], the authors have shown that the equilibrium (w + , g + ) is locally asymptotically stable if ξ 1 ≤ 2ξ 2 . Now we will discuss the stability of (w + , g + ) for the case ξ 1 > 2ξ 2 . Define (2.5) Note that ξ 2 − 2ξ 1 < 0. Then if a < a 0 or a > a 0 and 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 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.
Proof. It is useful to note that aN 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 ≈ 0 which implies the local stability of (w + , g + ) when b ≈ 0.
When 0 b <b, the discussion is divided into the following cases: For these cases, we have T (N + ) < 0 which implies that the equilibrium (w + , g + ) is locally asymptotically stable. (ii) N * <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 <b which implies that the equilibrium (w + , g + ) is locally asymptotically stable when 0 < b < b 1 and unstable when b 1 < b <b.
Then we have the following results about the stability of the equilibrium (w + , g + ).
Theorem 2.4. Assume all parameters a, b, µ 1 , µ 2 , ξ 1 , ξ 2 are positive and a > µ Then the equilibrium (w + , g + ) is locally asymptotically stable when ξ 1 ≤ 2ξ 2 . When ξ 1 > 2ξ 2 , defineb as in (2.3), a 0 , 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 1 > 0. When a > a 0 and 1 > 0, the following cases may occur: 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.

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 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 Here, we use the facts that Then we obtain the following result. 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 µ 1 = 1, µ 2 = 1, ξ 1 = 4, ξ 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 ∞ 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 µ 1 = 1, µ 2 = 1, ξ 1 = 4, ξ 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 ± , g ± ) 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 µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 ∞ 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 µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 ± , g ± ) 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 µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 ∞ 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 µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 ± , g ± ) 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 µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 ∞ 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 µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 ± , g ± ) 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.

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).
Then system (1.2) admits a unique positive equilibrium (w,ḡ) (at the saddle-node bifurcation point) defined byw and Det(J(w,ḡ)) = 0, Trace(J(w,ḡ)) = 0. Therefore, if conditions (i) and (ii) are satisfied then the Jacobian matrix at (w,ḡ) 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 µ 1 = 1, µ 2 = 1, ξ 1 = 4, ξ 2 = 0.51 in Figure 9 (a) and µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 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 Γ 1 and the Hopf curve Γ 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 Γ 1 and the Hopf curve Γ 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 Γ 1 and the Hopf curve Γ 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.  Figure 9. The blue curve represents the saddle-node bifurcation curve and the cyan curve represents the Hopf bifurcation curve. The "BT" mark indicates a Bogdanov-Takens bifurcation point and the "GH" mark indicates a generalized Hopf point where the first Lyapunov coefficient vanishes while the second Lyapunov coefficient does not vanish, which indicates that it is nondegenerate, i.e. Hopf bifurcation changes from subcritical to supercritical. Parameters: (a) µ 1 = 1, µ 2 = 1, ξ 1 = 4, ξ 2 = 0.51; (b) µ 1 = 1, µ 2 = 2, ξ 1 = 11, ξ 2 = 0.95.

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". where