Research article

Stage-structured discrete-time models for interacting wild and sterile mosquitoes with beverton-holt survivability

  • Received: 18 June 2018 Accepted: 03 October 2018 Published: 11 January 2019
  • The sterile insect technique (SIT) is an effective weapon to prevent transmission of mosquito-borne diseases, in which sterile mosquitoes are released to reduce or eradicate the wild mosquito population. To study the impact of the sterile insect technique on the disease transmission, we formulate stage-structured discrete-time models for the interactive dynamics of the wild and sterile mosquitoes using Beverton-Holt type of survivability, based on difference equations. We incorporate different strategies for releasing sterile mosquitoes, and investigate the model dynamics. Numerical simulations are also provided to demonstrate dynamical features of the models.

    Citation: Yang Li, Jia Li. Stage-structured discrete-time models for interacting wild and sterile mosquitoes with beverton-holt survivability[J]. Mathematical Biosciences and Engineering, 2019, 16(2): 572-602. doi: 10.3934/mbe.2019028

    Related Papers:

    [1] Chen Liang, Hai-Feng Huo, Hong Xiang . Modelling mosquito population suppression based on competition system with strong and weak Allee effect. Mathematical Biosciences and Engineering, 2024, 21(4): 5227-5249. doi: 10.3934/mbe.2024231
    [2] Shuyang Xue, Meili Li, Junling Ma, Jia Li . Sex-structured wild and sterile mosquito population models with different release strategies. Mathematical Biosciences and Engineering, 2019, 16(3): 1313-1333. doi: 10.3934/mbe.2019064
    [3] Liming Cai, Shangbing Ai, Guihong Fan . Dynamics of delayed mosquitoes populations models with two different strategies of releasing sterile mosquitoes. Mathematical Biosciences and Engineering, 2018, 15(5): 1181-1202. doi: 10.3934/mbe.2018054
    [4] Martin Bohner, Jaqueline Mesquita, Sabrina Streipert . The Beverton–Hold model on isolated time scales. Mathematical Biosciences and Engineering, 2022, 19(11): 11693-11716. doi: 10.3934/mbe.2022544
    [5] Luis Almeida, Michel Duprez, Yannick Privat, Nicolas Vauchelet . Mosquito population control strategies for fighting against arboviruses. Mathematical Biosciences and Engineering, 2019, 16(6): 6274-6297. doi: 10.3934/mbe.2019313
    [6] Martin Bohner, Sabrina Streipert . Optimal harvesting policy for the Beverton--Holt model. Mathematical Biosciences and Engineering, 2016, 13(4): 673-695. doi: 10.3934/mbe.2016014
    [7] Rajivganthi Chinnathambi, Fathalla A. Rihan . Analysis and control of Aedes Aegypti mosquitoes using sterile-insect techniques with Wolbachia. Mathematical Biosciences and Engineering, 2022, 19(11): 11154-11171. doi: 10.3934/mbe.2022520
    [8] Xiaoli Wang, Junping Shi, Guohong Zhang . Bifurcation analysis of a wild and sterile mosquito model. Mathematical Biosciences and Engineering, 2019, 16(5): 3215-3234. doi: 10.3934/mbe.2019160
    [9] Zhongcai Zhu, Xiaomei Feng, Xue He, Hongpeng Guo . Mirrored dynamics of a wild mosquito population suppression model with Ricker-type survival probability and time delay. Mathematical Biosciences and Engineering, 2024, 21(2): 1884-1898. doi: 10.3934/mbe.2024083
    [10] Yuanxian Hui, Genghong Lin, Qiwen Sun . Oscillation threshold for a mosquito population suppression model with time delay. Mathematical Biosciences and Engineering, 2019, 16(6): 7362-7374. doi: 10.3934/mbe.2019367
  • The sterile insect technique (SIT) is an effective weapon to prevent transmission of mosquito-borne diseases, in which sterile mosquitoes are released to reduce or eradicate the wild mosquito population. To study the impact of the sterile insect technique on the disease transmission, we formulate stage-structured discrete-time models for the interactive dynamics of the wild and sterile mosquitoes using Beverton-Holt type of survivability, based on difference equations. We incorporate different strategies for releasing sterile mosquitoes, and investigate the model dynamics. Numerical simulations are also provided to demonstrate dynamical features of the models.


    To prevent mosquito-borne diseases, the sterile insect technique (SIT) has been applied to reduce or eradicate the wild mosquitoes and has shown promising results in laboratory studies [1,8,38], but predicting the impact of releasing sterile mosquitoes into the field of wild mosquito populations is still challenging. Mathematical models have proven useful in gaining insights into interactive dynamics of wild and sterile mosquito populations, and there are models in the literature for such studies [4,5,6,7,13,17,18,28,29]. However, most of them assume homogeneous mosquito populations without distinguishing the metamorphic stages of mosquitoes.

    Mosquitoes undergo complete metamorphosis, going through four distinct stages of development during a lifetime, egg, pupae, larva, and adult. After a female mosquito drinking blood, she can lay from 100 to 300 eggs at a time in standing water or very slow-moving water. In her lifetime, she can produce from 1000 to 3000 eggs [34]. Within a week, the eggs hatch into larvae, which will use their tubes to breathe air by poking above the surface of the water. Larvae eat a bit of floating organic matter and each other. Larvae molt four times totally as they grow and after the fourth molt, they are called pupae. Pupae also live near the surface of water and breathe through two horn-like tubes (called siphons) on their back. But pupae do not eat. When the skin splits after a few days from a pupae, an adult mosquito emerges. The adults live for only a few weeks and a full life-cycle of a mosquito takes about a month [2,9].

    To have more realistic modeling of mosquitoes, we need to include stage structure since the different stages have different responses to environment and regulating factors to the population [36]. While the interspecific competition and predation are rare events and could be discounted as major causes of larval mortality, the intraspecific competition could represent a major density dependent source. Thus the effect of crowding could be an important factor in the population dynamics of mosquitoes [15,19,35].

    Moreover, since the first three stages in a mosquito's life time are aquatic and the major density dependent source comes from the larval stage, following the line in [24,26], we group the three aquatic stages of mosquitoes into one class and divide the whole mosquito population into only two classes to keep our mathematical modeling as simple as possible. We call the class consisting of the first three stages larvae and the other class adult. We assume that the density dependence is based on larvae not the adults. We still simplify our models such that no male and female individuals are distinguished.

    For the density-dependent mortality, most existing works in the literature, including our previous studies, have assumed the Ricker-type nonlinearity [24,25,27,29,31]. The dynamics of the Ricker-type nonlinearity are complex, causing, e.g. period-doubling bifurcations even without any other interactions. As the sterile mosquitoes are included, the model dynamics become more complex and it is not clear whether the complexity is from the baseline model without interaction already or from the interaction. Thus, we assume that the mosquito population follow the nonlinearity of Beverton-Holt type [11,12] in this paper.

    We first investigate the dynamics of the general stage-structured model with no releases of sterile mosquitoes in Section 2. We then introduce sterile mosquitoes into the model and formulate the interactive stage-structured models in Section 3. Similar to those in [4,5,13,29], we consider three strategies of releases. The case of constant releases is studied in Section 3.1. Complete mathematical analysis for the model dynamics is given. We then formulate a model for the case where the number of sterile mosquito releases is proportional to the wild mosquito population size in Section 3.2. Mathematical analysis and numerical simulations are provided to demonstrate the complexity of the model dynamics. Considering different sizes of wild mosquito population, we consider a different releasing strategy as in [30,31] in Section 3.3, where the releases of sterile mosquitoes are proportional to the wild mosquitoes size when the wild mosquitoes size is small but is saturated and approaches a constant as the wild mosquitoes size is sufficiently large. We provide complete mathematical analysis for the model dynamics. We finally provide a brief discussion on our findings, particularly on the impact of the three different strategies on the mosquito control measures in Section 4.

    We first consider a stage-structured model of wild mosquitoes in the absence of sterile mosquitoes. Let xn and yn be the numbers of mosquito larvae and adults at generation n, respectively, and assume that population dynamics of the mosquitoes are described by the following system:

    xn+1=f(xn,yn)yns1(xn,yn),yn+1=g(xn,yn)xns2(xn,yn),

    where f is the number of the offspring produced per adult, s1 is the survival probability of larvae or the fraction of larvae who survive, g is the progression rate of larvae or the adults emergence rate, and s2 is the survival probability of adults.

    We assume a constant birth rate and denote it as f:=a. Since the intraspecific competition mainly takes place within the aquatic stages of mosquitoes, we assume that the death and the progression rates of larvae are density-dependent only on the larvae size, with the Beverton-Holt type of nonlinearity, such that s1(xn,yn)=k11+η1xn and g(xn,yn)=γ1+η2xn, where k1 is the maximum survival probability, η1 and η2 are density-dependent factors, and γ is the maximum progression rate.

    We assume that food is abundant for mosquito adults so that the adults survival rate is constant, denoted as s2(xn,yn):=s2. Then the model equations become

    xn+1=ayn1+η1xn,yn+1=γxn1+η2xn, (2.1)

    where we merge k1 and s2 into a and γ, respectively, but we still use a and γ for those parameters without confusion.

    The origin (0,0) is a trivial fixed point of system (2.1). Define the intrinsic growth rate of the stage-structured mosquito population r0:=aγ. The trivial fixed point is locally asymptotically stable if r0<1 and is unstable if r0>1.

    The x component of a positive fixed point of system (2.1) satisfies

    η1η2x2+(η1+η2)x+1r0=0. (2.2)

    Clearly, the quadratic equation (2.2) has no positive root and hence system (2.1) has no positive fixed point if r01. Equation (2.2) has a unique positive root and hence system (2.1) has a unique positive fixed point ¯E:=(ˉx,ˉy) if r0>1, with

    ˉx=(η1+η2)+Δ2η1η2,ˉy=γ((η1+η2)+Δ)2η1η2+η2((η1+η2)+Δ), (2.3)

    where Δ=(η1η2)2+4η1η2aγ.

    The Jacobian matrix of system (2.1) at ¯E has the form

    J1:=(η1ˉx1+η1ˉxa1+η1ˉxγ(1+η2ˉx)20). (2.4)

    Fixed point ¯E is locally asymptotically stable if

    |trJ1|<1+detJ1<2

    [22,33], that is,

    η1ˉx1+η1ˉx<111+η2ˉx<2.

    Thus fixed point ¯E is locally asymptotically stable if

    η1<η2,

    and is unstable if

    η1>η2.

    System (2.1) may have periodic cycles of different periods. We first consider 2-cycles with xn+2=xn0 and yn+2=yn0, for all n0.

    It follows from system (2.1) that

    xn+2=ayn+11+η1xn+1=aγxn(1+η2xn)(1+η1xn+1),yn+2=γxn+11+η2xn+1=aγyn(1+η1xn)(1+η2xn+1). (2.5)

    Then there exists a positive nontrivial 2-cycle if and only if

    (1+η2xn)(1+η1xn+1)=(1+η1xn)(1+η2xn+1)=aγ,

    which implies

    (η2η1)(xnxn+1)=0,

    for all n0.

    If η1η2, there exist no positive 2-cycles. If η1=η2:=η, there may exist positive nontrivial 2-cycles. Before we investigate their existence and dynamics, we consider nonnegative synchronous 2-cycles which are not strictly positive, but are non-negative with alternating zero and positive components [14]. In such a situation, the mosquito larvae and adults are synchronized in such a way as to appear and vanish alternately in one time unit.

    Synchronous 2-cycles can be found by looking for nontrivial equilibria of system (2.5) which have one component equal to zero. It follows from

    xn+2=aγxn(1+η2xn)(1+η1xn+1)

    that

    x=aγx(1+η2x)(1+η10),

    which yields

    x=r01η2. (2.6)

    Then it follows from

    xn+2=ayn+11+η1xn+1,

    that

    x=ay1+η10,

    which leads to

    y=xa=aγ1aη2. (2.7)

    Thus the system will undergo a unique synchronous 2-cycle as

    (x0)(0y)(x0)(0y),

    where x is given in (2.6) and y is given in (2.7), for r0>1.

    At the synchronous 2-cycle, the Jacobian matrix is

    J2:=(aγ(1+η2x)2aη1x1+η1x0aγ1+η1x)=(1r0aη1x1+η1x0r01+η1x).

    Then, it follows from

    trJ2=1r0+r01+η1x,detJ2=11+η1x,

    that the unique synchronous 2-cycle is locally asymptotically stable if

    1r0+r01+η1x<1+11+η1x,

    that is,

    η1>η2.

    We now assume η1=η2:=η to analyze the case of positive 2-cycles. Then the unique positive fixed point has the components

    ˉx0=r01η,ˉy0=γ(r01)aη, (2.8)

    whose stability is determined by the eigenvalues of the Jacobian matrix in (2.4). The eigenvalues are

    λ1=1,λ2=1r0<1,

    which implies the possibility of bifurcated 2-cycles. We now explore the existence of such 2-cycles.

    It follows from system (2.1) that

    xn+2=aγxn(1+ηxn)(1+ηxn+1)=r0xn1+ηxn+aηyn,yn+2=aγyn(1+ηxn)(1+ηxn+1)=r0yn1+ηxn+aηyn. (2.9)

    For a positive 2-cycle with initial values (x0,y0), it follows from (2.9) that the two components satisfy the following linear equation:

    1+ηx0+aηy0=r0. (2.10)

    Thus, if we let x0(0,r01η) and y0=r01ηx0aη, then a solution of (2.1) with such initial values is a positive 2-cycle.

    Any positive 2-cycle has the form

    (x1y1)(x2y2)(x1y1)(x2y2)

    with (x1,y1) and (x2,y2) on the straight line given in (2.10). Hence there exists a continuum of positive 2-cycles of system (2.1). To determine the asymptotic behavior of these 2-cycles, we first show that for any initial point (x0>0,y0>0), the distance between the point (xk,yk) and the straight line given in (2.10) after k1 steps is smaller and smaller until converges to zero.

    If the initial point (x0,y0) is above the line given in (2.10), the distance between the line and the point (xk,yk) after k steps is

    dk=ηxk+aηyk+1r0η2+(aη)2.

    Thus

    dk+1dk=η(xk+1xk)+aη(yk+1yk)η2+(aη)2,

    where

    xk+1xk=aykxk(1+ηxk)1+ηxk,yk+1yk=γxkyk(1+ηxk)1+ηxk.

    Therefore,

    dk+1dk=ηxk(1+ηxk)η2+(aη)2(r0(1+ηxk+aηyk)), (2.11)

    where r0<1+ηxk+aηyk since the point (xk,yk) is above the straight line 1+ηx+aηy=r0 for any positive k.

    We use mathematical induction method to prove that (xk,yk) is above the straight line 1+ηx+aηy=r0 given that (x0,y0) is above the same line. If the point (xm,ym) is above the line, then 1+ηxm+aηym>r0. For the point (xm+1,ym+1), we have

    1+ηxm+1+aηym+1r0=1+ηaym1+ηxm+aηγxm1+ηxmr0=1+ηxm+aηymr01+ηxm>0,

    which implies that the point (xm+1,ym+1) is also above the same line. Thus, for any positive k, the point (xk,yk) is always above the straight line 1+ηx+aηy=r0.

    Then we have

    dk+1dk<0,

    which indicates that the distance {dk} is a nonnegative strictly decreasing sequence bounded below by zero. Thus, limkdk:=d exists. Taking the limit in (2.11), we then have limk(1+ηxk+aηyk)=r0, that is, (xk,yk) approaches the line given in (2.10) and thus d=0.

    Similarly, if the initial point is below the line, we can show that the distance {dk} is also a nonnegative strictly decreasing sequence with a limit equal to zero and (xk,yk) approaches the line given in (2.10) as well. Therefore, the line given in (2.10) is a global attractor and a continuum of positive 2-cycles of system (2.1) occurs, where each of the positive 2-cycle is locally stable. We give an example to demonstrate the existence of a positive 2-cycle in Example 1.

    Example 1. Given parameters

    a=5,γ=0.4,η2=0.3, (2.12)

    there exists a continuum of positive 2-cycles of system (2.1) if η1=0.3. The line given in (2.10) is a global attractor. With initial values on the line, the solutions of (2.1) are positive 2-cycles, which are locally stable. Initial value (x0,y0)=(0.9524,0.4762) is on the line given in (2.10) and creates a positive 2-cycle as shown in the left figure in Figure 1. If η1=0.5, there exists a unique synchronous 2-cycle with components x=3.3333 and y=0.6667 which is globally asymptotically stable as shown in the right figure in Figure 1.

    Figure 1.  The parameters are given in (2.12). The line given in (2.10) is globally attractive if η1=0.3. Initial value (x0,y0)=(0.9524,0.4762) is on the line and creates a positive 2-cycle as shown in the left figure. If η1=0.5, there exists a unique synchronous 2-cycle with components x=3.3333 and y=0.6667 which is globally asymptotically stable as shown in the right figure.

    We would like to point out that if we define parameter Γ:=η1η2, it can be used as a bifurcation parameter. As Γ<1, that is η1<η2, there exists a unique positive fixed point which is globally asymptotically stable. It loses its stability when Γ1, that is η1η2. At Γ=1, a continuum of positive 2-cycles is bifurcated and for Γ>1, a unique asymptotically stable synchronous 2-cycle appears. A simple bifurcation diagram is given in Figure 2.

    Figure 2.  When η1<η2, that is Γ<1, there exists a unique positive fixed point which is globally asymptotically stable. At Γ=1, a continuum of positive 2-cycles is bifurcated which is shown as the vertical line. When Γ>1, that is, η1>η2, the positive fixed point becomes unstable, and a unique asymptotically stable synchronous 2-cycle appears.

    We further prove, in Appendix 4, that other kcycles, for k3, do not exist, which indicates the chaotic phenomenon cannot occur [16]. We summarize our results as follows.

    Theorem 2.1. The trivial fixed point (0,0) for (2.1) is globally asymptotically stable if the intrinsic growth rate r0<1 and unstable if r0>1. System (2.1) has a unique positive fixed point ¯E=(ˉx,ˉy) given in (2.3) and a unique synchronous 2-cycle if r0>1. The positive fixed point is globally asymptotically stable when η1<η2 and is unstable when η1>η2. The unique synchronous 2-cycle is globally asymptotically stable if η1>η2 and unstable if η1<η2. A continuum of positive 2-cycles with initial values on the straight line given in (2.10) appears when η1=η2. The straight line is globally attractive and every positive 2-cycle with an initial value on the straight line in (2.10) is locally stable.

    Now suppose sterile mosquitoes are released into the field of wild mosquitoes. Since sterile mosquitoes do not reproduce, the birth input term will be their releases rate. Let Bn be the number of sterile mosquitoes released at generation n. After the sterile mosquitoes are released, the mating interaction between wild and sterile mosquitoes takes place. We assume harmonic means for matings such that the per capita birth rate is given by

    bn=C(Nn)aynyn+Bn=C(Nn)aynNn,

    where C(Nn) is the number of matings per mosquito, per unit of time, with Nn=yn+Bn and a is the number of wild larvae produced per wild mosquito. The interactive dynamics of wild and sterile mosquitoes are then described by the following system:

    xn+1=C(Nn)aynyn+Bnyn1+η1xn,yn+1=γxn1+η2xn. (3.1)

    We first consider the case where Bn:=b is a constant which means sterile mosquitoes are constantly released for each generation, and assume that the number of matings C(Nn) is a constant and is merged into the birth rate a with the same notation for convenience. Then the system (3.1) becomes

    xn+1=aynyn+byn1+η1xn=ay2n(yn+b)(1+η1xn),yn+1=γxn1+η2xn. (3.2)

    Clearly, the origin (0,0) is a fixed point and is always locally asymptotically stable. Let (x,y) be a positive fixed point. It satisfies the following equations:

    x=ay2(y+b)(1+η1x),y=γx1+η2x,

    which lead to

    ay(y+b)(1+η1x)γ1+η2x=1.

    Solving for b then yields

    b=aγy(1+η1x)(1+η2x)y=γx(1+η1x)(1+η2x)2(r0(1+η1x)(1+η2x)):=γH(x), (3.3)

    for (1+η1x)(1+η2x)r0, i.e. xˉx, where ˉx is given in (2.3).

    Clearly, there exists no positive fixed point if r01. For the existence of positive fixed point, function H(x) first needs to satisfy H(x)>0 for r0>1. We then only consider xΩ where set Ω is defined by

    Ω:={x:0<x<ˉx}. (3.4)

    Since

    H(x)=1(1+η1x)(1+η2x)3(r0(1η2xη1x(1+η2x)1+η1x)(1+η1x)(1+η2x)), (3.5)

    we define L(x):=r0(1η2xη1x(1+η2x)1+η1x) and F(x):=(1+η1x)(1+η2x). Then H(x)=0 for x0 if and only if L(x)=F(x) for x0.

    Since

    L(x)=r0(η2η1+2η1η2x+η21η2x2(1+η1x)2)<0,

    L(x) is decreasing in the set Ω and L(0)=r0>1. Notice that F(x) is increasing in Ω and F(0)=1. Then there exists a unique intersection point between the curves of F(x) and L(x), denoted by xc, with 0<xcˉx since F(xc)=L(xc)<r0.

    With this unique xc, we have H(x)|x=xc=0. Clearly, L(x)>F(x) for 0<x<xc, and L(x)<F(x) for xc<x<ˉx. Thus

    {H(x)>00<x<xc,H(x)<0xc<x<ˉx.

    We define

    bc:=γH(xc). (3.6)

    Then bc determines the threshold value for releases of sterile mosquitoes such that system (3.2) has no positive fixed point, one positive fixed point (x,y), or two positive fixed points (xi,yi), i=1,2 with x1<xc<x2, if b>bc, b=bc, or b<bc, respectively.

    To investigate the existence of synchronous 2-cycles, we have such a form

    (x0)(0y)(x0)(0y),

    where

    x=ay2b+y,y=γx1+η2x. (3.7)

    To solve for x and y, we have

    aη2y2+(1r0)y+b=0. (3.8)

    Since r0>1, the existence of positive solutions of (3.8) depends on the discriminant Δ=(r01)24aη2b. Define the threshold value

    b1:=(r01)24aη2. (3.9)

    Then there exists no, one synchronous 2-cycle with

    y=r012aη2,

    or two synchronous 2-cycles with

    y{1}=r01(r01)24aη2b2aη2,y{2}=r01+(r01)24aη2b2aη2, (3.10)

    if b>b1, b=b1, or b<b1, respectively.

    Now we claim bc<b1. Thus, if b>b1, there exists neither synchronous 2-cycle nor positive fixed point thus no positive 2-cycles, which makes the trivial fixed point globally asymptotically stable. To this end, we define a function

    P(x):=γH(x)b1=Q(x)4aη2(1+η1x)(1+η2x)2,

    where

    Q(x):=Ax3Bx2+Cx(r01)2

    with

    A=4r0η1η22+η1η22(r01)2>0,B=4r0η2(η1+η2)+η22(r01)2+2η1η2(r01)2>0,C=4r0η2(r01)2η2(r01)2η1(r01)2. (3.11)

    Then the function P(x) has the same sign as Q(x). If C<0, then the function Q(x) is negative for all xΩ since all of the coefficients are negative. Thus the function P(x) is negative for all xΩ, which implies bc<b1.

    Assume C>0. It follows from

    Q(x)=3Ax22Bx+C

    that Q(x) has a maximum value

    Q(x0)=Ax30Bx20+Cx0(r01)2,

    at

    x0=B2+3ACB3A>0,

    where, A, B and C are given in (3.11).

    Using η1 as a variable, we have

    Q(x0):=q(η1)=A(η1)x30(η1)B(η1)x20(η1)+C(η1)x0(η1)(r01)2,

    where η1(0,1].

    Taking the derivative of q(η1) with respect to η1, we have

    q(η1)=(3A(η1)x20(η1)2B(η1)x0(η1)+C(η1))x0(η1)+(A(η1)x30(η1)B(η1)x20(η1)+C(η1)x0(η1))=A(η1)x30(η1)B(η1)x20(η1)+C(η1)x0(η1)

    since

    3Ax202Bx0+C=Q(x0)=0.

    Notice that

    A(η1)=4r0η22+η22(r01)2>0,B(η1)=4r0η2+2η2(r01)2>0,C(η1)=(r01)2<0.

    Thus q(η1)<0 for η1(0,1], and then q(η1) is monotone decreasing.

    Moreover, it follows from A(0)=0, B(0)=η22(r0+1)2, and C(0)=2η2(r201) that

    x0(η1)|η10=limη10x0(η1)=C(0)2B(0),

    and thus

    limη10q(η1)=B(0)x20(0)+C(0)x0(0)(r01)2=B(0)C2(0)4B2(0)+C(0)C(0)2B(0)(r01)2=C2(0)4B(0)(r01)2=14B(0)(C2(0)4B(0)(r01)2)=0.

    Hence q(η1)<0 for η1(0,1], that is, Q(x0)<0, and then P(x)<0 for all xΩ. Therefore bc<b1.

    We next investigate the stability of the positive fixed points and the synchronous 2-cycles.

    The Jacobian matrix evaluated at a positive fixed point has the form

    J:=(η1x1+η1xxyy+2by+byx11+η2x0).

    Since

    trJ=η1x1+η1x,detJ=y+2b(y+b)(1+η2x),

    a positive fixed point (x,y) is locally asymptotically stable if

    η1x1+η1x+y+2b(y+b)(1+η2x)<1,

    which is equivalent to

    b(1+2η1xη2x)(η2η1)xy<0. (3.12)

    If η2<η1, then b(1+2η1xη2x)(η2η1)xy=b(1+(2η1η2)x+(η1η2)xy>0, which implies that all positive points are unstable if they exist.

    We then assume η2>η1 such that the wild mosquitoes maintain a locally steady state before the sterile mosquitoes are released. Since the component x of a positive fixed point (x,y) is a solution of b=γH(x) in (3.3), condition (3.12) is equivalent to

    γH(x)(1+2η1xη2x)(η2η1)xγx1+η2x=γx(1+η1x)(1+η2x)(r0(1+2η1xη2x)(1+η1x)2(1+η2x))<0.

    Define the following function h(x) to determine the stability of positive fixed points

    h(x):=r0(1+2η1xη2x)(1+η1x)2(1+η2x).

    Thus, the positive fixed point (x,y) is locally asymptotically stable if h(x)<0 and unstable if h(x)>0.

    For b=bc, there exists a unique positive fixed point (xc,yc) where H(xc)=0 from (3.5) and hence

    (1+η1xc)(1+η2xc)=r0(1η2xcη1xc(1+η2xc)1+η1xc).

    Then we have

    h(xc)=r0(1+2η1xcη2xc)(1+η1xc)r0(1η2xcη1xc(1+η2xc)1+η1xc)=r0(1+2η1xcη2xc)r0((1+η1xc)η2xc(1+η1xc)η1xc(1+η2xc))=r0(2η1xc+2η1η2x2c)=2r0η1xc(1+η2xc)>0,

    and thus this unique positive fixed point (xc,yc) is unstable.

    For b<bc, there exist two positive fixed points (xi,yi), i=1,2, with x1<x2, where H(x1)>0 and H(x2)<0.

    We first consider the positive fixed point (x1,y1), with H(x1)>0 which is equivalent to

    (1+η1x1)(1+η2x1)<r0(1η2x1η1x1(1+η2x1)1+η1x1).

    Then we have

    h(x1)>r0(1+2η1x1η2x1)(1+η1x1)r0(1η2x1η1x1(1+η2x1)1+η1x1)=r0(1+2η1x1η2x1)r0((1+η1x1)η2x1(1+η1x1)η1x1(1+η2x1))=r0(2η1x1+2η1η2x21)=2r0η1x1(1+η2x1)>0,

    and thus fixed point (x1,y1) is unstable.

    Next we consider the positive fixed point (x2,y2) with H(x2)<0. Simple calculation shows that h(x2) can be negative or positive. We then define a threshold value xs satisfying h(xs)=0, that is,

    h(xs)=η21η2x3s(η21+2η1η2)x2s+(r0(2η1η2)2η1η2)xs+(r01)=0.

    Notice that h(x)=0 is a cubic equation and r01>0. Then it follows from Descartes' rule of sign that there exists a unique positive solution to h(x)=0. Moreover, since

    h(ˉx)=r0(η1η2)ˉx<0

    and h(x1)>0, the threshold value xs satisfies x1<xs<ˉx, and h(x)>0 for x1<x<xs and h(x)<0 for xs<x<ˉx.

    Then, the positive fixed point (x2,y2) is locally asymptotically stable if x2>xs with h(x2)<0 and unstable if x2<xs with h(x2)>0. We define the corresponding threshold value of releases for the stability of (x2,y2) as

    bs:=γH(xs), (3.13)

    where H(x) is defined in (3.3). Then positive fixed point (x2,y2) is locally asymptotically stable if b<bs and unstable if b>bs. Notice that bs<bc since bc is the maximum value of the function γH(x).

    For the stability of the synchronous 2-cycles with components x and y given in (3.7), the corresponding Jacobian matrix has the form

    (λ000),

    where

    λ=2r0y(1+η2x)(b+y)2aη2y2(1+η2x)(b+y)r0y2(1+η2x)2(b+y)2.

    Since r0y=(1+η2x)(b+y) at the synchronous 2-cycles,

    λ=2r0yr0y2aη2y2r0yr0y2r20y2=22η2γy1r0.

    If b<b1, there exist two synchronous 2-cycles with components x{1}<x{2} and y{1}<y{2} in (3.10). The corresponding eigenvalues are

    λ{1}=2r0(r01)24baη2r0=r0+(r01)24baη2r0>1,

    at (x{1},y{1}), and

    λ{2}=2r0+(r01)24baη2r0=r0(r01)24baη2r0<1,

    at (x{2},y{2}), respectively. Thus synchronous 2-cycle (x{1},y{1}) is unstable and synchronous 2-cycle (x{2},y{2}) is locally asymptotically stable.

    We summarize our results in the following theorem and Table 1.

    Table 1.  Summary table for the existence of positive fixed points and synchronous 2-cycles with r0>1.
    (PFP stands for positive fixed point and STC stands for synchronous 2-cycle.)
    b<bs bs<b<bc bc<b<b1 b1<b
    η1<η2 Two PFP Two PFP No PFP
    One stable both unstable
    One unstable
    Two STC No STC
    One stable
    One unstable
    η2<η1 Two PFP No PFP
    both unstable
    Two STC No STC
    One stable
    One unstable

     | Show Table
    DownLoad: CSV

    Theorem 3.1. The trivial fixed point (0,0) is always locally asymptotically stable for system (3.2). Suppose the intrinsic growth rate of the stage-structured mosquito population r0>1. We define the threshold values for the releases of the sterile mosquitoes bc in (3.6), b1 in (3.9), and bs in (3.13), respectively, with bs<bc<b1. System (3.2) has no, one, or two positive fixed points if b>bc, b=bc, or b<bc, and has no, one, or two synchronous 2-cycles if b>b1, b=b1, or b<b1 respectively. If b>b1, the trivial fixed point is globally asymptotically stable, which makes the population of the wild mosquitoes to go extinct when sufficient sterile mosquitoes are released. If bc<b<b1, there exist no positive fixed points and two synchronous 2-cycles, where the one with larger components is locally asymptotically stable and the other is unstable. If b<bc, it depends on the relation of η1 and η2. If η2<η1, there exist two positive fixed points, where both of them are unstable, and two synchronous 2-cycles, where the one with larger components is locally asymptotically stable and the other is unstable. Suppose η1<η2. If bs<b<bc, there exist two positive fixed points, but both are unstable, and two synchronous 2-cycles, where the one with larger components is locally asymptotically stable and the other is unstable. If 0<b<bs, there exist two positive fixed points and two synchronous 2-cycles. The positive fixed point with larger components and the synchronous 2-cycle with larger components are both locally asymptotically stable, and the other positive fixed point and the other synchronous 2-cycle are both unstable.

    We then give the following example to demonstrate the results in Theorem 3.1, but only address the case of η1<η2.

    Example 2. Choosing the following parameters

    a=25,γ=0.8,η1=0.2,η2=0.7, (3.14)

    we have r0=aγ=20>1 and η1<η2 . The threshold values of releases are bc=4.1523, b1=5.1571 and bs=2.9252. For b=1, there exist two positive fixed points E1=(0.0741,0.0564) and E2=(5.1949,0.8964), where E1 is unstable and E2 is locally asymptotically stable since b<bs. For the same release value b=1, there also exist two synchronous 2-cycles with components (x{1},0)(0,y{1}) and (x{2},0)(0,y{2}), where the synchronous 2-cycle with bigger components x{2}=13.0700 and y{2}=1.0302 is locally asymptotically stable while the one with smaller components x{1}=0.0728 and y{1}=0.0554 is unstable. Notice that the origin (0,0) is always locally asymptotically stable. Therefore, for b=1, solutions approach the origin, positive fixed point E2, or synchronous 2-cycle with components x{2} and y{2}, depending on their initial values, as shown in the upper left, upper right, or the lower figure in Figure 3, respectively.

    Figure 3.  With the parameters given in (3.14), the threshold values of releases are bc=4.1523, b1=5.1571, and bs=2.9252. For b=1, there exist two positive fixed points E1=(0.0741,0.0564) and E2=(5.1949,0.8964), where E1 is unstable and E2 is locally asymptotically stable since b<bs. For the same release value b=1, there also exist two synchronous 2-cycles with components (x{1},0)(0,y{1}) and (x{2},0)(0,y{2}), where the synchronous 2-cycle with bigger components x{2}=13.0700 and y{2}=1.0302 is locally asymptotically stable while the one with smaller components x{1}=0.0728 and y{1}=0.0554 is unstable. Notice that the origin (0,0) is always locally asymptotically stable. Solutions approach the origin, positive fixed point E2, or synchronous 2-cycle (x{2},0)(0,y{2}), depending on their initial values, as shown in upper left, upper right, or the lower figure, respectively.

    Instead of constant releases of sterile mosquitoes, we assume that the releases are proportional to the population size of the wild mosquitoes such that the number of releases is B():=by where b is a constant [13].

    We assume that there is no mating difficulty for mosquitoes. Then the model dynamics are governed by the following system:

    xn+1=aynyn+bynyn1+η1xn=ayn(1+b)(1+η1xn)=ˉayn1+η1xn,yn+1=γxn1+η2xn, (3.15)

    where ˉa=a1+b.

    Mathematically, the system is the same as system (2.1). It follows from Theorem 2.1 that if ¯r0=ˉaγ>1, the trivial fixed point is globally asymptotically stable. We define the sterile mosquitoes release threshold value by

    bc:=aγ1=r01.

    Then the trivial fixed point is globally asymptotically stable if b>bc and unstable if b<bc.

    If b<bc, there exists a unique positive fixed point E:=(x,y) with

    x=(η1+η2)+Δ2η1η2,y=γx1+η2x=γ((η1+η2)+Δ)2η1η2+η2((η1+η2)+Δ), (3.16)

    where Δ=(η1η2)2+4η1η2r01+b, and a unique synchronous 2-cycle

    (x0)(0y)(x0)(0y) (3.17)

    with x=r0(1+b)η2(1+b)>0 and y=r0(1+b)r0η2>0. It follows from Theorem 2.1 that the positive fixed point E is globally asymptotically stable if η1<η2 and the synchronous 2-cycle is globally asymptotically stable if η1>η2. In summary, we have the following theorem.

    Theorem 3.2. The trivial fixed point (0,0) for system (3.15) is globally asymptotically stable if b>bc and is unstable if b<bc where the sterile mosquito release threshold bc:=r01. If b<bc, there exist a unique positive fixed point and a unique synchronous 2-cycle. The unique positive fixed point E, given in (3.16), is globally asymptotically stable if η1<η2 and unstable if η1>η2. The synchronous 2-cycle, given in (3.17), is globally asymptotically stable if η1>η2 and unstable if η1<η2. A continuum of locally stable positive 2-cycles exists if and only if η1=η2.

    It follows from Theorem 3.2 that if b>bc such that sufficient sterile mosquitoes are released, wild mosquitoes can be wiped out. On the other hand, if b<bc such that not enough sterile mosquitoes are released, the two types of mosquitoes coexist. In such a case, although the stability condition remains the same as in system (2.1) even we have released sterile mosquitoes, it follows from (3.16) that x becomes smaller for larger b. Moreover, since y is a strictly increasing function of x, y becomes smaller with smaller x. That is to say, increasing the releases of sterile mosquitoes can reduce the population size of the wild mosquitoes. We demonstrate our findings in Example 3.

    Example 3. Given the parameters

    a=5,γ=0.8, (3.18)

    we have threshold value bc=3. For b=4>3, the trivial fixed point is globally asymptotically stable as shown in the upper left figure in Figure 4. For b=0.5<bc=3, the trivial fixed point is unstable and there exist a unique positive fixed point E=(2.5519,1.1563) and a unique synchronous 2-cycle with components x=5.5556 and y=1.1111. The synchronous 2-cycle is unstable and the positive fixed point E is globally asymptotically stable when η1=0.2<η2=0.3 as shown in the upper right figure in Figure 4. With the same values η1 and η2, as the threshold value of releases of sterile mosquitoes is increased to b=1.5, the unique positive fixed point E=(1.0642,0.6453) is still globally asymptotically stable but has smaller magnitudes of x and y compared to those for b=0.5 as shown in the lower left figure in Figure 4. If we change the parameter values to η1=0.5>η2=0.3, the unique positive fixed point E=(0.6667,0.4444) is unstable and the unique synchronous 2-cycle with x=2 and y=1.25 becomes globally asymptotically stable as shown in the lower right figure in Figure 4.

    Figure 4.  The parameters are given in (3.18). We have threshold value bc=3. For b=4, the origin is globally asymptotically stable as shown in the upper left figure. For b=0.5, we choose η1=0.2<η2=0.3, the unique positive fixed point E=(2.5519,1.1563) is globally asymptotically stable as shown in the upper right figure. For b=1.5, the unique positive fixed point E=(1.0642,0.6453) is still globally asymptotically stable but has smaller magnitudes of x and y compared to those for b=0.5 as shown in the lower left figure. When η1 and η2 are changed to η1=0.5>η2=0.3, the unique positive fixed point E=(0.6667,0.4444) is unstable and the unique synchronous 2-cycle with x=2 and y=1.25 is globally asymptotically stable as shown in the lower right figure.

    Compared to the case of constant releases, the proportional releases may be a good strategy when the size of the wild mosquito population is small since the size of releases is also small. However, if the wild mosquito population size is significantly large, the release size would be large as well, which may exceed our affordability. Then, as in [13], we consider a different strategy where the number of releases is proportional to the wild adult mosquito population size when it is small, but it is saturated and approaches a constant when the wild adult mosquito population size is sufficiently large. To this end, we let the releases be of Holling-Ⅱ type [21] such that B():=by1+y. Then we consider the following system of equations:

    xn+1=aynyn+byn1+ynyn1+η1xn=ayn(1+yn)(1+b+yn)(1+η1xn),yn+1=γxn1+η2xn. (3.19)

    We assume r0=aγ>1 and define an initial sterile mosquitoes release threshold value b0:=r01 such that the origin (0,0) is locally asymptotically stable if b>b0 and unstable if b<b0.

    A positive fixed point E=(x,y) satisfies

    a(1+y)(1+b+y)(1+η1x)γ1+η2x=1,

    that is,

    b=(1+γx1+η2x)(r0(1+η1x)(1+η2x)1)=:G(x). (3.20)

    Let xΩ where set Ω is defined in (3.4). It follows from

    G(x)=A1x3+A2x2+A3x+A4(1+η1x)2(1+η2x)3, (3.21)

    where

    A1=η21η2γ<0,A2=2η1η22r02η1η2γr0η21γ2η1η2γ<0,A3=3η1η2r0η22r0η2γr02η1γη2γ<0,A4=γ((aγaη2)(1+aη1)),

    that if A4<0, that is, aγaη2<1+aη1, then G(x)<0 for all xΩ. Thus G(x) is decreasing with the maximum value G(0)=b0=r01, and hence there exists no, or one positive fixed point if b>b0, or b<b0.

    Next we assume A4>0, that is, aγaη2>1+aη1. Then G(0)=A4>0 and function G(x) is increasing for x positive and near 0. It is clear that the numerator of G(x) has only one positive root, denoted by xc, according to Descartes’ rule of sign [3,20]. Since the denominator of G(x) is positive, G(x)=0 has a unique positive solution at xc, and hence G(x) has a unique maximum value at x=xc. We write bc:=G(xc). Then function G(x) is increasing for x(0,xc), starting with the point (0,b0), and decreasing for x(xc,ˉx), ending with the point (ˉx,0) where G(ˉx)=0 as above. Write x1=G1(b0)>0. Then for each b(0,b0), there exists a unique x(x1,ˉx) such that G(x)=b. On the other hand, for each b(b0,bc), there exist two x(0,x1) such that G(x)=b. Therefore, with these two threshold values b0<bc, there exists no positive fixed point if b>bc, one positive fixed point if b=bc or 0bb0, or two positive fixed points if b0<b<bc.

    The existence of the positive fixed points, based on the release value of the sterile mosquitoes b, is illustrated in Figure 5 where the xaxis is b.

    Figure 5.  In the left figure, we have a=2.25, γ=.8, η1=.1, and η2=.3 such that A4=0.08<0. Thus there exists no positive fixed point, or one positive fixed point if b>b0, or b<b0, respectively. In the right figure, we have a=2.25, γ=.8, η1=.01, and η2=.1 such that A4=0.44>0. Then there exists no positive fixed point if b>bc, one positive fixed point if b=bc or 0bb0, or two positive fixed points if b0<b<bc.

    The system may also have positive cycles of different periods. We only consider synchronous 2-cycles with components (x,0) and (0,y). It follows from

    (x0)(0y)(x0)(0y)

    that

    x=ay(1+y)1+b+y,y=γx1+η2x,

    where

    aη2y2(b0aη2)y+bb0=0. (3.22)

    If b0aη2<0, that is, aγaη2<1, there exists no or one synchronous 2-cycle if bb0 or b<b0. If b0aη2>0, that is, aγaη2>1, the existence of positive solutions depends on the discriminant of the quadratic equation (3.22). We define the threshold value

    b1:=(b0+aη2)24aη2b0.

    Then there exists no synchronous 2-cycle if b>b1, one synchronous 2-cycle if b=b1 or 0bb0, or two synchronous 2-cycles if b0<b<b1, respectively. Similarly as in Section 3.1, the two threshold values bc and b1 have bc<b1.

    We illustrate our results in Table 2.

    Table 2.  Summary table for the existence of positive fixed points and synchronous 2-cycles.
    (PFP stands for positive fixed point and STC stands for synchronous 2-cycle.)
    b<b0 b0<b<bc bc<b<b1 b1<b
    a(γη2)<1 One PFP No PFP
    One STC No STC
    1<a(γη2)<1+aη1 One PFP No PFP
    One STC Two STC No STC
    1+aη1<a(γη2) One PFP Two PFP No PFP
    One STC Two STC No STC

     | Show Table
    DownLoad: CSV

    We next investigate the stability of the positive fixed points. The Jacobian matrix at a positive fixed point E=(x,y) has the form

    ˉJ:=(η1x1+η1x(b(1+2y)+(y+1)2)(1+η2x)γ(1+y)(1+y+b)γ(1+η2x)20).

    Since

    trˉJ=η1x1+η1x,detˉJ=b(1+2y)+(y+1)2(1+η2x)(1+y+b)(1+y),

    E is locally asymptotically stable if

    b(1+2y)+(y+1)2(1+η2x)(1+y+b)(1+y)<1η1x1+η1x,

    that is,

    [(1+η1x)(1+2y)(1+η2x)(1+y)]b<(η2η1)x(1+y)2. (3.23)

    Substituting y=γx1+η2x for y in (3.23), we then define function

    Φ(x):=[(1+η1x)(1+2y)(1+η2x)(1+y)]b(η2η1)x(1+y)2=γx(η2x+γx+1)(1+η1x)(1+η2x)3ϕ(x),

    where

    ϕ(x)=η21η2x3(η21+2η1η2)x2+(η1(2b0+aη2)η2(r0+1+aη2))x+b0aη2+aη1.

    Thus the positive fixed point is locally asymptotically stable if ϕ(x)<0 and unstable if ϕ(x)>0.

    It is clear from (3.23) that if η2<η1, a positive fixed point is unstable. We then assume η1<η2. If b0aη2+aη1<0, that is, aγaη2<1aη1, the coefficient of the linear term in ϕ(x) becomes

    η1(2b0+aη2)η2(r0+1+aη2)<η1(2b0+aη2)η2(r0+1+b0+aη1)=2b0η12(b0+1)η2=2b0(η1η2)2η2<0,

    which implies that ϕ(x)<0, for all xΩ, and there exists a unique positive fixed point for b<b0. Hence this unique positive fixed point is always locally asymptotically stable.

    If b0aη2+aη1>0, that is, aγaη2>1aη1, there exists a unique positive solution to ϕ(x)=0 from Descartes' rule of sign. Denote it as xs such that ϕ(xs)=0, and define

    bs:=G(xs),

    where G(x) is defined in (3.20). Then a positive fixed point (x,y) is locally asymptotically stable if x>xs and unstable if x<xs.

    For the case of 1aη1<aγaη2<1+aη1, it follows from table 2 that there exists a unique positive fixed point if b<b0, where b0 is the maximum value of function G(x) and bs<b0. Then this positive fixed point is locally asymptotically stable if b<bs and unstable if b>bs.

    For the case of aγaη2>1+aη1 with bs<b0, there exists one or two positive fixed points if b<b0 or b0<b<bc, respectively. Then the unique positive fixed point is locally asymptotically stable if b<bs and unstable if bs<b<b0. The two positive fixed points are unstable if b0<b<bc.

    For the case of aγaη2>1+aη1 with bs>b0. If b<b0, the unique positive fixed point is always locally asymptotically stable. If b0<b<bs, there exist two positive fixed points (xi,yi), i=1,2, with x1<xs<x2. Thus the positive fixed point with larger components, (x2,y2), is locally asymptotically stable while the one with smaller components, (x1,y1), is unstable. If bs<b<bc, we claim xc<xs. Then the two positive fixed points (xi,yi), i=1,2, are both unstable since x1<xc<x2<xs.

    To prove xc<xs we denote the numerator of G(x) in (3.21) by

    dG(x):=A1x3+A2x2+A3x+A4,

    where dG(xc)=0 since G(xc)=0. It follows from

    γϕ(x)dG(x)=(2η1η22r0+2η1η2γr0)x2+(4η1η2r0+2η1γr0)x+2r0η1>0,

    for all xΩ, that γϕ(xc)>dG(xc)=0, which implies ϕ(xc)>0. Thus xc<xs since ϕ(x) is monotone decreasing for positive x and ϕ(xs)=0.

    We use Table 3 to summarize our results, and give Example 4 to demonstrate the existence and stability results for the positive fixed point of system (3.19).

    Table 3.  Summary table for the stability of positive fixed points.
    (PFP stands for positive fixed point and L.A.S stands for locally asymptotically stable.)
    a(γη2)<1aη1 b<b0 b0<b
    One PFP No PFP
    L.A.S. -
    1aη1<a(γη2)<1+aη1 b<bs bs<b<b0 b0<b
    One PFP One PFP No PFP
    L.A.S. unstable -
    1+aη1<a(γη2) with bs<b0 b<bs bs<b<b0 b0<b<bc bc<b
    One PFP One PFP Two PFP No PFP
    L.A.S. unstable both unstable -
    1+aη1<a(γη2) with b0<bs b<b0 b0<b<bs bs<b<bc bc<b
    One PFP Two PFP Two PFP No PFP
    L.A.S. larger one L.A.S. both unstable -

     | Show Table
    DownLoad: CSV

    Example 4. With the following parameters,

    a=2.25,γ=0.8,η1=0.1,η2=0.3, (3.24)

    we have b0=0.8, bs=0.6957 and 1aη1=0.775<a(γη2)=1.125<1+aη1=1.225. There exists no positive fixed point and no synchronous 2-cycles if b>b0, which implies that the trivial fixed point is globally asymptotically stable. For b=1>b0, the trivial fixed point is globally asymptotically stable as shown in the upper left figure in Figure 6. For b=0.79<b0, there exists a unique positive fixed point (0.0822,0.0641), which is unstable since b>bs and a unique synchronous 2-cycle with components x=0.3380 and y=0.2455, which is locally asymptotically stable, as shown in the upper right figure in Figure 6. For b=0.5<bs, there exists a unique positive fixed point (0.8428,0.5382), which is locally asymptotically stable as shown in the lower figure in Figure 6.

    Figure 6.  With the parameters given in (3.24), we have b0=0.8, bs=0.6957 and 1aη1=0.775<a(γη2)=1.125<1+aη1=1.225. For b=1>b0, the trivial fixed point is globally asymptotically stable as shown in the upper left figure. For b=0.79<b0, there exists a unique positive fixed point (0.0822,0.0641), which is unstable since b>bs and a unique synchronous 2-cycle with components x=0.3380 and y=0.2455, which is locally asymptotically stable as shown in the upper right figure. For b=0.5<bs, there exists a unique positive fixed point (0.8428,0.5382), which is locally asymptotically stable as shown in the lower figure.

    Example 5. Given the following parameters,

    a=2.25,γ=0.8,η1=0.1,η2=0.2, (3.25)

    we have the threshold values b1=0.8681, bs=0.5505<b0=0.8000, bc=0.8061 and a(γη2)=1.350>1+aη1=1.225. For b=0.3<bs, there exists a unique positive fixed point (1.6967,1.0135), which is locally asymptotically stable as shown in the upper left figure in Figure 7. For b=0.7<b0, there exists a unique positive fixed point (0.7306,0.5100), which is unstable since b>bs and a unique synchronous 2-cycle with x=1.6667 and y=1.0000, which is locally asymptotically stable as shown in the upper right figure in Figure 7. For b0<b=0.803<bc, there exist two positive fixed points (0.0353,0.0281) and (0.2176,0.1668), which are both unstable, and two synchronous 2-cycles, where the one with larger components x=1.1902 and y=0.7691 is locally asymptotically stable as shown in the lower figure in Figure 7.

    Figure 7.  With the parameters given in (3.25), we have b1=0.8681, bs=0.5505<b0=0.8000, bc=0.8061 and a(γη2)=1.350>1+aη1=1.225. For b=0.3<bs, there exists a unique positive fixed point (1.6967,1.0135), which is locally asymptotically stable as shown in the upper left figure. For b=0.7<b0, there exists a unique positive fixed point (0.7306,0.5100), which is unstable since b>bs and a unique synchronous 2-cycle with x=1.6667 and y=1.0000, which is locally asymptotically stable as shown in the upper right figure. For b0<b=0.803<bc, there exist two positive fixed points (0.0353,0.0281) and (0.2176,0.1668), which are both unstable, and two synchronous 2-cycles, where the one with larger components x=1.1902 and y=0.7691 is locally asymptotically stable as shown in the lower figure.

    Example 6. With the following parameters,

    a=2.25,γ=0.8,η1=0.02,η2=0.1, (3.26)

    we have the threshold values b1=1.1674, bs=1.0034>b0=0.8 , bc=1.0625 and a(γη2)=1.575>1+aη1=1.045. For b=0.6<b0, there exists a unique positive fixed point (4.0936,2.3237), which is locally asymptotically stable as shown in the upper left figure in Figure 8. For b=0.9>b0, there exist two positive fixed points (0.2713,0.2113) and (2.8617,1.7800), where smaller one is unstable while the larger one (2.8617,1.7800) is locally asymptotically stable since b<bs as shown in the upper right figure in Figure 8. For b=1.03>bs, the two positive fixed points (0.8742,0.6431) and (2.0203,1.3446) are both unstable and there are two synchronous 2-cycles exist, where the one with larger components x=3.4660 and y=2.0591 is locally asymptotically stable as shown in the lower figure in Figure 8.

    Figure 8.  With the parameters given in (3.26), we have b1=1.1674, bs=1.0034>b0=0.8, bc=1.0625 and a(γη2)=1.575>1+aη1=1.045. For b=0.6<b0, there exists a unique positive fixed point (4.0936,2.3237), which is locally asymptotically stable as shown in the upper left figure. For b=0.9>b0, there exist two positive fixed points (0.2713,0.2113) and (2.8617,1.7800), where the one with smaller components is unstable and the one with larger components (2.8617,1.7800) is locally asymptotically stable since b<bs as shown in the upper right figure. For b=1.03>bs, the two positive fixed points (0.8742,0.6431) and (2.0203,1.3446) are both unstable, and there are two synchronous 2-cycles, where the one with larger components x=3.4660 and y=2.0591 is locally asymptotically stable as shown in the lower figure.

    In this paper, we first formulated a discrete-time stage-structured mosquito model where the mosquito population is divided into two groups, the larvae and the adults. We assume that the survivability and progression of larvae are both of Beverton-Holt type nonlinearity. We determined the existence and stability for the positive fixed points and the synchronous 2-cycles, respectively. When the intrinsic growth rate of the population r0<1, the trivial fixed point is the only nonnegative fixed point and is globally asymptotically stable. If r0>1, the trivial fixed point is unstable, and the model dynamics depend also on parameters η1 and η2. If η1<η2, there exists a unique positive fixed point, which is globally asymptotically stable; if η1=η2, there exists a continuum of positive 2-cycles, each of which is locally stable; and if η1>η2, the positive fixed point becomes unstable, and there exists a unique synchronous 2-cycle, which is globally asymptotically stable.

    We then introduced sterile mosquitoes in the stage-structured wild mosquito population and considered three different strategies for the releases of sterile mosquitoes in model system (3.2) where the sterile mosquitoes are released constantly, (3.15) where the releases are proportional to the size of the wild mosquitoes, and (3.19) where the releases are of Holling-Ⅱ type, respectively. We established threshold value bc and b1 for the existence of the positive fixed points or synchronous 2-cycles, and threshold value bs for the stability of the positive fixed points or synchronous 2-cycles for each of the model systems. If b>bc, there exists no positive fixed point for all of the three model systems. If b<bc, there exist two positive fixed points for system (3.2) and (3.19), and a unique positive fixed point for system (3.15). A positive fixed point is locally asymptotically stable if b<bs, and unstable if b>bs. We also defined threshold value b1 for the existence of synchronous 2-cycle with b1>bc for system (3.2) and (3.19). If b<b1, there exist two synchronous 2-cycles, where the one with bigger components is locally asymptotically stable and the one with smaller components is unstable. If b>b1, there exist no synchronous 2-cycles and no positive fixed points. Thus the trivial fixed point is globally asymptotically stable. Details can be found in Theorems 3.1 and 3.2 and Tables 1, 2, and 3.

    We note that, in the absence of sterile mosquitoes, if the density-dependent death has less effect than the density-dependent progression from the larvae, that is, η1xn<η2xn, the structured population tends to be more stable in the sense that the positive fixed point is globally asymptotically stable while the unique synchronous 2-cycle is unstable. On the other hand, if the density-dependent death has more effect than the density-dependent progression from the larvae, that is, η1xn>η2xn, the structured population tends to be less stable such that the positive fixed point becomes unstable and it tends asymptotically to the oscillatory state, that is, the synchronous 2-cycle.

    Such dynamical features are similarly carried out when the sterile mosquitoes are released constantly or proportionally. More specifically, with the constant release rate and in the case of b<bs, there are two positive fixed points. When the density-dependent death has less effect than the density-dependent progression from the larvae such that η1xn<η2xn, one of the two positive fixed points is asymptotically stable, whereas the two positive fixed points are both unstable when the density-dependent death has more effect than the density-dependent progression from the larvae such that η1xn>η2xn. The picture for the proportional releases with saturation is not as clear as for the other two release strategies. Other parameters play a role as well.

    For any of the three release strategies, it is not surprising that the amount of releases changes the model dynamics. As small amounts of sterile mosquitoes are released, there exist stable positive fixed points or synchronous 2-cycles. When the release amount is gradually increasing greater than the stability thresholds, first positive fixed points or synchronous 2-cycles become unstable, and then all disappear leading to the extinction of wild mosquitoes.

    We would also like to point out that the outcomes from the models studied in this paper seem to be similar to those with the Ricker-type nonlinearity in [27]. However, it is well known that the Beverton-Holt nonlinearity excludes the possibility of the period doubling bifurcation and chaotic feature for the models without stage structure. When stage structure is included, the relatively simple dynamics with no period doubling bifurcation and chaotic feature are carried out, which makes the analysis more tractable. Such model structure has been applied to the discrete-time malaria transmission models incorporating releases of sterile mosquitoes in [32].

    The authors thanks the two anonymous reviewers for their careful reading and valuable comments and suggestions.

    All authors declare no conflicts of interest in this paper.

    Proof. To prove that there exist no positive k cycles for system (2.9), we consider k=3 first. The model (2.9) then becomes

    xn+3=a2γyn1+aηyn+η(r0+1)xn,yn+3=aγ2xn1+aηyn+η(r0+1)xn. (4.1)

    Any point (x,y) on a 3-cycle should satisfy xr01η, 1+ηx+aηyr0 and y=γax. According to the first equation of (4.1), we have

    1+aηy+η(r0+1)x=a2γyx.

    Plugging y=γax into it, we have

    x=r01η,

    which is exactly the fixed point from (2.8). Thus 3-cycles do not exist and as a consequence, any (2n+1)-cycles do not exist with integer n>0.

    We next check for 4-cycles. If k=4, the model (2.9) becomes

    xn+4=r20xn1+η(r0+1)xn+aη(r0+1)yn,yn+4=r20yn1+η(r0+1)xn+aη(r0+1)yn. (4.2)

    For a point (x,y) on 4-cycle, it satisfies

    r20=1+η(r0+1)x+aη(r0+1)y.

    Equivalently, we can write it in this form

    r201=(r0+1)(ηx+aηy),

    that is,

    r0=1+ηx+aηy.

    It is exactly the positive 2-cycle from (2.10). Thus there exist no 4-cycles, and as a consequence, there exist no 2n-cycles for integer n>1.

    Therefore, there exist no positive kcycles for k3.



    [1] L. Alphey, M. Benedict, R. Bellini, G. G. Clark, D. A. Dame, M. W. Service and S. L. Dobson, Sterile-insect methods for control of mosquito-borne diseases: An analysis, Vector-Borne Zoonot., 10 (2010), 295–311.
    [2] AMCA, Life Cycle. Available from: http://www.mosquito.org/page/lifecycle.
    [3] B. Anderson, J. Jackson and M. Sitharam, Descartes's rule of signs revisited, Amer. Math. Monthly 105 (1998), 447–451.
    [4] H. J. Barclay, The sterile insect release method for species with two-stage life cycles, Res. Popul. Ecol., 21 (1980), 165–180.
    [5] H. J. Barclay, Pest population stability under sterile releases, Res. Popul. Ecol., 24 (1982), 405– 416.
    [6] H. J. Barclay, Mathematical models for the use of sterile insects, in Sterile Insect Technique, Principles and Practice in Area-Wide Integrated Pest Management, V. A. Dyck, J. Hendrichs, and A. S. Robinson , (Eds. Springer, Heidelberg), (2005), 147–174.
    [7] H. J. Barclay and M. Mackuer, The sterile insect release method for pest control: a density dependent model, Environ. Entomol., 9 (1980), 810–817.
    [8] A. C. Bartlett and R. T. Staten, Sterile Insect Release Method and other Genetic Control Strategies, Radcliffe's IPM World Textbook, 1996, available from: http://ipmworld.umn.edu/chapters/bartlett.htm.
    [9] N. Becker, Mosquitoes and Their Control, Kluwer Academic/Plenum, New York, 2003.
    [10] M. Q. Benedict and A. S. Robinson, The first releases of transgenic mosquitoes: an argument for the sterile insect technique, Trends Parasitol., 19 (2003), 349–55.
    [11] R. J. H. Beverton and S. J. Holt, On the Dynamics of Exploited Fish Populations, Springer Science and Business Media, Jun 30, 1993.
    [12] M. Bohner and H. Warth, The Beverton-Holt dynamic equation, Appl. Anal., 86 (2007), 1007– 1015.
    [13] L. M. Cai, S. B. Ai and J. Li, Dynamics of mosquitoes populations with different strategies for releasing sterile mosquitoes, SIAM, J. Appl. Math., 74 (2014), 1786–1809.
    [14] J. M. Cushing and J. Li, On Ebenman's model for the dynamics of a population with competing juveniles and adults, Bull. Math. Biol., 51 (1989), 687–713.
    [15] C. Dye, Intraspecific competition amongst larval aedes aegypti: Food exploitation or chemical interference, Ecol. Entomol., 7 (1982), 39–46.
    [16] S. Elaydi, Discrete Chaos: With Applications in Science and Engineering, 2nd edition, Boca Raton, CRC Press, 2007.
    [17] K. R. Fister, M. L. McCarthy, S. F. Oppenheimer and Craig Collins, Optimal control of insects through sterile insect release and habitat modification, Math. Biosci., 244 (2013), 201–212.
    [18] J. C. Floresa, A mathematical model for wild and sterile species in competition: immigration, Physica A, 328 (2003), 214–224.
    [19] R. M. Gleiser, J. Urrutia and D. E. Gorla, Effects of crowding on populations of aedes albifasciatus larvae under laboratory conditions, Entomol. Exp. Appl., 95 (2000), 135–140.
    [20] D. J. Grabiner, Descartes's rule of signs: Another construction, Amer. Math. Monthly, 106 (1999), 854–855.
    [21] C. S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Can., 45 (1965), 5–60.
    [22] E. I. Jury, The inners approach to some problems of system theory, IEEE, Trans. Automatic Contr., AC-16 (1971), 233–240.
    [23] E. S. Krafsur, Sterile insect technique for suppressing and eradicating insect populations: 55 years and counting, J. Agr. Entomol., 15 (1998), 303–317.
    [24] J. Li, Simple stage-structured models for wild and transgenic mosquito populations, J. Diff. Eqns. Appl., 17 (2009), 327–347.
    [25] J. Li, Modeling of mosquitoes with dominant or recessive transgenes and Allee effects, Math. Biosci. Eng., 7 (2010), 101–123.
    [26] J. Li, Malaria model with stage-structured mosquitoes, Math. Biol. Eng., 8 (2011), 753–768.
    [27] J. Li, Stage-structured models for interacting wild and sterile mosquitoes, Journal of Shanghai Normal University, 43 (2014), 511–522.
    [28] J. Li, New revised simple models for interactive wild and sterile mosquito populations and their dynamics, J. Biol. Dynam., 11(S2) (2017). 316–333.
    [29] J. Li and Z. L. Yuan, Modeling releases of sterile mosquitoes with different strategies, J. Biol. Dynam., 9 (2015), 1–14.
    [30] Y. Li, Discrete-time Structured Models and their Dynamics for Interactive Wild and Sterile Mosquitoes and Malaria Transmissions, Ph.D thesis, The University of Alabama in Huntsville, 2017.
    [31] Y. Li and J. Li, Discrete-time models for releases of sterile mosquitoes with Beverton-Holt-type of survivability, Ricerche di Matematica, 67 (2018), 141–162.
    [32] Y. Li and J. Li, Discrete-time model for malaria transmission with constant releases of sterile mosquitoes, preprint, 2018.
    [33] R. M. May, G. R. Conway, M. P. Hassell and T. R. E. Southwood, Time delays, density-dependence and single-species oscillations, J. Anim. Ecol., 43 (1974), 747–770.
    [34] Mosquito Facts, 33 Things You Didn't Know About Mosquitoes, https://www.megacatch. com/mosquito-faqs/mosquito-facts/.
    [35] M. Otero, H. G. Solari and N. Schweigmann, A stochastic population dynamics model for Aedes aegypti: formulation and application to a city with temperate climate, Bull. Math. Biol., 68 (2006), 1945–1974.
    [36] D. Ruiz, G. Poveda, I. D. Velez, M. L. Quinones, G. L. Rua, L. E.Velasquesz and J. S. Zuluage, Modeling entomological-climatic interactions of Plasmodium falciparum malaria transmission in two Colombian endemic-regions: contributions to a National Malaria Early Warning System, Malaria J., 5 (2006), 66.
    [37] Wikipedia, Mosquito-borne disease. Available from: https://en.wikipedia.org/wiki/ Mosquito-borne_disease.
    [38] Wikipedia, Sterile insect technique. Available from: https://en.wikipedia.org/wiki/ Sterile_insect_technique.
  • This article has been cited by:

    1. Manuel De la Sen, Asier Ibeas, Aitor J. Garrido, Stage-Dependent Structured Discrete-Time Models for Mosquito Population Evolution with Survivability: Solution Properties, Equilibrium Points, Oscillations, and Population Feedback Controls, 2019, 7, 2227-7390, 1181, 10.3390/math7121181
    2. Martin Bohner, Jaqueline Mesquita, Sabrina Streipert, The Beverton–Hold model on isolated time scales, 2022, 19, 1551-0018, 11693, 10.3934/mbe.2022544
    3. Jianshe Yu, Jia Li, Discrete-time models for interactive wild and sterile mosquitoes with general time steps, 2022, 346, 00255564, 108797, 10.1016/j.mbs.2022.108797
    4. Manuel De la Sen, Santiago Alonso-Quesada, Asier Ibeas, Aitor J. Garrido, Ewa Pawluszewicz, On an Extended Time-Varying Beverton–Holt Equation Subject to Harvesting Monitoring and Population Excess Penalty, 2023, 2023, 1607-887X, 1, 10.1155/2023/5052799
    5. Yijun Lou, Ruiwen Wu, Modeling insect growth regulators for pest management, 2024, 88, 0303-6812, 10.1007/s00285-024-02091-y
    6. Shouzong Liu, Yang Xu, Mingzhan Huang, Comparative analysis of sterile mosquito release strategies based on a population suppression model, 2024, 9, 2473-6988, 23344, 10.3934/math.20241135
    7. Mensah Folly-Gbetoula, On a family of higher order recurrence relations: symmetries, formula solutions, periodicity and stability analysis, 2023, 12, 2193-5343, 541, 10.1007/s40065-023-00438-9
    8. Manuel De la Sen, Santiago Alonso-Quesada, Asier Ibeas, Aitor J. Garrido, A. E. Matouk, On a Coupled Time-Varying Beverton–Holt Model with Two Habitats Subject to Harvesting, Repopulation, and Mixed Migratory Flows of Populations, 2023, 2023, 1607-887X, 1, 10.1155/2023/6050789
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5086) PDF downloads(742) Cited by(8)

Figures and Tables

Figures(8)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog