Rich and complex dynamics of a time-switched differential equation model for wild mosquito population suppression with Ricker-type density-dependent survival probability

1.   Introduction

Dengue, mainly transmitted by Aedes aegypti and Aedes albopictus, is one of the top ten global health threats, causing more than $390$ million people to be infected annually in the tropical and subtropical world ^[1,2,3]. In China, the dengue epidemic first occurred in 1978 Guangdong, a coastal province with a mild climate ^[4]. Since then, the Chinese people in the warm temperate zone have been suffering from the disease, and the people in Guangdong are the most afflicted. In $2014$, there were $45,053$ cases reported in Guangdong. From $2017$ to $2021$, $34,036$ cases were recorded, among which $10,594$ cases were identified in Guangdong ^[5,6]. As there are neither specific drugs nor safety vaccines available, the most appropriate method to prevent dengue outbreak is to suppress the density of vector mosquitoes such that the level of the dengue viruses present in mosquitoes is below the disease epidemic proportion ^[7,8,9].

For wild mosquitoes, figuring out their intrinsic population growth dynamics contributes a lot to the control of their density. Mosquitoes undergo four developmental stages during a lifetime: egg, larva, pupa and adult. The density dependence mainly occurs in the first three stages, called the aquatic stages ^[10]. Keeping these facts in mind and inspired by ^[11], we formulate the ordinary differential equation model

with Ricker-type density-dependent survival probability ^[12,13,14] to depict the population dynamics of wild mosquitoes. In Eq (1.1), $w = w(t)$ denotes the number of wild mosquitoes at time $t$, $a > 0$ represents the per capita daily egg production rate, $\left.{1}\middle/{b}\right. > 0$ estimates the size at which the population reproduces at its maximum rate ^[15,16], and $\mu > 0$ describes the density-independent death rate of wild mosquitoes. We assume $a > \mu$ in this work such that (1.1) has a unique positive equilibrium $e^* = \left.{\ln(\left.{a}\middle/{\mu}\right.)}\middle/{b}\right.$ besides the trivial equilibrium $e_0 = 0$, and a study of the direction field associated with (1.1) shows the instability of $e_0$ and the global asymptotic stability of $e^*$ ^[17,18].

Among the tools for combating wild mosquitoes, the incompatible insect technique (IIT) has proven to be an effective, eco-friendly and scalable tactic ^[19,20,21]. IIT relies on releases of the obligate intracellular bacteria Wolbachia-infected male mosquitoes (we refer to them as W-males hereafter) reared in laboratories or mosquito factories to sterilize wild females and thus to suppress wild mosquito population. The feasibility of IIT is due to Wolbachia in each released male mosquito, which not only induces the cytoplasmic incompatibility (CI) mechanism, as the eggs laid by a wild female mosquito that mated with an infected male mosquito will hatch partially, but also brings fitness alterations to the hosts ^[22,23]. Based on ^[24,25], we assume that Wolbachia induces complete CI, which means that no viable offspring will be produced when a wild female mosquito mates with an infected male mosquito, and it causes no apparent changes to the mating competitiveness of the infected male mosquito compared with the wild one.

Since the sole role that the W-males play in the suppression dynamics is to sterilize wild females through matings, Yu in ^[26] proposed an innovative modeling idea that only those W-males being sexually active are included. Therefore, due to the releases of W-males, the compatible mating probability between wild and W-males in Eq (1.1) decreases from $1$ to $\left.{w}\middle/{(w+g)}\right.$, where $g = g(t)$ is the number of sexually active W-males at time $t$, and the population dynamics of wild mosquitoes thus can be governed by

In addition, for W-males, their sexual lifespan is shorter than their lifespan ^[10]. Thus, the death of those W-males can be neglected. Consequently, $g(t)$ can be treated as a given function in advance, and it is determined by the release strategy implemented. Meanwhile, the author in ^[26] proposed a periodic and impulsive release strategy that a constant amount of $c$ of W-males is released after a constant waiting period $T$, that is, the W-males are released periodically and impulsively at discrete time points $T_k = kT, k = 0, 1, 2, \cdots$.

Following ^[26], the authors in ^[27] introduced $\bar{T}$ as the sexual lifespan of W-males. Regarding the release period $T$ and the sexual lifespan $\bar{T}$, the critical case is $T = \bar T$, under which model (1.2) becomes

