Modeling and analysis of release strategies of sterile mosquitoes incorporating stage and sex structure of wild ones

1.   Introduction

Mosquito-borne infectious diseases are spread by mosquitoes and are prevalent worldwide, including malaria, dengue fever, filariasis, yellow fever, and other highly harmful infectious diseases. Research shows that there are more than 80 such diseases, many of which have caused significant losses to human lives and health in different countries and regions. In recent years, the risk of mosquito-borne infectious disease outbreaks has increased due to climate change, as well as the rise of commercial trade and personnel exchanges.

Due to the lack of an effective vaccine for mosquito-borne diseases, controlling the population of mosquitoes has become crucial in curbing the spread of such diseases. Traditional physicochemical control methods have significant disadvantages such as negative environmental effects and drug resistance. In recent years, new biological control approaches such as insect sterility technology (SIT) and its derivatives have been extensively studied and practiced. These approaches involve the release of sterile or Wolbachia-infected mosquitoes into the wild to disrupt or suppress the development of the wild mosquito population. Certain studies specifically explored the use of sterile insect release for targeted population control ^[1,2,3]. Other research considered how the endosymbiotic bacterium Wolbachia can induce resistance to the dengue virus in Aedes aegypti ^[4], and there has been research into the control of particular mosquito-borne diseases ^[5]. Furthermore, in ^[6], a combination of incompatible and sterile insect techniques has been implemented to eradicate wild mosquitoes. During the exploration of these technologies, detailed analyses have been conducted on specific population characteristics such as population diffusion, periodic changes, the Allee effect, and others ^{[7,8,9,10,11]}. To ensure successful suppression of mosquito populations, continuous release over a long period is often necessary. Researchers have studied release strategies for different situations through various methods. For instance, uninterrupted continuous release strategies, based on ODE or DDE models, have undergone extensive examination ^{[12,13,14,15]}. More recently, impulsive release strategies, which better align with practical operations, have received considerable attention. Switched systems have been widely employed in simulating these release strategies ^{[16,17,18,19,20]}. In addition, many researchers have shown a preference for impulsive differential dynamical systems ^{[21,22,23,24]}. Some have integrated the duration of sexual activity of sterile mosquitoes with impulsive release behavior to investigate release strategies that effectively suppress mosquito populations ^[25,26,27].

Typically, male sterile mosquitoes are released into the wild to mate with female wild mosquitoes, rendering them unable to lay or hatch eggs. In ^{[18,26,27,28]}, the authors proposed a mathematical model that considers the sexual lifespan of sterile mosquitoes, including only those within their sexually active period. Due to the short sexual lifespan of sterile mosquitoes, natural deaths during the sexual activity period are often neglected. Instead of using a single equation to describe the change in the number of sterile mosquitoes, the number of sterile mosquitoes is used as a control function ^{[15,18,26,27,28,29]}. Based on these ideas, the authors provided a nearly complete characterization of the release strategy using the following model.

where $\omega(t)$ and $g(t)$ represent the number of wild mosquitoes and the number of sterile male mosquitoes in their sexually active period. The birth rate per wild mosquito, denoted by $a$, follows logistic growth. Additionally, $\xi$ and $\mu$ represent the density-dependent and independent death rates, respectively.

We know that mosquito growth includes three aquatic stages, namely egg, larva, and pupa, as well as an adult stage. Intraspecific competition significantly affects the aquatic stages, while the adult stage is rarely affected. Building on model (1.1), the authors in ^[30] proposed and studied a stage-structured model for suppressing mosquito populations as follows:

where $J$ and $A$ represent the number of wild mosquitoes in the aquatic and adult stages, respectively, while $g$ represents the number of sterile males. The birth rate of adults is denoted by $\beta$. $\mu_i$ and $\xi$ represent the natural mortality and intraspecific competition coefficient, respectively, where $i = 0, 1$. The parameter $\alpha$ represents the emergence rate from larvae to adults.

As mentioned earlier, in order to suppress the wild mosquito population, male sterile mosquitoes are released in most scenarios. In addition, the role of female and male mosquitoes in interacting with sterile mosquitoes is different. Motivated by the previous research works, this paper aims to further investigate and develop a wild mosquito suppression model with both stage and sex structure.

The paper is organized as follows: In Section 2, we develop a mosquito population suppression model with both stage and sex structure. The model is a time-switched system with the number of sexually active sterile males serving as a control function. In Section 3, we first define two important thresholds for the release amount and then determine the conditions for the stability of the wild mosquito extinction equilibrium and the existence of positive periodic solutions. In Section 4, we provide numerical results to validate our theoretical findings. Finally, we conclude with a brief summary in Section 5.

2.   Model formulation

L. Almeida et al. in ^[31] describes the dynamics of the mosquito population as follows

where $E(t)$, $L(t)$, $P(t)$, $F(t)$ and $M(t)$ stands for the numbers of eggs, larvae, pupa, adult females and adult males at time $t$, respectively. The parameters $\beta_E$, $K$, and $\tau_E$ represent the oviposition rate, environmental capacities and hatching rate for eggs, respectively. Pupas will develop into males or females, and $\nu\in(0, 1)$ reflects the corresponding proportion. $\tau_L$, $\tau_P$ are transition rates, $c$ is the intraspecific competition of larvae, and $\delta_E$, $\delta_L$, $\delta_P$, $\delta_F$ and $\delta_M$ are the corresponding death rates.

