Mirrored dynamics of a wild mosquito population suppression model with Ricker-type survival probability and time delay

1.   Introduction

Dengue fever is a mosquito-borne viral disease caused by any of the four dengue virus serotypes that spread between humans and Aedes mosquitoes, primarily Aedes aegypti. Its clinical symptoms include high fever, severe headache, vomiting, muscle and joint pain and rash ^[1,2]. In some cases, the disease can develop into life-threatening dengue shock syndrome, which may result in hypotension, marked thrombocytopenia, plasma leakage, leucopenia, hepatomegaly, hypoproteinemia, circulatory failure and mortality in 1–5% of infected individuals ^[3,4]. Nowadays, the emergency and global incidence of the disease are rapidly increasing, as shown in Figure 1. According to ^[5], the disease may cause an estimated 400 million infections and 100 million symptomatic cases worldwide every year, posing a serious threat to human health and life security ^[6].

The sterile insect technique (SIT) ^[7,8,9,10] and the incompatible insect technique (IIT) ^{[11,12,13,14]}, which aim at suppressing the indigenous mosquito population, are two effective and eco-friendly weapons for combating the disease. SIT (IIT) relies on massive production and remarkable release of the radiation-treated (the endosymbiotic bacterium Wolbachia-infected) male mosquitoes (we refer to these two kinds of male mosquitoes as sterile mosquitoes hereafter) reared in the labs or mosquito factories to sterilize wild female mosquitoes and, hence, to lower the density of the wild mosquito population ^{[15,16,17,18]}.

Mathematics, as a tool discipline, can provide in-depth insights into many problems that are challenging to be dealt with in many other disciplines ^{[19,20,21,22,23,24,25,26,27,28,29]}. Once sterile mosquitoes are released into the target environment, these two types of mosquitoes will spontaneously give rise to interactions. To study the interactive dynamics between wild and sterile mosquitoes and to aid in seeking optimal release strategies, a body of mathematical models have been formulated and dissected; we refer to ^{[30,31,32,33,34,35]} for ordinary differential models, ^{[36,37,38,39,40,41]} for delay differential models and the references cited therein.

Very recently, Z. Zhu and X. He ^[42] developed the next ordinary differential equation model

to explore the impact of the Ricker-type density-dependent survival probability $e^{-bw}$ of wild mosquitoes (in the aquatic stages) on the suppression effect. Here, $w = w(t), g = g(t)$ denote the numbers of wild and sexually active sterile mosquitoes at time $t$, respectively, $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 ^[43,44] and $\mu > 0$ describes the density-independent death rate of wild mosquitoes. The authors assumed that the energetic sterile mosquitoes are released in a periodic and impulsive manner, such that a constant amount $c$ of sterile mosquitoes is released after a constant waiting period $T$. That is, the sterile mosquitoes are released periodically and impulsively at discrete time points $T_i = iT, i = 0, 1, 2, \cdots$, and the sexual lifespan of the sterile mosquitoes $\bar T$ is less than the release period $T$, under which $g(t)$ has the following structure of the form

By injecting (1.2) into (1.1), the authors obtained the following switched ordinary differential equations

and

where $i = 0, 1, 2, \cdots$. After the identifications of an implicit threshold $c^*$, together with two explicit thresholds

the authors investigated the model dynamics and compared it with that of ^[34], then they asserted that, under some specific conditions, model (1.1) drives a more satisfactory suppression effect while possessing a much lower cost.

It is widely known that the growth of a mosquito in each stage requires time ^[45]. Thus, the population densities of wild and sterile mosquitoes at some previous time $t-\tau$ play significant roles in determining the population size of wild adults at time $t$. Here, $\tau$ is the average waiting time from parent mating to the emergence of the reproductive offspring. Therefore, models that take the maturation time of the mosquito species into account are more realistic and meaningful. Such models can exhibit much richer and more complex dynamics than those delay-free models ^{[34,36,37,38]}.

Based on the above consideration and by incorporating the maturation delay into (1.1), we obtain the following delay differential equation model

The initial condition for the model is

where $\varphi: [-\tau, 0]\to[0, +\infty)$ is continuous.

Once the eggs are laid, it takes 1–2 days for an egg to hatch into a larva. The larva enters the pupal stage within 7–10 days, then it takes 2–3 days for the pupa to develop into an adult ^[37]. It is worth noting that the maturation period of sterile mosquitoes is almost the same as that of wild mosquitoes. For the sake of saving manpower and cost, we further assume in this paper that the waiting period between two consecutive releases is exactly the maturation period; that is, $T = \tau$. This scenario results in Eq (1.6) being transformed into the following three sub-equations

and

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

The rest of the work is organized as follows. Section 2 mainly gives a lemma, which acts as a bridge to our main results. In Section 3, based on whether the solution of model (1.7)–(1.10) oscillates, we first adopt an ingenious transformation, which turns (1.6) into an equivalent equation. Next, our attention is focused on the exploration of the solution of that equation. Subsequently, by investigating the qualitative properties of the solution's limit superior and limit inferior, we derive the expected theoretical results. Finally, Section 4 offers some discussions, expansions and conjectures about the theoretical results of the delay model.

2.   Preliminaries

When $c > c^*$ and $T > T^*$, Theorem 3.1 of ^[42] manifests that the time-switched ordinary differential equations model (1.3)–(1.4) possesses a unique periodic solution denoted by $\bar w(t)$; it then follows that the associated delay differential equations model (1.7)–(1.10) also admits $\bar w(t)$ as the unique periodic solution.

For the attractivity of $\bar w(t)$, we introduce the change of variables

under which (1.6) can be transformed to

In the following, we will turn to explore some qualitative properties of $x(t)$. First, for the boundedness of $x(t)$, we have the next lemma.

Lemma 2.1. The function $x(t)$ is bounded.

Proof. We only prove that $x(t)$ is bounded from above (the proof of $x(t)$ bounded from below is similar and, hence, omitted). Otherwise, there exists a sequence of points $\{t_n\}$ such that

and

which, together with Eq (2.2), gives

that is,

In view of $\lim\limits_{n\to+\infty}x(t_n) = +\infty$, we have

which, in fact, is impossible since

where $\bar w_* = \min\limits_{0\le t\le T}\bar w(t)$. The proof is completed. □

3.   Global attractivity of the unique periodic solution

From Section 2, we have already known that model (1.7)–(1.10) has a unique periodic solution $\bar w(t)$ provided that $c > c^*$ and $T > T^*$. Furthermore, the boundedness of $x(t)$ implies that there exists two constants $x_*$ and $x^*$ such that $x_*\le x(t)\le x^*$ holds for all $t\ge0$.

In this section, we aim at proving the global attractivity of $\bar w(t)$. To this end, we first define

where $c^{**}$ is specified in (1.5). (Clearly, we have $l > 1$.) Under the following hypotheses

${\bf (H1)}$ $\bar w_*\ge\bar b = \left.{l}\middle/{b}\right.$;

${\bf (H2)}$ $\frac{a}{\mu}\le\frac{m_1le^l}{m_2(2l+|x^*|)}(m_1, m_2\in\mathbb{Z}_+$ and $m_1 < m_2)$;

${\bf (H3)}$ $\gamma: = \left.{a(1-e^{-\mu T})}\middle/{\mu}\right.\neq1$,

we derive the following theorem.

Theorem 3.1. Assume that the release period $T$ and release amount $c$ satisfy $c > c^*$ and $T > T^*$, respectively, then model (1.7)–(1.10) has a unique $T$-periodic solution $\bar w(t)$. Furthermore, if the release intensity $c$ satisfies $c^* < c < c^{**}$ and the above hypotheses (H1)–(H3) hold, then for any positive solution $w(t)$, we have

Proof. On the basis of the above, we only need to prove that under the conditions $c^* < c < c^{**}$ and $T > T^*$ and the hypotheses (H1)–(H3), $\bar w(t)$ is globally attractive, which in view of (2.1), is equivalent to show $\lim\limits_{t\to+\infty}x(t) = 0$.

To that effect, we distinguish the following two cases.

Case 1: Suppose that the function $x(t)$ is oscillatory about $0$, that is, the equation $x(t) = 0$ possesses arbitrarily large roots.

Let a sequence of points $\{\xi_n\}$ satisfying $\tau\le\xi_1 < \xi_2 < \cdots < \xi_n < \cdots$ and $\lim\limits_{n\to+\infty}\xi_n = +\infty$ be the roots of $x(t) = 0$, and $x(t)$ assumes sign-changing in each interval $(\xi_n, \xi_{n+1})$. For each $n\ge1$, we set $t_n$ and $s_n$ as the points belonging to $(\xi_n, \xi_{n+1})$, such that

Undoubtedly, the relations $x(t_n) > 0, x'(t_n) = 0$ and $x(s_n) < 0, x'(s_n) = 0$ hold for each $n\ge1$. We first demonstrate that the following statements

and

are valid for each $n\ge1$. See Figure 2 for illustration.

In fact, if (3.3) is not true, then we have $\xi_n < t_n-\tau < \xi_{n+1}$ and $x(t_n-\tau) > 0$. However, combining the fact $x'(t_n) = 0$ and (2.2), we obtain

Define

Clearly, $\Xi(0) = 0, \Xi(+\infty) = \left.{v^2}\middle/{(v+g)}\right. > 0$, and

where $g = g(t-\tau)$ is defined in (1.2).

Set

then it is easy to see that $\kappa(u) > 0$ holds for $u > 0$. Thus, $\Xi(u) > 0, u > 0$. Choosing $v = \bar w(t_n-\tau)\ge\bar w_*\ge\bar b$ (as the hypothesis (H1) holds) and $u = x(t_n-\tau) > 0$, we observe that equality (3.5) is impossible.

On the other hand, if (3.4) is not valid, then we obtain $\xi_n < s_n-\tau < \xi_{n+1}$ and $x(s_n-\tau) < 0$. Since $x'(s_n) = 0$, we have from (2.2) that

which gives

Further computations offer

since $x(s_n-\tau) < 0$ and the hypothesis (H2) holds. Hence, we get $x(s_n) > x(s_n-\tau)$, which is a contradiction to the second equality of (3.2).

Next, let

and by choosing $v = \bar w(t-\tau) = \bar w(t)(\ge\bar w_*\ge\bar b)$, then we have

and (2.2) becomes

The mean value theorem implies that the above equation can be written as

where

and $\eta(t)$ is located between $0$ and $x(t-\tau)$, which results in

In view of $x(t-\tau) > -b\bar w(t-\tau) = -b\bar w(t)$, we gain $\eta(t) > -b\bar w(t)$.

Note that

then by setting $u = \eta(t)(>-bv)$, we have

Multiplying Eq (3.6) by $e^{\mu t}$ and then integrating over $[\alpha_n, t_n]$ and $[\beta_n, s_n]$, we obtain

and

respectively.

Denote

Now, for any given $\varepsilon > 0$, there exists a sufficiently large $n_0$ such that for any $t\ge n_0$, we have $x(t-\tau) < \bar\lambda+\varepsilon$ and $x(t-\tau) > \underline\lambda-\varepsilon$, i.e., $-x(t-\tau) < -\underline\lambda+\varepsilon$. Thus, $|x(t-\tau)| < M+\varepsilon, t\ge n_0$, where $M = \max\{\bar\lambda, -\underline\lambda\}$. Injecting this into (3.8) and (3.9) together with the facts that $t_n-\tau\le\alpha_n, s_n-\tau\le\beta_n$ and (3.7), we get

and

that is,

and

Considering the fact that $\varepsilon > 0$ is arbitrary and taking $\limsup$ as $n\to+\infty$ in (3.11) and (3.12), we arrive at

and

Finally, combining (3.13) and (3.14), we get $0\le M\le\gamma M$. It then follows from the hypothesis (H3) that $M = 0$; thus, $\bar\lambda = \underline\lambda = 0$, so $\lim\limits_{t\to+\infty}x(t) = 0$.

Case 2: Suppose that the function $x(t)$ is not oscillatory about $0$. Without loss of generality, we assume that $x(t)$ is eventually positive; in other words, there exists $t_0\ge0$ such that

Now, assume by contradiction that (3.1) does not hold. Revisit (2.2) and we find that $\lim\limits_{t\to+\infty}x(t) = 0$ as long as the limit of $x(t)$ exists. Thus, the limit of $x(t)$ does not exist, so