The dynamics of model (1.3) are simple, as they are only determined by the release amount $c$. However, the release strategy with $T = \bar{T}$ is not practical, since this critical release case is not that easy to implement. A more reasonable release strategy is $T > \bar{T}$ when the density of wild mosquitoes is slightly below the threshold required for dengue prevalence, or $T < \bar{T}$ when the density is above the threshold.

In the current study, we focus on the case with $T > \bar{T}$, under which $g(t)$ takes the form

and then Eq (1.2) can be specified as

and

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

The rest of the paper is organized as follows. Section 2 gives release amount thresholds $c^*$ and $c^{**}$ and the release period threshold $T^*$. In addition, we define a Poincaré map for seeking periodic solutions of (1.5) and (1.6) and offer a lemma to build the connection between the monotonicity of the Poincaré map and the asymptotic stability of the origin. Section 3 shares and proves two theorems: One guarantees the existence of a unique globally asymptotically stable periodic solution, and the other provides sufficient conditions for the existence of at most two periodic solutions. We point out here that a fixed point of the Poincaré map corresponds to an initial value evolving a periodic solution, and the stability of the periodic solution is determined by the relation between the derivative of the Poincaré map at the fixed point and $1$. Different from the extant works ^[27,28,29], the derivative of the Poincaré map at a fixed point has no explicit expressions. Nevertheless, we skillfully find that the derivative of the Poincaré map can be generated with two crucial functions. Then, by investigating the qualitative properties of the two functions, we obtained the above theorems. Next, we offer two numerical examples to support and expand the theoretical results in Section 4. Finally, in Section 5, through comparing the conclusions of (1.5) and (1.6) and those derived from the existing models, we make clear the rationality of the choice of the Ricker-type survival probability in our model and provide a short discussion on our future works.

2.   Preliminaries

Obviously, (1.5) and (1.6) possess the origin, denoted by $E_0$, as the unique equilibrium. Assume that $w(t) = w(t; 0, u)$ is the solution of (1.5) and (1.6) which initiates at $w(0) = u > 0$. Then, $w(t)$ is a continuous and piecewise differentiable function defined on $[0, +\infty)$ which satisfies Eq (1.5) with initial value $w(0) = u$ on $[0, \bar{T})$ and Eq (1.6) with initial value $w(\bar{T})$ on $[\bar{T}, T)$, where $w(\bar{T})$ is defined as the right limit of $w(t)$ at $\bar{T}$. As to the other intervals $[iT, (i+1)T) (i = 1, 2, 3, \cdots)$, the solution can be depicted by a similar method. It then follows from the continuous dependence and differentiable theorems of solutions upon initial values ^[30] that $w(T) = w(T; 0, u)$ is continuously differentiable with respect to $u$, since $w(\bar{T})$ is continuously differentiable with respect to initial value $u$, and $w(T)$ is continuously differentiable with respect to initial value $w(\bar{T})$.

In biology, the sign of the per capita growth rate of a specific population is an indicator which shows whether the density of the population under consideration is increasing. In Eq (1.5), the per capita growth rate of the wild mosquito population satisfies

To determine the sign of $F(w, c)$, we first calculate the partial derivative of $F(w, c)$ with respect to $w$, which yields

Define $G(w, c) = -bw^2-bcw+c$. Then, $\frac{\partial F(w, c)}{\partial w} = \frac{ae^{-bw}}{(w+c)^2} G(w, c)$, and it is easy to see that $G(w, c) = 0$ has a unique positive real root, denoted by

So, $G(w, c) > 0$ holds for $w\in[0, w_+)$, and $G(w, c)\le0$ holds for $w\in[w_+, +\infty)$. Hence, $F(w, c)$ increases for $w\in[0, w_+)$ and decreases for $w\in [w_+, +\infty)$. Together with the fact $F(0, c) = -\mu < 0$, $F(w, c)$ attains its maximum at $w = w_+$ with

Moreover, with simple algebraic operations from (2.3), we know that the function $F(w_+, c)$ is strictly decreasing with respect to $c$. This, combined with the fact $F(w_+, 0) = a-\mu > 0$, shows that there exists a unique $c^*$ such that $F(w_+, c^*) = 0$, and

See Figure 1 for illustration. Thus, we obtain the following theorem.

Theorem 2.1. For Eq (1.5) with $t\ge 0$, there exists a positive implicit threshold $c^*$ such that the following statements are valid.