To simplify the mosquito population model, we make the following assumptions, similar to those in ^[31]: (i) the dynamics of larvae and pupae are fast, and intraspecific competition at the larvae stage is negligible $(i.e., \; c\ll1)$; (ii) the probability of a pupa developing into a female is the same as that of a male $(\nu = \frac{1}{2})$; (iii) male and female mosquitoes have the same death rate $(\delta_F = \delta_M)$.

With these assumptions and the third and fourth functions of system (2.1), we obtain the equilibrium equations for larvae and pupae, which are:

After substituting Eq (2.2) into system (2.1) and introducing a release function $u(\cdot)$ for sterile male mosquitoes, we can derive the following interactive dynamical system describing the population dynamics of two kinds of mosquitoes

where $M_s(t)$ is the number of sterile males in the field, $\gamma$ measures the mating competitiveness of sterile males, and $\beta_F = \frac{\tau_P\tau_L\tau_E}{(\tau_P+\delta_P)(\tau_L+\delta_L)}$.

In this paper, we also only consider sexually active sterile mosquitoes in the dynamics. Let $M_s(t)$ be the number of sexually active sterile males at time $t$ and ignore their natural death as in previous works ^{[15,18,26,27,28,29]}. By taking $M_s(t)$ as a control function, we convert the system (2.1) into

Let $T$ be the waiting time between two consecutive release, and sterile mosquitoes are released with an amount of $m$ at time $t = kT$, $k = 0, 1, 2, \cdots$. We use $\bar{T}$ to represent the sexual lifespan of sterile mosquitoes, ignoring their natural death as in ^{[15,18,26,27,28,29]}. There are three possible cases for $T$ and $\bar{T}$: $T = \bar{T}$, $T > \bar{T}$, and $T < \bar{T}$. Previous works ^[30] studied the effect of releasing sterile mosquitoes on wild mosquito populations in the first two cases.

In this work, we mainly study the last case: $T < \bar{T}$. In fact, experimental studies have indicated that the lifespan of sterile mosquitoes can be as long as 51 days ^[32] such as the $Wolbachia$-infected mosquitoes, and their sexual lifespan can reach 14 days ^[33]. Nonetheless, the release of sterile mosquitoes is relatively frequent, such as in the experimental field of Shazai Island in Guangzhou, China, where they were released three times a week ^[32]. Hence, it is meaningful to consider the situation where the sexual lifespan of sterile mosquitoes exceeds the release period, that is $T < \bar{T}$.

Similar to the discussion in ^[17], when $T < \bar{T}$, there must exist a unique positive integer $p$ and a non-negative number $q$ that satisfy

where $q\in[0, T)$ and $p = [\bar{T}/T]$, which is the integer closest to $\bar{T}/T$ and not greater than $\bar{T}/T$. Since $pT\leq\bar{T}$, we can obtain the first $p$ releases

where $k = 0, 1, 2, \cdots, p-1$. From the $(p+1)$th release, it is easy to get

where $k = p, p+1, \cdots$. If $q = 0$, then $M_s(t) = pm$, and the characteristics of the system are completely consistent with the simple system of constant value continuous release. In the present work, we assume that $0 < q < T$. As we are studying the asymptotic behavior of the system, for the sake of analysis convenience and without loss of generality, we can further assume that for $k = 0, 1, 2, \cdots$,

Based on the above assumption, the model (2.4) is transformed as follows:

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

Population dynamics of wild mosquitoes with sterile mosquito release then depend on the continuous switching between systems (2.7) and (2.8). In this paper, based on the switched systems (2.7) and (2.8), we mainly investigate the extinction and periodic change of wild mosquito population.

3.   Dynamics of the switched system

3.1.   preliminaries

First, we consider the dynamic behavior of the wild mosquito population without the release of sterile mosquitoes. Let $M_s(t)\equiv0$, then the system (2.4) or (2.7) and (2.8) can be rewritten as follows:

Let $\Omega = \{(E, F)|0\leq E\leq K, 0\leq F\}$. Obviously, $\Omega$ is a positive invariant set of systems (3.1), (2.7) and (2.8), and every solution of systems (3.1), (2.7) and (2.8) with nonnegative initial values is positive and bounded. Therefore, in this work, we will investigate the dynamics of these systems in $\Omega$. Additionally, it is easy to verify that system (3.1) is a monotonic system on $\Omega$, so there are no closed orbits. By defining the basic offspring number of the wild mosquito as

we list the following conclusions about the system (3.1):

Lemma 1. $(i)$ If $R_0 \leq1$, then the extinction equilibrium $A_0(0, 0)$ is the unique equilibrium of system (3.1) that is globally asymptotically stable in $\Omega$.

$(ii)$ If $R_0 > 1$, then $A_0(0, 0)$ is unstable, and the unique positive equilibrium $A^*(E^*, F^*)$ of system (3.1) is globally asymptotically stable in $\Omega$, where

According to the necessity of releasing sterile mosquitoes, we will only consider the case $R_0 > 1$ in this paper.

For the convenience of illustration, the following system is introduced to uniformly express the two subsystems (2.7) and (2.8)

where $\tilde{p} = p+1$ or $p$.

To compute the jacobian matrix of the system (3.2), we obtain

It has non-negative extra-diagonal coefficients on $\Omega = \{(E, F)|0\leq E\leq K, 0\leq F\}$, indicating that the system (3.2) is monotone on $\Omega$.

Clearly, the wild mosquito extinction equilibrium $A_0(0, 0)$ is an equilibrium of both systems (2.7) and (2.8). By using the general form (3.2) of systems (2.7) and (2.8), we can now calculate the positive equilibria by finding the positive roots of the following equations,

which can be further equivalently transformed into

It is easy to know that the discriminant of the above equation is

then we define the following release threshold for the system (3.2)

Lemma 2. Let $R_0 > 1$. The trivial equilibrium $A_0$ of system (3.2) is always locally asymptotically stable. Moreover,

(i) There exists $m^*(\tilde{p}) > 0$ such that the system (3.2) admits no positive equilibrium if $m > m^*(\tilde{p})$, one positive equilibrium if $m = m^*(\tilde{p})$, and two positive equilibria $A_{\cdot j}(E_{\cdot j}(\tilde{p}), F_{\cdot j}(\tilde{p})), j = 1, 2$ if $0 < m < m^*(\tilde{p})$, where 5

(ii) If $m > m^*(\tilde{p})$, then the trivial equilibrium $A_0$ is globally asymptotically stable in $\Omega$. If $0 < m < m^*(\tilde{p})$, then $A_{\cdot 1}(E_{\cdot 1}(\tilde{p}), F_{\cdot 1}(\tilde{p}))$ is unstable, while $A_{\cdot 2}(E_{\cdot 2}(\tilde{p}), F_{\cdot 2}(\tilde{p}))$ is locally asymptotically stable. If $m = m^*(\tilde{p})$, then the unique positive equilibrium is a saddle node.

Proof. To calculate the Jacobian matrix of the system (3.2) at the trivial equilibrium $A_0(0, 0)$, we have

It has two negative eigenvalues, therefore the trivial equilibrium $A_0$ is always locally asymptotically stable.

(i) By calculating the discriminant $\Delta$ in Eq (3.5), we can see that $\Delta = 0$ when $m = m^*(\tilde{p})$, $\Delta < 0$ when $m > m^*(\tilde{p})$, and $\Delta > 0$ when $0 < m < m^*(\tilde{p})$. Moreover, the root of the quadratic equation (3.4) must be positive roots if it exists. These observations lead to the conclusions in (i).

(ii) Since the system (3.2) is monotone on $\Omega$, it has no closed orbits. If $m > m^*(\tilde{p})$, then the extinction equilibrium $A_0$ is the unique equilibrium which is stable, so it is globally asymptotically stable in $\Omega$.

If $0 < m < m^*(\tilde{p})$, we calculate Jacobian matrix of the system (3.2) at $A_{\cdot j}, j = 1, 2$ by using jacobian matrix (3.3) and get

We can easily get $tr J_{A_{\cdot j}} = \frac{-\nu\beta_E\beta_FF_{\cdot j}}{\delta_F(m_1+F_{\cdot j})}-\delta_F < 0$ and

where $\Phi(F) = \frac{-G(F)}{\gamma \tilde{p}m+F}.$

Since $-\Phi'(F) = -\frac{G(F)}{(\gamma \tilde{p}m+F)^2}+\frac{G'(F)}{(\gamma \tilde{p}m+F)}$, $G(F_{\cdot 1}) = G(F_{\cdot 2}) = 0$ and $G'(F_{\cdot 1}) > 0, G(F_{\cdot 2}) < 0$, we can obtain $|J(A_{\cdot 1})| < 0$ and $|J(A_{\cdot 2})| > 0$. Therefore, $A_{\cdot 1}(E_{\cdot 1}, F_{\cdot 1})$ is an unstable saddle, while $A_{\cdot 2}(E_{\cdot 2}, F_{\cdot 2})$ is locally asymptotically stable.

If the release amount $m$ increases to $m = m^*(\tilde{p})$, then the two positive equilibria $A_{\cdot 1}$ and $A_{\cdot 2}$ merge into a single semi-stable saddle node. The proof is completed.

Remark 1. According to Lemma 2, the constant continuous release system (3.2) has a threshold $m(\tilde{p})$ for the release amount $m$. When the release amount $m > m(\tilde{p})$, the wild mosquito population in the field will be successfully suppressed.

Referring to the two subsystems (2.7) and (2.8), we denote their corresponding release thresholds as $m_1^*$ and $m_2^*$, respectively, given by:

When there are two positive equilibria, the positive equilibria corresponding to systems (2.7) and (2.8) are denoted by $A_{1j}(E_{1 j}, F_{1 j})$ and $A_{2j}(E_{2 j}, F_{2 j}), j = 1, 2$, respectively, where

Lemma 2 provides insights into the behavior of systems (2.7) and (2.8). We can now state the following result:

Lemma 3. Let $R_0 > 1$. The trivial equilibrium $A_0$ of system (2.7) (or (2.8)) is always locally asymptotically stable. In addition,

$(i)$ If $0 < m < m_1^*$, system (2.7) (or (2.8)) has two positive equilibria $A_{11}(E_{11}, F_{11})$ and $A_{12}(E_{12}, F_{12})$ (or $A_{21}(E_{21}, F_{21})$ and $A_{22}(E_{22}, F_{22})$), where $A_{11}(E_{11}, F_{11})$ (or $A_{21}(E_{21}, F_{21})$) is unstable, and $A_{12}(E_{12}, F_{12})$ (or $A_{22}(E_{22}, F_{22})$) is locally asymptotically stable.

$(ii)$ If $m_1^* < m < m_2^*$, system (2.8) has two positive equilibria $A_{21}(E_{21}, F_{21})$ and $A_{22}(E_{22}, F_{22})$, where $A_{21}(E_{21}, F_{21})$ is unstable, and $A_{22}(E_{22}, F_{22})$ is locally asymptotically stable. However, system (2.7) has no positive equilibrium.

$(iii)$ If $m > m_2^*$, both systems (2.7) and (2.8) have no positive equilibrium, and the unique equilibrium $A_0(0, 0)$ is globally asymptotically stable in $\Omega$.

3.2.   Periodic and impulsive releases

When the release period of sterile mosquitoes is shorter than their sexual lifespan, that is, $T < \bar{T}$, the number of sterile mosquitoes that are sexual active will continually switch between two levels, as shown in function (2.6). The population dynamics of wild ones will then depend on the switched systems (2.7) and (2.8).

With the appearance of continual switching, the extinction equilibrium $A_0(0, 0)$ becomes the unique equilibrium of the systems (2.7) and (2.8). We will first discuss its stability.

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

then for any given $0 < \sigma < 1$, denote

Theorem 1. Assuming $R_0 > 1$ and $\bar{T} = pT+q$, for any solution $(E(t), F(t))$ of systems (2.7) and (2.8) with $0 < E(0)+F(0) < \ell_1$, the inequality

holds, where $r = (1-\sigma)\alpha > 0$. Thus, the extinction equilibrium $A_0(0, 0)$ of systems (2.7) and (2.8) is exponentially asymptotically stable.

Proof. Suppose $(E(t), F(t))$ is a solution of the systems (2.7) and (2.8) in $\Omega$ with $0 < E(0)+F(0) < \ell_1$ for fixed $\sigma\in(0, 1)$. Let $P(t): = E(t)+F(t)$ for $t\geq0$. According to the first subsystem (2.7), we have for $t\in[kT, kT+q]$, $k = 0, 1, 2, \cdots$,

We claim that if $P(kT) < \ell_1 < \ell$, then $P(t) < \ell_1$ for all $t\in [kT, kT+q ]$, $k = 0, 1, 2, \cdots$. Let's prove this by contradiction. Suppose it does not holds, then there must exist $\bar{t}\in (kT, kT+q ]$ such that $P(t) < \ell_1 < \ell$ for $t\in[kT, \bar{t})$ and $P(\bar{t}) = \ell_1$, then $P'(\bar{t})\geq 0$. However, from inequality (Eq 3.10), we deduce

which leads a contradiction. Therefore, we have $P(t) < \ell_1$ for all $t\in [kT, kT+q ]$ and $P'(t)\leq (\frac{\beta_E}{\gamma(p+1)m}\ell_1-\alpha)P(t) = -(1-\frac{p}{p+1}\sigma)\alpha P(t) < 0$. Then we can further deduce

Similarly, for the second subsystem (2.8), we can obtain

Then if $P(kT+q) < \ell_1$, then $P(t) < \ell_1$ for all $t\in [kT+q, (k+1)T]$, $k = 0, 1, 2, \cdots$. In addition,

Combining with inequalities (Eq 3.11) and (Eq 3.12), we know that

where $\xi: = e^{-(1-\sigma)\alpha T} < 1$.

Thus,

If $t\in[kT, kT+q]$, then from inequality (Eq 3.11) we get

If $t\in[kT+\bar{T}, (k+1)T]$, then from inequality (Eq 3.12) we get

By combining inequalities (Eq 3.13) and (Eq 3.14), we can see that the inequality (Eq 3.9) holds and $A_0(0, 0)$ is exponentially asymptotically stable. The proof is completed.

In Theorem 1, we show the exponential asymptotic stability of $A_0(0, 0)$. However, this stability is only a local feature. Due to the periodic release of sterile mosquitoes, we are concerned about the possibility of a periodic change in the wild mosquito population, specifically, whether the systems (2.7) and (2.8) has a positive periodic solution. In the following discussion, we focus on the existence of positive periodic solutions.

Based on Lemma 3, when $R_0 > 1$, $T < \bar{T}$, and $0 < m < m_1^*$, both subsystems (2.7) and (2.8) have a pair of positive equilibria, namely $A_{11}(E_{11}, F_{11})$, $A_{12}(E_{12}, F_{12})$ and $A_{21}(E_{21}, F_{21})$, $A_{22}(E_{22}, F_{22})$. It is easy to verify that $F_{22} > F_{12} > F_{11} > F_{21} > 0$ and $K > E_{22} > E_{12} > E_{11} > E_{21} > 0$. Let $\Omega_0: = [E_{12}, E_{22}]\times[F_{12}, F_{22}]$ be a rectangle in the $EF$-phase plane. Then we can get the following result.

