Dynamic analysis of a mosquito population model with a stage structure and periodic releases of sterile males

1.   Introduction

Mosquito-borne diseases, such as malaria, dengue fever and yellow fever, are highly detrimental infectious diseases transmitted by mosquitoes. These diseases are prevalent worldwide and pose a significant threat to human life and health. Effectively controlling mosquito-borne diseases has become a primary focus for public health departments and researchers. Since there are no effective vaccines for these diseases, suppressing mosquito populations has emerged as a crucial strategy for controlling their spread. In recent times, apart from conventional physical and chemical approaches, the emergence of biological techniques has been observed, which include the release of sterilized mosquitoes or those carrying the bacterium Wolbachia. This approach has gained widespread recognition and proven to be highly effective in suppressing wild mosquito populations. Numerous studies delved specifically into the use of sterile insect release as a means for targeted population control ^[1,2]. Concurrently, some research investigated the role of the endosymbiotic bacterium, Wolbachia, in inducing resistance to the dengue virus within Aedes aegypti ^[3,4]. Efforts have also been made to devise methods for controlling specific mosquito-borne diseases ^[5,6,7]. Throughout these explorations, comprehensive analyses of distinct population characteristics, including population diffusion and periodic changes, among others, have been undertaken ^[8,9]. In addition, a study ^[10] has successfully implemented a blend of incompatible and sterile insect techniques to eradicate wild mosquitoes.

Researchers have employed a diverse range of mathematical models to simulate various release modes of sterile mosquitoes, thereby exploring effective strategies under distinct scenarios. For instance, uninterrupted continuous release strategies, which are grounded in ordinary differential equation or delay differential equation models, have been thoroughly investigated ^[11,12,13]. Zheng, along with other researchers, developed several switched systems to scrutinize the impact of releasing Wolbachia-infected mosquitoes ^[14,15,16]. Moving beyond one-dimensional models, the authors of ^[17] devised and examined a suite of high-dimensional ODE models with a stage structure. Furthermore, the release strategies for both extended and finite periods were probed from an optimal control perspective ^[18,19], while ^[20] paid specific attention to impulsive release strategies within the framework of open-loop and closed-loop control. Yu and his colleague introduced the notion of the sexual lifespan of sterile mosquitoes during the construction of mosquito population suppression models ^[21,22,23]. Building upon this foundation, scholars initiated a series of dialogues on the interplay between sexual lifespan and pulse release strategies. They further constructed some switched dynamic systems incorporating delay effects ^[24,25] and the stage structure ^[26,27,28].

Frequently, sterile male mosquitoes are introduced into the environment to vie for copulation chances alongside their wild counterparts. When females engage in reproduction with these sterile males, they either cannot lay eggs or their eggs do not hatch. Consequently, repeated release of sterile male mosquitoes can progressively diminish the wild mosquito population. However, male mosquitoes have a relatively short lifespan and remain sexually active for even shorter periods. Only sexually active sterile mosquitoes can effectively compete with wild male mosquitoes for mating opportunities ^{[12,14,22,23,25]}. Considering these factors, Yu and his collaborators suggested that only sexually active sterile mosquitoes should be taken into account in the dynamic system of interaction between the two types of mosquitoes. Moreover, since sterile mosquitoes remain sexually active for a brief period, the impact of their natural death can be neglected when constructing mathematical models ^{[14,22,23,25]}. For instance, Yu and his partners gave nearly comprehensive analysis of the population suppression model based on the following equation

where $w(t)$ and $g(t)$ represent the numbers of wild mosquitoes and sexually active sterile ones in the field, respectively. The birth rate per wild mosquito, denoted by $a$, adheres to logistic growth. $\xi$ and $\mu$ correspond to the density-dependent and density-independent mortality rates, respectively.

Lin and Hui in ^[27] investigated a similar model, grounded in the same assumptions and research methodology. Meanwhile, Ai et al. in ^[26] explored a two-dimensional stage structured model of interacting wild and sterile mosquitoes, wherein the wild mosquito population comprises both larvae and adults.

The mosquito life cycle is known to consist of several stages, each with unique environmental adaptation necessities. For example, the development of eggs and larvae during the aquatic phase is substantially impacted by competition within the species, while adult mosquitoes are barely affected. Moreover, the released mosquitoes are primarily sterile males, which compete with wild male mosquitoes for mating opportunities with wild female mosquitoes ^[29]. Consequently, it is essential for mathematical models to factor in the various growth stages of mosquitoes and distinguish between male and female groups. We observe that most existing mosquito population suppression models overlook the stage structure of wild mosquitoes when characterizing their growth. Although a few models consider the stage structure of mosquito growth, such as ^[26], the sex structure is ignored.

In this work, we will develop a mosquito population model incorporating both stage and sex structure, as well as periodic releases of sterile individuals, while considering only sterile mosquitoes with mating ability in the dynamic system.

We organize this paper as follows. We will first develop a novel mosquito population model incorporating stage and sex structure in Section 2, in which sterile males are assumed to be periodically released and possess mating capability only for a limited duration. Then we will examine the stability of the extinction equilibrium and the presence of positive periodic solutions of the proposed model in Section 3. Subsequently, some numerical simulations are presented in Section 4 to validate the theoretical findings. Finally, a brief conclusion will be offered in Section 5.

2.   Model formulation

Almeida et al. in ^[30] considered a population model for mosquitoes as follows:

where $E(t), L(t), P(t), F(t)$ and $M(t)$ represent the number of eggs, larvae, pupae, adult females and adult males in the field, respectively. $\beta_E F (1-\frac{E}{K})$ is interpreted as a "skip oviposition" behavior whereby females are capable of selecting their breeding sites and seeking oviposition sites with high food content and low intraspecific competition pressure ^[6]. $\beta_E$ denotes the oviposition rate, while $K$ and $\tau_E$ are the environmental capacity and the hatching rate of eggs, respectively. $\nu$ and $1-\nu$ respectively represent the proportion of pupae developing into females and males. $\tau_L$ and $\tau_P$ are transition rates, $c$ represents the intraspecific competition of larvae and $\delta_E, \delta_L, \delta_P, \delta_F, \delta_M$ are the respective death rates.

To simplify this system, the authors of ^[30] proposed the following assumptions:

(Ⅰ) Assume that the time dynamics in the larvae and pupae compartments are rapid.

(Ⅱ) Assume that the competition at larval stage is negligible (i.e., $c \ll 1$).

(Ⅲ) Assume that males and females have the same death rate ($\delta_F = \delta_M$) and the probability of pupae emerging as females or males is equal ($\nu = \frac{1}{2}$) such that male and female consequently satisfy the same equation. Without loss of generality, we can let $F = M$.

Under these assumptions, the equations of $P(\cdot)$ and $L(\cdot)$ are at equilibrium. That is

By incorporating (2.2) into system (2.1) and introducing a release function $u(\cdot)$ for sterile male mosquitoes, we can derive the following simplified model for the interaction between wild and sterile mosquitoes, as presented in ^[30]:

where $M_s(t)$ represents the number of sterile males at time $t$ and $\gamma$ describes the mating competitiveness of sterile males. A new notation, $\beta_F = \frac{\tau_P\tau_L\tau_E}{(\tau_P+\delta_P)(\tau_L+\delta_L)}$, is introduced.

In this paper, akin to ^{[12,14,22,23,25]}, we present the concept of the sexual lifespan of sterile mosquitoes. We denote the number of sterile mosquitoes with mating ability at time $t$ as $M_s(t)$ and disregard the natural mortality of sexually active sterile mosquitoes. Consequently, system (2.3) can be rewritten in the following form:

Assume that sterile males are impulsively introduced into the field in quantities of $m$ at discrete time $t = kT, k = 0, 1, 2, \cdots$, with $T$ representing the period of releases. Furthermore, let $\bar{T}$ denote the sexual lifespan of sterile males; upon release, they will lose their mating capabilities after this time span. There are three cases concerning the relationship between $T$ and $\bar{T}$: (1) $T = \bar{T}$; (2) $T > \bar{T}$; (3) $T < \bar{T}$.

For the first case, we can easily deduce $M_s(t) \equiv m$ for $t\geq 0$. As a result, system (2.4) transforms into:

A number of experimental investigations have shown that the sexual lifespan of sterile males is generally short ^[14]. Consequently, in some prior research, it has been assumed that the sexual lifespan of sterile males is shorter than the period during which they are released. In the present study, we also examine the case where $T > \bar{T}$. Under such circumstances, the sexually viable sterile males in the environment are not replenished at a sufficient rate, which implies that when the formerly released sterile ones have lost their capacity to mate, new sterile mosquitoes have not been introduced yet. In this scenario, the quantity of sterile males possessing mating capabilities, $M_s(t)$, is a stepwise constant function with period $T$ as follows

Subsequently, system (2.4) undergoes a transformation

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

Clearly, system (2.7) is a switched system, consisting of two continuous-time subsystems and a governing rule that directs the transition between these two subsystems.

In this work, we will study the dynamic behaviors of system (2.4) for the cases $T = \bar{T}$ and $T > \bar{T}$, that is, the dynamic behaviors of systems (2.5) and (2.7). We mainly investigate the stability of the extinction equilibrium and the existence of positive periodic solutions.

3.   Dynamic analysis

3.1.   Preliminaries

First, we will examine the scenario in which no sterile males are released into the environment, meaning that $M_s(t) = 0$. In this case, system (2.4) becomes

Through direct calculation, we can get that when $\nu\beta_E\beta_F > \delta_F(\tau_E+\delta_E)$, system (3.1) possesses a positive equilibrium $A^*(E^*, F^*)$, where

When $\nu\beta_E\beta_F\leq\delta_F(\tau_E+\delta_E)$, system (3.1) has no positive equilibrium.

Denote

then, system (3.1) possesses a unique positive equilibrium $A^*(E^*, F^*)$ when $\aleph_0 > 1$, where $E^* = (1-\frac{1}{\aleph_0})K, \quad F^* = \frac{\nu\beta_F}{\delta_F}(1-\frac{1}{\aleph_0})K.$ Obviously, $\aleph_0$ is the classical basic offspring number of wild mosquitoes ^[6].

According to system (2.4), we can easily get that

is a globally attractive positive invariant set of (2.4). In this paper, we mainly study system (2.4) in the region $\Omega$.

In addition, by selecting Dulac function $B(E, F) = 1$, we can verify that system (2.4) has no closed orbits in $R^2_+$. Then we can obtain the subsequent result related to (3.1) when $M_s(t) = 0$.

Theorem 1. For system (3.1), if $\aleph_0\leq1$, the extinction equilibrium $A_0(0, 0)$ is a globally asymptotically stable node and no positive equilibrium exists. Conversely, if $\aleph_0 > 1$, the system has a unique positive equilibrium $A^*(E^*, F^*)$ that is globally asymptotically stable, while the extinction equilibrium $A_0$ becomes unstable.

Proof. To discuss the stability of the equilibria, we first calculate the Jacobian matrix of (2.5) as follows

Then we can deduce

and

Obviously, $tr(J_{A_0}) < 0$ and $tr(J_{A^*}) < 0$. □

When $\aleph_0 < 1$, we can easily get $Det(J_{A_0}) > 0$ and that the unique equilibrium $A_0(0, 0)$ is an asymptotically stable node. Since system (2.4) has no closed orbits in $R^2_+$ and there is no other equilibrium, we can claim that $A_0(0, 0)$ is the positive limit set of all trajectories in $R^2_+$. That is, $A_0(0, 0)$ is globally asymptotically stable. When $\aleph_0 > 1$, we have that $Det(J_{A_0}) < 0, Det(J_{A^*}) = \frac{\delta_F\beta_EF^*}{K} = \delta_F\beta_E\beta_F(\aleph_0-1)$. Then $A_0(0, 0)$ is unstable and $A^*(E^*, F^*)$ is asymptotically stable. Since $A^*(E^*, F^*)$ is the unique equilibrium which is stable, we can claim that it is the positive limit set of all trajectories in $R^2_+$. Thus $A^*(E^*, F^*)$ is globally asymptotically stable if $\aleph_0 > 1$. The proof is completed.

Given that our focus is on cases where the wild mosquito population in the field is significantly large and requires control measures, we assume $\aleph_0 > 1$ in this study.

In what follows, we investigate the dynamic properties of systems (2.5) and (2.7) within the region $\Omega$ when $\aleph_0 > 1$.

3.2.   Constant releases

We will now examine the case of constant releases, where $M_s(t)\equiv m > 0$. In this case, system (2.4) is transformed into system (2.5).

It is easy to obtain that the extinction equilibrium $A_0$ of (2.5) is locally asymptotically stable. Next, we explore the existence of positive equilibria. For convenience, let

In order to confirm the existence of a positive steady state $(E, F)$, we must solve the subsequent algebraic equations

By direct calculation, we get $E = \frac{F\delta_F}{\nu\beta_F}$, and $F$ is a positive root of

Let

then the discriminant of this quadratic equation is

Define a threshold value for the release amount $m$ as

then we have the following conclusions.

Lemma 1. (i) If $0 < m < m^*$, then $\Lambda(F) > 0$ and system (2.5) has two positive equilibria $A_1(E_1, F_1)$ $A_2(E_2, F_2)$, where

(ii) If $m = m^*$, then $\Lambda(F) = 0$ and system (2.5) has a unique positive equilibrium $A_{m^*}(E_{m^*}, F_{m^*})$, where

(iii) If $m > m^*$, then $\Lambda(F) < 0$ and system (2.5) has no positive equilibrium point.

To determine the stability of positive equilibria, we compute the Jacobian matrix of (2.5) as follows

We can easily get that $tr J = \frac{-\nu\beta_E\beta_FF}{\delta_F(m_1+F)}-\delta_F < 0$ and

where

Since $-H'(F) = -\frac{\Theta(F)}{\aleph_0(m_1+F)^2}+\frac{\Theta'(F)}{\aleph_0(m_1+F)}$ and $\Theta(F_1) = \Theta(F_2) = 0$, we can obtain that $|J(A_1)| < 0$ and $|J(A_2)| > 0$. Therefore, $A_1(E_1, F_1)$ is an unstable saddle, while $A_2(E_2, F_2)$ is locally asymptotically stable.

To sum up, we get the subsequent conclusions.

Theorem 2. For system (2.5), the extinction equilibrium $A_0$ is always locally asymptotically stable. In addition, when $\aleph_0 > 1$,

(i) if $0 < m < m^*$, then $A_1(E_1, F_1)$ is an unstable saddle, while $A_2(E_2, F_2)$ is locally asymptotically stable;

(ii) if $m = m^*$, then the unique positive equilibrium $A_{m^*}(E_{m^*}, F_{m^*})$ is a saddle node;

