1. Introduction

For more than a century, human beings have attempted to control blood-feeding mosquitoes. This is because of the significant mortality and morbidity burden associated with mosquito-borne diseases (e.g., malaria, dengue fever, and West Nile virus), which are transmitted between humans via blood-feeding mosquitoes [7,38]. Various control approaches have been explored, which include the development of more effective drug treatments, vaccines, and vector (mosquito) suppression. Vector control measures include the elimination or reduction of their nesting places by draining stagnant water deposits. This measure has been very effective in places where the prevalence of the disease was not very high. Indoor and outdoor insecticide spraying has been applied for many years for controlling mosquito populations. Despite some mosquito-borne diseases have been successfully controlled in many regions through vector-targeted intervention such as insecticide-treated bed nets (ITNs) and indoor residual sprays (IRS), massive and long time spraying of adulticide is not recommended. Since they have commonly been chemically-based, the effectiveness of this measure has been hampered by the appearance of insecticide resistant vector strains [1,2,8]. The genetically-altered or transgenic mosquitoes may provide a new and potentially effective weapon to fight these mosquito-borne diseases [14,20,29]. Sterile Insect Technique (SIT) is an environmentally friendly alternative strategy that is gaining renewed interest for the control of mosquito populations. The technique involves the massive release of male mosquitoes (sterilized through radiological or chemical means) into the environment to mate with wild mosquitoes that are present in the environment. Female mosquitoes mating successfully with a sterile male will either not reproduce, or produce eggs which will not hatch. Repetitive releases of sterile mosquitoes or releasing a significantly large number of sterile mosquitoes may eventually wipe out or suppress a wild mosquito population [5,37,3].

Mathematical models (see[11,13,17,22,18,34,25,9,26,40,41] and the references therein) have been formulated to explore the interactive dynamics of wild and sterile mosquito populations, and potentially estimate the effectiveness of the control strategy of mosquitoes. In particular, the models in [9,26] are ODE systems under the homogeneous assumption for the wild and sterile mosquitoes. However, these models hardly take into account the stage structured life-history of the mosquitoes, which can have significant effects on their dynamics (see[17,10,27,39]). In fact, all mosquito species go through four distinct stages (egg, pupa, larva, and adult) during a lifetime. The first three stages occur in water, but the adult is an active flying insect. Only the female mosquito bites and feeds on the blood of humans or other animals. Murdoch et al. [31] have shown short-period population oscillations in abundance occur from the developmental lags between mosquito life-history stages. In recent years, several interesting mathematical models have been developed to investigate the dynamics of mosquitoes with stage structure, for instance, the discrete models in [28], the continuous time models in [16,27] and the delayed models in [15,24].

Based on the models in [9,26,32,15], in this paper, we assume that the wild mosquito population growth is stage-structured. Let $w(a,t)$ denote the density of wild mosquitoes at time $t$ of age $a$. Let $\tau$ be the maturation time for all mosquitoes (the total time from egg to adult). During the larvae maturation, they are subject only to the possibility of natural death. Let $\mu_0$ be the per-capita natural mortality. We have

$\label{larave} \dfrac{\partial w(a,t)}{\partial t}+ \dfrac{\partial w(a,t)}{\partial a} = -\mu_0 w(a,t), \ \ \ \ 0<a<\tau, t>0.$

(1)

Let $W(t) = \int_\tau ^{\infty} w(a,t)da$ be the wild adult mosquitoes population at time $t\geq 0$, $\mu_1$ the per-capita natural death rate of adult mosquitoes, and $\xi_1$ the density-dependent death rate of adult mosquitoes. Let $g(t)$ be the sterile mosquitoes population at time $t\geq 0$, $\mu_2$ and $\xi_2$ the density-independent and dependent death rates of sterile mosquitoes. The interactive dynamics of wild adult mosquitoes and sterile mosquitoes are governed by the following equations:

$\begin{split} \dfrac{\partial w(a,t)}{\partial t}+ \dfrac{\partial w(a,t)}{\partial a} & = - \left[ \mu_1 + \xi_1 (W(t)+g(t)) \right]w(a,t), \ \ \ \ a>\tau, t>0, \\ \dfrac{dg}{dt} & = B(\cdot) - \bigl[\mu_2 + \xi_2 (W(t) +g(t))\bigr]g(t), \ \ \ \ t>0, \end{split}$

where $B(\cdot)$ is the release rate of the sterile mosquitoes.

Assume that $w(\infty,t) = 0$. Integrating the first equation of the above system over the interval $[\tau,\infty)$ with respect to $a$ gives

$\label{interact} \begin{split} \dfrac{ dW(t)}{d t} & = w(\tau, t) - \left( \mu_1 + \xi_1 (W(t)+g(t)) \right)W(t), \ \ \ \ t>0, \\ \dfrac{dg}{dt} = &\ B(\cdot) - \bigl(\mu_2 + \xi_2 ( W(t) +g(t))\bigr)g(t), \ \ \ \ t>0. \end{split}$

(2)

Use the method of characteristics to solve (1) and get $w(a,t) = w (0,t-a ) e^{-\mu_0 a}$ for $0\leq a\leq \tau$ and $t\geq a$. In particular, $w(\tau, t) = w (0,t-\tau) e^{-\mu_0 \tau}$ for $t\geq \tau$. Supposed that the interactions between the two types of mosquitoes lead to egg-laying rate given by

$w(0,t) = \dfrac{a_0 W (t )}{W(t )+g(t)}W(t), \ \ \ \ \ \ t \geq 0,$

where $a_0$ is the number of wild offspring produced by per mating between the wild mosquitoes. Then it follows that

$w(\tau, t ) = \dfrac{a_0W (t-\tau )}{W(t-\tau )+g(t-\tau)} W(t-\tau)e^{-\mu_0\tau}, \ \ \ \ \ \ t\geq \tau.$

(3)

Substituting this formula into (2) and replacing $W(t)$ by $w(t)$ (for easing the notation) leads to the delayed model of the adult wild and sterile mosquitoes:

$\label{basis} \begin{split} \dfrac{dw(t)}{dt} = &\ \dfrac{a_0 e^{-\mu_0\tau} w^2 (t-\tau)}{w(t-\tau)+g(t-\tau)} - \left( \mu_1 + \xi_1 (w+g) \right)w, \ \ \ \ t\geq \tau, \\ \dfrac{dg(t)}{dt} = &\ B(\cdot) - \bigl(\mu_2 + \xi_2 (w+g)\bigr)g, \ \ \ \ t\geq \tau. \end{split}$

(4)

The factor $e^{-\mu_0\tau}$ in (4) stands for the survival rate of the immature mosquitoes who were born at time $t-\tau$ and still remain alive at the time $t$. The initial conditions for system (4) take the following form

$\label{initial} w(t_0+\theta) = \phi(\theta)>0, \ \ \ \ g(t_0+\theta) = \psi(\theta)>0,\ \ \ \ \ \theta\in [-\tau,0],$

(5)

where $t_0\geq \tau$, $\phi(\theta)$ and $\psi(\theta)$ are positive continuous functions for $\theta\in [-\tau,0]$.

The main aim of this paper is to investigate how the delay $\tau$ affects the dynamics of system (4) with different strategies of releasing rate of sterile mosquitoes. In Section 2, we first explore the model (4) with the constant release rate of sterile mosquitoes and analyze the effect of the delay $\tau$ on the dynamics of mosquito populations. Our results show that the delay does not affect the dynamics of the model. In Section 3, we investigate (4) with the release rate of sterile mosquitoes proportional to the wild mosquito population, and obtain the conditions on the stability of the positive equilibria and the occurrence of Hopf bifurcations for certain values of the delay. In Section 4, numerical examinations are provided to demonstrate the complexity of the model dynamics in the latter case. In Section 5, we provide a brief discussions on our findings, in particular on the impact of the time delay on the mosquito control measures.