Theorem 2. If $R_0 > 1$ and $0 < m < m_1^*$, then the systems (2.7) and (2.8) has a $T$-periodic solution in $\Omega_0$. Suppose $(E(t), F(t))$ is a solution of the systems (2.7) and (2.8) with $E(t_0)\geq E_{11}$ and $F(t_0)\geq F_{11}$ for some $t_0\geq0$, then $E(t)\geq E_{11}$ and $F(t)\geq F_{11}$ 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 $T$-periodic, then $(E(t), F(t))\in\Omega_0$ for all $t\geq0$.

Proof. For the conditions $R_0 > 1$, $T < \bar{T}$, and $0 < m < m_1^*$, we obtain the existence of $A_{ij}$, $i, j = 1, 2$. Let $\eta(E): = \frac{KE(\tau_E+\tau_E)}{K-E}$, $E\in[0, K]$. The vertical isoclines $E'(t) = 0$ of systems (2.7) and (2.8) in the first quadrant of the $EF$-phase plane are given by

respectively (see Figure 1). Since $L_1(E) > L_2(E)$, the curve of $F = L_1(E)$ is strictly above the curve of $F = L_2(E)$ for $E > 0$. Additionally, $F = G(E): = \frac{\nu\beta_EE}{\delta_F}$ is the shared horizontal isocline of systems (2.7) and (2.8).

Select any two positive numbers $E_m$ and $E_M$ satisfying $E_{11} < E_m < E_{12} < E_{22} < E_M < K$, then $F_m = G(E_m)$ and $F_M = G(E_M)$. Due to the strict monotonic increase property with respect to $E$ of $G(E)$, it follows that $F_m < F_M$. 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)$.

Clearly, $F_{11} < F_m < F_{12} < F_{22} < F_M$ and $\Omega_0\subset\Omega_1$. Since $G(E)$ is strictly monotone increasing with respect to $E$, the segment $\overline{P_1P_2}$ lies strictly below the isocline $F = G(E)$ except for the point $P_1$, while the segment $\overline{Q_1Q_2}$ lies strictly above the isocline $F = G(E)$ except for the point $Q_2$. According to the second equation of systems (2.7) and (2.8), we have

Thus, the vector fields of (2.7) and (2.8) point into the closed rectangle $\Omega_1$ on $\overline{P_1P_2}$ and $\overline{Q_1Q_2}$.

In addition, since $L_1(E)$ and $L_2(E)$ are both strictly monotone increasing with respect to $E$ on $[0, K]$, we can see that the entire segment $\overline{P_1Q_1}$ lies strictly above the isocline $F = L_1(E)$, and the entire segment $\overline{P_2Q_2}$ lies strictly below the isocline $F = L_2(E)$. From system (2.7), we obtain

where $L_1^{-}(E) = \frac{1}{2\beta_E}[-\eta(E)+\sqrt{\eta(E)^{2}+4\gamma (p+1)m\beta_E\eta(E)}] > 0$.

From the system (2.8), we can obtain

where $L_2^{-}(E) = \frac{1}{2\beta_E}[-\eta(E)+\sqrt{\eta(E)^{2}+4\gamma pm\beta_{E}\eta(E)}] > 0$.

Clearly, on the segment $\overline{P_1Q_1}$, we have $E\equiv E_m$, so $F > L_1(E_m)$ and $F > L_2(E_m)$. As a result, $\frac{dE}{dt} > 0$ holds for both systems (2.7) and (2.8) on $\overline{P_1Q_1}$. Similarly, on the segment $\overline{P_2Q_2}$, we have $F < L_1(E_M)$ and $F < L_2(E_M)$, leading to $\frac{dE}{dt} < 0$ for both systems (2.7) and (2.8) on this segment. Therefore, the vector fields of systems (2.7) and (2.8) point towards the inside of $\Omega_1$ on both $\overline{P_1Q_1}$ and $\overline{P_2Q_2}$.

In summary, $\Omega_1$ represents a positive invariant set of the switched systems (2.7) and (2.8), and any point $A\in\Omega_1$ will have the trajectory of the systems (2.7) and (2.8) starting from it stay within $\Omega_1$.

We define a map $\Phi:\Omega_1\rightarrow \Omega_1$ as given in ^[30]:

Here, $\phi_t(A)$ and $\psi_t(A)$ are the solutions of systems (2.7) and (2.8) starting from the point $A$, respectively. Since $\Phi$ is continuous and maps $\Omega_1$ to $\Omega_1$, by Brouwer's fixed point theorem, we know that $\Phi$ has a fixed point $\bar{A}\in\Omega_1$. Then, the solution of systems (2.7) and (2.8) passing through the point $\bar{A}$ is a continuous $T$-periodic solution whose trajectory lies entirely in $\Omega_1$.

It is worth noting the arbitrariness of $E_m$ and $E_M$. If $E_m\rightarrow E_{12}$ and $E_M\rightarrow E_{22}$, then $\Omega_1\rightarrow \Omega_0$, and the continuous $T$-periodic solution of systems (2.7) and (2.8) in $\Omega_1$ lies in $\Omega_0$. If $E_m\rightarrow E_{11}$ and $E_M\rightarrow \infty$, then $\Omega_1\rightarrow [E_{11}, \infty)\times[F_{11}, \infty)$. Since $\Omega_1$ is an attraction domain, $[E_{11}, \infty)\times[F_{11}, \infty)$ is also a positively invariant set of systems (2.7) and (2.8). Then, for any solution of systems (2.7) and (2.8) with $E(t_0)\geq E_{11}$ and $F(t_0)\geq F_{11}$, it follows that $E(t)\geq E_{11}$ and $F(t)\geq F_{11}$ for all $t > t_0$. Moreover, $dist((E(t), F(t)), \Omega_0)\rightarrow 0$ as $t\rightarrow \infty$. Furthermore, if this solution is periodic, then $(E(t), F(t))\in\Omega_0$. Thus, the proof is completed.

When we increase the release amount $m$ to exceed the smaller threshold $m_1^*$, that is, $m_1^*\leq m < m_2^*$, we obtain the following result:

Theorem 3. Suppose $R_0 > 1$ and $m_1^* < m < m_2^*$. Then there exist positive constants $T_0 > 0$ and $\delta_0 > 0$ such that the systems (2.7) and (2.8) has a continuous and positive $T$-periodic solution when

Proof. According to Lemma 3, when $R_0 > 1$, $T < \bar{T}$ and $m_1^* < m < m_2^*$, the subsystem (2.8) has two positive equilibria $A_{21}$ and $A_{22}$, while (2.7) has only the extinction equilibrium.

Similar to Theorem 2, we first select two arbitrary positive numbers $E_m$ and $E_M$ such that $0 < E_m < E_{22} < E_M$ and $G(E_M) < L_1(E_m)$ ($G(E)$ and $L_i(E), i = 1, 2$ are defined in the proof of Theorem 2). Let $F_m = G(E_m)$ and $F_M = G(E_M)$. Then we construct a closed rectangle $\Omega_1$ with vertices at $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, let $P_3$ be the intersection point of the segment $P_1P_2$ and the isocline $F = L_2(E)$, and for any $E_0\in(E_m, E_3)$, denote points $P_0: = (E_0, F_m)$ and $Q_0: = (E_0, F_M)$ (see Figure 2). Then $\Omega_1$ can be regarded as a rectangle composed of two parts, namely, $\Omega_1 = \Omega_1^1\bigcup \Omega_1^2$, where $\Omega_1^1$ is the rectangle with $P_1$, $P_0$, $Q_1$, and $Q_0$ as its vertices, and $\Omega_1^2$ is the rectangle with $P_0$, $P_2$, $Q_0$, and $Q_2$ as its vertices.

We define $\omega_1$ and $\omega_2$ as follows:

where $G_1(E, F)$ and $G_2(E, F)$ are the right-hand side functions of the first equations of systems (2.7) and (2.8), respectively.

We assume that

where $T_0: = \frac{E_0-E_m}{-\omega_1}$, $\delta_0: = \frac{E_0-E_m}{\omega_2}$.

In the following, we will prove that the Poincare map $\Phi$ defined in Theorem 2 is continuous and maps $\Omega_1^2$ to itself. Then, similarly to the proof of Theorem 2, we can deduce that the systems (2.7) and (2.8) has a continuous and positive $T$-periodic solution in $\Omega_1^2$.

The vector fields of system (2.7) point towards the interior of $\Omega_1$ along the segments $\overline{P_1P_2}$, $\overline{Q_1Q_2}$, and $\overline{P_2Q_2}$. Therefore, for any solution of system (2.7) with initial value $(E(0), F(0))\in\Omega_1^2$, it can only cross the segment $\overline{P_0Q_0}$ if it leaves $\Omega_1$. We will now show that this is impossible by contradiction.

Suppose there exists a time $t_0\in[0, q)$ such that $(E(t_0), F(t_0))\in\overline{P_0Q_0}$ and $(E(t), F(t))\in\Omega_1^2$ for $t\in[0, t_0)$. Then, we have $E'(t_0) < 0$, which means that the solution immediately enters $\Omega_1^1$ after $t_0$. When the solution leaves $\Omega_1^1$, it can only go through the segment $\overline{P_1Q_1}$.

Now, assume there exists a time $\tilde{t}\in(t_0, q]$ such that $(E(\tilde{t}), F(\tilde{t}))\in\overline{P_1Q_1}$ and $(E(t), F(t))\in\Omega_1^1$ for $t\in[t_0, \tilde{t})$. Then, we have $E'(t) = G_1(E(t), F(t))\geq \omega_1$ for $t\in[t_0, \tilde{t}]$, and $E(\tilde{t}) = E_m\geq E(t_0)+\omega_1(\tilde{t}-t_0)\geq E_0+\omega_1q$. Since $\omega_1 < 0$ and $q < T_0$, we have $E_m < E_0+\omega_1T_0$, which is a contradiction. Therefore, any solution $(E(t), F(t))$ of system (2.7) with $(E(0), F(0))\in\Omega_1^2$ will remain in $\Omega_1$ for $t\in[0, q]$.