(iii) if $m > m^*$, then the extinction equilibrium $A_0$ is globally asymptotically stable.

3.3.   Periodic releases

In case where the release period of sterile males is longer than the sexual lifespan, specifically, $\bar{T} < T$, the number of sterile males with mating ability $M_s(t)$ switches between two levels, as illustrated in (2.6). Consequently, system (2.4) is converted to system (2.7).

Evidently, systems (2.5) and (3.1) are two subsystems of (2.7). The continual switching between these two subsystems is crucial to the dynamics of (2.7). With the occurrence of switching, the extinction equilibrium $A_0(0, 0)$ is the unique equilibrium of system (2.7). We first study the stability of $A_0$.

If $\aleph_0 = \frac{\nu\beta_E\beta_F}{\delta_F{(\tau_E+\delta_E)}} > 1$, then we have that $\nu\beta_E\beta_F > \delta_F(\tau_E+\delta_E)$. Since $\beta_F = \frac{\tau_P\tau_L\tau_E}{(\tau_P+\delta_P)(\tau_L+\delta_L)} < \tau_E$ and $\nu\in(0, 1)$, then we get that $\tau_E+\delta_E > \nu\beta_F$ and $\beta_E > \delta_F$.

Let

Obviously, $\lambda > \lambda_0\geq \lambda_1$. Given any $0 < \sigma < 1$, we further define

Theorem 3. If $\aleph_0 > 1$ and $\bar{T} < T < \omega\bar{T}$, then for any $E(0)+F(0) < \theta$, the inequality

holds, where $\kappa = \lambda_1 (1-\sigma)\frac{\bar{T}}{T}-(\lambda-\lambda_0)(1-\frac{\bar{T}}{T}) > 0$. Thus the extinction equilibrium $A_0(0, 0)$ of system (2.7) is locally asymptotically stable.

Proof. Suppose that $(E(t), F(t))$ is a solution of system (2.7) with $E(0)+F(0) < \theta$ for a fixed $\sigma\in(0, 1)$. Let $Q(t) = E(t)+F(t)$, $t\in[kT, kT+\bar{T}]$, $k = 0, 1, 2, \cdots$.

Based on the first two equations of system (2.7), we can deduce that

If $Q(kT) < \theta$, we can assert that $Q(t) < \theta$ for $t \in [kT, kT+\bar{T}]$. We will prove this statement via contradiction. Suppose that there exists a $\overline{t} > 0$ for which $Q(t) < \theta$ when $t\in [kT, \overline{t})$, but that $Q(\overline{t}) = \theta$. Consequently, it follows that $Q'(\overline{t})\geq 0$. However, from (3.7), we have

which leads to a contradiction.

Thus we can get

which can yield $Q'(t) < 0$ and $Q(t) < Q(kT) < \theta$, $t\in[kT, kT+\bar{T}]$. Therefore we obtain

Similarly, according to the last two equations of system (2.7) for $t\in[kT+\bar{T}, (k+1)T]$, we can easily get

and then

By (3.8) and (3.9), we obtain

where $\eta: = e^{-(1-\sigma)\lambda_1\bar{T}+(\lambda-\lambda_0)(T-\bar{T})} < 1$.

Obviously, $Q((k+1)T)\leq \eta Q(kT)$ and $Q(kT) < \eta^k Q(0), \; k = 0, 1, 2, \cdots.$

If $t\in[kT, kT+\bar{T}]$, then $k > \frac{(t-\bar{T})}{T}$, and from (3.8) we know that

If $t\in[kT+\bar{T}, (k+1)T]$, then $k+1 > \frac{t}{T}$, and from $(3.10)$ we know that

According to (3.11) and (3.12), we can easily obtain (3.6), which indicates the local asymptotic stability of the extinction equilibrium $A_0(0, 0)$ of system (2.7). The proof is completed. □

Clearly, Theorem 3.3 only provides sufficient conditions for the local stability of $A_0(0, 0)$. Owing to the limitations of the theorem's conditions, the applicability of this result is somewhat restricted. Apart from examining $A_0(0, 0)$, we are also intrigued by the periodic fluctuations in mosquito population levels caused by the recurring release of sterile males. We will now investigate the presence of positive periodic solutions for system (2.7) in scenarios where $0 < m < m^*$ and $m > m^*$, respectively.

Switched system (2.7) is composed of two subsystems (2.5) and (3.1). According to the previous discussion, if $\aleph_0 > 1$ and $0 < m < m^*$, system (2.5) has two positive points $A_i(E_i, F_i)$, $i = 1, 2$ (see Lemma 1), while system (3.1) has a unique positive point $A^*(E^*, F^*)$ given in (3.2). It is easy to verify that $F^* > F_2 > F_1$ and $K > E^* > E_2 > E_1$.

Let

be a rectangle in the $E F$-phase plane; then, we obtain the subsequent result.