Hence, there is a sequence of points $\{t_n\}$ satisfying $t_n\ge t_0+\tau$ and

and

Clearly, we have the following inequality:

Combining (2.2) and (3.15), we get

and further computations provide

Let

where $g = g(t-\tau)$ is defined in (1.2), then

Clearly, $L'(u) < 0$ holds for $u\ge\bar b$ and $p > 0$, which implies that the function $L(u)$ is strictly decreasing for $u\ge\bar b$ and $p > 0$. Letting $u = \bar w(t_n-\tau) = \bar w(t_n)(\ge\bar w_*\ge\bar b)$ and $p = \frac{1}{b}x(t_n-\tau)(>0)$ in $L(u)$, we gain

which, together with the hypothesis (H2), reduces (3.16) to $x(t_n) < \frac{m_1}{m_2}x(t_n-\tau)$. Thus

This contradicts the facts $\bar\lambda > 0, m_1, m_2\in\mathbb{Z}_+$ and $m_1 < m_2$ and establishes the proof. □

4.   Discussion

Based on ^[42], in which a time-switched wild mosquito population suppression model with Ricker-type density-dependent survival probability was proposed and investigated, focusing more on the biological sense, here we introduced the maturation delay of the mosquito species into that model and obtained the core model (1.6). Meanwhile, by adopting the specific and critical release strategy that the release period $T$ of sterile males equals the maturation period $\tau$ of wild mosquitoes ^[37], we got the resulting time-switched delay differential equations model (1.7)–(1.10). Note that, regarding the core model, its growth term lacked monotonicity, which considerably obstructed the mathematical exploration of its dynamical features as there was no existing broad-spectrum route to investigate. Nevertheless, via an ingenious transformation presented in (2.1), we gained an equivalent equation (2.2) of the core model. Subsequently, through concentrating our attention on (2.2) and studying the qualitative properties of the limit superior and limit inferior of its solution, we derived the mirrored dynamics elucidated in Theorem 3.1. The findings in the current work can not only act as a supplement to the theoretical study of the impact of time delay on the dynamical features harbored by the mosquito population suppression models, but also aid in designing more cost-effective strategies for preventing and controlling dengue fever and other vector-borne diseases.

It should be stressed here that it was the periodically release, not the time lag, that made the delay model admit $\bar w(t)$ as the unique periodic solution. This assertion can be inferred from the fact that the corresponding delay-free model also admitted $\bar w(t)$ as the unique periodic solution. In fact, the mirrored dynamics must be driven by the specific release strategy of $T = \tau$, and it will disappear as long as this release strategy is slightly changed such that $T\neq\tau$. Furthermore, inspired by the findings in ^[39] and ^[46] together with the accumulated knowledge about the dynamical features of the model, we guessed that there were two thresholds $\tau^*$ and $\tau^{**}$ that played a central role in determining the changing trend of the size of the amplitude of the periodic solution. More precisely, with $\tau$ increasing, whether the amplitude was being increased depended on the relation between $\tau$ and the two thresholds: It was increasing when $\tau$ was less than $\tau^*$ and it was decreasing if $\tau$ was between $\tau^*$ and $\tau^{**}$, while it exhibited no clear-cut regularity once $\tau$ exceeded $\tau^{**}$.

Last but not least, we'd like to point out here that we only offered some theoretical insights into the Ricker-type time-switched delay differential equations model (1.7)–(1.10), and there was neither suitable data nor numerical examples provided to support and illustrate the analytical results. A possible reason is that model (1.7)–(1.10) contained the unimodal and non-monotonic population growth term, the time delay and the periodic switches, each of which can induce extremely rich and rather complicated model dynamics together with the resulting effects of the combination of the hypotheses (H1)–(H3). Nevertheless, with a strong belief, we are sure that in the near future, the expected data will be found in either a specific ecological environment of wild mosquitoes or some other disciplines.

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 National Natural Science Foundation of Guangxi Province (No: 2020JJG110003), National Natural Science Foundation of China (Nos: 12331017, 12371484, 12301621), Shanxi Scholarship Council of China (No: 2022-176), National Natural Science Foundation of Shanxi Province (No: 202203021211115) and China Postdoctoral Science Foundation (No: 2023M740802).

Conflict of interest

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