Next, we show that for any solution of system (2.8) with $(E(q), F(q))\in\Omega_1$, the solution will stay in $\Omega_1$ for $t\in[q, T]$ with $(E(T), F(T))\in\Omega_1^2$. Noting that the vector fields of system (2.8) point to the inside of both $\Omega_1$ and $\Omega_1^2$, we see that both the rectangles $\Omega_1$ and $\Omega_1^2$ are positive invariant sets of system (2.8). Therefore, $(E(q), F(q))\in\Omega_1$ implies $(E(t), F(t))\in\Omega_1$ for all $t\in[q, T]$. Furthermore, we can claim that $(E(T), F(T))\in\Omega_1^2$. If this does not hold, then $(E(t), F(t))\in\Omega_1^1$ for all $t\in[q, T]$ due to the properties of the invariant set $\Omega_1^2$. This implies that $E'(t)\geq\omega_2$ for $t\in [q, T]$ and $E(T)\geq E(q)+\omega_2(T-q)$. We further deduce that $E_0-E_m\geq E(T)- E(q)\geq\omega_2(T-q)$, which contradicts the inequality (3.18).

In summary, we have shown that for any solution $(E(t), F(t))$ of systems (2.7) and (2.8) starting from the initial point $(E(0), F(0))\in\Omega_1^2$, the solution will stay in $\Omega_1$ for $t\in[0, T]$ and $(E(T), F(T))\in\Omega_1^2$. Thus, the Poincare map $\Phi:\Omega_1^2\rightarrow \Omega_1^2$ is well-defined and has at least one fixed point $A^*$ in $\Omega_1^2$. Then the solution of systems (2.7) and (2.8) through the point $A^*$ is a continuous and positive $T$-periodic solution. The proof is completed.

Remark 2. If $m = m_1^*$, then the two equilibria $A_{11}(E_{11}, F_{11})$ and $A_{12}(E_{12}, F_{12})$ of (2.7) coincide into one, denoted by $A_{1c}^*(E_{1c}^*, F_{1c}^*)$. By slightly modifying the proof of Theorem 3.3 and selecting values of $E_m$ and $E_M$ such that $E_{1c}^* < E_m < E_{22} < E_M$ and $G(E_M) < L_1(E_m)$, we can similarly prove that the systems (2.7) and (2.8) has a positive $T$-periodic solution.

Next, we consider the case where the release amount $m$ is further increased until $m\geq m_2^*$, and we present the following result.

Theorem 4. If $R_0 > 1$ and $m > m_2^*$, then the unique extinction equilibrium $A_0$ of systems (2.7) and (2.8) is globally asymptotically stable.

Proof. According to Theorem 1, if $R_0 > 1$ and $m > m_2^*$, the extinction equilibrium $A_0(0, 0)$ of systems (2.7) and (2.8) is locally asymptotically stable. We will now prove that it is also globally attractive.

To this end, we construct the following system:

Similar to the discussion about system (3.1), we can show that system (3.19) is monotone on $\Omega$. Furthermore, we have

Thus, system (3.19) can be regarded as a comparison system for systems (2.7) and (2.8). With the same initial values $\tilde{E}(0) = E(0)$ and $\tilde{F}(0) = F(0)$, we obtain

From Lemma 2, we know that if $R_0 > 1$ and $m > m_2^*$, the trivial equilibrium $(0, 0)$ of system (3.19) is globally asymptotically stable. That is, for any non-negative initial value $(\tilde{E}(0), \tilde{F}(0))$, the solution of system (3.19) from it satisfies

By the comparison theorem, we can easily see that for any non-negative initial value $(E(0), F(0))$, the solution of systems (2.7) and (2.8) from it satisfies

Thus, the extinction equilibrium $A_0(0, 0)$ of systems (2.7) and (2.8) is also globally attractive. The proof is completed.

Remark 3. If $m = m_2^*$, then the extinction equilibrium $A_0(0, 0)$ of system (2.7) is globally asymptotically stable, meaning that every solution of system (2.7) will tend to the extinction equilibrium $A_0(0, 0)$. Moreover, the two equilibria $A_{21}(E_{21}, F_{21})$ and $A_{22}(E_{22}, F_{22})$ of system (2.8) coincide into one, denoted by $A_{2c}^*(E_{2c}^*, F_{2c}^*)$. According to Lemma 2, this case presents a possible semi-stability, where the basin of attraction of $A_0(0, 0)$ contains the interval $[0, A_{2c}^*)$, denoted as $\{(E, F)\in R_2^+:0\leq E < E_{2c}^*, 0\leq F < F_{2c}^*\}$, and the basin of attraction of $A_{2c}^*(E_{2c}^*, F_{2c}^*)$ contains the interval $(A_{2c}^*, \infty)$, denoted as $\{(E, F)\in R_2^+: E > E_{2c}^*, F > F_{2c}^*\}$. Due to the continual switching between two subsystems (2.7) and (2.8), any solution of systems (2.7) and (2.8) must enter $[0, A_{2c}^*)$ after some time and will eventually tend to the extinction equilibrium $A_0(0, 0)$. Therefore, the extinction equilibrium $A_0(0, 0)$ of systems (2.7) and (2.8) is globally asymptotically stable.

Remark 4. According to Theorems 1–4 and Remark 3, the periodic and impulsive release systems (2.7) and (2.8) has a threshold $m_2^*$ for the release amount $m$, and only when the release amount $m\geq m_2^*$, can the wild mosquito population in the field be successfully suppressed.

4.   Numerical simulation