---

2. Constant release rate of the sterile mosquitoes

In this section, we investigate the model (4) with constant release rate of the sterile mosquitoes, i.e., we take $B(\cdot) = b$ with $b>0$. Then (4) becomes

$\label{const} \begin{split} \dfrac{dw(t)}{dt} = &\ \dfrac{a w^2 (t-\tau)}{w(t-\tau)+g(t-\tau)} - \left( \mu_1 + \xi_1 (w+g) \right)w, \ \ \ \ \ \ a = a_0e^{-\mu_0\tau}, \\ \dfrac{dg(t)}{dt} = &\ b - \bigl(\mu_2 + \xi_2 (w+g)\bigr)g. \end{split}$

(6)

The initial conditions of system (6) satisfy (5). We first establish the following result:

Lemma 2.1. For any given initial data in (5) that are positive on $[t_0-\tau, t_0]$ for some $t_0\ge 0$, system (6) has a unique solution $(w(t), g(t))$, which are defined, positive and bounded for all $t\ge t_0-\tau$.

Proof. The local existence and uniqueness of the solution $(w(t),g(t))$ of (6) subject to a prescribed initial condition (5) follows from the standard theory for delay differential equations (e.g., Theorem 3.1 in [19]). The global existence of this solution in $[t_0-\tau, +\infty)$ follows from the positiveness and boundedness of this solution that we prove below.

Let $[t_0-\tau,T)$ be the maximal interval of existence of $(w(t),g(t))$, and let

$K_1> \max\Big \{\dfrac{a-\mu_1}{\xi_1}, \max\limits_{\theta\in [-\tau, 0]}\phi(\theta) \Big \}, \ \ \ \ \ \ K_2> \max\Big\{\dfrac{b}{\mu_2}, g(t_0) \Big\}.$

We claim that $( w(t),g(t)) \in\Omega_1 : = \left\{(w,g): \,\, 0<w< K_1, \,\,\,0<g <K_2 \right\}$ for $[t_0-\tau,T)$. First, since $w(t_0)>0$ and $g(t_0)>0$, and for $t\in [t_0,T)$

$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ w(t)\\ = e^{-\int_{t_0}^t [\mu_1+\xi_1(w+g)]\,ds } \left(w(t_0) + \int_{t_0}^t \frac{a w^2(s-\tau)}{w(s-\tau)+g(s-\tau)} e^{-\int_{t_0}^s [\mu_1+\xi_1(w+g)]\,d\eta}\,ds \right),\\ \ \ \ \ \ \ \ \ \ \ g(t) = e^{-\int_{t_0}^t [\mu_2+\xi_2(w+g)]\,ds } \left(g(t_0) +b \int_{t_0}^t e^{-\int_{t_0}^s [\mu_2+\xi_2(w+g)]\,d\eta}\,ds \right),$

it follows that $w(t)>0$ and $g(t)>0$ for $t\in [t_0,T)$.

Now we prove that $w(t)<K_1$. If not, there would exist the smallest $t_1>t_0$ such that $w(t_1) = K_1$, $w'(t_1)\geq 0$ and $w(t) < K_1$ for any time $t_0-\tau \leq t<t_1$. However, from the first equation of system (6) we have

$w'(t_1) < a w (t_1-\tau) - \left( \mu_1 + \xi_1 w(t_1) \right)w(t_1) < (a-\mu_1-\xi_1K_1)K_1 < 0,$

a contradiction. This shows that $w(t) < K_1$ for $t \in [ t_0,T)$. Similarly using the second equation of (6) we show that $g(t)<K_2$ for $t\in [t_0,T)$.

We thus conclude the above claim. Since $\Omega_1$ is bounded, it follows that $T = \infty$ and $(w(t),g(t))\in \Omega_1$ for $t\geq t_0$.

We note that the above proof implies that for any $K_1>\dfrac{a-\mu_1}{\xi_1}$ and $K_2>\dfrac{b}{\mu_2}$, the set $\Omega$ defined in the above proof is positively invariant for the flows of (6).

Since the presence of the delay does not change the number of the equilibrium solutions in system (6), the existence of the equilibria of (6) follows from the same argument as for ODE system in [9]. Let

$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \bar N = \bar N(\tau):\\ = \frac{1}{3\xi_1\xi _2}\left(\sqrt{\bigl(\xi_2(a-\mu_1) - \xi_1\mu_2\bigr)^2+ 3\xi_1\xi _2 (a-\mu_1)\mu_2} + \xi_2(a-\mu_1) - \xi_1\mu_2\right),$

and define the threshold release value of the sterile mosquitoes as

$b_0(\tau): = \dfrac{1}{2a} \left(\xi_1\xi _2 \bar N^2 + (a-\mu_1)\mu_2\right) \bar N.$

We have the following result similar to the one in the reference [9].

Theorem 2.2. Let $\tau\geq 0$ and $b>0$. Then the system (6) has a unique boundary nonnegative equilibrium $E_0(0, g^0)$ where $g^0$ is given in (7). Furthermore, regarding the existence of positive equilibria, we have:

(i) For $b>b_0(\tau)$, system (6) does not have any positive equilibrium;

(i) For $b = b_0(\tau)$, system (6) has a unique positive equilibrium $E^*(w^*, g^*)$, where

$g^0 = \dfrac{\sqrt{\mu_2^2 + 4b\xi_2} - \mu_2}{2\xi_2},\ \ w^* = \dfrac{ \bar N(\mu_1 + \xi_1 \bar N)}{a}, \ \ \ g^* = \dfrac{b}{\mu_2 + \xi_2 \bar N};$

(7)

(iii) For $0<b<b_0(\tau)$, system (6) has exactly two positive equilibria $E_1^*(w_1^*, g_1^*)$ and $E_2^*(w_2^*, g_2^*)$ given by:

$w_1^* = \dfrac{N_1^*(\mu_1 + \xi_1 N_1^*)}{a} < w_2^* = \dfrac{N_2^*(\mu_1 + \xi_1 N_2^*)}{a}, \\ \ \ \ \ \ \ \ \ \ g_2^* = \dfrac{b}{\mu_2 + \xi_2 N_2^*} < g_1^* = \dfrac{b}{\mu_2 + \xi_2 N_1^*},$

where $0<N_1^*<N_2^*$ are the positive roots of

$P(N) : = \xi_1\xi_2 N^3- [(a-\mu_1)\xi_2- \xi_1\mu_2]N^2 -(a-\mu_1)\mu_2 N +ab = 0.$

We now investigate the stability of the equilibria of system (6). Upon linearizing system (6) at the equilibrium $E_0$, $E_1^*$, or $E_2^*$, we obtain the characteristic equation

$\label{charac} det(\lambda I- A - e^{-\lambda \tau} B ) = 0,$

(8)

where

$A = \begin{pmatrix} -(\mu_1 +\xi_1(2w+ g ))&-\xi_1w \\ -\xi_2g&- (\mu_2+ \xi_2(w+2 g)) \end{pmatrix},\\ \ \ \ \ \ \ \ \ B = \begin{pmatrix}\dfrac{aw(w+2g)}{(w+g)^2} &-\dfrac{aw^2 }{(w+g)^2} \\ 0& 0 \end{pmatrix}.$

Direct calculations give

$\lambda I- A - e^{-\lambda \tau} B = \begin{pmatrix} \lambda+ \mu_1 +\xi_1(2w+ g )-e^{- \lambda \tau}\dfrac{aw(w+ 2g)}{(w+g)^2} & \xi_1w +e^{- \lambda \tau}\dfrac{aw^2}{(w+g)^2}\\ \xi_2g&\lambda+ \mu_2+ \xi_2(w+2 g) \end{pmatrix}.$

At the boundary equilibrium $E_0$, it is easy to obtain that $B$ is zero matrix and the characteristic equation (8) has only two roots $\lambda_1 = - (\mu_1 +\xi_1 g^0)<0,\ \lambda_2 = - (\mu_2 +\xi_2 g^0)<0$. So $E_0$ is locally asymptotically stable for any $\tau>0$.

At the positive equilibria $E_1^*$ and $E_2^*$ (when $0<b<b_0$), by direct calculations we reduce the characteristic equation (8) as:

$\label{characposit1} L(\lambda): = \lambda^2+ a(\tau) \lambda + b(\tau) \lambda e^{-\lambda \tau} + c(\tau) +d(\tau) e^{-\lambda \tau} = 0,$

(9)

where, with $N = N_1^*$ or $N_2^*$,

$\label{characcoeff} \begin{split} a(\tau) & = \mu_1+\xi_1 N + \mu_2+\xi_2 N + \dfrac{ \xi_1 N(\mu_1+\xi_1N)}{a}+\dfrac{b\xi_2}{ \mu_2+\xi_2N },\\ b(\tau) & = -\bigg[ \mu_1+\xi_1N + \dfrac{b(\mu_1+\xi_1N)}{N(\mu_2+\xi_2N)} \bigg],\\ c(\tau)& = ( \mu_1+\xi_1N)(\mu_2+\xi_2 N)+\dfrac{\xi_1 N(\mu_1+\xi_1N)(\mu_2+\xi_2N)}{a}+ \dfrac{b\xi_2 ( \mu_1+\xi_1 N)}{\mu_2+\xi_2N}, \\ d(\tau)& = -\bigg[( \mu_1+\xi_1N)(\mu_2+\xi_2 N)+\dfrac{b(\mu_1+\xi_1N) }{N}+ \dfrac{2b\xi_2 ( \mu_1+\xi_1 N)}{\mu_2+\xi_2N }\bigg ]. \end{split}$

(10)

Recall that for $\tau = 0$, the system (6) becomes the system (2.2) in[9]. From Theorem 2.2 in [9], we know that the positive equilibrium $E_1^*$ is a saddle, and $E_2^*$ is a stable node or focus. The following results show that the stabilities of these equilibria remain the same for $\tau>0$.

Theorem 2.3. Let $\tau \geq 0$.

(i) For any $b>0$, $E_0$ is locally asymptotically stable.

(ii) For $0< b<b_0(\tau)$, $E_1^*$ is unstable and $E_2^*$ is locally asymptotically stable.

Proof. The stability of $E_0$ is already proved above. So we only need to prove (ⅱ). We shall complete the proof by the following two steps.

Step 1. Show that $E_1^*$ is unstable. It follows from (9) and (10) with $N = N_1^*$ that

$\label{characsign1} \begin{split} L(0)& = c (\tau)+d (\tau)\\ & = ( \mu_1+\xi_1N_1^*)(\mu_2+\xi_2 N_1^*)+\dfrac{\xi_1 N_1^*(\mu_1+\xi_1N_1^*)(\mu_2+\xi_2N_1^*)}{a}+ \dfrac{b\xi_2 ( \mu_1+\xi_1 N_1^*)}{\mu_2+\xi_2N_1^*}\\ &-\bigg[( \mu_1+\xi_1N_1^*)(\mu_2+\xi_2 N_1^*)+\dfrac{b(\mu_1+\xi_1N_1^*) }{N_1^*}+ \dfrac{2b\xi_2 ( \mu_1+\xi_1 N_1^*)}{\mu_2+\xi_2N_1^* }\bigg]\\ & = \dfrac{\xi_1 N_1^*(\mu_1+\xi_1N_1^*)(\mu_2+\xi_2N_1^*)}{a}-\dfrac{b(\mu_1+\xi_1N_1^*) }{N_1^*}-\dfrac{b\xi_2 ( \mu_1+\xi_1 N_1^*)}{\mu_2+\xi_2N_1^* }\\ & = \dfrac{\xi_1 N_1^*(\mu_1+\xi_1N_1^*)(\mu_2+\xi_2N_1^*)}{a}-\dfrac{b(\mu_1+\xi_1N)(\mu_2+2\xi_2N ) } {N(\mu_2+\xi_2N_1^*) }\\ & = \dfrac{ \mu_1+\xi_1N_1^*}{aN_1^* }\bigg[\xi_1 (\mu_2+ \xi_2N_1^* ) (N_1^*)^2 - \dfrac{ab(\mu_2+2\xi_2N_1^*) } { \mu_2+\xi_2N_1^* }\bigg ]\\ & = \dfrac{ \mu_1+\xi_1N_1^*}{aN_1^* }P'(N_1^*). \end{split}$

(11)

The above last expression is obtained by the same calculations for $det J_1$ on the page 1790 of [9] (indeed, $L(0)$ is equal to $det J_1(E_1^*)$ in [9] with $a: = a_0e^{-\mu_0\tau}$). Since $P'(N_1^*)<0$, it follows that $L(0) = c(\tau)+d(\tau)<0$. Note that, for any $\tau> 0$, $L(\lambda)>0$ for sufficiently large $\lambda>0$. Thus, there is at least one positive root for $L(\lambda) = 0$. This yields that $E_1^*$ is unstable for $\tau> 0$.

Step 2. $E_2^*$ is stable when $\tau>0$. To do so, we need the following Proposition (for the clearness of the proof we write $L_{\tau}(\lambda)$ for $L(\lambda)$ for $\tau \geq 0$):

Proposition. For each $\tau>0$, $L_{\tau}(\lambda) = 0$ has countably many roots, all lying to the left of a vertical line $Re \lambda = \beta$ for some real number $\beta$; moreover, for any $\alpha<\beta$ there are at most finitely many of these roots lying in the vertical strip $\alpha< Re \lambda <\beta$.

Using this proposition, the Rouché's theorem [12], and the fact that for $\tau = 0$, $L_{0}(\lambda) = 0$ has only two roots with negative real parts, we conclude that for some sufficiently $\delta>0$ and any $0<\tau<\delta$, all the roots of $L_{\tau}(\lambda) = 0$ has negative real parts (NRP). For simplicity, we say $L_{\tau}$ has the NRP-property if all its roots have negative real parts. Recall we are trying to show that for any $\tau>0$, $L_{\tau}$ has the NRP-property. Assume by the contradiction that this is false. Then, letting

$\tau_0 = {{\text{inf}}\{\tau>0: L_{\tau} \ {\text{does not have the NRP-property}}}\},$

we have that $\tau_0\geq \delta >0$, and $L_{\tau_0}$ has at least one root, say, $\lambda(\tau_0)$, such that $Re \lambda(\tau_0) \geq 0$. Using the above property, and the Rouché's theorem we can exclude the case that $Re \lambda(\tau_0) > 0$, yielding $Re \lambda(\tau_0) = 0$. Hence we may assume that $\lambda(\tau_0) = i\omega$ for some $\omega\geq 0$ and then insert it into (9) with $\tau = \tau_0$ to get

$\label{imcharc} -\omega^2+ a(\tau_0) \omega i + b(\tau_0) \omega e^{-i\omega\tau_0} + c(\tau_0) +d(\tau_0) e^{-i\omega \tau_0} = 0.$

(12)

Taking the absolute values of both sides yields

