Modeling environmental transmission of MAP infection in dairy cows

1.   Introduction

As a rapidly spreading mosquito-borne disease, dengue fever poses severe economic burden and public health threat in tropical and subtropical areas ^[1]. In recent years, dengue virus has infected over $50-100$ million people each year, and almost half the world's population has been threatened ^[1]. Brazil has experienced a serious outbreak of dengue, with over $6.3$ million dengue cases and $4483$ death cases since the start of 2024. In the absence of licensed vaccines and effective therapeutic methods, mosquito control remains the main method to control dengue fever, which focuses on larval source reduction and adult control based on chemical pesticides ^[1,2]. The growing problems in pesticide resistance and environmental pollution give rise to failure of the current control measures. Innovative methods are sought to prevent and control dengue fever.

The incompatible insect technique (IIT) based on Wolbachia has been proven to be a promising technology for dengue control ^[3,4]. Some Wolbachia strains can block the reproduction of dengue virus in Aedes mosquitoes, and cause a mechanism called cytoplasmic incompatibility (CI) that causes infected male mosquitoes to be effectively sterile when they mate with uninfected females ^[4,5]. By releasing millions of factory-reared Wolbachia infected males in Shazai island in Guangzhou since 2015, the IIT approach has been applied successfully to suppress Aedes albopictus populations, with over $97\%$ reduction of adults in the control areas ^[4].

As a frequently used tools in theoretical research, mathematical models have played a critical role in designing release strategies and optimizing release programmes ^[4,6,7]. Typically, differential equations or difference equations are often used to define the target mosquito populations, identify the threshold release level of infected males for population suppression, and quantify the impact of various factors on the threshold. These factors include environmental stochasticity ^[8,9], spatial diffusion ^[10,11], mosquito development period ^[12,13,14], incomplete CI and impaired mating competitiveness of infected males ^[15,16], imperfect maternal transmission and fecundity cost ^[17,18], density dependence ^[19], and so on. Interestingly, by using impulsive and periodic release strategy, Yu et al gave the conditions to guarantee the periodical oscillation of wild mosquito populations ^[14,20,21]. Since the density-dependent competition primarily occurs in the larval stage, stage-structured models, including the aquatic and terrestrial stages are more suitable to describe the mosquito dynamics of wild and suppression populations ^[12,16,22].

We consider a wild mosquito population distributing uniformly in terms of space and gender. Let $L(t)$ and $A(t)$ be the total number of larvae and adults at time $t$, respectively. The target mosquito population is suppressed by releasing a total number $R(t)$ of Wolbachia infected male adults, which induce complete CI in wild females and have equal mating competitiveness with wild males. The field trails of Aedes albopictus population suppression experiment in shazai island, Guangzhou suggested that human activities facilitate mosquito immigration into release sites and compromise the efficiency of Aedes albopictus elimination ^[4]. Besides, theoretical studies also identified that the immigration of fertile female mosquitoes from surrounding uncontrolled areas has potential to seriously hinder the suppression efficiency ^[23,24,25]. In this paper, we use a stage-structured model incorporating larval density dependence to analyze how the immigration of wild fertile female and male mosquitoes from surrounding uncontrolled areas hinders the IIT suppression efficiency by omitting the migration of native mosquitoes from the control area to surrounding areas. Let $D(t)$ be the immigration number of wild female adults from surrounding areas into the target area. Note that there is also number $D(t)$ of wild males immigrating into the control area. Under random mating, the incompatible mating probability of a wild female with an infected male at time $t$ is the ratio of the number $R(t)$ of infected males over the total number $A(t)/2+D(t)+R(t)$ of males in the control area $R(t)/[A(t)/2+D(t)+R(t)]$. Prompting by field experiment studies ^[26,27,28], we assume that wild females immigrate into the control area after mating with wild males. Let $\beta > 0$ be the average number of first instar larvae produced by a female from compatible mating with wild male, and $\tau_1 > 0$ be the average development period from the eclosion of female adults to the hatching of eggs of next generation. Hence, the production rate of larvae at time $t$ is given by

The competition for limited breeding sites and food supply, mostly in larval stage, has shown to be a major factor impairing mosquito growth by delaying development time and elevating mortality rate ^{[29,30,31,32,33]}. We follow the classical logistic model to describe the larval density-dependent competition

where $m > 0$ is the natural mortality rate of larvae, and $K_L > 0$ is a constant characterizing the intensity of density dependence ^{[12,34,35,36]}. As usual, we assume first order stage transition with $\mu\in(0, 1]$ be the pupation rate of larvae and $\alpha\in(0, 1]$ be the eclosion rate of pupae, and first order natural death in adults with $\delta > 0$ be the natural mortality rate of adults. Let $\tau_2 > 0$ be the average development period of larvae and pupae. By combining (1.1) and (1.2), we derive the following delay differential equations:

By letting $s = (m+\mu)t$, and using the following change of variables, we get

supplemented with the conversions of parameters

and we transform (1.3) into the following equations

We consider the compensation release strategy such that the infected males in the control area are maintained almost a constant with $R(t)\equiv R\ge 0$ by replacing the loss of infected males with new release ^[22,37]. Besides, we assume that the immigration number of fertilized females from surrounding areas keeps as a constant with $D(t)\equiv D\ge 0$. By replacing $s$ with $t$, and omitting the overlines in $\overline{D}$, $\overline{R}$, $\overline{\alpha}$, $\overline{\tau}_1$, and $\overline{\tau}_2$ derive

We study the dynamics of (1.7) under the initial conditions

for some fixed time $t_0\ge 0$ and continuous functions $\phi(t)$ and $\psi(t)$ on $[t_0-\tau, t_0]$.

We analyze the global stabilities of (1.7) and (1.8) in Section 2, which are summarized in Theorems 2.1–2.3, and interpret them in terms of the original system parameters in Section 3. Theorem 2.1 reveals that the immigration of fertilized females makes it impossible to completely eliminate the target mosquito populations. Furthermore, Theorems 2.1 and 2.2 identify the threshold immigration number

over which (1.3) has a unique positive equilibrium point and displays globally asymptotically stable dynamics. When $0 < D < D^*$, Theorem 2.1 identifies two threshold release numbers

and

Theorem 2.2 shows that (1.3) has a globally asymptotically stable equilibrium point $E^*(L^*, A^*)$ when $R\ge R^*$, or $0 < D < D^*$ and $0 < R\le R_D$. Otherwise, Theorems 2.2 and 2.3 verify that (1.3) may display bi-stability or global asymptotical stability when $0 < D < D^*$ and $R_D < R < R^*$, depending on the number of positive equilibria. Our simulations show that the combination of small immigration number and moderate release intensity leads to bistable dynamics with one of the stable equilibrium near $A = 0$. Furthermore, we identify the maximum possible suppression efficiency by identifying the infima

of $L^*(R)$ and $A^*(R)$, respectively. $A^*_\infty$ defines the maximum suppression efficiency for wild adults with $A(t) > A^*_\infty$. We use the suppression rate index defined in (3.12) to assess the permitted most immigration number $D_{p_0}$ of wild females to reach a given suppression target $p_0\in (0, 1]$. Besides, we estimate the least release number $R_m(D)$ of infected males to reduce up to $90\%$ of wild adults in the peak season within two months. Our simulations show that $R_m(D)$ increases in the immigration number $D$, and increases near-vertically as $D$ approaches to the most immigration number $D_{0.1}$ of wild females, about $0.38\%$ of the carrying capacity of wild adults in the target area.

2.   Stability analysis of system (1.7)

In this section, we enumerate the non-negative equilibria of (1.7) and study their stabilities. We note that the solution $(x(t), y(t))$ of the initial-value problem (1.7) and (1.8) remains positive and bounded in $[t_0, \infty)$. In fact, if the positivity fails, then there is $t_1 > t_0$ or $t_2 > t_0$ such that $x, y > 0$ in $[t_0, t_1)$ and $x(t_1) = 0$, or $x, y > 0$ in $[t_0, t_2)$ and $y(t_2) = 0$. If the first case occurs, then $x'(t_1)\le 0$. By substituting $t = t_1$ in the first equation of (1.7), and using the nonnegativity of $D$ and $R$, we obtain an obvious contradiction with

Similarly, if the second case occurs, then $y'(t_2)\le 0$ and the second equation of (1.7) gives

The contradictions verify the positivity of $(x(t), y(t))$ for all $t\ge t_0$. The boundedness of $(x(t), y(t))$ can be proved by a similar method as the proof of Lemma 2.1 in ^[22], and we omit it.

2.1.   The enumeration of the equilibria

We note that the model (1.7) degenerates to model $(6)$ in ^[22] in an isolated area without mosquito immigration with $D = 0$. To study the impact of fertile mosquitoes immigrating from surrounding areas on the suppression efficiency, we assume $D > 0$ in the rest of our discussion. As in ^[22], we maintain the following basic assumption

By the definition of $b$ and $\overline{\alpha}$ in (1.5), under the parameters in origin system (1.3), the inequality in the condition (2.1) holds if and only if

which gives a threshold condition ensuring the persistence of isolated populations.

Let $E(x, y)$ be an equilibrium of equations (1.7), which satisfies $y = \alpha x+\rho D$ and $x(1+x)-bD = by(y+D)/(y+D+R)$. Hence, $x$ is a positive root of function

Note that the complete suppression state $E_0(0, 0)$ is no longer an equilibrium of (1.7) when $D > 0$. The number of positive roots of $g(x)$ depends on the signs of its extreme points, satisfying

The discriminant of $g'(x)$ is $4\Delta_g$ with

If $\Delta_g\le 0$, then $g'(x)\ge 0$, and $g(x)$ increases strictly in $x\ge 0$. By combining