$(1)$ If $0 < c < c^*$, then (1.5) with $t\ge 0$ has three equilibria: the origin, denoted by $e_0$, and two positive equilibria: $e_1, e_2$ with $e_1 < w_+ < e_2$, where $w_+$ is defined in (2.2). Furthermore, the equilibria $e_0$ and $e_2$ are asymptotically stable, while $e_1$ is unstable.

$(2)$ If $c = c^*$, then $e_1$ and $e_2$ coincide to a unique positive equilibrium $e^*$, and the whole equilibria of (1.5) with $t\ge 0$ are $e_0$ and $e^*$. Moreover, $e_0$ is also asymptotically stable and $e^*$ is semi-stable: stable from the right-hand side and unstable from the left-hand side.

$(3)$ If $c > c^*$, then (1.5) with $t\ge 0$ admits a unique equilibrium $e_0$, which is globally asymptotically stable.

Theorem 2.1 and the dynamical analysis of (1.1) show that the dynamics of single Eq (1.5) or (1.6) are simple. However, with the occurrence of the switches, (1.5) and (1.6) can generate very complicated dynamics, which may not be dealt with mathematically. In consideration of making the relevant dynamical analyses tractable, we discuss the situation $c > c^*$ hereafter.

The elimination of the wild mosquito population eventually, which mathematically is shown as the attractivity of $E_0$, can prevent the outbreak of dengue fever. While achieving the goal of eradicating the wild mosquito population is not that easy, it is unclear whether the lack of the ecological niche of the wild mosquitoes is safe for the local ecosystem. Thus, when applying IIT to combat dengue, suppressing wild mosquitoes to ensure their density be kept below the risk threshold required for the prevalence of the disease, which biologically means the coexistence state of the wild and W-males, is a more realistic and economic strategy than wiping them out.

The coexistence state of two species, a common phenomenon in ecosystems ^[18], can be indicated dynamically by the existence of a periodic solution, and the robustness of this state can be determined by the stability of the periodic solution. To investigate the number of periodic solutions of (1.5) and (1.6) and their corresponding stabilities, we define

Then, we have $h_0(u) = u$, and by induction, we obtain

For convenience, we label

The relation (2.4), along with the existence and uniqueness theorem of solutions of initial value problems, implies that $h_{i+1}(u) > h_i(u)$ provided $h(u) > u$, where $i = 0, 1, 2, \cdots$. More precisely, we have the following lemma.

Lemma 2.2 (^[27]). Let $H(u) = h(u)-u$. Then, for any given initial value $u > 0$, the following statements are valid.

$(1)$ If $H(u) > 0$, then sequences $\{\bar{h}_n(u)\}$ and $\{h_n(u)\}$ are both strictly increasing.

$(2)$ If $H(u) = 0$, then $h_n(u)\equiv u, n = 0, 1, 2, \cdots$. Moreover, $w(t)$ is a $T$-periodic solution of (1.5) and (1.6).

$(3)$ If $H(u) < 0$, then sequences $\{\bar{h}_n(u)\}$ and $\{h_n(u)\}$ are both strictly decreasing.

Lemma 2.2 $(2)$ says that the fixed points of the Poincaré map are precisely the initial values of the periodic solutions of (1.5) and (1.6). Thus, for obtaining periodic solutions, in the following, we are going to seek points $u$ solving $h(u) = u$.

To this end, we first set $\mathbb{A} = \frac{1}{b}\ln\frac{a}{\mu}$, and then Eq (1.6) becomes

which shows that $\left.{dw}\middle/{dt}\right.\le0$ holds for Eq (1.6) when $w\ge\mathbb{A}$. Moreover, it follows from Theorem 2.1 $(3)$ that $\left.{dw}\middle/{dt}\right. < 0$ holds for Eq (1.5) for any initial value, that is, we reach

which implies that any potential initial value evolving some periodic solution lies in the interval $(0, \mathbb{A})$ and gives us an inspiration that we only need to discuss initial values contained in this interval for seeking periodic solutions.

Furthermore, Lemma 2.2 also implies that the sign of $H'(u)$, or, equivalently, the relation between $h'(u)$ and $1$, plays a crucial role in determining the number of periodic solutions and their corresponding asymptotic stabilities ^[31]. Hence, subsequently, we intend to derive the expression for $h'(u)$.