$\label{characim} \omega^4- (b^2(\tau_0) + 2 c(\tau_0) -a^2(\tau_0) )\omega^2+ c^2(\tau_0)-d^2(\tau_0) = 0.$

(13)

From (10), we have $c(\tau_0) - d(\tau_0) >0$. Note that by the same argument used in (11) we have $c(\tau_0) + d(\tau_0) = \dfrac{ \mu_1+\xi_1N_2^*}{aN_2^* }P'(N_2^*)>0$. Thus we conclude that $c^2(\tau_0)-d^2(\tau_0)>0$.

Now we show $b^2(\tau_0) + 2 c(\tau_0) -a^2(\tau_0) <0$. To simplify the expressions in the following calculations, we let $A_1 = \mu_1+\xi_1N$ and $A_2 = \mu_2+\xi_2N.$ It follows from direct calculations that

$\begin{split} B_1&\stackrel{\it def}{ = }b^2(\tau_0) + 2 c(\tau_0) -a^2(\tau_0) = (b(\tau) -a(\tau_0))(b(\tau_0) +a(\tau_0)) + 2 c(\tau_0)\\ & = -\bigg[2 A_1+\dfrac{bA_1}{NA_2}+ A_2 + \dfrac{\xi_1NA_1}{a}+\dfrac{b\xi_2}{ A_2 }\bigg]\bigg[ A_2 + \dfrac{\xi_1 NA_1}{a}+\dfrac{b\xi_2}{ A_2} -\dfrac{bA_1}{NA_2}\bigg] \\& +2A_1A_2 +\dfrac{2\xi_1 NA_1A_2}{a}+ \dfrac{2b\xi_2 A_1}{A_2} \\ & = -\bigg[ \dfrac{bA_1}{NA_2}+ A_2 + \dfrac{\xi_1NA_1}{a}+\dfrac{b\xi_2}{ A_2 }\bigg] \cdot \bigg[ A_2+ \dfrac{ \xi_1 NA_1}{a}+\dfrac{b\xi_2}{ A_2} -\dfrac{bA_1}{NA_2}\bigg] \\& +\dfrac{2bA_1^2}{ NA_2 }- \dfrac{2\xi_1NA_1^2}{a}+ \dfrac{2\xi_1NA_1A_2}{a}, \end{split}$

and so

$\begin{split} B_1& = \dfrac{2bA_1^2}{ NA_2 }- \dfrac{2\xi_1NA_1^2}{a}+ \dfrac{2\xi_1NA_1A_2}{a} - \dfrac{bA_1}{N}-\dfrac{b\xi_1A_1^2}{aA_2}-\dfrac{b^2\xi_2A_1}{NA_2^2}+\dfrac{b^2A_1^2}{(NA_2)^2}\\ &- A_2^2-\dfrac{\xi_1NA_1A_2}{a}-b\xi_2 + \dfrac{bA_1}{N}- \dfrac{\xi_1NA_1A_2}{a}-\dfrac{(\xi_1NA_1)^2}{ a^2}- \dfrac{b\xi_2\xi_1NA_1}{aA_2}\\ &+\dfrac{b\xi_1NA_1^2}{aNA_2} - b\xi_2- \dfrac{b\xi_2\xi_1NA_1}{ aA_2}-\dfrac{(b\xi_2)^2}{A_2^2}+\dfrac{b^2\xi_2A_1}{NA_2^2}\\ & = \dfrac{2bA_1^2}{ NA_2 }- \dfrac{2\xi_1NA_1^2}{a} -\dfrac{b\xi_1A_1^2}{aA_2}+\dfrac{b^2A_1^2}{(NA_2)^2}- A_2^2-b\xi_2 -\dfrac{(\xi_1NA_1)^2}{ a^2}- \dfrac{b\xi_2\xi_1NA_1}{aA_2}\\ &+\dfrac{b\xi_1NA_1^2}{aNA_2} - b\xi_2- \dfrac{b\xi_2\xi_1NA_1}{ aA_2}-\dfrac{(b\xi_2)^2}{A_2^2}\\ & = \bigg[\dfrac{ bA_1}{ N }+\dfrac{b\xi_2A_1 }{ A_2 }-\dfrac{\xi_1NA_1A_2}{ a } \bigg]\bigg[\dfrac{2A_1}{ A_2 } + \dfrac{bA_1}{N A_2^2 } + \dfrac{N\xi_1A_1}{ aA_2 } \bigg] \\ & -\dfrac{2b\xi_2A_1^2}{ A_2^2 } -\dfrac{b^2\xi_2A_1^2}{ N A_2^3 }-\dfrac{ bN\xi_2\xi_1A_1^2}{ a A_2^2 } - A_2^2-2b\xi_2 - \dfrac{b\xi_2\xi_1NA_1}{aA_2}-\dfrac{(b\xi_2)^2}{A_2^2}. \end{split}$

(14)

From $c(\tau_0) +d(\tau_0) >0$, we have $\dfrac{ bA_1}{ N }+\dfrac{b\xi_2A_1 }{ A_2 }-\dfrac{\xi_1NA_1A_2}{ a } <0$. It follows from (14) that $B_1 = b(\tau_0)^2+2c(\tau_0)-a^2(\tau_0)<0$.

Now using the above estimates for $B_1$ and $c^2(\tau_0)-d^2(\tau_0)$, we conclude that (13) cannot hold for any real $\omega$. This contradiction shows $\tau_0 = \infty$, thereby completing the proof of Step 2 and the proof of Theorem 2.3.

In Figure 1 and Figure 2, we perform simulations to show that in system (6), the bistable phenomenon still occur and the time delay $\tau$ of stage-structure growth has only effect on the level of the positive equilibria with $\tau$ being varied.

---

3. Release rate proportional to wild mosquito population

In this section, we consider the system (4) with the release rate $B(\cdot)$ proportional to the wild mosquito population. As in the reference [9], we also incorporate the Allee effect [35] to account for the difficulty and stochasticity of finding mates when the populations of mosquitoes are small. Then the system (4) becomes:

$\label{proportion} \begin{split} \dfrac{dw}{dt} = &\ \dfrac{a_0e^{-\mu_0\tau} w^2 (t-\tau)}{1+w(t-\tau)+g(t-\tau)}- \left( \mu_1 + \xi_1 (w+g) \right)w,\\ \dfrac{dg}{dt} = &\ bw - \bigl(\mu_2 + \xi_2 (w+g)\bigr)g. \end{split}$

(15)

The initial conditions of system (15) satisfy (5). Note that when $\tau = 0$, the system reduces to the system (3.1) in [9]. As in the previous section, we investigate how the delay $\tau$ affects the dynamics of (15) when varying $\tau$.

First, by a similar proof for Theorem 2.1 we have the following result:

Theorem 3.1. For any given initial data in (5) that are positive on $[t_0-\tau, t_0]$ for some $t_0\ge 0$, system (15) has a unique solution $(w(t), g(t))$, which are defined, positive and bounded for all $t\ge t_0-\tau$.

Again use the fact that the number of the steady state solutions of the system (15) with $\tau>0$ are the same as that with $\tau = 0$. Hence we can apply Theorem 3.1 in [9] (with $a = a_0e^{-\mu_0 \tau}$) to obtain Theorem 3.2 below.

Let $a_0e^{-\mu_0 \tau} > (\sqrt \mu_1 +\sqrt\xi_1)^2$. Set

$G(N): = \frac{\bigl(a_0e^{-\mu_0 \tau} N - (1+N)(\mu_1+\xi_1N)\bigr)(\mu_2+\xi_2N)}{(1+N)(\mu_1+\xi_1N)},$

and