for $R\ge 0$ and $D > 0$, we derive that $g\left(x \right)$ has a unique positive root $x^*$, and (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$. If $\Delta_g > 0$, then $g'(x)$ has two roots

In this case, the number of positive roots of $g(x)$ is determined by the signs of the extremums $g(x_1)$ and $g(x_2)$. We classify the equilibria of (1.7) in the following theorem.

Theorem 2.1. Let (2.1) hold. Denote

and $g(x)$ defined in (2.2), $x_1$ and $x_2$ defined in (2.6). Then (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$ when one of the following conditions holds:

(i) $D\ge D^*$;

(ii) $R\ge R^*$;

(iii) $0 < D < D^*$ and $0 < R\le R_D$;

(iv) $0 < D < D^*$, $R_D < R < R^*$, and $g(x_1)g(x_2) > 0$.

Moreover, if $0 < D < D^*$, $R_D < R < R^*$, and $g(x_2) < 0 < g(x_1)$, then (1.7) has three positive equilibria $E^*_1(x^*_1, y^*_1)$, $E^*_2(x^*_2, y^*_2)$ and $E^*_3(x^*_3, y^*_3)$, satisfying

Proof. By using the discriminant of $\Delta_g$

we obtain that $\Delta_D > 0$ if and only if

For the case $D > D_1$, we derive $\Delta_D < 0$, and $\Delta_g > 0$ for any $R\ge 0$. Hence the two roots $x_1$ and $x_2$ of $g'(x)$ satisfy

where

If $D > D_1$ and $R\ge R_D$, then $x_1x_2\ge 0$, and

Hence $x_1+x_2 < 0$, and $x_1 < x_2\le 0$, which implies that $g'(x) > 0$ and $g(x)$ increases strictly in $x > 0$. It follows from (2.5) that $g(x)$ has a unique positive root $x^*$. If $D > D_1$ and $0 < R < R_D$, then $x_1x_2 < 0$, and $g'(x)$ has two opposite-sign roots $x_1 < 0 < x_2$. By combining (2.5), $g'(x) < 0$ for $0\le x < x_2$, and $g'(x) > 0$ for $x > x_2$, we obtain that $g(x)$ has a unique positive root $x^* > x_2$. In summary, (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$ when $D > D_1$ and $R > 0$.

For the case $D = D_1$, we have $\Delta_D = 0$, which implies that $\Delta_g\ge 0$ for all $R\ge 0$, and $\Delta_g$ has a unique positive root

In this case, we have

If $D = D_1$ and $R = R_1$, then $\Delta_g = 0$, which gives $g'(x)\ge 0$ for all $x$, and $g'(x) = 0$ if and only if

Hence $g'(x) > 0$ for all $x\ge 0$. If $D = D_1$, $R\ge R_D$, and $R\neq R_1$, then $\Delta_g > 0$, and the two roots $x_1$ and $x_2$ of $g'(x)$ defined in (2.6) satisfy (2.10). Hence, $x_1x_2 = (R-R_D)/(3\alpha)\ge 0$, and

which lead to $x_1 < x_2\le 0$, and $g'(x) > 0$ for all $x > 0$. Therefore, $g'(x) > 0$, and $g(x)$ increases strictly in $x > 0$, which implies $g(x)$ has a unique positive root $x^* > 0$, for $D = D_1$ and $R\ge R_D$. If $D = D_1$ and $0 < R < R_D$, then $x_1x_2 = (R-R_D)/(3\alpha) < 0$, and $x_1 < 0 < x_2$. Similarly to case $D > D_1$ and $0 < R < R_D$, we have $g(x)$ has a unique positive root $x^* > x_2$. Summarily, (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$ when $D = D_1$ and $R > 0$.

For the case $0 < D < D_1$, we have $\Delta_D > 0$, and $\Delta_g$ has two roots

satisfying

and

for all $D > 0$. Hence $0 < R^* < R_2$. If $0 < D < D_1$ and $R^* < R < R_2$, then $\Delta_g < 0$, and $g'(x) > 0$ for all $x$. If $0 < D < D_1$ and $R = R^*$ or $R = R_2$, then $\Delta_g = 0$, $g'(x)\ge 0$, and $g'(x) = 0$ has a unique solution $\tilde{x} = (\alpha b^*-R-D/\delta)/(3\alpha)$. If $0 < D < D_1$ and $R > R_2$, then $\Delta_g > 0$, and

and

which derive that the two roots of $g'(x)$ satisfy $x_1 < x_2 < 0$, and $g'(x) > 0$ for all $x\ge 0$. Hence, $g(x)$ increases strictly in $x\ge 0$ when $0 < D < D_1$ and $R\ge R^*$. In this case, it follows from (2.5) that $g(x)$ has a unique positive root $x^*$, and (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$.

If $0 < D < D_1$ and $0 < R < R^*$, then $\alpha(2b^*+3)-4\alpha bD/\delta > 3\alpha/2 > 0$, and

Hence $R^*\ge R_D$, and $R^* = R_D$ if and only if

Note that $D_1 > D^*$ by $D_1 = D^*+3\delta/(8b)$.

If $D = D^*$ and $0 < R < R^*$, then $R^* = R_D$, and the two roots $x_1$ and $x_2$ of $g'(x)$ satisfy $x_1x_2 = (R-R_D)/(3\alpha) = (R-R^*)/(3\alpha) < 0$, which implies $x_1 < 0 < x_2$. Similarly, if $0 < D < D_1$, $D\neq D^*$, and $0 < R < R_D$, then $x_1x_2 = (R-R_D)/(3\alpha) < 0$, and $x_1 < 0 < x_2$. If $0 < D < D_1$, $D\neq D^*$, and $R = R_D$, then

Obviously, $x_1 = 0 < x_2$ when $0 < D < D^*$ and $R = R_D$, and $x_1 < 0 = x_2$ when $D^* < D < D_1$ and $R = R_D$. If $D^* < D < D_1$, and $R_D < R < R^*$, then $x_1x_2 = (R-R_D)/(3\alpha) > 0$ and

which imply $x_1 < x_2 < 0$, and $g'(x) > 0$ for all $x\ge 0$. It follows from (2.5) that $g(x)$ has a unique positive root $x^* > 0$, and (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$, when $0 < D < D_1$ and $0 < R\le R_D$, or $D^*\le D < D_1$ and $R_D\le R < R^*$.

For the case $0 < D < D^*$ and $R_D < R < R^*$, we obtain $x_1x_2 = (R-R_D)/(3\alpha) > 0$, and

which imply that $g'(x)$ has two positive roots $0 < x_1 < x_2$. Hence $g'(x) > 0$ when $x < x_1$ or $x > x_2$, and $g'(x) < 0$ when $x\in (x_1, x_2)$, which lead to $g(x_1) > g(x_2)$. By using $R+D/\delta-\alpha b^* < 0$ and $R-(2b^*+1)D/\delta = R-R_D > 0$, the parameter values of $g(x)$ change signs three times. By Descart's Rule of Signs ^[38], $g(x)$ has 3 or $3-2 = 1$ positive roots, which implies $g(x_1)g(x_2)\neq 0$. If $g(x_1)g(x_2) > 0$, then $g(x_1) > g(x_2) > 0$ or $0 > g(x_1) > g(x_2)$. In both cases, using (2.5), we have $g(x)$ has a unique positive root $x^*$, and (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ with $y^* = \alpha x^*+\rho D$. Similarly, if $g(x_1)g(x_2) < 0$, then $g(x_1) > 0 > g(x_2)$, and $g(x)$ has three positive roots $x^*_1 < x^*_2 < x^*_3$ satisfying $0 < x^*_1 < x_1 < x^*_2 < x_2 < x^*_3$. In this case, (1.7) has three positive equilibria $E^*_1(x^*_1, y^*_1)$, $E^*_2(x^*_2, y^*_2)$ and $E^*_3(x^*_3, y^*_3)$, with $y^*_i = \alpha x^*_i+\rho D$ for $i = 1, 2, 3$. □

We use numerical example to show that the polynomial $g(x)$ has one or three positive roots when $0 < D < D^*$ and $R_D < R < R^*$. Let the parameter values $\beta = 2$, $\alpha = 0.95$, $m = 0.1$, $\mu = 0.1$, $\delta = 0.1$, $K_L = 20000$, and $D = 100$. By using the change of variables in (1.4), and the conversions of parameters in (1.5), we derive

If $R = 27000$, then

which satisfy $0 < \overline{D} < \overline{D}^*$, $\overline{R}_D < \overline{R} < \overline{R}^*$, and $g(x_2) < 0 < g(x_1)$. As shown in Figure 1, $g(x)$ has three positive roots $x^*_1 = 0.5755$, $x^*_2 = 0.4895$, and $x^*_3 = 1.7295$. Similarly, if $R = 31000$, then

Hence $0 < \overline{D} < \overline{D}^*$, $\overline{R}_D < \overline{R} < \overline{R}^*$, and $g(x_1)g(x_2) > 0$. As shown in Figure 1, $g(x)$ has a unique positive equilibrium.

Note that (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ when one of the conditions (i)–(iv) of Theorem 2.1 holds. We show that $E^*(x^*, y^*)$ is globally asymptotically stable.

Theorem 2.2. Let (2.1) hold. Then $E^*(x^*, y^*)$ is globally asymptotically stable when one of the following conditions holds:

(i) $D\ge D^*$;

(ii) $R\ge R^*$;

(iii) $0 < D < D^*$ and $0 < R\le R_D$;

(iv) $0 < D < D^*$, $R_D < R < R^*$, and $g(x_1)g(x_2) > 0$,

where $D^*$, $R^*$ and $R_D$ defined in (2.7), $g(x)$ defined in (2.2), $x_1$ and $x_2$ defined in (2.6).

Proof. Theorem 2.1 shows that (1.7) has a unique positive equilibrium $E^*(x^*, y^*)$ when one of the conditions (i)–(iv) holds. In this case, $g(x)$ switches signs from negative in $[0, x^*)$ to positive in $(x^*, \infty)$ with

For any positive constants $c_1$ and $c_2$ satisfying $c_1 < x^* < c_2$, we first claim, if the initial data $\phi(t)$ and $\psi(t)$ satisfy $c_1 < \phi(x) < c_2$ and $\alpha c_1+\rho D < \psi(t) < \alpha c_2+\rho D$ in $[t_0-\tau, t_0]$, then the solution $(x(t), y(t))$ of the initial value problem (1.7) and (1.8) satisfies

Otherwise, let $\bar{t} > t_0$ be the least time such that the solution $(x(\bar{t}), y(\bar{t}))$ reaches the boundary of the rectangular area $[c_1, c_2]\times [\alpha c_1+\rho D, \alpha c_2+\rho D]$ with $x(\bar{t}) = c_1$ or $x(\bar{t}) = c_2$, $c_1 < x(t) < c_2$, and $\alpha c_1+\rho D < y(t) < \alpha c_2+\rho D$ for $t\in[t_0, \bar{t})$, or $y(\bar{t}) = \alpha c_1+\rho D$ or $y(\bar{t}) = \alpha c_2+\rho D$, $c_1 < x(t) < c_2$, and $\alpha c_1+\rho D < y(t) < \alpha c_2+\rho D$ for $t\in[t_0, \bar{t})$. If the first case occurs with $x(\bar{t}) = c_1$, then

which gives

Hence

and $g(c_1) > 0$, which contradict to the fact $c_1 < x^*$ and $g(x) < 0$ in $[0, x^*)$ by (2.15). A similar argument derives that the second case $x(\bar{t}) = c_2$ would not occur. For the third case $y(\bar{t}) = \alpha c_1+\rho D$, it is easy to derive an obvious contradiction

A similar contradiction can be obtained for the fourth case $y(\bar{t}) = \alpha c_2+\rho D$. Hence, the claim (2.16) holds, which verifies the global stability of $E^*(x^*, y^*)$.

By using the fluctuation lemma (Appendix A.5 in ^[36]), there are two increasing and divergent sequences $\{s_n\}$ and $\{t_n\}$ along which

as $n\to \infty$. By substituting $t = t_n$ in the second equation of (1.7) leads to

Taking the limit by letting $n\to \infty$ derives

and $\underline{y}\ge \alpha\underline{x}+\rho D$. Taking the limit of (1.7) along with the sequence $\{s_n\}$ implies

The inequality $\underline{y}\ge \alpha\underline{x}+\rho D$ leads to

and $g(\underline{x})\ge 0$. It follows from (2.15) that $\underline{x}\ge x^*$ and $\underline{y}\ge \alpha x^*+\rho D = y^*$. By repeating the same argument for two divergent sequences, denoted by $\{s_n\}$ and $\{t_n\}$ again, along which $x(s_n)\to \overline{x} = \limsup_{t\to\infty}x(t)$, $x'(s_n)\to 0$, $y(t_n)\to \overline{y} = \limsup_{t\to\infty}y(t)$, and $y'(t_n)\to 0$, as $n\to \infty$, we can prove $\overline{x}\le x^*$ and $\overline{y}\le \alpha x^*+\rho D = y^*$. Taken together, we obtain

which imply that $E^*(x^*, y^*)$ is globally asymptotically stable with $\lim_{t\to\infty}(x(t), y(t)) = E^*(x^*, y^*)$. □

Theorem 2.1 shows that (1.7) has three positive equilibria $E^*_1(x^*_1, y^*_1)$, $E^*_2(x^*_2, y^*_2)$ and $E^*_3(x^*_3, y^*_3)$ when $0 < D < D^*$, $R_D < R < R^*$, and $g(x_1)g(x_2) < 0$. In this case, we show that (1.7) displays bistable dynamics with $E^*_1(x^*_1, y^*_1)$ and $E^*_3(x^*_3, y^*_3)$ being stable and $E^*_2(x^*_2, y^*_2)$ being unstable.

Theorem 2.3. Let (2.1) hold, and $(x(t), y(t))$ be the solution of (1.7) and (1.8). If $0 < D < D^*$, $R_D < R < R^*$, and $g(x_1)g(x_2) < 0$, then the following conclusions hold:

(i) $E^*_1(x^*_1, y^*_1)$ is asymptotically stable, and $\lim_{t\to\infty}(x(t), y(t)) = E^*_1(x^*_1, y^*_1)$ when $0 < \phi(t) < x^*_2$ and $0 < \psi(t) < y^*_2$ on $[t_0-\tau, t_0]$.

(ii) $E^*_2(x^*_2, y^*_2)$ is unstable.

(iii) $E^*_3(x^*_3, y^*_3)$ is asymptotically stable, and $\lim_{t\to\infty}(x(t), y(t)) = E^*_3(x^*_3, y^*_3)$ when $\phi(t) > x^*_2$ and $\psi(t) > y^*_2$ on $[t_0-\tau, t_0]$.

Proof. (i) Since $x^*_i$ for $i = 1, 2, 3$ are the roots of $g(x)$ defined in (2.2), using (2.5), we have

By a similar argument as in the proof of Theorem 2.2, for any positive constants $c_1$ and $c_2$ satsifying $c_1 < x^*_1 < c_2 < x^*_2$, we obtain that the claim (2.16) holds when the initial data $\phi(t)$ and $\psi(t)$ satisfy $c_1 < \phi(t) < c_2$ and $\alpha c_1+\rho D < \psi(t) < \alpha c_2+\rho D$ in $[t_0-\tau, t_0]$. Hence $E^*_1(x^*_1, y^*_1)$ is stable. Furthermore, the claim (2.16) verifies that

and

By using (2.17), and the same argument as in the proof of the second part in Theorem 2.2, we derive $\underline{x} = \overline{x} = x^*_1$, and $\underline{y} = \overline{y} = y^*_1$, which verifies the first part (i).

(ii) The instability of $E^*_2(x^*_2, y^*_2)$ follows directly from the second part of (i).

(iii) For any positive constants $c_3$ and $c_4$ satisfying $x^*_2 < c_3 < x^*_3$ and $c_4 > x^*_3$, if $c_3 < \phi(x) < c_4$ and $\alpha c_3+\rho D < \psi(t) < \alpha c_4+\rho D$ in $[t_0-\tau, t_0]$, we claim

If the claim (2.18) is not true, let $\bar{t}$ be the least time such that the solution $(x(\bar{t}), y(\bar{t}))$ reaches the boundary of the rectangular area $[c_3, c_4]\times [\alpha c_3+\rho D, \alpha c_4+\rho D]$. If $x(\bar{t}) = c_3$, then $c_3 < x(t) < c_4$ and $c_3+\rho D < y(t) < \alpha c_4+\rho D$ for $t\in[t_0, \bar{t})$, and $x'(\bar{t})\le 0$. By substituting $t = \bar{t}$ in the first equation of (1.7) implies

and $g(c_3) > 0$. Since $g(x) < 0$ in $(x^*_2, x^*_3)$ and $x^*_2 < c_3 < x^*_3$, we derive $g(c_3) < 0$, contradicting to the positiveness of $g(c_3)$. A similar contradiction can be obtained for the case $x(\bar{t}) = c_4$. If $y(\bar{t}) = c_3+\rho D$, then $\alpha c_1+\rho D < y(t) < \alpha c_2+\rho D$, $c_1 < x(t) < c_2$ for $t\in[t_0, \bar{t})$, and $y'(\bar{t})\le 0$. By letting $t = \bar{t}$ in the second equation in (1.7) gives

For the case $y(\bar{t}) = \alpha c_4+\rho D$, we can derive a similar contradiction. Hence the claim (2.18) holds, which verifies the stability of $E^*_3(x^*_3, y^*_3)$.

To complete the verification of (iii), for any $x^*_2 < c_3 < x^*_3$ and $c_4 > x^*_3$, the claim (2.18) shows that

and

By using a similar argument as in the proof of the second part in Theorem 2.2, let $\{s_n\}$ and $\{t_n\}$ be two divergent sequences along which $x(s_n)\to \underline{x}$, $x'(s_n)\to 0$, $y(t_n)\to \underline{y}$, and $y'(t_n)\to 0$, as $n\to \infty$. By substituting $t = t_n$ in the second equation of (1.7) gives

Hence, $\underline{y}\ge \alpha\underline{x}+\rho D$. Taking the limit in the first equation of (1.7) along the sequence $\{s_n\}$ leads to

Hence, $g(\underline{x})\ge 0$, which implies $\underline{x}\ge x^*_3$ by (2.17) and $\underline{x} > x^*_2$. Let $\{s_n\}$ and $\{t_n\}$ again be the two divergent sequences such that $x(s_n)\to \overline{x}$, $x'(s_n)\to 0$, $y(t_n)\to \overline{y}$, and $y'(t_n)\to 0$, as $n\to \infty$. Taking the limits along these sequences gives $\overline{y}\le \alpha\overline{x}+\rho D$, and

which imply $g(\overline{x})\le 0$. The fact $\overline{x} > x^*_2$ and $g(x)$ switching signs from negative in $(x^*_2, x^*_3)$ to positive in $(x^*_3, +\infty)$ verify that $\overline{x}\le x^*_3$. The fact $x^*_3\le \underline{x}\le \overline{x}\le x^*_3$ implies $\underline{x} = \overline{x} = x^*_3$. By using the following inequality

we obtain $\underline{y} = \overline{y} = \alpha x^*_3+\rho D$, which verifies (iii). □

3.   Implications in mosquito population suppression

3.1.   Interpreting the theorems in original system parameters

By reversing the change of variables defined in (1.4) with

and $s = (m+\mu)t$, the system (1.7) converts back to the original system (1.3). The initial conditions (1.8) change to

for $t\in [(t_0-\tau)/(m+\mu), t_0/(m+\mu)]$, where $\tau = \max\{\tau_1, \tau_2\}$. The parameters $\overline{D}^*$, $\overline{R}_D$, and $\overline{R}^*$ defined in (2.7) are converted to

and

If $D\ge D^*$, Theorems 2.1 and 2.2 show that the unique positive equilibrium point $E^*(L^*, A^*)$ of (1.3) is globally asymptotically stable. Otherwise, if $0 < D < D^*$, Theorem 2.1 identifies two threshold release numbers $R^*$ and $R_D$. When $R\ge R^*$, or $0 < D < D^*$ and $0 < R\le R_D$, (1.3) has a globally asymptotically stable equilibrium point $E^*(L^*, A^*)$. For the last case $0 < D < D^*$ and $R_D < R < R^*$, Theorem 2.1 shows that (1.3) may has one or three positive equilibria, depending on the signs of the extremums $g(x_1)$ and $g(x_2)$, which displays globally asymptotically stable or bistable dynamics behaviour. Specifically, if $0 < D < D^*$, $R_D < R < D^*$, and $g(x_2) < 0 < g(x_1)$, then (1.3) has three positive equilibria $E^*_1(L^*_1, A^*_1)$, $E^*_2(L^*_2, A^*_2)$ and $E^*_3(L^*_3, A^*_3)$, satisfying

In this case, (1.3) displays bistable dynamics, with $E^*_1$ and $E^*_3$ being asymptotically stable and $E^*_2$ being unstable.

For a population in an isolated area with $D = 0$, (1.3) degenerates to system $(6)$ in ^[22]. As proved in ^[22], there is a threshold release number $R_0$, over which the complete suppression equilibrium $E_0(0, 0)$ is globally asymptotically stable, which suggests that the native mosquito populations in control area can be completely eliminated by releasing many Wolbachia-infected males. Specially, if the mosquito population has not been interfered by releasing with $R = D = 0$, the Theorem 2.2 in ^[22] obtained the carrying capacities of larvae and adults in target area

If $D > 0$, then $E_0(0, 0)$ is no longer an equilibrium of (1.3), which indicates that the immigration of fertilized females from surrounding areas into the control area rules out the possibility of complete eradication. Under this case, to prevent and control dengue fever effectively, we should reduce the wild mosquitoes to a low level below which dengue fever could not cause an outbreak.

In this section, we use Aedes albopictus populations as an example to further discuss the impact of wild mosquito immigration on the suppression efficiency. Based on laboratory data and field data, we have estimated the life table parameters of Aedes albopictus population in Guangzhou ^[12,22], which are listed in Table 1. Note that the life table parameters are sensitive to climatic conditions such as temperature and precipitation. The hot and rainy summer in Guangzhou is very suitable for the breeding of Aedes albopictus mosquitoes. After the hatching of diapause eggs from early March, the mosquito density peaks in late September and early October, which overlaps with the high-risk period of dengue fever ^[27,29,30].

To display our discussion more specific and transparent, we fix the system parameters in our simulations as follows:

Note that $K_L$ scales the size of control area and is determined mainly by environmental conditions, such as the availability of breeding sites, and resource competition ^[43,44]. We take $K_L = 20000$ as an example, since $K_L$ has little impact on the suppression dynamics. All other parameter values are within the range of parameter values listed in Table 1. By substituting the parameters specified in (3.6) into (3.3) and (3.5), we derive

To display the rich dynamics of the system (1.3), we take $D = 100 < D^*$ as an example. By using the definition in (3.3) and (3.4), and the estimated parameters in (3.6), we obtain

In Figure 2A, the curves correspond to the number $A(t)$ of wild adult mosquitoes for different release numbers and initial data. As we expected, the curves with identical $R$ and differential initial data verify the global asymptotical stability of the unique positive equilibrium $E^*(L^*, A^*)$. For the case $R = 40,000 > R^*$, the system (1.3) has a unique positive equilibrium $E^*(1650, 3368)$. For both initial data $\phi(t) = \psi(t) = 1000$ (red dotted curve) and $\phi(t) = \psi(t) = 220,000$ (black dotted curve) on $t\in [-12, 0]$, the correspond number $A(t)$ of adults converge quickly to the equilibrium $A^*(R) = 3368$, about $2.4\%$ of the carrying capacity $A^*_0$ of adults. The suppression dynamics for the case $R < R_D$ is similar to the above case $R > R^*$. For instance, if $R = 8000 < R_D$, our simulations show that (1.3) has a unique positive equilibrium $E^*(131460,126,679)$. Both of the adult mosquito numbers converge steeply to $A^* = 126,679$ for $\phi(t) = \psi(t) = 1000$ (greet dotted curve) and $\phi(t) = \psi(t) = 220,000$ (blue dotted curve) on $t\in [-12, 0]$. Furthermore, for the case $R_D < R < R^*$, the system (1.3) may has one or three positive equilibria. For instance, for a relatively small release number $R = 20,000\in (R_D, R^*)$, we simulations show that $x_1 = 0.123$, $x_2 = 1.6401$, $g(x_1) = -0.006$, and $g(x_2) = -3.3234$, satisfying $g(x_1)g(x_2) > 0$. In this case, (1.3) has a unique positive equilibrium $E^*(95966, 92975)$. Both of the curves corresponding to the number $A(t)$ of wild adults for $\phi(t) = \psi(t) = 1000$ (cyan) and $\phi(t) = \psi(t) = 220,000$ (pink) on $t\in [-12, 0]$ converge to $A^* = 92,975$.

On the other hand, for a relatively large release number $R = 27,000\in (R_D, R^*)$, our simulations show that $x_1 = 0.2576$, $x_2 = 1.2599$, $g(x_1) = 0.1295$, and $g(x_2) = -0.8272$, which satisfy $g(x_2) < 0 < g(x_1)$. Theorems 2.1 and 2.3 verify that (1.3) has three positive equilibria $E^*_1(2301, 3987)$, $E^*_2(20512, 20150)$, and $E^*_3(69191, 67520)$, which display bistability dynamics with $E^*_1$ and $E^*_3$ being asymptotically stable. Moreover, $E^*_2$ is unstable. In Figure 2B, the curves correspond to the number $A(t)$ of wild adults for the identical initial number of larvae and adult on $t\in [-12, 0]$ with $\phi(t) = \psi(t) = 100$ (red), $20,000$ (blue), $22,000$ (green), and $120,000$ (black), respectively. As we expected, all the four curves show that the number $A(t)$ of wild adults converges to $A^*_3 = 67,520$ for the initial data $(\phi(t), \psi(t))$ on $t\in [-12, 0]$ above the instable equilibrium $E^*_2(20512, 20150)$, while $A(t)$ converges to $A^*_1 = 3987$ for $(\phi(t), \psi(t))$ below $E^*_2$. Our simulations show that a moderate release number $R$ with small immigration number $D$ can lead to bistable dynamics in adult abundance. The separate increase of release number $R$ or immigration number $D$, or concurrent decrease of $R$ and $D$ can change the bistable dynamics to globally asymptotically stable dynamics, by reducing the positive equilibria from three to one.

3.2.   Maximum possible suppression efficiency

We consider the maximum possible suppression efficiency due to the release of Wolbachia-infected male mosquitoes, which is obtained formally by letting $R\to\infty$ ^[25]. In this case, the native females in the control area are rendered sterile. For the case $D > 0$, we use $L^*(R)$ and $A^*(R)$ to denote the dependence of the unique positive equilibrium states of larval and adult mosquitoes on $R\ge 0$, respectively. In Figure 3A, the curves correspond to the unique positive equilibrium number $A^*(R)$ of adults for $D = 400$ (black), $700$ (red), and $900$ (blue), as $R$ increases from $0$ to $80,000$. It is seen that, as we expected, $A^*(R)$ decreases strictly in the release number $R$ for fixed $D$. For the case $D = 900$, about $0.63\%$ of $A^*_0$, our simulations show that $A^*(R)$ decreases almost linearly in $R$. For instance, $A^*(R)$ reduces to $1.035\times 10^5$ (about $72.63\%$ of $A^*_0$) when $R = 4\times 10^4$, and reduces to $7.299\times 10^4$ (about $51.22\%$ of $A^*_0$) when $R = 6\times 10^4$. A similar simulation shows that $L^*(R)$ also decreases strictly in $R$ for fixed $D$. On the other hand, the immigration of fertilized females contributes to rapid population recovery and hinders control efforts. It is seen that both of $L^*(R)$ and $A^*(R)$ increases strictly in the immigration number $D$ of fertilized females for fixed $R$, as shown in Figure 3B. We take $R = 1.2\times 10^5$, about $84.2\%$ of $A^*_0$, as an example. Our simulations show that $A^*(R) = 6098$ (about $4.3\%$ of $A^*_0$) when $D = 200$, which increases to $13576$ (about $9.5\%$ of $A^*_0$) when $D = 400$, and increases to $17867$ (about $12.5\%$ of $A^*_0$) when $D = 600$.

Since both of $L^*(R)$ and $A^*(R)$ have lower bounds for fixed $D > 0$, by monotone and boundedness theorem, their infima $L_*$ and $A_*$ are the limits of $L^*(R)$ and $A^*(R)$ as $R\to+\infty$, respectively, with

In fact, if the release number $R$ is much greater than $A(t)+2D$, then the solution $(L(t), A(t))$ of system (1.3) can be approximated by the following simple system

with the same initial data (3.2). Note that the limit system (3.8) has a unique positive equilibrium $E^*_{\infty}(L^*_{\infty}, A^*_{\infty})$ with

It is easy to show that $E^*_{\infty}(L^*_{\infty}, A^*_{\infty})$ is globally asymptotically stable by a similar argument as in the proof of Theorem 2.2. Denote the solutions of (1.3) and (3.8) by $(L(t), A(t))$ and $(L_{\infty}(t), A_{\infty}(t))$ for $t\ge t_0$, respectively. Hence

Combining (3.7) and (3.10) gives $(L_*, A_*) = (L^*_{\infty}, A^*_{\infty})$, and

Hence, $L^*_{\infty}$ and $A^*_{\infty}$ define the maximum possible suppression efficiencies of larval and adult mosquitoes, respectively. The limit in (3.11) allows us use $L^*_{\infty}$ and $A^*_{\infty}$ to approximate the sizes of $L^*(R)$ and $A^*(R)$ for large $R$. The above analysis shows that it is impossible to bring down the number of wild adults in the peak season to a level below the maximum possible suppression efficiency $A^*_\infty$.

Figure 4 displays the estimation of $(L^*(R), A^*(R))$ by $(L^*_\infty, A^*_\infty)$. The curves in Figure 4 correspond to the maximum possible suppression efficiency $A^*_{\infty}$ (black), and the equilibrium number $A^*(R)$ of native adult mosquitoes for $R = 10^7$ (red), $5\times 10^6$ (blue), and $10^6$ (pink), as $D$ increases from $0$ to $5000$. It is seen that $A^*(R)$ can be approximated perfectly by $A^*_\infty$ for fixed $D > 0$, and a better estimation is obtained under larger value of the release number $R$. Let the parameters are specified in (3.6). We take $D = 500$, about $0.35\%$ of the carrying capacity $A^*_0$, as an example. Our simulations show that $A^*_\infty = 13,270$ and $A^*(R) = 13,650$ when $R = 10^6$, about $6.7$ times of $A^*_0$. If the release number $R$ increases $5$ times, then the relative error of $A^*(R)$ decreases from $2.78\%$ to $0.56\%$ compared with the theoretical value $A^*_\infty$. Furthermore, as $R$ increases $10$ times to $10^7$, the relative error decreases to $0.29\%$. Similarly, the positive equilibrium number $L^*(R)$ of larvae can be estimated perfectly by $L^*_\infty$ defined in (3.9) for a sufficiently large release number $R$.

3.3.   Further quantification of the hindrance by immigration

To further quantify the hindrance of mosquito immigration on suppression efficiency, we define the suppression rate $p(t)$ as the ratio of the wild adult number $A(t)$ of suppressed population with immigration over the carrying capacity $A^*_0$ of adults in an isolated area

If $D = 0$, The Theorem 2.5 in ^[22] showed that there is a threshold release number

such that

when $R\ge R_0$. By using Theorems 2.2 and 2.3, we derive that $p(t) > 0$ when $D > 0$. Furthermore, the limits in (3.7) and (3.11) derive

The decrease of $A^*(R)$ in $R$ indicates that $p^*$ is the infimum of the suppression rate $p(t)$. Note that $p^*$ increases strictly in $D$ and

For a given suppression rate target $p_0\in (0, 1]$, there is a unique threshold immigration number, denoted by $D_{p_0}$, such that $A^*_{\infty}/A^*_0\le p_0$ if and only if $D\le D_{p_0}$. We call $D_{p_0}$ the permitted most immigration number $D$ such that $A^*_{\infty}\le p_0 A^*_0$.

We take $p_0 = 0.1$ as an example. By the parameters specified in (3.6), our simulations derive $D_{0.1} = 537$, about $0.377\%$ of the carrying capacity $A^*_0$ of wild adults in target area. If the immigration number $D$ of wild females is larger than the permitted most immigration number $537$, the wild adults can not be reduced to a level less than $A^*_0/10 = 14,250$. In Figure 5A, the curves correspond to the number $A(t)$ of wild adult mosquitoes for $D = 550$ and $R = 10^{10}$, under the initial data $\phi(t) = \psi(t) = 12,000$ (blue) and $\phi(t) = \psi(t) = 20,000$ (red) on $[-12, 0]$. It is seen that, for both cases with the initial data less than and larger than the suppression target $14,250$, the number of wild adults are always larger than this target ultimately, even though the release number is as large as $R = 10^{10}$. For instance, under the initial values $\phi(t) = \psi(t) = 12,000$ of larvae and adults, the number of wild adults raise quickly, and stabilize around $A^*(10^{10}) = 1.456\times 10^4$ in $60$ days, about $10.2\%$ of $A^*_0$, which is larger than the suppression target of up to $90\%$ reduction in wild adults in the peak season. Our analysis suggests that the immigration number $D$ of wild females hinder largely the population suppression, and restricts the most suppression efficiency. To reduce the wild adult mosquitoes up to $90\%$ in the peak season, the total immigration number of the fertile males and females from surrounding areas should be less than $2D_{0.1} = 1074$ per day, about $0.754\%$ of $A^*_0$.

A desirable and feasible dengue control strategy is to reduce rapidly the number of wild adult mosquitoes to a low level within finite time to prevent disease transmission. For a given suppression rate target $p_0 = 0.1$ and immigration number $0 < D < D_{0.1}$, we further estimate the required least release number, denoted by $R_m(D)$, to reduce the number of wild adults at peak season to a level $A^*_0/10$ within 60 days. Let the parameters be specified in (3.6). The high-incidence season of dengue fever in Guangzhou coinciding with the peak season of Aedes albopictus populations promotes us to consider the mosquito peak season as the beginning of suppression such that the wild mosquitoes stabilize around these carrying capacities ^[4,27]. We take the initial data $\phi(t) = L^*_0 = 150,000$ and $\psi(t) = A^*_0 = 142,500$ on $[-12, 0]$. In Figure 5B, the curves correspond to the number of wild adult mosquitoes $A(t)$ for $D = 0$ (blue), $150$ (green), $300$ (red), and $500$ (black), and estimated least release number $R_m(D) = 499,424$, $R_m(D) = 723,397$, $R_m(D) = 1,276,938$, and $R_m(D) = 108,999,526$, respectively. It is interesting to see that the four curves display similar suppression dynamics. The number of wild adults $A(t)$ decreases sharply to the target level of $10\%$ of $A^*_0$ after increasing to its peak larger than $A^*_0$ in about $10$ days when $D > 0$. For instance, if $D = 300$, our simulations show that the estimated least release number is $R_m(D) = 1,276,938$, about $9$ times of $A^*_0$. The number of wild adults firstly increases to its peak $145,500$ in $8$ days, about $102\%$ of $A^*_0$, then decreases steeply to the target of $14,250$ in the rest of $52$ days. Furthermore, the least release number $R_m(D)$ increases strictly in the immigration number $D$. The dependence of $R_m(D)$ on $D$ is characterized by the curve in Figure 5C. The moderate increase of $R_m(D)$ in $D$ is almost linear when $D\le D_{0.1}/2$. For instance, the least release number $R_m(D) = 584,990$, about $4.1$ times of $A^*_0$, when $D = 70$, which increases about $1.2$ times to $702,814$ when $D = 140$, and increases about $2$ times to $1,160,186$ when $D = 280$. However, $R_m(D)$ increases sharply when $D > D_{0.1}/2$, and increases near-vertically as $D$ approaches to the most immigration number $D_{0.1} = 537$. For instance, $R_m(D) = 6,731,344$, about $47.2$ times of $A^*_0$, when $D = 462$, which increases about $1.5$ times to $10,321,783$ when $D = 476$, and increases about $3.3$ times to $22,360,512$ when $D = 490$, about $156.9$ times of $A^*_0$. To reduce up to $90\%$ of wild adults in the peak season within two months, an economically viable strategy is to reduce the immigration number of wild females less than $300$, about $0.21\%$ of the carrying capacity of adults $A^*$ in the control area. Setting buffer zones around the control area may be a viable strategy to reduce the immigration number into the control area, and raise the suppression efficiency.

3.4.   Discussion

The incompatible insect technique based on the endosymbiotic bacterium Wolbachia has been proved to be a promising avenue to control mosquitoes, the vectors for mosquito-borne diseases, such as dengue fever, malaria, and Zika. However, both of theoretical studies and field experiments verify that the immigration of fertilized females from surrounding areas rules out the possibility of complete mosquito eradication, and compromises the suppression efficiency ^[4,23,24,25]. It has been observed that the density induced intra-specific competition occurs mostly in the larval stage, which is found to be the major determinant for the mosquito population growth that elevates mortality rates, delays development period, and influences the female size and fecundity ^{[29,30,31,32,33]}. Stage-structured models including the aquatic and terrestrial stages with their corresponding development periods are considered to be more suitable to model the mosquito dynamics. In this paper, we modestly attempt to discuss the joint impacts of the release of Wolbachia-infected male adults and the immigration of fertile females and males on the suppression efficiency using a framework of delay differential equations integrated larval density-dependent competition. We classify the release number of infected males and immigration number of fertilized females, to ensure that the system of delay equations display globally asymptotically stable or bistable dynamics.

Our theoretical results show that the suppression efficiency is more susceptible to the immigration of fertilized females than the release of infected males. The immigration of fertilized females makes the complete eradication becomes impossible, since the immigration of fertilized females are immune to the releases of Wolbachia-infected males. The suppression efficiency, characterized by the equilibrium state $A^*(R)$ of wild adults, decreases strictly in the release number $R$, and increases strictly in the immigration number $D$. The immigration of fertile females restricts the maximum possible suppression efficiency such that the wild adult mosquitoes in target area can not be reduced to a level below $A^*_\infty$ defined in (3.9). Furthermore, our theoretical results show that a small number of immigration and moderate release can lead to bistable dynamics with one of the three positive equilibria $E_1(L^*_1, A^*_1)$ being asymptotically stable and $A^*_1$ being near $A = 0$. An economically viable suppression strategy is to set buffer zones around the control area to minimize the entry of external fertile adults, and treat first the target mosquito populations by combining other control tools, including insecticide spraying, to bring down the number of wild mosquitoes to a lower level below $A^*_1$. In the buffer zones, the sustained release of Wolbachia-infected males may be a viable approach to reduce the immigration number $D$ of fertile adult mosquitoes into the control area. To bring down the wild adults in the target area up to $90\%$ in the peak season within two months, our simulations show that the immigration number of fertilized females should be less than $0.38\%$ of the carrying capacity of wild adults $A^*_0$. To reach the above suppression target in two months, the required least release number increases near-vertically as $D$ approaches to the permitted most immigration number $0.38\% A^*_0$.

The mathematical modeling framework presented in this paper is based on a series of simplifying assumptions by ignoring the impact of environment factors such as temperature and precipitation, and the heterogeneity of mosquito spatial distribution. In most previous studies, the authors considered the impact of immigration fertilized females on suppression dynamics by ignoring the immigration of fertile males ^[23,25]. In our model (1.3), we consider the joint impact of the immigration of fertile females and males from surrounding areas with one-to-one sex ratio on the suppression dynamics. The mathematical framework to model the larval density dependence, which is important in characterizing the population dynamics by limiting the carrying capacity of a population in the target area, relies on the accuracy of the data set and fitting methods ^[25,34]. The development period and survival rate of the aquatic and terrestrial stages are sensitive to climatic factors, especially temperature ^[27,29,30]. In this paper, different from ^[23,25], we follow the idea of the classical logistic model to describe the density-dependent mortality of larvae ^[35,36]. Although the larval density-dependent competition is considered to primarily prolong the larval development period, theoretical analyses have proved that both of the development periods of the aquatic and terrestrial stages have little impacts on the mosquito suppression dynamics ^{[12,15,22,37]}. Furthermore, to make the model mathematically analyzable, we ignore the impact of climatic factors such as temperature on the life-table parameters and simplify them as constants, which could partly give an instructive insights on population suppression. To deeply understand the mosquito suppression dynamics, experimental and theoretical research are required to further understand the density dependence mechanism and the impact of climatic factors on the life-table parameters of mosquito populations. Furthermore, experiment studies on Aedes albopictus populations have showed that some Wolbachia strains may cause incomplete CI ^[4,33]. Researchers should focus on the joint impact of incomplete CI and immigration of fertile females and males from surrounding areas on the suppression dynamics by utilising the present framework.

Our estimation for the least release number $R_m(D)$ increases strictly in the immigration number $D$ of fertile females in nonlinear form to reduce up to $90\%$ of wild adult mosquitoes in the peak season within two months. It is important in the practice of mosquito suppression to estimate exactly the immigration number of wild females from surrounding areas into the target area, which may not be straightforward. The immigration number is impacted by environmental factors and human activities such as frequent traffic and the flow of people ^[4]. Our model framework provides useful suggestions in deciding the least release number of infected males to obtain a given suppression target within finite time period, by reducing the migration number to a lower level. Setting buffer zones around the control area and releasing infected males in the buffer zones may be a viable strategy. Our model focuses on the impact of the immigration of wild mosquitoes from surrounding areas and ignores the migration of wild mosquitoes from the target area to surrounding areas. Further studies should focus on the joint impact of the immigration and migration of wild mosquitoes on the suppression dynamics, since the migration of wild mosquitoes from the target area may impact the suppression dynamics, especially in the early stage of suppression.

Use of AI tools declaration

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

Conflict of interest

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