For fulfilling this goal, we ought to solve initial value problem (1.5) with $w(0) = u$ to obtain the expression for $\bar h(u)$ and then solve initial value problem (1.6) with $w(\bar{T}) = \bar{h}(u)$ to gain the relation between $h(u)$ and $\bar h(u)$. We begin with separating variables on both sides of (1.5) and then integrating from $0$ to $\bar{T}$, which provides

Note that the function $\left.{(w+c)}\middle/{\left\{w\left[\left(ae^{-bw}-\mu\right)w-\mu c\right]\right\}}\right.$ is not integrable, which indicates that the method used in ^[27] is not applicable to the analyses of our model dynamics anymore. Nevertheless, by deliberating ^[27], we find that the expression for $h'(u)$ can be obtained even in the sense that the expression for $\bar h(u)$ is implicit.

Keeping this in mind, we go further by taking the derivative of both sides of (2.6) with respect to $u$, which yields

Then, by the same token as above, with the integral upper and lower limits being $\bar{T}$ and $T$, respectively, we attain

Finally, substituting (2.7) into (2.8), we arrive at

Let

and then we achieve

Clearly, under the assumptions $u\in(0, \mathbb{A})$ and $c > c^*$, we get

respectively. Set $\Gamma(u) = \left.{\Omega_1(u)}\middle/{\Omega_2(u)}\right.$. Then, we have $\Gamma(u) < 0$, and (2.11) becomes

The facts

show that the relation between $h(u)$ and $u$ when $u > 0$ and sufficiently approaches zero is vital to the existence of periodic solutions of (1.5) and (1.6). In other words, the relation between $h'(0)$ and $1$ directly determines the existence of periodic solutions of the model. We thus focus on computing $h'(0)$ in the following.

First, noting that $\bar{h}(u)\to0$ and $h(u)\to0$ as $u\to0$, we achieve $w(t)\to0$ as $u\to0$, where $t\in[0, T]$. When $t\in[0, \bar{T})$, (1.5) gives

and then integrating (2.15) from $0$ to $\bar{T}$, we obtain

Thus, we get

Furthermore, when $t\in[\bar{T}, T)$, (1.6) yields

and then integrating (2.17) from $\bar{T}$ to $T$, we have

which yields

Then, (2.16) and (2.18) imply

We then define the release period threshold $T^*$ and the other release amount threshold $c^{**}$ as follows:

Thus, we arrive at

Moreover, simple calculations give $c^{**} > c^*$. In fact, from (2.3), we have

Note that

and we obtain $F(w_+, c^{**}) < 0$. As $F(w_+, c)$ is decreasing with respect to $c$ and $F(w_+, c^*) = 0$, we derive $c^{**} > c^*$.

On the basis of the above, we have been primed to share our main results below.

3.   Model dynamics

According to (2.21), the condition $T > T^*$ guarantees $h'(0) > 1$, which, together with (2.14), implies that (1.5) and (1.6) have at least one periodic solution. In fact, (1.5) and (1.6) have a unique globally asymptotically stable periodic solution provided $T > T^*$, that is, the following theorem.

Theorem 3.1. Assume $T > T^*$. Then, (1.5) and (1.6) have a unique periodic solution, which is globally asymptotically stable.

Proof. By (2.21), $T > T^*$ signals that there exists $\delta > 0$ sufficiently small such that

The above inequality, together with (2.14), implies the existence of periodic solutions of (1.5) and (1.6).

We now turn to prove the uniqueness of periodic solutions. For convenience, we divide the proof into two cases, $c\ge c^{**}$ and $c^* < c < c^{**}$.

We first consider the case $c\ge c^{**}$. Assume by contradiction that the model has at least two periodic solutions. Then, there must exist $\bar{u} > 0$ such that

Inserting $h(u) = u$ into (2.13), we have

Since $\bar{h}(u) < u$, to reveal the relation between $h'(u)$ and $1$ with $u$ satisfying $h(u) = u$, we need to check the monotonicity of $\Gamma(u)$.

Taking the derivative of $\Gamma(u)$, we have

Note that the function $\left.{(ae^{-bu}-\mu-b\mu u)}\middle/{(b\mu)}\right.$ is strictly decreasing with respect to $u$ when $u\in(0, \mathbb{A})$, and its maximum is less than $\left.{(a-\mu)}\middle/{(b\mu)}\right.( = c^{**})$, which indicates $\Gamma'(u) > 0$ when $c\ge c^{**}$, that is, $\Gamma(u)$ is strictly increasing in this case. This, combined with the facts $\bar{h}(u) < u$ and $\Gamma(u) < 0$, gives $0 < h'(u) < 1$, a contradiction to (3.1).

Subsequently, we analyze the case $c^* < c < c^{**}$. Assume by contradiction that (1.5) and (1.6) possess at least two periodic solutions, and then there exist $v_1, v_2$ contained in $(0, \mathbb{A})$ with $v_1 < v_2$, such that

Furthermore, revisit the expression for $\Gamma'(u)$, and we know that there is $u^* = u^*(c)\in(0, \mathbb{A})$ such that

which means that $\Gamma(u)$ is strictly decreasing in $(0, u^*)$ and strictly increasing in $(u^*, \mathbb{A})$. Then, the facts $\Gamma(u) < 0$ and $\bar{h}(u) < u$, coupled with (3.2) and (3.3), imply that there is $\hat{u}\in(u^*, \mathbb{A})$ such that

Obviously, the critical conclusion (3.4) excludes the possibility that $h(u)$ has three or more fixed points. Next, we focus on obviating the situation that $h(u)$ has exactly two fixed points. To this end, we need to consider the following two cases:

See Figure 2 for illustration.

Moreover, (3.4) can also show that case $(a)$ of (3.5) is not true, and we thus only need to exclude case $(b)$ to finish the proof of the uniqueness of periodic solutions of (1.5) and (1.6).

Now, we turn to a perturbation method for obtaining a contradiction. Define $H_k(u) = h(u)-ku$. Let $k-1 > 0$ be sufficiently small such that $H_k(u) = 0$ has exactly three roots, denoted by $\bar{v}_1, \bar{v}_2, \bar{v}_3$, satisfying $\bar{v}_1 < \bar{v}_3 < \bar{v}_2$. Then, we have

Substituting $h(u) = ku$ into (2.9), we get

and (3.6) becomes

which, together with the facts $\bar h(\bar v_1) < \bar v_1, \bar h(\bar v_3) < \bar v_3$, and (3.4), gives $\bar{v}_3\le\hat u\le\bar{v}_1$. This contradicts the fact $\bar{v}_3 > \bar{v}_1$ and completes the proof of the uniqueness of periodic solutions of (1.5) and (1.6).

Denote the unique periodic solution of the model by $\bar{w}(t)$, and assume that its initial value is $\bar{u}_0\in(0, \mathbb{A})$. From the above analyses, we can deduce that

Then by ^[31], we know that $\bar{w}(t)$ is stable. Moreover, (3.8) and Lemma 2.2 imply that the basin of attraction of $\bar{w}(t)$ is $\{u|u > 0, u\neq\bar{u}_0\}$, which means that $\bar{w}(t)$ is globally attractive, and the proof is finished. □

Theorem 3.1 signals the potential relation between the stability of the origin and the exact number of periodic solutions of (1.5) and (1.6): The instability of the origin is sufficient to ensure the existence of periodic solutions. However, if the origin is stable, then the analyses on the exact number of periodic solutions become rather complicated, and we may also need to restrict the release amount $c$ to seek the exact number of periodic solutions, which is reflected in the following theorem.

Theorem 3.2. Assume that $T < T^*$. Then, (1.5) and (1.6) have at most two periodic solutions, and the origin is asymptotically stable. Moreover, the following two conclusions hold.

(ⅰ) If $c\ge c^{**}$, then the model has no periodic solutions, and the origin is globally asymptotically stable.

(ⅱ) If $c^* < c < c^{**}$, and the model admits one or two periodic solutions, then there exist two positive numbers $\widetilde{\delta}$ and $\widetilde{u}$, with $\widetilde{\delta}$ sufficiently small and $\widetilde{u}\in(0, \mathbb{A})$, such that at least one initial value of the periodic solution(s) is located in the interval $(\widetilde{\delta}, \widetilde{u})$.

Proof. Since $T < T^*$, from (2.21), we have $h'(0) < 1$, which shows that there exists $\widetilde{\delta} > 0$ sufficiently small such that

This implies the asymptotic stability of the origin.

(ⅰ) By utilizing a similar method to that of Theorem 3.1, we can also derive $0 < h'(u)|_{u = h(u)} < 1$ in this situation.