$N_{1,2} = \frac{1}{2\xi_1}\left((a_0e^{-\mu_0 \tau}-\mu_1-\xi_1)\pm \sqrt{(a_0e^{-\mu_0 \tau}-\mu_1-\xi_1)^2-4 \mu_1\xi_1}\right).$

(16)

Let $\bar N$ be the point in $(N_1,N_2)$ such that $G'(\bar N) = 0$. We then define the threshold release value of the sterile mosquitoes as

$b_0(\tau): = G(\bar N).$

Theorem 3.2. Let $\tau\geq 0$, $b>0$, and $a_0e^{-\mu_0 \tau}> (\sqrt \mu_1+\sqrt\xi_1)^2$. Then the system (15) has a unique nonnegative boundary equilibrium $E_0(0, 0)$. Furthermore, regarding the existence of positive equilibria, we have:

(i) For $b>b_0(\tau)$, system (15) does not have any positive equilibrium;

(ii) For $b = b_0(\tau)$, system (15) has a unique positive equilibrium $E^*(w^*, g^*)$, where

$g^0 = \dfrac{\sqrt{\mu_2^2 + 4b\xi_2} - \mu_2}{2\xi_2},\ \ w^* = \dfrac{\bar N(\mu_1 + \xi_1 \bar N)}{a}, \ \ \ \ g^* = \dfrac{b}{\mu_2 + \xi_2 \bar N};$

(iii) For $0<b<b_0(\tau)$, system (15) has exactly two positive equilibria $E_1^*(w_1^*,g_1^*)$ and $E_2^*(w_2^*,g_2^*)$ given by, for $i = 1, 2$,

$\begin{split} w_i^* = \dfrac{(1+N_i^*)(\mu_1 + \xi_1 N_i^*)}{a_0e^{-\mu_0 \tau} } ,\ \ \ \ \ g_i^* = \dfrac{bw_i^*}{\mu_i +\xi_2 N_i^*}, \end{split}$

(17)

where $N_1^*<N_2^*$ are the positive roots of

$\begin{split} F(N): = & \xi_1\xi _2 N^3 +\bigl(\xi_1(b+\mu_2)+\xi_2(\mu_1+\xi_1-a)\bigr)N^2\\ & +(\mu_1\xi_2+b(\mu_1+\xi_1)+ ( \mu_1+\xi_1-a) \mu_2) N +\mu_1(b+\mu_2), \end{split}$

(18)

where $a = a_0e^{-\mu_0 \tau}$. (Note that the dependence of $N_i^*$ and $(w_i^*,g_i^*)$ on $\tau$ are suppressed.)

To determine the stability of (15) at its equilibria, we again calculate the characteristic equation $H(\lambda): = det(\lambda I- A - e^{-\lambda \tau} B) = 0$ of (15) at these equilibria $(w,g) = (0,0), (w_1^*,g_1^*)$ or $(w_2^*,g_2^*)$, where

$\begin{array}{l} A = \left( {\begin{array}{*{20}{c}} { - ({\mu _1} + {\xi _1}(2w + g))}&{ - {\xi _1}w}\\ {b - {\xi _2}g}&{ - ({\mu _2} + {\xi _2}(w + 2g))} \end{array}} \right),\;\\ \;\;\;\;\;\;B = \left( {\begin{array}{*{20}{c}} {\frac{{aw[w + 2(1 + g)]}}{{{{(1 + w + g)}^2}}}}&{ - \frac{{a{w^2}}}{{{{(1 + w + g)}^2}}}}\\ 0&0 \end{array}} \right). \end{array}$

Since $B$ is zero matrix at $E_0$, it follows that the eigenvalues are $-\mu_1$ and $-\mu_2$, and so $E_0$ is locally asymptotically stable for any $\tau\geq 0$.

To analyze the stability at $E_1^*$ and $E_2^*$, we need to compute $H(\lambda) = 0$, which becomes after some computations:

$\label{charac2} H(\lambda): = \lambda^2+ a_1(\tau) \lambda + b_1(\tau) \lambda e^{-\lambda \tau} + c_1(\tau) +d_1(\tau) e^{-\lambda \tau} = 0,$

(19)

where, with $(w,g) = (w_1^*,g_1^*)$ or $(w_2^*,g_2^*)$,

$\begin{split} a_1(\tau) & = \mu_2+\xi_2 (w+2g) + \mu_1+\xi_1 (2w+ g),\\ b_1(\tau)& = -\dfrac{aw}{1+ w+g }-\dfrac{aw(1+g)}{(1+w+g)^2},\\ c_1(\tau)& = (\mu_2+\xi_2 (w+2g))( \mu_1+\xi_1 (2w+ g)) +\xi_1 w( b-\xi_2g) \\& = (\mu_2+\xi_2 (w+2g))( \mu_1+\xi_1 (2w+ g)) +\xi_1 ( \mu_2+\xi_2g) g ,\\ d_1(\tau)& = -\bigg[\dfrac{aw^2+2aw(1+g)}{(1+w+g)^2}(\mu_2+\xi_2 (w+2g))+(\xi_2g-b) \dfrac{aw^2 }{(1+w+g)^2}\bigg]. \end{split}$

(20)

Theorem 3.3. Let $\tau\geq 0$, $a_0e^{-\mu_0 \tau}> (\sqrt \mu_1 +\sqrt\xi_1)^2$, and $b<b_0(\tau)$. Then the positive equilibrium $E_1^*$ of system (15) is unstable.

Proof. Let $N = N_1^*$ and $(w,g) = (w_1^*, g_1^*)$ in (19) and (20). From (19), we see $H(0) = c_1(\tau) +d_1(\tau)$. From the calculations on pages 1972 and 1806-1807 in [9], we know that $H(0) = c_1(\tau) +d_1(\tau) = \dfrac{w_1^*}{1+N_1^*}F'(N_1^*) <0$. Noting that $\lim_{\lambda \to \infty}H(\lambda) = \infty$, we conclude that $H(\lambda) = 0$ has at least one positive root, and so $E_1^*$ is unstable.

Next we investigate the stability of $E_2^*$. In the rest of this section we let $N = N_2^*$ and $(w,g) = (w_2^*, g_2^*)$ in (19) and (20). We are interested in the situation that when $\tau = 0$ and $\tau>0$ is small, $E_2^*$ is stable, and as $\tau$ increases, $E_2^*$ may lose the stability via a Hopf bifurcation. When $\tau = 0$, the characteristic equation (19) becomes

$\label{charac3} \lambda^2+ (a_1(0) +b_1(0)) \lambda + + c_1(0) +d_1(0) = 0.$

(21)

From Theorem 3.2 in [9], we know in this case that $E_2^*$ is asymptotically stable if and only if $a_1(0)+b_1(0) = trace(A(0)+B(0))<0$ (since $c_1(0)+d_1(0)>0$), which is equivalent to (see [9]):

$({\bf{H_0}}): \ \ \ \ \ \ \mu_1-\mu_2+(\xi_1-2\xi_2)N_2^* - (\xi_1-\xi_2)w_2^* - \frac{a_0 (w_2^*)^2}{(1+N_2^*)^2} < 0 \mbox{ for } \tau = 0.$

(Note that a sufficient condition for ($H_0$) is $\mu_1 \le \mu_2$ and $\xi_1\le 2\xi_2$ from [9].) When ($H_0$) holds, all roots of (21) have negative real parts.