Theorem 4. If $\aleph_0 > 1$ and $0 < m < m^*$, then system (2.7) exhibits a continuous $T$-periodic solution within $\Omega_0$. Furthermore, if a solution $(E(t), F(t))$ of system (2.7) fulfills the conditions $E(t_0)\geq E_1$ and $F(t_0)\geq F_1$ for a certain $t_0\geq 0$, then $E(t)\geq E_1$ and $F(t)\geq F_1$ hold for all $t > t_0$ and dist$((E(T), F(t)), \Omega_0)\rightarrow 0$ as $t\rightarrow \infty$. In addition, if such a solution is periodic, then $(E(t), F(t))\in\Omega_0$ for all $t\geq0$.

Proof. Let $\zeta(E): = \frac{K(\tau_E+\delta_E)E}{K-E}, E\in[0, K]$; then, the vertical isoclinics $E' = 0$ of (2.5) and (3.1) are

respectively. Since $H_1(E) > H_2(E)$, the positions of these two isoclinics are as shown in Figure 1. In addition, $F = G(E): = \frac{\nu\beta_E {E}}{\delta_F}$ is the shared horizontal isoclinic of (2.5) and (3.1).

Take any two positive numbers $E_m, E_M$ such that $E_1 < E_m < E_2 < E^* < E_M < K$ and denote $F_m: = G(E_m)$ and $F_M: = G(E_M)$. Since $G(E)$ is strictly monotonically increasing with respect to $E$, we have that $F_m < F_M$. Now let $\Omega_1$ be the closed rectangle with the vertices $P_1, P_2, Q_1$ and $Q_2$, where $P_1 = (E_m, F_m), P_2 = (E_M, F_m), Q_1 = (E_m, F_M)$ and $Q_2 = (E_M, F_M)$ (see Figure 1).

Obviously, $F_1 < F_m < F_2 < F^* < F_M$ and $\Omega_0$ lies in the inside of $\Omega_1$. Given the strictly monotonic increasing of $G(E)$, the segment $\overline{P_1P_2}$ lies strictly below the isoclinic $F = G(E)$, except the point $P_1$. Similarly, the whole segment $\overline{Q_1Q_2}$ lies strictly above the isoclinic $F = G(E)$ except the point $Q_2$. In addition, from the second equations of (2.5) and (3.1), we can derive

and

Therefore, the vector fields of (2.7) point to the interior of $\Omega_1$ on $\overline{P_1P_2}$ and $\overline{Q_1Q_2}$.

Furthermore, because $H_1(E)$ and $H_2(E)$ are both strictly monotonically increasing for $K > E > 0$, we can easily obtain that the whole segment $\overline{P_1Q_1}$ lies strictly above the isoclinic $F = H_1(E)$ and the whole segment $\overline{P_2Q_2}$ lies strictly below the isoclinic of $F = H_2(E)$. From system (2.5), we can obtain

where $H^-_{1}(E) = \frac{1}{2\beta_E}[-\zeta(E)+\sqrt{\zeta^2(E)+4\beta_E m_1 \zeta(E)}] > 0$.

From system (3.1), we have

Then we can easily get that $F > H_1(E_m)$ and $F > H_2(E_m)$ hold on $\overline{P_1Q_1}$; therefore, $\frac{dE}{dt} > 0$ holds on $\overline{P_1Q_1}$ for both (2.5) and (3.1). Similarly, $F < H_1(E_M)$ and $F < H_2(E_M)$ hold on $\overline{P_2Q_2}$ so that $\frac{dE}{dt} < 0$ hold on $\overline{P_2Q_2}$ for both systems. Then the vector fields of (2.7) point to the inside of $\Omega_1$ on segments $\overline{P_1Q_1}$ and $\overline{P_2Q_2}$.

To sum up, $\Omega_1$ is a positively invariant set of system (2.7). Then for any point $A\in \Omega_1$, the trajectory of (2.7) from point $A$ will always stay in $\Omega_1$. Assume that $\phi_t (A)$ and $\psi_t (A)$ are the solutions of (2.5) and (3.1) starting from point $A$, respectively. Then both $\phi_t(A), t\in[0, \bar{T}]$ and $\psi_{t-\bar{T}}(\phi_{\bar{T}}(A)), t\in[\bar{T}, T]$ remain in $\Omega_1$.

Similar to the discussion in ^[26], define a Poincare map $\Phi:\Omega_1\rightarrow \Omega_1$ as follows

Since $\Phi$ is continuous and maps $\Omega_1$ into $\Omega_1$, by applying Brouwer's fixed point theorem, we can conclude that $\Phi$ possesses a fixed point, denoted by $\tilde{A}\in\Omega_1$. Then the solution of (2.7) through $\tilde{A}$ is a continuous $T$-periodic solution which lies entirely in $\Omega_1$.

We note that when $E_m\rightarrow E_2, E_M\rightarrow E^*$, it follows that $\Omega_1\rightarrow \Omega_0$. According to the arbitrariness of $E_m, E_M$ and the above discussion, we can deduce that the continuous $T$-periodic solution of (2.7) in $\Omega_1$ lies entirely in $\Omega_0$.