Assume by contradiction that the model has a periodic solution with initial value $\underline{u}\in(0, \mathbb{A})$. Thus, combining (3.9) and the fact $h(\mathbb{A}) < \mathbb{A}$, we can also obtain $h'(\underline{u})\ge1$, a contradiction to $0 < h'(u)|_{u = h(u)} < 1$, which finishes the proof of the nonexistence of periodic solutions for the model.

Since the model has no periodic solutions when $u\in(0, \mathbb{A})$, we get

which gives, together with (2.5),

This illustrates the global asymptotic stability of the origin and completes the proof of case (ⅰ).

(ⅱ) Suppose by contradiction that the model has at least three periodic solutions. We first analyze the case that the model has exactly three periodic solutions and denote their initial values by $u_1, u_2$ and $u_3$. Then, inevitably, we need to discuss the following four cases:

See Figure 3 for illustration.

Similar to the proof of Theorem 3.1 for the case when $c\in(c^*, c^{**})$, we can see that (3.4) is also true in the current situation, which excludes the possibilities of cases $(1)$–$(3)$. Next, we concentrate on excluding case $(4)$.

Similarly, let $k-1 > 0$ be small enough such that $H_k(u) = 0$ has exactly four roots, denoted by $\bar{u}_1, \bar{u}_{21}, \bar{u}_{22}, \bar{u}_3$, with $\bar{u}_1 < \bar{u}_{21} < \bar{u}_{22} < \bar{u}_3$. Then, at point $u$ satisfying $h(u) = ku$, we have

Since (3.7) is also valid, the above four inequalities can be transformed to

and

respectively. Then, the facts $\bar{h}(\bar{u}_{21}) < \bar{u}_{21}, \bar{h}(\bar{u}_{22}) < \bar{u}_{22}$ and (3.4) imply $\bar{u}_{22}\le\hat u\le\bar{u}_{21}$, which contradicts the fact $\bar{u}_{22} > \bar{u}_{21}$.

Obviously, (3.4) can also exclude the possibilities of the existence of four or more periodic solutions for the model. Thus, the model admits at most two periodic solutions in this case.

Next, we prove that there exists a positive number $\widetilde{u}\in(\widetilde{\delta}, \mathbb{A})$, such that at least one initial value of the potential periodic solutions is located in $(\widetilde{\delta}, \widetilde{u})$.

To this end, we need to investigate the monotonicities of $\Omega_1(u)$ and $\Omega_2(u)$, which are defined in (2.10). We first analyze the monotonicity of $\Omega_1(u)$. Since $\Omega'_1(u) = ae^{-bu}(1-bu)-\mu$, we get

and $\Omega''_1(u) = abe^{-bu}(bu-2)$. Thus, $\Omega''_1(u)$ has a unique zero $\breve{u} = \left.{2}\middle/{b}\right. > 0$, and then we consider the following two cases to compare $\breve{u}$ and $\mathbb{A}$ and judge the sign of $\Omega''_1(u)$.

Case 1: $\mu < a\le\mu e^2$. In this case, $\breve{u}\ge\mathbb{A}$, so $\Omega''_1(u) < 0$ holds for all $u\in(0, \mathbb{A})$, which means that $\Omega'_1(u)$ is decreasing in $(0, \mathbb{A})$. Hence, there is $\widetilde{u}_1\in(0, \mathbb{A})$, such that $\Omega'_1(\widetilde{u}_1) = 0$ and

Case 2: $a > \mu e^2$. We have $0 < \breve{u} < \mathbb{A}$ in this case. Then,

or, equivalently, $\Omega'_1(u)$ is decreasing in $(0, \breve{u})$ and increasing in $(\breve{u}, \mathbb{A})$. Then

where $\Omega'_1(u)_{min}$ denotes the minimum of $\Omega'_1(u)$. Thus, there exists $\widetilde{u}_2\in(0, \mathbb{A})$, such that $\Omega'_1(\widetilde{u}_2) = 0$ and

Based on the above analyses, we know that $\Omega'_1(u)$ has a unique zero, that is, we have $\widetilde{u}_1 = \widetilde{u}_2$. We denote the unique zero by $\widetilde{u}$. Then, we obtain

In the following, we discuss the monotonicity of $\Omega_2(u)$. Since