As $\tau$ increases, the only way that $E_2^*$ may lose its stability is that some characteristic roots of (19) cross the imaginary axis and move to the right-half complex plane ([6], P.1146(ⅱ)). To see the possibility of this happening, we use the similar analysis as we did for the characteristic equation (9). We assume that (19) has a pair of purely imaginary roots $\lambda = \pm i\omega$ $(\omega>0)$ for some $\tau>0$. Substituting $\lambda = i\omega$ into (19) and separating the real and imaginary parts yield

$\begin {split} -\omega^2 +b_1(\tau)\omega \sin (\omega \tau) + c_1(\tau) + d_1(\tau) \cos(\omega \tau) = 0,\\ a_1(\tau) \omega +b_1(\tau)\omega \cos( \omega \tau ) -d_1 (\tau) \sin (\omega \tau) = 0. \end {split}$

Solving for $\cos(\omega\tau)$ and $\sin(\omega\tau)$, we obtain

$\label{sincos} \begin {split} \sin(\omega \tau) = \dfrac{d_1(\tau)a_1(\tau)\omega+ b_1(\tau)\omega (\omega^2- c_1(\tau))}{b_1^2(\tau)\omega^2 + d_1^2(\tau)},\\ \cos(\omega \tau) = \dfrac{d_1(\tau)(\omega^2- c_1(\tau))- b_1(\tau) a_1(\tau) \omega^2}{b_1^2(\tau)\omega^2 + d_1^2(\tau)}.\\ \end {split}$

(22)

Squaring and adding both equations of (22), we see that $\omega$ must satisfy

$\label{xuroot2} H(\omega) = \omega^4- (b_1^2(\tau) + 2 c_1(\tau) -a_1^2(\tau) )\omega^2+ c_1^2(\tau)-d_1^2(\tau) = 0.$

(23)

From (20), we have that $c_1(\tau) - d_1(\tau) >0$. Similar to that in the proof of Theorem 3.3, we have

$c_1(\tau) + d_1(\tau) = \dfrac{w_2^*}{1+N_2^*}F'(N_2^*) >0,$

and hence $c_1^2(\tau)-d_1^2(\tau)>0$. Let $B_2(\tau): = b_1^2(\tau) + 2 c_1(\tau) -a_1^2(\tau)$. Then using the quadratic formula, we conclude that if and only if

$({\bf{ H_1}}): \ \ \ \ \ \ B_2(\tau) \geq 2\sqrt{c_1^2(\tau)-d_1^2(\tau)},$

equation (23) has two positive roots $\omega_1$ and $\omega_2$ given by

$\begin{split} \omega_{\pm}(\tau) = \sqrt{\dfrac{ B_2(\tau)\pm \sqrt{\Delta(\tau)}}{2}}, \ \ \ \ \ \ \Delta(\tau): = (B_2(\tau))^2-4(c_1^2(\tau)-d_1^2(\tau)). \end{split}$

(24)

Consequently, if $B_2(\tau) < 2\sqrt{c_1^2(\tau)-d_1^2(\tau)}$ for all $\tau \in [0,\tau_0)$ with some $\tau_0>0$, then (23) does not hold for any real $\omega$, which yields that $E_2^*$ is stable; If $(H_1)$ holds for some $\tau^*\in (0,\tau_0)$ and,