This section aims to illustrate the theoretical results developed in the previous section through several numerical examples. We use the values of most model parameters from ^[31,34], which are summarized in Table 1.

Another important parameter to be determined is the carrying capacity, denoted by $K$. We consider an island with an area of 74 hectares, as in ^[31,34], with an estimated male population of about 69 males per hectare. If there are no sterile mosquitoes in the field and the wild population reaches a stable state, the number of adult wild males is assumed to be $M^* = F^* = 69\times 74 = 5106$. Then, the number of eggs at the stable state is estimated as $E^* = \frac{\delta_F F^*}{\nu \beta_F} = 40848$, and the carrying capacity can be calculated as

For the sexual lifespan of sterile mosquitoes, we take the value from ^[17], and set $\bar{T} = 14$. Through direct calculation, we can get

When no sterile ones are released, the unique positive equilibrium $A^*(E^*, F^*)$ is globally and uniformly asymptotically stable (see Figure 3(a)), which is consistent with the result in Lemma 1.

For the constant release system (3.2), we can take $\tilde{p} = 4$ without loss of generality, yielding $m^*(4)\approx 4667$. When $0 < m < m^*(4)$, the system (3.2) displays bistability, as shown in Figure 3(b), where the extinction equilibrium $A_0$ and a positive equilibrium are locally stable. Increasing the release amount to $m = m^*(4)$ results in a unique positive equilibrium that is semi-stable (see Figure 3(c)). However, if $m > m^*(4)$, then the unique equilibrium $A_0$ becomes globally stable (see Figure 3(d)).

For the switched systems (2.7) and (2.8), textcolor{red}{we choose the release period $T = 3$ }and set $\bar{T} = pT+q$ with $p = 4$ and $q = 2$. After a simple calculation, we obtain $m_1^* = m^*(p+1)\approx 3735$ and $m_2^* = m^*(p)\approx 4667$.

Firstly, let $m = 3250 < m_1^*$ and we note that the systems (2.7) and (2.8) displays bistability. In addition to the locally stable extinction equilibrium, there also exists a stable positive periodic solution (see Figure 4 (a)). This is consistent with the conclusions of Theorems 1 and 2. As shown in Figures 4(b, c), when the release amount increases to $m_1^* < m < m_2^*$, the systems (2.7) and (2.8) continues to exhibit bistability, which is consistent with the conclusions of Theorems 1 and 3. Furthermore, we can see that the attraction region of the extinction equilibrium increases with the release amount $m$ until the release amount is not less than $m_2^*$. If $m < m_2^*$, the extinction equilibrium becomes globally attractive (see Figure 4(d)).

5.   Conclusions

In this work, we constructed and analyzed a two-dimensional switched system of the wild mosquito population with stage structure and periodic releases of sterile mosquitoes. To consider the role of sterile mosquitoes in the interaction between the two mosquito populations, we referred to some previous studies and only included sexually active sterile males in the system. Existing research data suggests that the sexual lifespan $\bar{T}$ of sterile mosquitoes can be shorter or longer than their release period $T$, and in ^[30], the authors studied the case $\bar{T}\leq T$. Therefore, we focused on investigating the case $\bar{T} > T$ based on the switched systems (2.7) and (2.8).

We first investigated the non-release system (3.1) and the general form (3.2) of the two subsystems. We gave conditions for the existence and stability of their equilibria (Lemmas 2 and 3), which prepared for the later study of the switched system. Two important release amount thresholds ($m_1^*$ and $m_2^*$) were determined. Then, we proved that the extinction equilibrium $A_0(0, 0)$ of systems (2.7) and (2.8) must be exponentially asymptotically stable (Theorem 1). This differs from the conclusion in ^[30] when $\bar{T} < T$ since, in the latter case, the extinction equilibrium may be unstable. As systems (2.7) and (2.8) has no positive equilibrium, we further discussed the existence of positive periodic solutions. In Theorems 2 and 3, through the analysis of vector fields and using the fixed-point theorem, we show that systems (2.7) and (2.8) has a positive periodic solution when $0 < m < m_1^*$ and $m_1^*\leq m < m_2^*$. Further increasing the release of sterile mosquitoes, and when $m \geq m_2^*$, we confirmed that the extinction equilibrium $A_0(0, 0)$ of systems (2.7) and (2.8) becomes globally asymptotically stable. This is a previously unconfirmed dynamic phenomenon when studying the case of $\bar{T} < T$.

Regarding the periodic switched systems (2.7) and (2.8), due to the complexity of its model construction, we only prove the existence of positive periodic solutions using the fixed-point theorem. However, neither the number nor the stability of the positive periodic solutions are determined. In future research, we will try more methods to solve these problems, such as constructing Lyapunov functions to study the stability of periodic solutions ^[35,36,37]. In addition, we mainly focus on the theoretical study and the validation of the related results in this work. For the numerical analysis, we used data from a mosquito population on a 74-hectare island. Due to the relatively small number of wild mosquitoes on the island, the amplitude of the periodic fluctuation in the mosquito population is indeed small. In future research, we will collect and organize relevant data from other scenarios with large mosquito populations, and apply our research results to the development of management strategies for wild mosquito populations.

Acknowledgements

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

Conflict of interest

The authors declare no conflicts of interest in this paper.