the fact $\Omega_2(u) < 0$ and (3.11) yield that $\Omega'_2(u) < 0$ holds when $u\in(\widetilde{u}, \mathbb{A})$. This, together with (2.11) and (3.9), provides $h'(u) < 1$ when $u\in(0, \widetilde{\delta})\cup(\widetilde{u}, \mathbb{A})$. Consequently, if the model admits a unique periodic solution, then its initial value must be contained in $(\widetilde{\delta}, \widetilde{u})$. Otherwise, if the model admits exactly two periodic solutions, then the smaller initial value must be limited in $(\widetilde{\delta}, \widetilde{u})$. This completes the proof. □

Theorem 3.2 signals that when the release becomes more frequent such that $T < T^*$, the origin is always asymptotically stable, which biologically means that the local wild mosquito population can always be wiped out eventually, provided we employ the integrated mosquito management strategy to control the wild mosquitoes such that their density is being maintained within the basin of attraction of the origin. Furthermore, under this scenario, Theorem 3.2 implies that the model possesses at most two periodic solutions. Nevertheless, we have not determined the corresponding conditions for the model to admit zero, a unique or exactly two periodic solutions in the current work. As in ^[29], we perceive that there exists another release period threshold, $T^{**}$ say, such that the exact number of periodic solutions is governed by the relation between $T$ and $T^{**}$, and we intend to tackle the relevant theoretical analyses with further research efforts in future works.

4.   Numerical examples

In this section, we offer two numerical examples to support and expand Theorems 3.1 and 3.2. We first give a numerical example to support Theorem 3.1.

Example 1. Given parameters

Then, we have $\mathbb{A}\approx7.2759$, and

Fix $T = 15 > T^*$ and limit $c > c^*$, and then the conditions of Theorem 3.1 are satisfied. To numerically illustrate Theorem 3.1, we first select $c = 100 > c^{**}$ or $c = 70\in(c^*, c^{**})$ in panels (A) and (B) of Figure 4, respectively. Then, we plot the graph of $w(t)$ in panel (A) and find that (1.5) and (1.6) admit a unique periodic solution, which is globally asymptotically stable. In panel (B), from a different viewpoint, we draw the image of $H(u)$ when $u\in(0, 8)\supset(0, \mathbb{A})$, where the vertical blue dotted line represents $u = \mathbb{A}$. From panel (B), we observe that $H(u) = 0$ has a unique positive root in $(0, \mathbb{A})$, which corresponds to the initial value of the unique periodic solution in panel (A).

Then, we provide a numerical example to expand and extend Theorem 3.2.

Example 2. Suppose the parameters are given the same as that of (4.1) and fix $T = 14.2 < T^*$. If we choose $c = 100 > c^{**}$, then the conditions for Theorem 3.2 (ⅰ) are satisfied, and the origin $E_0$ is thus globally asymptotically stable as shown in panel (A) of Figure 5. Moreover, if we select $c = 70\in(c^*, c^{**})$, then the conditions for Theorem 3.2 (ⅱ) are satisfied. To numerically explore the exact number of periodic solutions in this situation, we plot the image of $H(u)$ in panel (B) of Figure 5, where the two vertical dotted lines in black and blue represent $u = \widetilde{u}$ and $u = \mathbb{A}$, respectively. This image implies that (1.5) and (1.6) have no periodic solutions in this case, which is consistent with the theoretical results of Theorem 3.2 (ⅱ).

5.   Discussion and conclusions

Nowadays, dengue is endemic in many countries, such as Pakistan ^[32] and Nepal ^[33]. From January to $28$ September $2022$, a total of $25,932$ confirmed cases and $62$ deaths (CFR $0.25%$) and $28,109$ detected and suspected cases and $38$ identified deaths (overall CFR $0.13%$) were recorded in Pakistan and Nepal, respectively. Since there are no specific treatments, and the vaccine lacks reliability due to the antibody-dependent enhancement effect ^[34], the methods for the prevention and control of dengue mainly focus on suppressing or eliminating wild Aedes mosquitoes: the primary vectors of the disease. There are various strategies for reducing the density of wild Aedes mosquitoes, including larval source reduction, the use of curtains and water container covers treated with insecticide and indoor residual spraying ^[35], to name a few. Nevertheless, disturbingly, the geographic range of dengue is expanding, and the severity of its outbreaks is increasing. There is an urgent need for appropriate alternatives to the prevention and control of dengue ^[22].