$({\bf{H_2}}): \ \ \ \ \ \ \left\{ \begin{array}{ll} \mbox{for }\ i = + \mbox{or } \ -, \\ \displaystyle \sin (\tau^* \omega_i(\tau^*)) = h_1(\omega_i(\tau^*),\tau^*),\ \ \ \ \ \ \cos (\tau^* \omega_i(\tau^*)) = h_2(\omega_i(\tau^*),\tau^*), \\ \displaystyle \frac{d}{d\tau} Re \lambda(\tau^*)\neq 0 \end{array}\right.$

where $h_1(\omega,\tau)$ and $h_2(\omega,\tau)$ are the right-hand sides of (22) respectively, and $\lambda(\tau)$ is a root of (23) (which is a simple root for $\tau$ in a neighborhood of $\tau^*$ with $Re \lambda(\tau^*) = \omega_i(\tau^*)$), then by the Hopf Bifurcation Theorem (see Kuang [23]), system (15) has a Hopf bifurcation at $\tau = \tau^*$. In summary, we have the following:

Theorem 3.1. Let $\tau_0>0$. Assume that

$0< b<b_0(\tau), \ \ \ \ \ \ (\sqrt \mu_1 +\sqrt\xi_1)^2< a_0 e^{-\mu_0 \tau}, \ \ \ \ \ \ \forall\, \tau \in [0,\tau_0),$

and ($H_0$) holds. We have:

(i) If $B_2(\tau) < 2\sqrt{c_1^2(\tau)-d_1^2(\tau)}$ for $\tau \in [0,\tau_0)$, then the positive equilibrium $E_2^*$ of system (15) is locally asymptotically stable for every $\tau \in [0,\tau_0)$.

(ii) If ($H_1$) and ($H_2$) hold for some $\tau^*\in (0,\tau_0)$, then the system (15) undergoes a Hopf bifurcations at $E_2^*$ when $\tau = \tau^*$.

In theorem 3.1 (ⅱ), the conditions $(H_0)-(H_2)$ guarantee the occurrence of the Hopf bifurcation at some $\tau^*$. However, since we do not have the explicit formulas of $E_2^* = (w_2^*, g_2^*)$ for the general parameters in (15), these conditions are hard to be verified analytically. Below we show that the Hopf bifurcations do occur for some parameter ranges.

To this end, we consider the system (15) with the parameters satisfying:

$a_0 = 1, \ \ \ \ \ \ \mu_1 = \frac14 e^{-\mu_0 \tau} -\xi_1, \ \ \ \ \ \ \mu_2 = b- \xi_2.$

Then (15) becomes

$\label{hopfequation} \begin{split} \dfrac{dw}{dt} = &\ \dfrac{ e^{-\mu_0 \tau} w^2 (t-\tau)}{1+w(t-\tau)+g(t-\tau)}- \left( \dfrac14 e^{-\mu_0 \tau} -\xi_1 + \xi_1 (w+g) \right)w,\\ \dfrac{dg}{dt} = &\ bw - \bigl(b- \xi_2 + \xi_2(w+g)\bigr)g. \end{split}$

(25)

It is easy to verify that $\left(\dfrac12,\dfrac12\right)$ is always a positive equilibrium of system (25). In order to present our results in the following theorem, we introduce the notations below:

$\begin{split} a_1^0 = b+ \dfrac{1}{4}+\dfrac12(\xi_1+\xi_2),\ \ \ b_1^0 = - \dfrac{7}{16}, \ \ \ c_1^0 = \dfrac{1}{4}b+b\xi_1 +\dfrac18\xi_2,\ \ \ d_1^0 = -\dfrac{1}{8} (3b+ 2\xi_2), \end{split}$

(26)

and

$\begin{align}\label{coeffh''} \begin{split}&\Delta_1^0: = -(a_1^0+b_1^0) = \dfrac{3}{16}- b - \dfrac 12 (\xi_1+ \xi_2), \ \ \ \ \ \ \Delta_2^0: = c_1^0+d_1^0 = (\xi_1-\frac18) b - \dfrac18 \xi_2, \\ & B_2^0: = 2c_1^0+( b_1^0)^2 -(a_1^0)^2 = 2 (\dfrac{1}{4}b+b\xi_1 +\dfrac18\xi_2) + (\frac{11}{16}+b+\frac12(\xi_1+\xi_2)) \Delta_1^0, \\ &\Delta^0: = (B_2^0)^2 -4[( c_1^0)^2 -(d_1^0)^2] = (B_2^0)^2 -\frac12 [(5+8\xi_1) b+ 3\xi_2] \Delta_2^0. \end{split} \end{align}$

(27)

Theorem 3.2. Assume that

$\Delta_1^0>0, \ \ \ \ \ \ \Delta_2^0>0, \ \ \ \ \ \ \Delta^0 >0.$

(28)

(i) Assume that $\mu_0 = 0$. Then the system (25) undergoes the Hopf bifurcations at $(\frac12, \frac12)$ for $\tau = \tau_k^{\pm}$ with $\lambda(\tau_k^{\pm}) = i\omega_{\pm}$, $\lambda'(\tau_k^{\pm})\neq 0$ and

$\frac{d}{d\tau} Re \lambda(\tau_{k}^{\pm}) = Re \frac{d}{d\tau}\lambda(\tau_{k}^{\pm}) = \frac{ \pm \sqrt{\Delta^0} |\lambda'(\tau_k^{\pm})|^2 }{(b_1^0)^2 \omega_{\pm}^2 +(d_1^0)^2 } \ne 0,$

(29)

where $k = 0, 1, 2, \cdots$,

$\begin{split} \omega_{\pm} = \sqrt{\dfrac{ B_2^0\pm \sqrt{\Delta^0}}{2}}, \end{split}$

(30)

$\tau_k^{\pm} = \left\{\begin{split}& \frac{1}{\omega_{\pm}} \Big[ \cos^{-1} \Big(\dfrac{d_1^0(\omega_{\pm}^2- c_1^0)- b_1^0 a_1^0 \omega_{\pm}^2}{(b_1^0)^2 \omega_{\pm}^2 + (d_1^0)^2} \Big)+ 2k\pi\Big], \ \ \ \ \mathit{if} \ \ \ \omega_{\pm}^2 \le \dfrac{b_1^0c_1^0-a_1^0d_1^0}{b_1^0},\\ & \dfrac{1}{\omega_{\pm}}\Big[ 2\pi- \cos^{-1} \Big(\dfrac{d_1^0(\omega_{\pm}^2- c_1^0)- b_1^0 a_1^0 \omega_{\pm}^2}{(b_1^0)^2 \omega_{\pm}^2 + (d_1^0)^2} \Big)+ 2k\pi\Big], \ \ \ \ \mathit{if} \ \ \ \omega_{\pm}^2 > \dfrac{b_1^0c_1^0-a_1^0d_1^0}{b_1^0}. \end{split}\right.$

(31)

Furthermore, $\tau_0^+<\tau_0^-$ and the equilibrium $(\frac12, \frac12)$ is locally asymptotic stable for $\tau\in [0,\tau_0^+)$.

(ii) Given $\xi_1>0$, $\xi_2>0$ and $b>0$ that satisfy (28), for any given $\tau_0> \tau_0^{\pm}$, if $\mu_0>0$ is sufficiently small, then the system (25) undergoes the Hopf bifurcations at $(\frac12, \frac12)$ for $\tau = \tilde \tau_k^{\pm} = \tau_k^{\pm}+O(\mu_0)$ with $k = 0, 1, \cdots, N$, where $N$ is the largest integer such that $\tau_N^{\pm}<\tau_0$. Furthermore, $\tilde \tau_0^+< \tilde \tau_0^-$ and the equilibrium $(\frac12, \frac12)$ is locally asymptotic stable for $\tau\in [0,\tilde \tau_0^+)$.

Proof. First note that the characteristic equation (25) at $(\dfrac12, \dfrac 12)$ is the equation (19), whose coefficients (20) become

$\begin{split} a_1(\tau) & = b+ \dfrac{1}{4}e^{-\mu_0 \tau} +\dfrac12(\xi_1+\xi_2), \ \ \ \ \ \ b_1(\tau) = - \dfrac{7}{16}e^{-\mu_0 \tau},\\ c_1(\tau)& = \dfrac{1}{4}b+b\xi_1 +\dfrac18\xi_2, \ \ \ \ \ \ d_1(\tau) = -\dfrac{1}{8} (3b+ 2\xi_2) e^{-\mu_0 \tau}. \end{split}$

(32)

When $\tau = 0$, we have $a_1 = a_1^0$, $b_1 = b_1^0$, $c_1 = c_1^0$ and $d_1 = d_1^0$. The condition $\Delta_1^0>0$ and $\Delta_2^0>0$ guarantee the stability of $(\dfrac12, \dfrac 12)$ when $\tau = 0$. This also implies that $E_2^* = (\dfrac12, \dfrac 12)$.

We prove (ⅰ) now. Since $\mu_0 = 0$, it follows that $a_1 = a_1^0$, $b_1 = b_1^0$, $c_1 = c_1^0$ and $d_1 = d_1^0$, and so $B_2 = B_2^0$. It is clear from the expression of $B_2^0$ in (27) that $B_2^0>0$. The condition $\Delta^0 >0$ guarantees the assumption ($H_1$), and so the positive roots of (23) are given by $\omega_{\pm}^0$ in (30). Since $\omega_{\pm}^0$ does not depend on $\tau$, one can compute $\tau_k^{\pm}$ explicitly from the systems (22) and get their formulas given in (31). Since the coefficients of (19) do not depend on $\tau$, we get the explicit formulas for $\frac{d}{d\tau} Re \lambda(\tau_{k}^{\pm})$ in (29) (see, e.g., [6]). It follows from the Hopf bifurcation theorem that the Hopf bifurcation occurs at $(\frac12, \frac12)$ for each $\tau = \tau_k^{\pm}$. Since $(\frac12, \frac12)$ is asumptotically stable for sufficiently small $\tau\geq 0$, and $\frac{d}{d\tau} Re \lambda(\tau_{0}^{+}) >0$ and $\frac{d}{d\tau} Re \lambda(\tau_{0}^{-}) <0$ in (29), it follows that the first bifurcation occurs at $\tau = \tau_0^+$. We thus conclude that $\tau_0^+<\tau_0^-$ and the equilibrium $(\frac12, \frac12)$ is locally asymptotic stable for $\tau\in [0,\tau_0^+)$. This shows (ⅰ). The assertions (ⅱ) follows from (ⅰ) and the implicit function theorem.

Corollary 1. Assume that $\xi_1$, $\xi_2$ and $b$ in (25) satisfy either

$\frac18< \xi_1 <\frac38, \ \ \ \mathit{and \ both}\ b>0\ and \ \xi_2>0 \ are \ sufficiently \ small \ with \ b>\dfrac{\xi_2}{8\xi_1 - 1} ,$

or

$0<b<\frac{3}{16}, \ \ \ and \ both \ \xi_1-\dfrac18 \ and \ \xi_2>0 \ are \ sufficiently \ small \ with \ \xi_1-\dfrac18> \dfrac{\xi_2}{8b} .$

Then the assumptions in (28) are satisfied. Consequently the assertions in Theorem 3.2 hold.

Proof. We verify the assumptions in (28). Under the either assumption of the corollary we have $\Delta_2^0 >0$. Since $\lim_{\xi_2\to 0}\Delta_1^0 = \frac{3}{16}- b-\frac12 \xi_1$, we see that $\Delta_2^0>0$ under the either assumption. Finally, if the first assumption holds, then

$\lim\limits_{(b,\xi_2)\to (0,0)} \Delta^0 = (\frac{11}{16}+\frac12 \xi_1) (\frac{3}{16}-\frac12 \xi_1) >0;$

if the second assumption holds, then

$\lim\limits_{(\xi_1,\xi_2)\to (\frac18,0)} \Delta^0 = (\frac{11}{16}+b) (\frac{3}{16}-b) >0.$

This shows that $\Delta^0>0$ under the either assumption. Applying Theorem 3.2 completes the proof of the corollary.

---

4. Numerical simulation of the model (15)

In this section, we will focus on the numerical simulation of the model (15). First, let the parameter values $a_0 = 30, \mu_0 = 0.1,\ \mu_1 = 0.5,\ \mu_2 = 1.5,\ b = 2,\ \xi_1 = 4, \xi_2 = 0.51$ in system (15). In Figure 3, we show that the time delay $\tau$ have effect on stability of the positive equilibrium $E_2^*$ and there is a stable periodic solution for $\tau = 4.9$. Then, we choose to use the same set of parameters as in Figure 3 except $\tau$ and $b$. We let $\tau$ vary and use it as a bifurcation parameter. Other parameter values are fixed except parameter $b$ in Figure 9, where $b$ is varying as a secondary free parameter to generate the two dimensional bifurcation diagram. Using Matlab, we obtain the bifurcation diagram (see Figure 4).

From the bifurcation diagram in Figure 4, we can see that there is a stable steady state for delay $\tau$ from 0 up to around 4.5. There exists a stable periodic solution for delay between 4.5 to 5.2. Interestingly, there is a discontinuity at $\tau\approx 5.2$. Actually after further exploration, we find the discontinuity occurs since there exist two stable periodic solutions for delay $\tau$ close to 5.2. One of them has a smaller amplitude and the other one has a larger amplitude. Due to the change of their attraction basin as delay $\tau$ varies, a solution which originally approaches to the periodic solution with a smaller amplitude, switches to the periodic solution with a larger amplitude. When this switch occurs, we see a jump in our bifurcation diagram in Figure 4. There is another interesting observation on the bifurcation diagram around $\tau\approx 6.8$ where the bifurcation diagram stop abruptly. We did further exploration there too and found the amplitude of the periodic solutions shrink abruptly within a very small interval of delay and our Matlab code fails to work appropriately due to a precision issue.

To demonstrate the existence of two stable periodic solutions, we choose delay $\tau = 5.18$ and choose two different sets of constant initial values $(w,g) = (10,5)$ and $(w,g) = (0.5,1.5)$. The solution associated with the initial values $(w,g) = (10,5)$, approaches to the periodic solution with a larger amplitude and is plotted in dotted line in Figure 5. The solution associated with the initial values $(w,g) = (0.5,1.5)$, approaches to the periodic solution with a smaller amplitude and is plotted in solid line in Figure 5.

Similarly there is a secondary discontinuity at delay $\tau\approx 6.5$. There also exist two stable periodic solutions for delay close to 6.5. To illustrate that, we choose $\tau = 6.55$ and two different set of constant initial values $(w,g) = (0.71,0.11)$ and $(w,g) = (0.5,1.5)$. Accordingly we obtain two stable periodic solutions which are plotted in Figure 6.

To have a better understanding about system (15), we use DDEbiftool to investigate the system further and obtain some bifurcation diagram. First we obtain a one dimensional bifurcation diagram using delay $\tau$ as a bifurcating parameter (see, Figure 7). In this figure we can see, the system has two positive equilibria which merge and disappear through a limit point bifurcation at around $\tau = 6.2$. For $\tau$ between 4.7 and 7.2, there are three periodic solutions. Some of them is stable and some of them is unstable. We calculate the stability of periodic solutions along these curves and place the result in Figure 8.

In Figure 8, it is easy to see that there exists two stable periodic solutions for delay between $\tau\in(4.8, 5.8)$. This observation is consistent with our simulation result in Figure 4. The first branch of periodic solution is stable for $\tau\in(4.55, 5.5)$. The second branch periodic solution is stable for delay $\tau\in(4.8, 5.7)$, and the third branch of periodic solution is stable from delay $\tau\in(5.55, 7.2)$. This bifurcation diagram confirms the existence of two stable periodic solutions and clearly indicate delay values where one can expect the two stable periodic solutions coexist.

In addition, via tracing the periodic branches, we generate a two dimensional bifurcation diagram (see, Figure 9) in the parameter space of $(\tau, b)$. From the bifurcation diagram, we can see the system undergoes supercritical Hopf bifurcation (the red curve), subcritical Hopf bifurcation (the cyan curve which is very close to red curves and very hard to see), the fold-Hopf bifurcation (blue curve) and torus bifurcation (the black curve). In addition, along Hopf bifurcation curves, other more degenerated points are detected including Hopf-Hopf bifurcation (six black circles) and Generalized Hopf bifurcation (one black square). At those degenerated points, Hopf bifurcation curves, Fold-Hopf bifurcation curves or torus bifurcation curves meet each other. We believe more interesting dynamics can be expected around these Hopf-Hopf points or Generalized Hopf point. For torus bifurcation, we also plot the profile of orbits on the torus at the bifurcation point and placed them in Figure 10. In summary, Figure 9 shows that very complicated dynamics can be expected from the model. In another follow up work, we will investigate the model further around the Hopf-Hopf bifurcation, Generalized-Hopf bifurcation points.

---

5. Discussion

In this work, we studied (both analytically and numerically) the dynamics of interactive systems of the wild and sterile mosquitoes with different releasing sterile mosquitoes strategies by incorporating the delay in the growth stage of wild mosquito populations. Our obtained results in Theorem 2.3 have shown that the growth time delay of the mosquito population has no destabilizing effect on the solution's behavior in the constant release rate of sterile mosquitoes. That is, in this case, for any size of the growth delay, the solutions of the system approach to either the boundary equilibrium or the stable positive equilibrium, depending on the size of the release rate. This implies that the control of the wild mosquito population is highly dependent on the rate at which the sterile mosquito are released, with only high release rates giving sufficient control. However, when the number of the sterile releases is s proportional to the wild mosquito population, analytical results in Theorems 3.1 and 3.2 and numerical results in Section 4 suggest that the solutions of the model exhibit complex dynamical behavior. Some phenomena are observed on fluctuating interactions between the wild and sterile mosquitoes. As the time delay varies, we have found Hopf bifurcations, Bautin bifurcations and Hopf-Hopf bifurcations in the system. The existence and stability of periodic solutions created in these three types of bifurcations are investigated by numerical analysis. These obtained results describe the equilibrium of system process. In particular, when a stable periodic orbit exists, it can be understood that the wild and sterile mosquitoes system can coexist for a long term although the wild mosquito is not eliminated. The conditions for the existence of the bifurcations indicate that the parameters of the system are important in controlling the development and progression of wild mosquitoes.

In summary, by investigating two different delay differential equations models, we have provided convincing evidences that a simple delay that accounts for such processes as the growth stages of mosquito populations is not alone responsible for sustained oscillations between the interactive dynamics of the wild and sterile mosquitoes. If such oscillations occur, then they must be the consequence of delayed dependence on other processes, or of a more complex dependence on the past population density and Allee effect [4,21]. Therefore, our research in this paper shows that the choice of the releasing sterile mosquito strategies is an important determinant of overall mosquito population dynamics.

Finally, we would like to stress that for many mosquito populations, seasonal environmental effects are an important feature of their life-history, giving rise to seasonal cyclic dynamics [36,30]. For example, seasonal rainfall can increase the abundance of mosquitoes, where the reproduction depends on the availability of suitable breeding sites such as water-filled containers. Moreover, many species of mosquitoes in temperate zones overwinter in a diapausal state [33]. Hence, in order to fully understand the influences of the releasing sterile mosquitoes on control wild mosquitoes, further investigation on seasonal environmental effects must be addressed in the modeling, which is planned in our future research.

---

Acknowledgments

The authors are very grateful to the Editor-in-Chief, Professor Yang Kuang, and the anonymous referees for their valuable comments and suggestions, which helped us to improve the presentation of this work significantly.

---