In addition, when $E_m\rightarrow E_1$ and $E_M\rightarrow \infty$, it follows that $\Omega_1\rightarrow [E_1, \infty)\times[F_1, \infty)$. Therefore, $[E_1, \infty)\times[F_1, \infty)$ is also a positively invariant set of system (2.7) for which $\Omega_0$ is a small attraction domain. Then for any solution $(E(t), F(t))$ of system (2.7) that satisfies $E(t_0)\geq E_1$ and $F(t_0)\geq F_1$, $t_0\geq 0$, it follows that $E(t)\geq E_1$ and $F(t)\geq F_1$ for all $t > t_0$ and dist$((E(T), F(t)), \Omega_0)\rightarrow 0$ as $t\rightarrow \infty$. Moreover, if such a solution is periodic, then $(E(t), F(t))\in\Omega_0$ for all $t\geq0$. The proof is completed. □

In the case where the release amount $m$ exceeds the release threshold $m^*$, meaning that $m > m^*$, we obtain the subsequent result concerning the existence of a positive periodic solution.

Theorem 5. If $\aleph_0 > 1$ and $m > m^*$, then there must exist $T_0 > 0$ and $\delta_0 > 0$ such that system (2.7) exhibits a positive and continuous $T$-periodic solution when $\bar{T} < T_0$ and $T-\bar{T} > \delta_0$.

Proof. It is easy to obtain that for system (2.5), there is no positive equilibrium, whereas system (3.1) possesses a unique positive equilibrium when $\aleph_0 > 1$ and $m > m^*$.

Similar to the proof in Theorem 4, we select arbitrarily two positive numbers $E_m$ and $E_M$ such that $E_m < E^* < E_M$ and $G(E_M) < H_1 (E_m)$. Using the same donations as in Theorem 4, let $F_m: = G(E_m)$, $F_M: = G(E_M)$ and $\Omega_1$ be the closed rectangle with vertices $P_1 : = (E_m, F_m)$, $P_2: = (E_M, F_m)$, $Q_1: = (E_m, F_M)$ and $Q_2: = (E_M, F_M)$. In addition, we denote the intersection point of the isoclinic $F = H_2 (E)$ and the segment $\overline{P_1P_2}$ by $P_3(E_3, F_m)$, and then we arbitrarily choose $E_0\in (E_m, E_3)$ and denote $P_0 : = (E_0, F_m)$ and $Q_0: = (E_0, F_M)$ (see Figure 2).

Obviously, the rectangle $\Omega_1$ can be divided into two parts which are denoted by $\Omega_1^1$ with vertices $P_1$, $P_0$, $Q_1$ and $Q_0$ and $\Omega_1^2$ with vertices $P_0$, $P_2$, $Q_0$ and $Q_2$, respectively. Denote

where $G_i(E, F), i = 1, 2$ correspond to the right-hand side expressions of the first equations in systems (2.5) and (2.7), respectively. Then let $T_0: = \frac{E_0 -E_m}{|\alpha_1|}$ and $\delta_0 : = \frac{E_0 -E_m}{\alpha_2}$ and assume that

then in what follows we will prove that the Poincare map $\Phi$ defined in (3.13) is continuous and maps $\Omega_1^2$ into $\Omega_1^2$. Thus by the same discussion method in Theorem 4, it can be deduced that system (2.7) possesses a positive and continuous $T$-periodic solution.

First, we demonstrate that for any given solution $(E(t), F(t))$ of (2.5) with $(E(0), F(0))\in \Omega_1^2$, it follows that $(E(t), F(t))\in \Omega_1$ for $t\in[0, \bar{T}]$. Analogous to the analysis in Theorem 4, it becomes evident that the vector field of (2.5) points to the inside of $\Omega_1$ on the segments $\overline{P_1P_2}$, $\overline{Q_1Q_2}$ and $\overline{P_2Q_2}$. Thus for any given solution $(E(t), F(t))$ of (2.5) with $(E(0), F(0))\in \Omega_1^2$, if it leaves $\Omega_1^2$, it can only cross $\overline{P_0Q_0}$. Suppose that there exists some $t_0\in(0, \bar{T})$ such that $(E(t_0), F(t_0))\in\overline{P_0 Q_0}$ and $(E(t), F(t))\in\Omega_1^2, t\in[0, t_0)$; then, $E'(t_0) < 0$, and the solution enters $\Omega_1^1$ after $t_0$. We can further prove that $(E(t), F(t))\in\Omega_1^1$ for $t\in(t_0, \bar{T}]$. If it leaves $\Omega_1^1$, it can only cross $\overline{P_1Q_1}$. Suppose that there exists some $\tilde{t}\in(t_0, \bar{T}]$ such that $(E(\tilde{t}), F(\tilde{t}))\in\overline{P_1 Q_1}$ and $(E(t), F(t))\in\Omega_1^1, t\in(t_0, \tilde{t}]$; then, $E'(t) = G_1 (E(t), F(t))\geq\alpha_1$ for $t\in[t_0, \tilde{t}]$, and $E(\tilde{t}) = E_m\geq E(t_0)+\alpha_1 (\tilde{t}-t_0)\geq E_0+\alpha_1 \bar{T}$. Since $\alpha_1 < 0$ and $\bar{T} < T_0$, we get that $E_m > E_0+\alpha_1T_0$, which contradicts with the definition of $T_0$. Thus, for any given solution $(E(t), F(t))$ of (2.5) with $(E(0), F(0))\in \Omega_1^2$, we can deduce $(E(t), F(t))\in \Omega_1$, $t\in[0, \bar{T}]$.