The incompatible insect technique (IIT) is a suitable and promising weapon for controlling the density of wild mosquitoes, and it relies on the massive production and sufficient releases of W-males to induce a considerable high sterility in the wild females, which causes the density of the target population to decline ^[19,36]. To guarantee the success of IIT at a relatively acceptable cost, the release strategy, which in general is prior to the actual releases, needs to be carefully designed, and there are numerous studies dedicated to revealing the relation between the release strategy and the suppression effect ^{[37,38,39,40,41,42,43]}.

Clearly, the survival probability of wild mosquitoes in the aquatic stage could also have a non-negligible impact on the suppression effect ^[28,44,45]. In order to assess the impact of the survival probability on the suppression dynamics, we developed a time-switched wild mosquito population suppression model with Ricker-type survival probability in this work, that is, (1.5) and (1.6). It is worth mentioning that the Ricker-type survival probability is commonly adopted in discrete models, for it shares the advantages such as nonnegativity and high accuracy. However, this type of survival probability is rarely applied in continuous models, especially in ordinary differential equation or equations models, since, inevitably, a transcendental equation must be analytically solved for obtaining equilibria of the concerned model in this case. To the best of our knowledge, it is the first time that the Ricker-type survival probability is incorporated into an ordinary differential equations model, and the resulting dynamics are very rich and rather complicated.

As we all know, thresholds, which can be identified via mathematical manners, play a crucial role in mathematical biology. In this work, we found three release thresholds, i.e., two amount thresholds $c^*$ and $c^{**}$ satisfying $c^{**} > c^*$ and a period threshold $T^*$. Comparing the associated thresholds defined in (2.20) with that of ^[28], we observed that the two release period thresholds coincide with each other. This is not surprising at all, since the two models share the following unified form:

where

Simple calculations, together with the fact $w\to0$ as $u\to0$, give $h'(0) = e^{(a-\mu)\left(T-\frac{a}{a-\mu}\bar{T}\right)}$, which implies that the release period threshold $T^*$ equals $\frac{a}{a-\mu}\bar T$, regardless of whether the form of $Ł$ is $Ł_1$ or $Ł_2$. However, the release amount thresholds specified in (2.20) and ^[28] are generally different: Their ratio is equal to $\left.{\xi}\middle/{b}\right.$. According to ^[45,46], there exists some habitat of wild mosquitoes in which the parameter $b$ is larger than the parameter $\xi$, so, to ensure the global asymptotic stability of the origin, the smallest amount of W-males needed from (1.5) and (1.6) is less than that of ^[28]. In other words, to eradicate the wild mosquitoes, the model which employs a Ricker-type survival probability yields a relatively lower cost than that of ^[28] in some settings. Hence, the results could stimulate the development of more efficient and cost-effective strategies for controlling mosquito-borne diseases including dengue fever.

In biology, the fluctuation alteration of two types of species in the system under consideration, which mathematically can be reflected by a specific periodic solution ^[47], is a common phenomenon, and it is vital to the evolution of the two species as well as the persistent stability of the relevant system ^[48]. Thus, for (1.5) and (1.6), it is valuable to investigate the number of periodic solutions and their corresponding stabilities. However, mathematically, it is very challenging to tackle these topics. To overcome these difficulties, in this study, via defining a Poincaré map $h(u)$ and performing rigorous dynamical analyses, we found that the initial value of a specific periodic solution is exactly the fixed point of $h(u)$, and the stability of the periodic solution is determined by the relation between $h'(u)|_{u = h(u)}$ and $1$. Then, by seeking the fixed points of $h(u)$ and computing $h'(u)|_{u = h(u)}$, we observed that $h'(u)|_{u = h(u)}$ could be expressed with two functions whose signs are opposite. Moreover, the analyses on the monotonicities of the two functions generate our main results: Theorems 3.1 and 3.2.

Yet, for the dynamical analyses of models with Ricker-type survival probability, we still have a long way to go. In particular, when the release parameter combination $(T, c)$ lies in the set $\{(T, c)|\bar{T} < T < T^*, c^* < c < c^{**}\}$ and the initial value $u\in[\widetilde{\delta}, \widetilde{u}]$, the theoretical proof for the nonexistence of periodic solutions has not been explored. In addition, the model dynamics under the critical case $T = T^*$ also have not been evaluated since it is difficult to get the explicit expression for $h(u)$. Finally, we have not touched the investigation of the model dynamics for the case when $0 < c < c^*$. Nevertheless, we will try our best to tackle these topics in our future efforts.

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 (11971127, 12071095).

Conflict of interest

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