Next, we prove that for any solution $(E(t), F(t))$ of (3.1) with $(E(\bar{T}), F(\bar{T}))\in \Omega_1$, it follows that $(E(t), F(t))\in \Omega_1$, $t\in[\bar{T}, T]$ and $(E(T), F(T))\in\Omega_1^2$. Obviously, the vector field of (3.1) points to the inside of $\Omega_1$ and $\Omega_1^2$; then, they are both positively invariant sets and for any given solution $(E(t), F(t))$ of (3.1) with $(E(\bar{T}), F(\bar{T}))\in \Omega_1$, it follows that $(E(t), F(t))\in \Omega_1$, $t\in[\bar{T}, T]$. Since $\Omega_1^2$ is also a positively invariant set, if $(E(T), F(T))\notin\Omega_1^2$, then $(E(t), F(t))\in\Omega_1^1, t\in [\bar{T}, T]$. However, from $E'(t) = G_2 (E(t), F(t))\geq \alpha_2$, we can deduce $E(T)\geq E(\bar{T})+\alpha_2 (T-\bar{T})$ and $E_0-E_m \geq E(T)- E(\bar{T}) \geq \alpha_2 (T-\bar{T})$, which contradicts with the definition of $\delta_0$ and (3.15).

On the basis of the above discussion, it implies that for any solution of (2.7) with $(E(0), F(0))\in\Omega_1^2$, it holds that $(E(t), F(t))\in\Omega_1$ for all $t\in[0, T]$ and $(E(T), F(T))\in\Omega_1^2$. Clearly, the Poincare map $\Phi$ defined in (3.13) maps $\Omega_1^2$ into $\Omega_1^2$ and has at least one fixed point in $\Omega_1^2$ which is denoted by $\tilde{A}_0\in\Omega_1^2$. Then the solution of (2.7) through $\tilde{A}_0$ is a positive continuous $T$-periodic solution. The proof is completed. □

4.   Numerical simulation

In this section, we present several numerical simulations to substantiate the theoretical findings presented in Section 3. The majority of parameter values have been sourced from ^[30,31] as follows

We utilize the method for determining the environmental capacity $K$ found in ^[30,31], which takes into account an island of 74 hectares. After estimation, the density of male mosquitoes on this island is approximately 69 per hectare. The assumption is that in the absence of sterile mosquitoes, when the mosquito population stabilizes, the quantity of wild male and female wild mosquitoes on the island is approximately $M^* = F^* = 5106$. Consequently, we can determine the number of eggs $E^* = \frac{\delta_F F^*}{\nu \beta_F} = 40848$ at this time, as well as the approximate value of the environmental capacity

For the selection of parameter $\bar{T}$, we adopt the value from ^[13] and set

Through direct calculation, we can get

For the constant releases, i.e., $T = \bar{T}$, we select varying values for the release amount $m$ and observe the mosquito population's changes. From Figure 3(a), it is evident that when no sterile males are present in the field, meaning that $M_s(t) = m = 0$, the unique positive equilibrium $A^*(E^*, F^*)$ is globally uniformly asymptotically stable, aligning with the conclusion in Theorem 1. When sterile mosquitoes are consistently released into the field, namely $M_s(t) = m > 0$, we find that system (2.5) exhibits a bistable phenomenon when $0 < m < m^*$, and that both the extinction equilibrium $A_0$ and a positive equilibrium are locally stable (refer to Figure 3(b)). As the release amount increases to $m = m^*$, these two positive equilibrium points merge into a unique positive equilibrium, which is semi-stable (see Figure 3(c)). Conversely, if $m > m^*$, the extinction equilibrium $A_0$ becomes the only equilibrium which is globally stable (see Figure 3(d)).

For the case $T > \bar{T}$, we first let $m = 1.3\times 10^4 < m^*$; then, after simple calculation, we can get that $\lambda_0 = \lambda_1 = 0.04$ and $\lambda = 10$. Without losing generality, we take the constant $\sigma = 0.1\in(0, 1)$ and then deduce $\omega = 1.0375$ and $\omega\bar{T} = 2.5938$. Choose $T = 2.59$; then, $\bar{T} < T < \omega\bar{T}$ and the extinction equilibrium $A_0$ of the switched system (2.7) is locally stable as shown in Figure 4 (a), aligning with the conclusion in Theorem 3.3. It is worth noting that besides $A_0$, system (2.7) also exhibits a positive periodic solution which is locally stable. In other words, system (2.7) demonstrates a bistable phenomenon as well. In addition, we find that the domain of attraction of the extinction equilibrium shown in Figure 4(a) is much larger than the region determined in Theorem 3.3, which further proves that the stability conditions given in Theorem 3.3 are only sufficient and unnecessary. Under the same conditions, it can have a larger domain of attraction.

Furthermore, we fix the release amount $m = 1.3\times 10^4 < m^*$ and extend the release period to $T = 3 > \omega \bar{T}$; then, we find that the extinction equilibrium loses its stability while the positive continuous periodic solution becomes globally stable (see Figure 4(b)). According to the expressions of equilibria $A^*(E^*, F^*)$ and $A_2(E_2, F_2)$ in (3.2) and Lemma 1, we can calculate that when $m = 1.3\times 10^4$,

Then the rectangle region $\Omega_0 = [E_2, E^*]\times[F_2, F^*]$ defined in Theorem 4 is $\Omega_0 = [31676, 40848]\times[3960, 5106]$. We observe that the stable periodic solutions in Figure 4(a) and (b) do fall into the region $\Omega_0$. This further confirms the conclusion in Theorem 4.

Then we increase the release amount to $m = 2\times 10^4 > m^*$ and choose $T = 2.55$ and $T = 3$, respectively. From Figure 5(a), we can see that when $\bar{T} < T = 2.55 < \omega \bar{T}$, system (2.7) also shows the bistable phenomenon, and that there is a locally stable positive periodic solution besides $A_0$ which is locally stable. When $T = 3$, the positive periodic solution becomes globally stable (see Figure 5(b)). Compared with the case $m = 1.3\times 10^4 < m^*$ in Figure 4, we can see that the population level at the positive stable state decreases significantly due to the increase of the release amount.

In addition, if we further increase the release amount to $m = 3\times 10^4$ with $T = 2.55$, then the positive periodic solution vanishes and $A_0$ attains global stability. (see Figure 6).

Upon analyzing the numerical simulation results, we observed that the system displays diverse equilibrium stabilities, including bistability, semistability and global stability. For a mosquito population suppression model, deciphering the phenomena of bistability, semi-stability and global stability is crucial as they are directly related to the dynamic variations of the population over time. These dynamics are especially significant when considering responses to interventions targeted at curtailing or eradicating the population. The comprehension of these phenomena aids in forecasting the long-term consequences of varied intervention strategies and assessing the potential for the population to either rebound or be extinguished.

5.   Conclusions

In this paper, we developed and analyzed a new mosquito population model featuring a stage structure and periodic releases of sterile males. Certain modeling concepts from prior research were incorporated, such as only considering sexually active sterile mosquitoes in the interactive dynamic system and periodically and impulsively releasing sterile mosquitoes.

In light of the relationship of two important parameters $T$ and $\bar{T}$, our primary focus is on examining the dynamical behaviors in two distinct scenarios: $T = \bar{T}$ and $T > \bar{T}$. Regarding the first scenario, the nonautonomous system (2.4) was transformed into the two-dimensional ordinary differential equation system (2.5). The dynamical behaviors of this system were analyzed in detail (Theorem 2) and an important threshold $m^*$ for the release amount was proposed. For the scenario $T > \bar{T}$, the nonautonomous system (2.4) was transformed into the switched system (2.7). We first explored sufficient conditions under which the extinction equilibrium of system (2.7) is locally stable in Theorem 3.3, and then we discussed the existence of positive periodic solutions of (2.7) for cases $0 < m < m^*$ and $m > m^*$, respectively (Theorems 4 and 5). In addition, we conducted a range of numerical simulations that not only corroborate our theoretical findings but also supplement the results of theoretical investigations (see Figures 3–6).

A clear demonstration of bistability in system (2.7) is evident from Figures 3–5. This crucial observation underlines the system's sensitivity to initial conditions. As a result, delving deeper into the fractal structure of the basins of attraction becomes immensely important. However, considering the intricacy of the switched system that we have formulated, we are presently facing challenges in identifying a suitable method to probe the stability of the positive periodic solutions. Identifying the precise attractor regions for each stable state also remains a complex task using our current methodology. Nevertheless, in future research work, we will further attempt new methods, such as constructing Lyapunov functions to determine the stability of periodic solutions ^[32,33]. Moreover, the developmental process of mosquitoes is dictated by a variety of factors. Therefore, incorporating elements such as random effects and time delay effects into our existing system can yield a more accurate representation of mosquito population growth. Concurrently, we are continuously refining our understanding of such mathematical modeling methods ^{[34,35,36,37]}, preparing us for more extensive research in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (11901502), Program for Science $&$ Technology Innovation Talents in Universities of Henan Province (21HASTIT026) and Graduate Research Innovation Fundation of XYNU (2022KYJJ055).

Conflict of interest

The authors declare no conflict of interest regarding the publication of this paper.