Monotonicity properties arising in a simple model of Wolbachia invasion for wild mosquito populations

1.   Introduction

It is widely known that female mosquitoes of the species Aedes aegypti are major transmitters of dengue and other vector-borne infections. When deliberately infected with Wolbachia, they lose their vector competence by becoming far less capable of developing a viral load sufficient for transmission of the virus to humans. Due to this remarkable feature, Wolbachia-based biocontrol of mosquito populations has recently emerged as a novel method for the prevention and control of vector-borne infections and is accepted as an ecologically friendly and potentially cost-effective method ^[1,2,3,4,5].

The goal of Wolbachia-based biocontrol is the eventual elimination of wild insects (capable of transmitting the virus to humans) by performing single or periodic releases of Wolbachia-carrying mosquitoes in some determined localities initially populated by wild mosquitoes. The practical implementation of this method requires the mass-rearing of a large quantity of Wolbachia-infected insects for posterior releases, and the desired result depends on the progressive Wolbachia invasion and its durable establishment in wild mosquito populations. The final outcome of this process is usually referred to as "population replacement" ^[6,7,8].

In the literature from the last decade, one can find a variety of Wolbachia infestation models created with the purpose, among others, of evaluating biocontrol strategies that seek suppression of wild mosquito populations. Previous studies include some frequency-based models formulated as a single equation ^[9,10,11], models considering only female insects ^{[6,12,13,14,15]}, models accounting for sex structure ^{[16,17,18,19]} or stage structure ^{[20,21,22,23]}, and more sophisticated models that take into account the mosquito maturation delay or their age structure ^[18,19,24].

The main features of all these models describing the natural dynamics of Wolbachia are related to mimicking vertical transmission and the interference of the reproductive outcomes induced by cytoplasmic incompatibility (or CI reproductive phenotype). In this context, it is meaningful to note that CI occurs when a female, uninfected by Wolbachia, is inseminated by an infected male, producing inviable eggs. Thus, the CI reproductive phenotype suppresses the growth of the uninfected population and facilitates Wolbachia spread.

Different strains of Wolbachia may induce either perfect (100%) or imperfect (less than 100 %) maternal transmission and CI reproductive phenotype. Notably, some of the existent models are designed to address imperfect maternal transmission and/or imperfect CI ^{[6,9,17,18,20,21,23,24]}. Nevertheless, laboratory trials evince that Aedes aegypti mosquitoes that were deliberately infected by the wMelPop Wolbachia strain (regarded as the best blocker of arboviral infections ^{[3,5,25,26,27,28]}) exhibit almost perfect CI and maternal transmission ^[25,29]. Table 1 summarizes the results of matings between infected and uninfected mosquitoes when the maternal transmission and CI are perfect (100%). On the other hand, many scholars also indicate that the wMelPop strain is associated with high "fitness costs" since it reduces female fecundity, the viability of eggs, and the lifespan of infected mosquitoes ^[7,8,10,25]. The latter makes the spread of wMelPop Wolbachia infection a rather challenging task.

The group of models assuming perfect CI and maternal transmission of Wolbachia exhibit a bistable nature that makes them fully compliant with the so-called principle of competitive exclusion ^[30]. This phenomenon implies the existence of a certain (dynamic) threshold in the current level (or frequency) of infection above which Wolbachia is capable of invading and persisting in the uninfected population and below which the Wolbachia-carrying population is driven toward extinction.

For each infection frequency, this threshold can be expressed in terms of the so-called minimal viable sizes of each population (infected and uninfected) that are tightly related to the frequency-dependent Allee effect. Anticipated knowledge of the minimal viable size of the Wolbachia-infected population corresponding to the current size of the wild population is the key issue for determining the appropriate size of release(s) for implementation of Wolbachia-based biocontrol, and one of the main goals of the present work consists of assessing the release size(s) through the use of a simplified bidimensional model of Wolbachia invasion.

In effect, many scholars have intuitively detected the presence of the aforementioned threshold in the infection frequency (see, e.g., ^{[6,10,11,12,14,15,16,17,18,19,22,23]}). However, none of these works were focused on explicit identification of the minimal viable sizes of infected and uninfected populations nor their useful significance for the practical implementation of Wolbachia-based biocontrol. The present paper intends to contribute to this strand of research and corroborate the result obtained in ^[31] based on a different model.

All models describing Wolbachia invasion are competitive; however, only a few of them rigorously exhibit the property of monotonicity, namely:

$–$ the four-dimensional stage-structured model proposed in ^[22];

$–$ the bidimensional model that is further reduced to a single equation in terms of the infection frequency ^[6];

$–$ the reduced bidimensional version ^[31] of the sex-structured model initially proposed in ^[16]*;

^*Notably, the four-dimensional sex-structured model developed in ^[16] is not monotone, while its reduced bidimensional version studied in ^[31] clearly exhibits the property of monotonicity.

$–$ the simplified version without maturation delay considered in ^[19].

This important property makes the theoretical analysis of these models simpler by making use of the variety of research results developed for monotone dynamical systems and assembled in ^[32]. In particular, the application of the theory of monotone systems to competitive dynamics allows identifying a partial order under which a competitive monotone system exhibits the so-called saddle-point behavior ^[33], and the latter bears a strong relationship to the principle of competitive exclusion ^[30]. In general terms, a system exhibiting the saddle-point behavior possesses two locally stable equilibria on the boundary and one unstable (saddle-point) equilibrium in the interior of the state domain. Moreover, the state domain is divided into three disjointed and invariant parts: two attraction basins of boundary equilibria and the so-called "threshold manifold" (or separatrix) containing the unstable equilibrium that separates the attraction basins. One of the goals of the present work is to identify the threshold manifold and study its properties in light of the practical implementation of Wolbachia-based biocontrol.

For that purpose, we have developed a simplified bidimensional variant of the four-dimensional model presented in ^[22] that retains all the key properties of the original four-dimensional model, including the property of monotonicity. Moreover, the model studied in ^[22] does not require the anticipated knowledge of the carrying capacity parameter (the maximum number of individuals supported by the environment), in contrast to the model proposed in ^[6].

Notably, our reduced model bears a certain degree of similarity to the models studied in ^[19] and ^[31], but they are not the same. The bidimensional model obtained in ^[19] by omitting the maturation delay assumes only the density-dependent mortality of both mosquito populations and ignores their natural mortalities, which differ for infected and uninfected insects. In contrast, our bidimensional model accounts for both mortality types (natural and density-dependent). Furthermore, the authors of ^[19] have, in effect, identified the threshold manifold for their bidimensional model, but they did not discuss its underlying properties in light of the practical implementation of Wolbachia-based biocontrol.

In its turn, the bidimensional model studied in ^[31] is the fruit of the reduction of a four-dimensional sex-structured Ricker-type model encompassing the evolution of only adult insects. Therefore, the Ricker-type model presented in ^[31] is structurally different from our model and the one presented in ^[19]. Namely, the density-dependence in ^[31] is included in the recruitment term (which is natural for Ricker-type models) and expresses the insect survival to the adult stage. In contrast, our model and the one studied in ^[19] exhibit the density-dependent mortality of the insect population attributed to inter- and intra-species competition at the aquatic stage.

Our main motivation in proposing a bidimensional (or planar) dynamical system is that, in general, such dynamical systems have several important advantages compared to higher-dimensional systems^†. First, they concede a comprehensive visualization of the system's behavior in the phase plane that not only facilitates the conceptual theoretical analysis of the model but also provides meaningful interpretations of the potential outcomes of the model. Second, for planar dynamical systems, there are numerous optimal control tools ^{[34,35,36,37]} that can be applied to the model proposed in this paper to obtain an analytically optimal synthesis of various optimization problems related to biocontrol purposes. Such analytical solutions can be later tested in more complex models. Although we do not analyze optimal control problems in this work, the proposed model has been designed to serve these purposes in the future and to provide major insights regarding the evolution of both mosquito populations under the action of biocontrol. In addition, as we discuss in the final section, the properties of the proposed model enable the impulsive releases of Wolbachia-infected mosquitoes in order to achieve the population replacement. Notably, in ^[38] the author proposes to employ feedback controls for reaching the steady state representing the population replacement. Motivated by that work and our exercise carried out with impulsive releases (Section 5) we believe that our proposed model is well-suited for determining optimal feedback in analytical form.

^†In this context, it is worth pointing out that planar systems are less precise in describing the evolution of populations than higher-dimensional models.

The paper is organized as follows. In Section 2, we introduce a bidimensional Wolbachia infestation model for populations of Aedes aegypti mosquitoes. In Section 3, we establish some basic properties of the introduced model. In Section 4, we prove that the reduced model retains the property of monotonicity inherent from the original four-dimensional model developed in ^[22]. Using this remarkable property and other tools from the theory of monotone dynamical systems, we establish the existence of an invariant threshold manifold for the proposed bidimensional model. Finally, in Section 5, we discuss the core properties of the points located on the threshold manifold and propose their practical interpretations for performing single and periodic releases of Wolbachia-infected mosquitoes to reach an eventual elimination of wild insects and thus achieve the population replacement.

2.   Simplified model of Wolbachia infestation

Let us consider two populations of mosquitoes, $P_N(t): = F_N(t) + M_N(t)$ and $P_W(t): = F_W(t) + M_W(t)$ present at the day $t \geq 0$ in some locality, where $F_N(t)$ and $M_N(t)$ stand, respectively, for the density of female and male insects that are free of Wolbachia symbiotic bacterium while $F_W(t)$ and $M_W(t)$ denote, respectively, the density of female and male insects infected with Wolbachia bacterium.

Scientific evidence ^[39,54] suggests that wild female and male mosquitoes are often evenly distributed; therefore, let us suppose that $N(t) : = F_N(t) = M_N(t)$ for all $t \geq 0$. A similar assumption can be introduced, for the sake of simplicity, regarding Wolbachia-carrying mosquitoes, that is, $W(t) : = F_W(t) = M_W(t)$ for all $t \geq 0$. Then, the frequency of Wolbachia infection in the total mosquito population, $\dfrac{P_W(t)}{P_N(t) + P_W(t)}$ can be determined by $N(t)$ and $W(t)$ as $\dfrac{W(t)}{N(t) + W(t)}$.

To propose a simplified model of Wolbachia invasion, we have chosen as a starting point the stage-structured model of Wolbachia infestation with four state variables developed by Bliman et al ^[22], and our final goal is to design a bidimensional (reduced) version of this model with similar characteristics. The original four-dimensional model ^[22] accounts for the larval (or aquatic) stage; therefore, it is worth recalling that immature mosquitoes undergo inter- and intra-species competition and exhibit density-dependent mortality at this stage. Thus, our reduced bidimensional version of the model should bear density-dependence in the mortality term, in contrast to another bidimensional model ^[31] where density-dependence is included in the recruitment term.

Let us err on the side of caution while reducing the model's dimension and recall that a "good" model for describing Wolbachia invasion with perfect maternal transmission and CI must necessarily account for the following features:

(ⅰ) Maternal transmission of the bacterium Wolbachia to the next generation. This feature implies that Wolbachia-infected mosquitoes are progenies of Wolbachia-carrying females. A wild female cannot produce Wolbachia-infected offspring (see Table 1).

(ⅱ) The reproductive phenotype of cytoplasmic incompatibility (CI). This feature implies that Wolbachia-carrying females are capable of producing viable and Wolbachia-infected offspring after mating with either wild or Wolbachia-carrying males. On the other hand, wild females produce inviable offspring after mating with Wolbachia-carrying males.

(ⅲ) Positive invariance, well-posedness, and the existence of the absorbing set that attracts all the system trajectories. Any model is "biologically viable" if its state variables are nonnegative for all $t \geq 0$ and their underlying trajectories are bounded.

(ⅳ) Bistable nature. This feature expresses the so-called principle of competitive exclusion ^[30] according to which only one mosquito population (either with or without Wolbachia) should ultimately survive and persist.

In mathematical terms, feature (ⅳ) implies the presence of two local attractors (boundary equilibria) and an unstable strictly positive coexistence equilibrium in between. From the biological standpoint, the last feature (ⅳ) is also directly related to the so-called frequency-dependent Allee effect and implies the existence of a certain threshold in the frequency of Wolbachia infection above which the wild population is eventually driven toward extinction, and Wolbachia successfully invade the mosquito population. Such a threshold is usually referred to as the "minimal viable population size" of wild mosquitoes.

In addition to key features given by (ⅰ)–(ⅳ), the four-dimensional model developed in ^[22] possesses another important property related to its monotonicity. Namely, the semiflow associated with that four-dimensional model is monotone and strongly order-preserving for the partial order induced by the cone $\mathcal{K}: = \mathbb{R}_{-} \times \mathbb{R}_{-} \times \mathbb{R}_{+} \times \mathbb{R}_{+}$. Notably, there are several models of Wolbachia invasion with perfect maternal transmission and CI, both four-dimensional ^[14,16,17] and bidimensional ^[12,13], that do not exhibit such an important property. Therefore, it is highly desirable to conserve this property in the reduced bidimensional version of the original model developed in ^[22].

To reduce the four-dimensional model developed in ^[22] to a bidimensional scenario, the aquatic and aerial stages of each mosquito class (with and without Wolbachia) can be "merged" into one state variable denoted by $N(t)$ and $W(t)$ for noninfected and Wolbachia-infected insects, respectively. As a result, we obtain the following model without the stage structure:

that is bidimensional and describes the time evolution of adult mosquito populations, $N(t)$ and $W(t)$.

Let us provide a brief description of positive constant parameters included in the model (2.1):

● The parameters $\rho_N$ and $\rho_W$ represent the fecundity rates of uninfected and Wolbachia-infected insects, respectively, in the absence of competition (i.e., a mean number of adult mosquitoes produced by one female on average per day during her lifespan).

● The parameters $\alpha_N$ and $\alpha_W$ refer to a natural mortality rate of uninfected and Wolbachia-infected insects, respectively (note that $1/\alpha_N$ and $1/\alpha_W$ express the average lifespan of noninfected and Wolbachia-infected mosquitoes).^‡

^‡Notably, system (2.1) bears a certain degree of similarity to the one studied in ^[19]. However, system (2.1) explicitly accounts for the natural mortalities, $\alpha_N$ and $\alpha_W$, of wild and Wolbachia-infected mosquitoes while the model in ^[19] ignores them and only considers the density-dependent mortalities (expressed by the last terms in both equations of (2.1)).

● The parameters $\beta_N$ and $\beta_W$ are associated with the competition between two mosquito populations for food resources, breeding sites, and mating opportunities, including larvae development into adults under density dependence.

To guarantee survival and persistence of each population in the absence of another, the following conditions are introduced:

It is worth pointing out that conditions (2.2) assert that each mosquito population is viable and persistent in the absence of another. In other words, the recruitment or birth rate $\rho_N$ (resp. $\rho_W$) is greater than the natural density-independent mortality rate $\alpha_N$ (resp. $\alpha_W$). Therefore, conditions (2.2) are very realistic, whereas their violation would imply the inviability of each mosquito population.

It is easy to verify that condition (ⅰ) on p. 1152 is fulfilled by Eq (2.1b) according to which Wolbachia-carriers are progenies of $W$ only (cf. the positive term in the right-hand side of (2.1b) that expresses the recruitment of $W$). Additionally, Eq (2.1a) captures the condition (ⅱ) on p. 1152 that refers to the CI-phenotype in the sense that recruitment of wild mosquitoes (cf. the positive term in the right-hand side of (2.1b)) is proportional to the number of matings between wild males and females. These outcomes agree with the description of maternal transmission and the reproductive phenotype of cytoplasmic incompatibility presented in Table 1.

To illustrate the remaining features, a deeper analysis of the model (2.1) is required, and the latter is presented in the next section.

3.   Primary properties of the reduced model

This section is focused on establishing the key features (ⅲ) and (ⅳ) on p. 1152 for the reduced bidimensional model (2.1). The well-posedness of the model is attested in subsection 3.1, and the stability analysis of the system (2.1) is carried out in subsection 3.2.

3.1.   Well-posedness of the model

To verify the condition (ⅲ) on p. 1152, we establish and prove the following result.

Proposition 1. With condition (2.2) in force, dynamical system (2.1) is invariant in the positive cone $\mathbb{R}_+^2$. Additionally, there exists a compact set $\mathcal{X} \subset \mathbb{R}_+^2$ such that the trajectories of (2.1) engendered by initial conditions $\big(N(0), W(0) \big) \in \mathcal{X}$ remain in $\mathcal{X}$ for all $t \geq 0$.

Proof. It is immediately noted that

which plainly indicates that system (2.1) is invariant in the positive cone $\mathbb{R}^2_+$, i.e., its trajectories $N(t)$ and $W(t)$ engendered by $\big(N(0), W(0) \big) \in \mathbb{R}_+^2$ remain in $\mathbb{R}_+^2$ for all $t \geq 0$, and thus

Furthermore, let us observe that for $F(N, W)$ and $G(N, W)$ the following hold:

Thus, one has

where these quantities are strictly negative when $N > \frac{\rho_N - \alpha_N}{\beta_N}$ and $W > \frac{\rho_W - \alpha_W}{\beta_W}$. Hence, it follows that

where

In other words, the compact set

is invariant in the sense that all trajectories $(N(t), W(t))$ of (2.1) engendered by the initial conditions $\big(N(0), W(0) \big) \in \mathcal{X}$ remain in $\mathcal{X}$ for all $t \geq 0$.

It stems from the proof of Proposition 1 that the reduced dynamical system (2.1) is dissipative ^[32], and $\mathcal{X}$ is referred to as an absorbing set. This means that all trajectories $N(t), W(t)$ of (2.1) engendered by $\big(N(0), W(0) \big) \in \mathbb{R}_+^2 \setminus \mathcal{X}$ are attracted to $\mathcal{X}$ and there is a finite $\bar{t} > 0$ such that

Thus, we can conclude that all solutions of (2.1) engendered by nonnegative initial conditions $\big(N(0), W(0) \big) \in {\mathbb R}_+^2$ are uniformly ultimately bounded.

Remark 1. It is easy to check that $\mathbb{R}_+^2$ contains three subsets that are invariant with respect to solutions of the system (2.1):

1. Set $\Omega_0 : = \big\{ (0, 0) \big\}$ containing the origin is invariant since for the initial condition $\big(N(0), W(0)\big) \in \Omega_0$ we have that $\big(N(t), W(t)\big) \in \Omega_0$ for all $t \geq 0$.

2. Set $\Omega_N : = \big\{ (N, W) \in \mathbb{R}_+^2: \ N > 0, W = 0 \big\}$ that contains the $N$-axis is invariant. In effect, $W(0) = 0$ implies the absence of Wolbachia-carriers for all $t \geq 0$ and the dynamical system (2.1) turns into the logistic equation for the wild population:

whose solutions engendered by $N(0) > 0$ tend to the carrying capacity $N_{\sharp}$ as $t \to \infty$. Therefore, if $\big(N(0), W(0) \big) \in \Omega_N$, we have that $\big(N(t), W(t) \big) \in \Omega_N$ for all $t \geq 0$.

3. Set $\Omega_W : = \big\{ (N, W) \in \mathbb{R}_+^2: \ N = 0, W > 0 \big\}$ that contains the $W$-axis is also invariant because $N(0) = 0$ implies the absence of wild mosquitoes for all $t \geq 0$ and the dynamical system (2.1) turns into the logistic equation for the Wolbachia-carrying population

whose solutions engendered by $W(0) > 0$ tend to the carrying capacity $W_{\sharp}$ as $t \to \infty$. Therefore, if $\big(N(0), W(0) \big) \in \Omega_W$, we have that $\big(N(t), W(t) \big) \in \Omega_W$ for all $t \geq 0$.

3.2.   Stability analysis of the model

To verify the condition (ⅳ) on p. 1152, it is necessary to determine all possible equilibria of the system (2.1) that are solutions of the algebraic system

After some manipulations, we determine that the dynamical system (2.1) has four steady states (or equilibria), all of them are located in $\mathcal{X} \subset \mathbb{R}_+^2$, defined in (3.2), and their coordinates can be explicitly expressed in terms of the model's parameters (the detailed derivation of the model's equilibria is provided in Appendix A).

● Trivial steady state ${\bf{E}}_0 = (0, 0)$ that corresponds to the extinction of both mosquito populations.

● Boundary steady state ${\bf{E}}_N = (N_{\sharp}, 0)$ that corresponds to the survival and persistence of the wild mosquito population and eventual extinction of the Wolbachia-carrying population.

● Boundary steady state ${\bf{E}}_W = (0, W_{\sharp})$ that corresponds to the survival and persistence of the Wolbachia-carrying population and eventual extinction of the wild mosquito population.

● Strictly positive steady state ${\bf{E}}_c = (N_c, W_c)$ that corresponds to the coexistence of both mosquito populations and where

and $N_c + W_c = W_{\sharp}$.

Figure 1 displays the abovementioned equilibria on the phase portrait of the system (2.1). This figure also shows two nullclines or zero-growth isoclines of the system: the $N$-nullcline (i.e., the curve $F(N, W) = 0$) is plotted in blue color, and the $W$-nullcline (i.e., the curve $G(N, W) = 0$) is plotted in red color.

Recall that the violation of conditions (2.2) (i.e., a situation with $\rho_N < \alpha_N$ or $\rho_W < \alpha_W$) represents the case in which, for each population, the mortality rate is higher than the recruitment rate. The only outcome in such a case would be the extinction of both populations, and the trivial steady state ${\bf{E}}_0 = (0, 0)$ would be globally asymptotically stable because the other (nontrivial) equilibria ${\bf{E}}_N, {\bf{E}}_W$, and ${\bf{E}}_c$ would become unfeasible (that is, with negative coordinates, cf. (3.1) and (3.6)). In other words, the absorbing set $\mathcal{X}$ defined by (3.2) would contain only the trivial equilibrium, and all the system trajectories would be attracted to ${\bf{E}}_0$. We will not take into consideration this case, and from now on, we will assume that both conditions in (2.2) are always in force.

On the other hand, the existence of a strictly positive equilibrium ${\bf{E}}_c = (N_c, W_c)$ requires that

From the biological standpoint, condition (3.7) is rather credible. Let us recall that Wolbachia infection negatively affects the individual fitness of its carriers by reducing the females' fecundity and increasing the natural mortality of mosquitoes (see the exhaustive review by Dorigatti et al., 2018 ^[25] and more precise references therein). Therefore, it holds that

Furthermore, several recent studies have determined that, at high levels of intraspecific competition, Wolbachia-infected larvae experience reduced survival ^[45,46,55]. In other words, the Wolbachia-carrying insects exhibit higher mortality due to intraspecific competition than wild insects, meaning that $\beta_N \leq \beta_W$. From the above rationale, it follows that

which is equivalent to (3.7).

The following result is claimed in order to complete the proof of the condition (ⅳ) on p. 1152.

Proposition 2. With conditions (2.2) and (3.7) in force, the steady states ${\bf{E}}_N = (N_{\sharp}, 0)$ and ${\bf{E}}_W = (0, W_{\sharp})$ are asymptotically stable while the coexistence equilibrium ${\bf{E}}_c = (N_c, W_c)$ is unstable. Furthermore, the trivial steady state ${\bf{E}}_0$ is a source (nodal repeller).

Proof. Local stability of each equilibrium can be determined by the signs of eigenvalues of the Jacobian matrix associated with the system (2.1):

Direct evaluation of (3.8) in the boundary steady state ${\bf{E}}_N = (N_{\sharp}, 0)$ renders that

which is an upper-triangular matrix and its eigenvalues $\lambda_i^N, i = 1, 2$ are located on the main diagonal. In effect,

are both negative due to the conditions (2.2), (3.7). Therefore, ${\bf{E}}_N = (N_{\sharp}, 0)$ is locally asymptotically stable (nodal attractor).

Furthermore, by direct substitution of the boundary steady state ${\bf{E}}_W = (0, W_{\sharp})$ in (3.8), we obtain

which is a lower-triangular matrix and its eigenvalues $\lambda_i^W, i = 1, 2$ are located on the main diagonal. Consequently,

are both negative by virtue of (2.2). Therefore, ${\bf{E}}_W = (0, W_{\sharp})$ is locally asymptotically stable (nodal attractor).

On the other hand, direct evaluation of the Jacobian matrix (3.8) in the coexistence steady state ${\bf{E}}_c = (N_c, W_c)$ defined by (3.6) results in a rather cumbersome approach. However, let us recall that $\det \Big(\mathbb{J}(N_c, W_c) \Big) = \lambda_1^c \lambda_2^c$ where $\lambda_i^c, i = 1, 2$ denote two eigenvalues of $\mathbb{J}(N_c, W_c)$. Therefore, to prove that the coexistence equilibrium ${\bf{E}}_c = (N_c, W_c)$ is unstable, it is sufficient to show that $\det \Big(\mathbb{J} (N_c, W_c) \Big)$ is strictly negative. In effect,

Thus, we have $\det \Big({\mathbb J} (N_c, W_c) \Big) = \lambda_1^c \lambda_2^c < 0$ meaning that $\lambda_1^c$ and $\lambda_2^c$ have opposite signs. Therefore, the coexistence equilibrium ${\bf{E}}_c = (N_c, W_c)$ is unstable (saddle point).

To show that ${\bf{E}}_0 = (0, 0)$ is repelling, we cannot merely substitute its coordinates in the Jacobian matrix (3.8) since $\mathbb{J}(0, 0)$ cannot be directly computed at this steady state. Let us now recall (see Remark 1) that the origin is a self-contained invariant set and that the system trajectories engendered by $W(0) = 0$ and $N(0) > 0$ (which can be arbitrarily small) are attracted to ${\bf{E}}_N$, while the trajectories engendered by $N(0) = 0$ and $W(0) > 0$ (which can also be arbitrarily small) are attracted to ${\bf{E}}_W$. Thus, every vicinity of ${\bf{E}}_0$ contains initial conditions from which the system trajectories move away from the origin ${\bf{E}}_0 = (0, 0)$ in the direction of either ${\bf{E}}_N$ or ${\bf{E}}_W$. This means that ${\bf{E}}_0$ is a source. This completes the proof of Proposition 2 and also establishes the validity of the condition (ⅳ) presented on p. 1152.

In the next section, we explore monotonicity and other important properties of the reduced system (2.1).

4.   Monotonicity and other important properties of the model

To show that a certain dynamical system is monotone, it should be recalled that a dynamical system is called monotone when the flow generated by this system preserves some partial order ^[32]. Therefore, first, we define the partial order induced by a convex cone $\mathcal{K} \subset \mathbb{R}^2$ which is typically a quadrant of $\mathbb{R}^2$ for bidimensional dynamical systems ^[30,53]. Given two elements ${\bf{X}} = (x_1, x_2) \in \mathcal{K}$ and ${\bf{Y}} = (y_1, y_2) \in \mathcal{K}$, we write:

● ${\bf{X}} \leq_{\mathcal{K}} {\bf{Y}}$ if ${\bf{Y}} - {\bf{X}} \in \mathcal{K}$;

● ${\bf{X}} < _{\mathcal{K}} {\bf{Y}}$ if ${\bf{Y}} - {\bf{X}} \in \mathcal{K}\setminus {\bf{0}}$;

● ${\bf{X}} \ll_{\mathcal{K}} {\bf{Y}}$ if ${\bf{Y}} - {\bf{X}} \in {\rm{Int }}\; \mathcal{K}$.

Partial order defined by the last two items is also referred to as "strict order" and "strong order", respectively.

The sets

are referred to as order intervals (open and closed, respectively) induced by the cone $\mathcal{K}$.

Let ${\bf{X}} \in \mathbb{R}_+^2$ such that ${\bf{X}} = (x_1, x_2): = (N, W)$ defines the state vector of the system (2.1). Denote as ${\bf{F}}: \mathbb{R}_+^2 \mapsto \mathbb{R}_+^2, {\bf{F}}: = \big(F({\bf{X}}), G({\bf{X}}) \big)$ the vector field whose components $F$ and $G$ are defined in (2.1) as scalar functions of ${\bf{X}} = (N, W)$. Using these notations, our system (2.1) can be written as

and its solution, engendered by an initial condition ${\bf{X}}_0 = \big(N(0), W(0) \big) \in \mathbb{R}_+^2$, can be denoted as $\Phi(t, {\bf{X}}_0)$.

The positive semiflow of the system (2.1) or (4.1) generated by the vector field ${\bf{F}}({\bf{X}})$ is then the continuous mapping defined by $\Phi: \ {\mathbb R}_+ \times {\mathbb R}_+^2 \mapsto \times {\mathbb R}_+^2$, where $\Phi(t, {\bf{X}})$ denotes the solution of (4.1) that satisfies ${\bf{X}}(0) = {\bf{X}}$. Here, we consider only the positive semiflow of the system (4.1), that is, with $t \in [0, \infty)$, since we are interested in the system's behavior in forward time.

The semiflow $\Phi(t, \cdot)$ is called monotone (resp. strongly monotone) on a subset $A \subset {\mathbb R}^2_+$, with respect to partial order induced by the cone $\mathcal{K}$ if

Furthermore, the semiflow $\Phi(t, \cdot)$ is called strongly order-preserving on $A \subset {\mathbb R}^2_+$ (or SOP, for briefness) if $\Phi (t, \cdot)$ is monotone on $A$ and, whenever ${\bf{X}} < _{\mathcal{K}} {\bf{Y}}$, there exist open neighborhoods $\mathcal{V}_X$ of ${\bf{X}}$ and $\mathcal{V}_Y$ of ${\bf{Y}}$ such that

In other words, the partial order induced by the cone $\mathcal{K}$ is preserved for every ordered pair $\hat{{\bf{X}}} \leq_{\mathcal{K}} \hat{{\bf{Y}}}$, where $\hat{{\bf{X}}} \in \mathcal{V}_X, \hat{{\bf{Y}}} \in \mathcal{V}_Y$ within these open vicinities $\mathcal{V}_X$ and $\mathcal{V}_Y$ for all $t > 0.$

Proposition 3. System (2.1), (4.1) is monotone in $\mathbb{R}_+^2$ and strongly order-preserving (SOP) in the interior of $\mathbb{R}_+^2$ with respect to the partial order induced by the cone $\mathcal{K}: = \mathbb{R}_+ \times \mathbb{R}_-$.

Proof. It is easy to see that the partial order induced by $\mathcal{K}: = \mathbb{R}_+ \times \mathbb{R}_-$ is related to the "standard" order (induced by $\mathbb{R}_+^2$) in the following sense. For any two elements ${\bf{X}} = (x_1, x_2) \in \mathbb{R}_+^2$ and ${\bf{Y}} = (y_1, y_2) \in \mathbb{R}_+^2$, it is fulfilled that

A similar statement also holds with "$\ll_{\mathcal{K}}$" replacing "$\leq_{\mathcal{K}}$" and "$<$" replacing "$\leq$".

Using the idea of order isomorphism (for more details, see [32,Section 3.5]) induced by the diagonal matrix

it is immediate to conclude that

According to ^[32,40], a simple change of variables

leads to the cooperative dynamical system

which is monotone with respect to the "standard" order (induced by $\mathbb{R}_+^2$). Furthermore, the Jacobian of this system $\mathbb{J}(\bf{Z})$ is linked to the Jacobian $\mathbb{J} (N, W)$ of (2.1), (4.1) by the relationship

and is a Metzler matrix

for all $(N, W) \in \mathbb{R}_+^2$. The latter implies that the semiflow generated by (4.2) preserves the "standard" order (induced by $\mathbb{R}_+^2$), while the semiflow $\Phi(t, \cdot)$ generated by the system (2.1), (4.1) preserves the partial order induced by the cone $\mathcal{K}: = \mathbb{R}_+ \times \mathbb{R}_-$ for all $(N, W) \in \mathbb{R}_+^2$. Moreover, the Jacobian matrix (4.3) is irreducible whenever $N \neq 0, W \neq 0.$ Therefore, the system (2.1), (4.1) is strongly order-preserving in the interior of $\mathbb{R}_+^2$ in accordance with the results from [32,Chapter 4].

Let us recall the basin of attraction of a locally asymptotically stable equilibrium ${\bf{X}}^{*}$ for the dynamical system (4.1) is the set of initial conditions ${\bf{X}}_0$ such that the solutions $\Phi\big(t; {\bf{X}}_0 \big)$ engendered by ${\bf{X}}_0$ converge to ${\bf{X}}^{*}$ as $t \to \infty$. According to Proposition 2, our system (2.1) possesses two local attractors, ${\bf{E}}_N = \big(N_{\sharp}, 0 \big)$ and ${\bf{E}}_W = \big(0, W_{\sharp} \big)$, whose respective basins of attraction can be written as

where the limits are understood in the "componentwise" sense.

To determine sets included in the above basins of attraction, we introduce the following lemma obtained from Proposition 3.

Lemma 1. Systems (2.1) and (4.1) are strongly order-preserving (SOP) on $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ with respect to the partial order induced by the cone $\mathcal{K}: = \mathbb{R}_+ \times \mathbb{R}_-$.

Proof. We follow the proof scheme used in [22,Theorem 5], using the fact that strong monotonicity implies the SOP property, as indicated in [32,Proposition 1.1.1]. First, from Proposition 1, it can be easily deduced that systems (2.1) and (4.1) are positively invariant in the interior of $\mathbb{R}_+^2$. Furthermore, from Remark 1, the sets $\mathbb{R}_+\times \{0\}$, $\{0\} \times \mathbb{R}_+$ and $\{0\}\times\{0\}$ are invariant as well. Proposition 3 proves that the flow is SOP in the interior of $\mathbb{R}_+^2$ which is invariant. In each invariant set $\mathbb{R}_+\times \{0\}$ and $\{0\} \times \mathbb{R}_+$, systems (2.1) and (4.1) are reduced to a quadratic equation, equivalent to a logistic equation growth in each case. In each of these sets, in the order relation restricted to each invariant set, the flow is strongly monotone, and thus, the flow is also SOP in these sets. The condition holds trivially in $\{0\}\times\{0\}$ given that this point is a steady state.

Finally, having proven that the flow is SOP for initial conditions in each of the invariant sets, from [32,Remark 5.1.1], we can deduce that given an initial condition in one of these sets (i.e., in the border of ${\mathbb R}^2_+$) and another initial condition in the interior of $\mathbb{R}_+^2$, strong monotonicity is also preserved, which concludes the proof that the flow is SOP on $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ because this interval is a subset of ${\mathbb R}^2_+$.

The combination of results established by Propositions 2 and 3 and Lemma 1 together provide an essential basis for applying the fundamental Order Interval Trichotomy Theorem (see, e.g., [32,Theorem 2.2.2]) to the order interval $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ and two subintervals it contains. Let us recall the formulation of this theorem.

Theorem 1 (^[32], Theorem 2.2.2). Let the semiflow $\Phi(t, \cdot)$ of the dynamical system given in general form (4.1) be SOP with respect to the partial order induced by some cone $\mathcal{C}$ on the order interval $\left[\!\left[ \right.\!\right. {\bf{X}}^{*}, {\bf{Y}}^{*} \left.\!\left. \right]\!\right] _{\mathcal{C}}$ where ${\bf{X}}^{*} < _{\mathcal{C}} {\bf{Y}}^{*}$ and ${\bf{X}}^{*}, {\bf{Y}}^{*}$ are equilibria of (4.1). If $\overline{\Phi \big(t, \left[\!\left[ \right.\!\right. {\bf{X}}^{*}, {\bf{Y}}^{*} \left.\!\left. \right]\!\right] _{\mathcal{C}} \big)}$ is compact for each $t > 0$, then one of the following holds:

(ⅰ) There exists another equilibrium $\bf{Z}^{*} \in \left.\!\left. \right]\!\right] {\bf{X}}^{*}, {\bf{Y}}^{*} \left[\!\left[ \right.\!\right. _{\mathcal{C}}$ of the system (4.1).

(ⅱ) For any ${\bf{X}}_0 \in \left[\!\left[ \right.\!\right. {\bf{X}}^{*}, {\bf{Y}}^{*} \left.\!\left. \right]\!\right] _{\mathcal{C}} \setminus \{ {\bf{Y}}^{*} \}$, all solutions $\Phi \big(t; {\bf{X}}_0 \big)$ are attracted to ${\bf{X}}^{*}$, that is, $\lim \limits_{t \to \infty} \Phi \big(t; {\bf{X}}_0 \big) = {\bf{X}}^{*}.$

(ⅲ) For any ${\bf{X}}_0 \in \left[\!\left[ \right.\!\right. {\bf{X}}^{*}, {\bf{Y}}^{*} \left.\!\left. \right]\!\right] _{\mathcal{C}} \setminus \{ {\bf{X}}^{*} \}$, all solutions $\Phi \big(t; {\bf{X}}_0 \big)$ are attracted to ${\bf{Y}}^{*}$, that is, $\lim \limits_{t \to \infty} \Phi \big(t; {\bf{X}}_0 \big) = {\bf{Y}}^{*}.$

To adapt the hypotheses of this theorem to the system (2.1) with $\mathcal{C} = \mathcal{K}$, ${\bf{X}}^{*} = {\bf{E}}_N$, and ${\bf{Y}}^{*} = {\bf{E}}_W$, it must be shown that $\Phi \big(t, \left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}} \big)$ has compact closure in $\mathbb{R}_+^2$. In this context, it is instructive to note that any orbit

of (2.1) has a compact closure due to the existence of the absorbing set $\mathcal{X} \subset \mathbb{R}_+^2$ (given by (3.2)) which, in effect, coincides with the order interval $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}} : = [0, N_{\sharp}] \times [0, W_{\sharp}]$.

It is clear that direct application of Theorem 1 to the order interval $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ reaffirms, via item (ⅰ), the existence of ${\bf{E}}_c = (N_c, W_c)$ such that ${\bf{E}}_N < _{\mathcal{K}} {\bf{E}}_c < _{\mathcal{K}} {\bf{E}}_W$. On the other hand, the following proposition relates the two order subintervals of $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ with the two basins of attraction $\mathcal{B}_{{\bf{E}}_N}$ and $\mathcal{B}_{{\bf{E}}_W}$ of the boundary equilibria ${\bf{E}}_N$ and ${\bf{E}}_W$ defined by (4.4).

Proposition 4. For the system (2.1), it is fulfilled that

where

Proof. Dynamical system (2.1) fulfills the hypotheses of Theorem 1, which will be applied to the order intervals $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_c \left.\!\left. \right]\!\right] _{\mathcal{K}}$ and $\left[\!\left[ \right.\!\right. {\bf{E}}_c, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$, sets included in $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ where the flow is SOP (see Lemma 1). It is immediate to check that there is no equilibrium point inside the order interval $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_c \left.\!\left. \right]\!\right] _{\mathcal{K}}$, and ${\bf{E}}_N$ is an attractor. By virtue of item (ii) of Theorem 1, all orbits $\mathcal{O} (N, W)$ started in $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_c \left[\!\left[ \right.\!\right. _{\mathcal{K}}$ are attracted to ${\bf{E}}_N$. Therefore, $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_c \left[\!\left[ \right.\!\right. _{\mathcal{K}}$ belongs to $\mathcal{B}_{{\bf{E}}_N}$. A similar rationale applies to the order interval $\left.\!\left. \right]\!\right] {\bf{E}}_c, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ making use of item (iii) of Theorem 1.

Figure 1 displays order intervals $\left[\!\left[ \right.\!\right. {\bf{E}}_N, {\bf{E}}_c \left.\!\left. \right]\!\right] _{\mathcal{K}}$ and $\left[\!\left[ \right.\!\right. {\bf{E}}_c, {\bf{E}}_W \left.\!\left. \right]\!\right] _{\mathcal{K}}$ as pink and green rectangles, respectively. The pink rectangle belongs to $\mathcal{B}_{{\bf{E}}_N} \cap \mathcal{X}$ while the green one lies inside $\mathcal{B}_{{\bf{E}}_W} \cap \mathcal{X}$, that is, each order subinterval is included in the intersection of the underlying basin of attraction and the absorbing set $\mathcal{X}$. The unmarked areas of $\mathcal{X}$ contain points ${\bf{X}} = (N, W)$ that may belong to either $\mathcal{B}_{{\bf{E}}_N}$ or $\mathcal{B}_{{\bf{E}}_W}$.

The following result establishes further properties of the reduced system (2.1) related to its bistability.

Proposition 5. Dynamical system (2.1) exhibits the saddle-point behavior and there exists an invariant threshold manifold that passes through the positive steady state ${\bf{E}}_c$ separating the attraction basins $\mathcal{B}_{{\bf{E}}_N}$ and $\mathcal{B}_{{\bf{E}}_W}$.

Proof. Notably, according to ^[33], a dynamical system is said to admit a "saddle-point behavior" if it possesses two locally stable equilibria on the boundary of the state domain and one unstable (saddle-point) equilibrium in the interior of the state domain. Furthermore, the state domain of the system can be divided into three disjoint and invariant parts: two attraction basins (each containing one stable equilibrium on the boundary) and the so-called "threshold" manifold containing the unstable equilibrium that separates the attraction basins of two locally stable equilibria. Such a manifold is also referred to as the separatrix of two attraction basins.

The saddle-point behavior of the system (2.1) will be shown by applying the result summarized by H. Smith [53,Theorem 3.2]. For that purpose, we establish the cogency of four necessary hypotheses:

(H1) The semiflow $\Phi(t, \cdot)$ generated by the system (2.1) is strictly order-preserving on $\mathbb{R}^2_+$ with respect to $< _{\mathcal{K}}$ and order-compact for each $t > 0$.

(H2) The origin ${\bf{E}}_0 = (0, 0)$ is a repelling equilibrium.

(H3) All orbits originated on the boundaries of $\mathbb{R}_+^2$ are confined to these boundaries.

(H4) If ${\bf{X}}, {\bf{Y}} \in \mathbb{R}_+^2$ satisfy ${\bf{X}} < _{\mathcal{K}} {\bf{Y}}$ and either ${\bf{X}}$ or ${\bf{Y}}$ belongs to Int $\mathbb{R}_+^2$, then $\Phi(t, {\bf{X}}) \ll_{\mathcal{K}} \Phi(t, {\bf{Y}})$ for $t > 0$. If ${\bf{X}} = (X_1, X_2) \in \mathbb{R}_+^2$ satisfies $X_i \neq 0, i = 1, 2,$ then $\Phi(t, {\bf{X}}) \in$ Int $\mathbb{R}_+^2$ for $t > 0$.

To show the validity of (H1), we recall the statement of Proposition 3 according to which the semiflow $\Phi(t, \cdot)$ generated by the system (2.1) is SOP in the interior of $\mathbb{R}_+^2$ and, therefore, it is also strictly order-preserving in Int $\mathbb{R}_+^2$. Strict monotonicity of $\Phi(t, \cdot)$ on the borders of $\mathbb{R}_+^2$ (which are, in fact, invariant sets $\Omega_N$ and $\Omega_W$, see Remark 1) follows from strict monotonicity of solutions of logistic equations (3.3) and (3.4). Furthermore, the semiflow $\Phi(t, \cdot)$ generated by the system (2.1) is order-compact since any orbit (4.5) of (2.1) has a compact closure due to the existence of the absorbing set (3.2).

The hypothesis (H2) is cogent by virtue of Proposition 2, and the hypothesis (H3) is justified by invariance of the boundaries (sets $\Omega_N$ and $\Omega_W$, see Remark 1).

Finally, the hypothesis (H4) is corroborated by the positiveness of the system trajectories engendered by positive initial conditions (see Proposition 1) along with the SOP property of the semiflow $\Phi(t, \cdot)$ on Int $\mathbb{R}_+^2$ (see Proposition 3).

With the hypotheses (H1)-(H4) in force, we can now apply the result recaptured by H. Smith [53,Theorem 3.2] which basically affirms the following. If there is a unique equilibrium in ${\bf{E}}_c \in \left.\!\left. \right]\!\right] {\bf{E}}_N, {\bf{E}}_W \left[\!\left[ \right.\!\right. _{\mathcal{K}} \: \cap \; {\rm{Int }}\; \mathbb{R}_+^2$ and it is a saddle point, then there exists an unordered positively invariant set

that contains the unstable equilibria ${\bf{E}}_0, {\bf{E}}_c$ and consists of points ${\bf{X}}_s \in {\rm{Int }}\; \mathbb{R}_+^2$ such that

This completes the proof of Proposition 5 and establishes the saddle-point behavior of the system (2.1).

5.   Practical applications and final remarks

The thorough analysis of the reduced model (2.1) performed in Sections 3 and 4 provides very useful insights for practical applications. By recalling the Stable Manifold Theorem (see, e.g., [50,p. 107]), it can be concluded that the invariant "threshold" manifold $\mathcal{S}_{{\bf{E}}_c}$ is, in effect, the stable invariant manifold of the saddle point ${\bf{E}}_c$ that is tangent to the eigenvector generated by the negative eigenvalue of $\mathbb{J}({\bf{E}}_c)$. Furthermore, there also exists the unstable invariant manifold of the saddle point ${\bf{E}}_c$ that is tangent to the eigenvector generated by the positive eigenvalue of $\mathbb{J}({\bf{E}}_c)$. The unstable manifold contains two monotone heteroclinic orbits ^[32] that connect the saddle point ${\bf{E}}_c$ with two local attractors ${\bf{E}}_N$ and ${\bf{E}}_W$.

Figure 2 displays the plots of the stable and unstable manifolds of ${\bf{E}}_c$ (green and purple curves, respectively). It is clearly shown that the green curve plays the role of separatrix $\mathcal{S}_{{\bf{E}}_c}$ that divides $\mathbb{R}_+^2$ into two attraction basins $\mathcal{B}_{{\bf{E}}_N}$ (unshaded area) and $\mathcal{B}_{{\bf{E}}_W}$ (green-shaded region). The situation displayed in Figure 2 fully agrees with previous findings obtained for the case of perfect maternal transmission and CI ^{[19,24,31,43]}.

From the biological standpoint, points in $\mathcal{S}_{{\bf{E}}_c}$ indicate minimal viable population sizes of wild and Wolbachia-carrying mosquito populations that are related to the frequency-dependent Allee effect. In fact, if the initial wild mosquito population size is $N_0 \geq 0$, one can compute the minimal viable population size $\hat W_0 = \hat W_0(N_0)$, such that $(N_0, \hat W_0) \in \mathcal{S}_{{\bf{E}}_c}$. Thus, if the initial Wolbachia-carrying mosquito population $W_0$ is lower than the minimal viable population size $\hat W_0$, then $(N_0, W_0) \in \mathcal{B}_{{\bf{E}}_N}$ and the underlying solution $\big(N(t), W(t) \big)$ of the system (2.1) is attracted to ${\bf{E}}_N = (N_{\sharp}, 0)$, meaning the ultimate persistence of the wild population and progressive extinction of the Wolbachia-infected population.

Similarly, if the initial condition $(N_0, W_0)$ assigned to the system (2.1) lies above $\mathcal{S}_{{\bf{E}}_c}$, that is, $W_0 > \hat W_0$, it implies the initial size of the Wolbachia-infected population $W_0$ exceeds its minimal viable population size $\hat{W}_0$. In such a case, we have that $(N_0, W_0) \in \mathcal{B}_{{\bf{E}}_W}$ and the underlying solution $\big(N(t), W(t) \big)$ of the system (2.1) is attracted to the desired equilibrium ${\bf{E}}_W = (0, W_{\sharp})$, meaning the ultimate persistence of the Wolbachia-infected population and progressive extinction of the wild mosquito population.

In this context, it is worthwhile to highlight that while staying in line with the results presented in [8], our work goes further. In particular, the authors of ^[22] have identified only the "secured parts" of the attraction basins $\mathcal{B}_{{\bf{E}}_N}$ and $\mathcal{B}_{{\bf{E}}_W}$ in the form of four-dimensional order intervals^§. In contrast, our model has allowed us not only to identify the entire basins of attraction $\mathcal{B}_{{\bf{E}}_N}$ and $\mathcal{B}_{{\bf{E}}_W}$ corresponding to both boundary equilibria but also to estimate the minimum viable population sizes. The latter was not addressed in the work ^[22].

^§Note that their bidimensional analogs corresponding to our model are displayed in Figure 1 as pink and green rectangles, respectively.

To plot all figures presented in this paper, we have used the parameter values given in Table 2. These values correspond to the wMelPop strain of Wolbachia, which is regarded as the best one for controlling dengue infections among human individuals since it confers the most profound resistance to the replication of dengue virus in mosquitoes ^[25,28]. However, many scholars point out that the wMelPop strain is associated with high "fitness cost" since it reduces female fecundity, the viability of eggs, and the lifespan of infected mosquitoes ^[7,8,25]. Therefore, the successful invasion of mosquitoes carrying the wMelPop strain of Wolbachia and their durable persistence even in small detached localities appears to be a challenging task.

It should be recalled that the ultimate goal of Wolbachia-based biocontrol consists of seeking the eventual elimination of wild insects (capable of transmitting the virus to human individuals) by performing periodic releases of Wolbachia-carrying mosquitoes in some determined localities initially populated by wild mosquitoes. The practical implementation of this method requires mass-rear a massive quantity of Wolbachia-infected insects for posterior releases, and the desired result is propelled by the progressive Wolbachia invasion and its durable establishment in wild mosquito populations. The final outcome of this process is usually referred to as "population replacement".

Let us now provide some useful insights and practical interpretations derived from Figure 2 while keeping in mind the primary goal of Wolbachia-based biocontrol. First, we recall that the attraction basin $\mathcal{B}_{{\bf{E}}_W}$ contains the initial conditions $\big(N(0), W(0) \big)$ starting from which the trajectories of the system (2.1) converge to the desired boundary equilibrium ${\bf{E}}_W = (0, W_{\sharp})$ and the population replacement will be eventually achieved.

Suppose now that the population replacement is sought to be achieved with a single (or inundative) initial release of Wolbachia-carrying insects. To determine the size of such an abundant release, the current size of the wild mosquito population should be first assessed by some known technique ^[44,49]. Once the abundance of wild mosquitoes $N_0$ is fairly estimated, the information regarding the minimal viable sizes $\hat W_0 = \hat W_0(N_0)$ such that $\big(N_0, \hat{W}_0 \big) \in \mathcal{S}_{{\bf{E}}_c}$ of Wolbachia-carrying mosquito populations will be of the utmost importance, and this information is explicitly supplied by the model (2.1) and its underlying parameters.

Upon closer inspection of Figure 2, we observe that the area inside $\mathcal{B}_{{\bf{E}}_W}$ surrounding $W_{\sharp}$ (portion of the green-shaded region around $W_{\sharp}$) is much smaller than the area inside $\mathcal{B}_{{\bf{E}}_N}$ surrounding $N_{\sharp}$ (portion of the unshaded region around $N_{\sharp}$). Therefore, when the size of the wild mosquito population is close to saturation or its carrying capacity $N_{\sharp}$, an extreme amount of Wolbachia-carriers will be necessary for a single inundative release.

On the other hand, it is instructive to recall that wild mosquito populations may exhibit seasonal size variations ^[47,49]. In this context, the best timing for Wolbachia-based biocontrol by a single inundative release will be the period of relatively low mosquito abundance. Such periods usually correlate with cooler and windier seasons in tropical and subtropical regions ^[49] or arise after carefully planned and thoroughly implemented vector control measures ^[41,52].

For different initial sizes $N_0 = \lambda N_{\sharp}$ of wild mosquito population expressed as fractions $\lambda \in \{0.25, 0.5, 0.75, 1\}$ of $N_{\sharp}$, one can estimate the corresponding minimal release sizes of Wolbachia-carrying insects $\hat{W}_0 = \hat{W}_0(\lambda N_{\sharp}) = \hat{\lambda} N_{\sharp}$ (also expressed as the multiplicatives of $N_{\sharp}$) that ensure the population replacement by a single inundative release. The corresponding values of $\lambda$ and $\hat{\lambda}$ are presented in Table 3 (considering the parameters of Table 2), and the points $\Big(\lambda N_{\sharp}, \hat{\lambda} N_{\sharp} \Big) \in \mathcal{S}_{{\bf{E}}_c}$ are displayed in Figure 3 in red color.

From Table 3 and Figure 3, we observe that the minimal Wolbachia-infected population sizes needed for single inundative releases always exceed the initial sizes of the wild mosquito population. When the initial size of the wild population is higher than $0.5 N_{\sharp}$, the population replacement becomes substantially more difficult to reach with a single inundative release, as it requires mass-rearing vast numbers of Wolbachia-carriers.

As an alternative to a single inundative release, one may perform several periodic (or inoculative) releases. This strategy may seem reasonable if the mass-rearing facility cannot produce the vast quantity $\hat{\lambda} N_{\sharp}$ of Wolbachia-carriers en masse but is capable of producing a smaller quantity of Wolbachia-carrying insects every $\tau$ days. However, under such a setting, the total amount of Wolbachia-infected mosquitoes necessary to reach the population replacement will be larger than in the case of a single inundative release.

From the mathematical standpoint, periodic inoculative releases of Wolbachia-carrying insects can be modeled by the following impulsive dynamical system:

where $\tau$ denotes the period of releases, $\Lambda = \rm{const}$ stands the release size (i.e., the number of Wolbachia-carriers to be released), and $n$ defines the number of releases. Notably, $W(i\tau^{\pm})$ denote the right and left limits of the function $W(t)$ at $t = i\tau$. Formal analysis of the impulsive system (5.1) is a challenging task and may be proposed as an object for further studies. Therefore, in the context of this paper, we limit ourselves to revising its numerical solutions in order to assess the practical value of periodic releases and to compare their overall performance with an outcome of a single inundative release.

By performing a series of numerical simulations, we have estimated the minimal release sizes $\Lambda: = \hat{\lambda} N_{\sharp}$ (also expressed as the multiplicatives of $N_{\sharp}$) for $\tau = 1$ day and $\tau = 3$ days and taking different initial sizes of wild mosquito population $N_0$, expressed as fractions $\lambda \in \{0.25, 0.5, 0.75, 1\}$ of $N_{\sharp}$.

In Table 4, we present minimum release sizes $\Lambda$ that ensure the population replacement by $n$ periodic releases even when the initial condition $\big(N_0, \Lambda \big)$ lies inside the attraction basin $\mathcal{B}_{{\bf{E}}_N}$ of the boundary equilibrium ${\bf{E}}_N = (N_{\sharp}, 0)$ (that is, strictly below the separatrix $\mathcal{S}_{{\bf{E}}_c}$). Revising the entries of Table 4, it is easy to detect several patterns or "tradeoffs" between the frequency of releases $\tau$, constant release size $\Lambda$, and the overall number of releases $n$ needed to ensure the population replacement. Namely, more frequent releases ($\tau = 1$) require smaller release sizes $\Lambda$ and shorter overall time of the release program but a greater number of releases $n$. The latter is quite logical, and not only aligns with common sense, but also bears similarities with other works dealing with periodic releases of mosquitoes ^[42]. In this context, the anticipated knowledge of the production costs related to the mass-rearing of Wolbachia-carrying insects and the logistics costs for performing field releases are important for choosing the release frequency. Although we have no reliable information regarding such costs, the impulsive system (5.1) may serve to be of potential utility in the future when healthcare entities eventually decide to evaluate this method of biological vector control. In particular, the model (5.1) may help adjust the frequency and the sizes of releases to the production capacity of artificial mass-rearing facilities of the Wolbachia-carrying insects.

On the other hand, the outcomes of numerical simulations performed on the original model (2.1) can be also compared with those obtained for the impulsive system (5.1). Contrasting the values of $\hat{\lambda}$ from Tables 3 and 4 for the same values of $N_0 = \lambda N_{\sharp}$, we observe that they bear a more striking difference for $\tau = 1$ than for $\tau = 3$. Moreover, the mentioned difference is smaller for the smaller values of $N_0$ (such as $0.25 N_{\sharp}$ and $0.5 N_{\sharp}$) and becomes more noticeable for the larger values of $N_0$ (such as $0.75 N_{\sharp}$ and $N_{\sharp}$). Thus, the release programs based on periodic inoculative releases seem more practicable when the initial size of the wild mosquito population is close to its saturation level $N_{\sharp}$.

Figure 4 displays simulation results for the impulsive system (5.1) when the wild mosquito population is at saturation $N_0 = N_{\sharp}$. The left column of Figure 4 corresponds to daily releases ($\tau = 1$ day) and the right one corresponds to inoculative releases performed every three days ($\tau = 3$). The upper charts of Figure 4 present the system's trajectories $N(t)$ and $W(t)$ drawn by blue and red-colored curves, respectively, and also bear two dashed lines marking the coordinates of the coexistence equilibrium ${\bf{E}}_c = \big(N_c, W_c)$. The lower charts exhibit the underlying parts of orbits $\mathcal{O}(N, W) = \Big\{ \big(N(t), W(t) \big) \in \mathbb{R}_+^{2}: \ t \geq 0 \Big\}$ in the phase space that start in $(N_0, W_0) \in \mathcal{B}_{{\bf{E}}_N}$ (red-colored point below the separatrix $\mathcal{S}_{{\bf{E}}_c}$) and move the system states to the attraction basin $\mathcal{B}_{{\bf{E}}_W}$ of the desired boundary equilibrium ${\bf{E}}_W = (0, W_{\sharp})$ (red-colored point above the separatrix $\mathcal{S}_{{\bf{E}}_c}$).

The periodic inoculative releases are suspended when the orbit $\mathcal{O}(N, W)$ of (5.1) crosses the separatrix $\mathcal{S}_{{\bf{E}}_c}$ and enters the attraction basin $\mathcal{B}_{{\bf{E}}_W}$ (green-shaded region in the lower charts of Figure 4). The latter is also clearly visible in the upper charts of Figure 4: the system trajectory $N(t)$ decays and crosses the blue-colored dashed line, while the trajectory $W(t)$ remains strictly above the red-colored dashed line.

Browsing once again the simulation results given in Tables 3 and 4 and contrasting them for each particular value of $N_0 = \lambda N_{\sharp}, \lambda \in \{ 0.25, 0.5, 0.75, 1 \}$, we may conclude that, from the practical standpoint, an implementation of a single inundative release seems more rational and operative than several periodic inoculative releases. Effectively, under the "worst scenario", i.e., when $N_0 = N_{\sharp}$ (this situation is illustrated in Figure 4) it is necessary to mass-rear at least $5.16 N_{\sharp}$ of Wolbachia-carriers during 12 days (with $\tau = 1$) or at least $11.12 N_{\sharp}$ of Wolbachia-carriers during 18 days (when $\tau = 3$), while a single inundative release only requires to mass-rear $1.85 N_{\sharp}$ of Wolbachia-infected insects, albeit all at once.

Thus, the reduced bidimensional model (2.1) has resulted in a quite handy and easily interpretable tool for determining the appropriate size of a single inundative release or periodic releases of Wolbachia-carrying insects since it explicitly yields the dependence between minimal viable population sizes of wild and Wolbachia-infected mosquito populations.

To conclude this study, let us mention that the proposed bidimensional model (2.1) is adaptable to other Wolbachia strains (e.g., wMEl, wAlb) that are being tested for vector control. Furthermore, this model can be adjusted to other insect species in which Wolbachia induces the virus inhibition (e.g., rice planthoppers and fruit flies ^[48,51]). Therefore, an anticipated knowledge of the separatrix $\mathcal{S}_{{\bf{E}}_c}$ may constitute the key element of Wolbachia-based biocontrol for human and crop diseases.

Acknowledgements

Diego Vicencio was supported by the program CONICYT PFCHA/ Doctorado Becas Chile/2017-21171813 and FONDECYT grant N 1200355 ANID-Chile program. Pedro Gajardo was partially supported by FONDECYT grant N 1200355 ANID-Chile program, and ANID/BASAL FB210005. Olga Vasilieva acknowledges financial support from the National Fund for Science, Technology, and Innovation (Autonomous Heritage Fund Francisco José de Caldas) by way of the Research Program No. 1106-852-69523 (Principal Investigator: Hector J. Martinez), Contract: CT FP 80740-439-2020 (Colombian Ministry of Science, Technology, and Innovation – Minciencias), Grant ID: CI-71241 (Universidad del Valle, Colombia). Olga Vasilieva also appreciates the endorsement obtained from the STIC AmSud Program for regional cooperation (20-STIC-05 NEMBICA project, international coordinator: Pierre-Alexandre Bliman, INRIA – France).

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix A.   Determination of the four equilibria of the system (2.1)

Let us recall that to find the steady states for our model, we must solve the following algebraic system:

Note that equation (A.1b) has a solution in $W = 0$. By replacing that value of $W$ in (A.1a), we can easily obtain the first possible steady state, ${\bf{E}}_0 = (0, 0)$. If we assume that $N\not = 0$, we obtain the following equation:

Taking this expression, it is evident that the value of $N$ must be $N_\sharp = \dfrac{\rho_N-\alpha_N}{\beta_N}$, thus obtaining the second steady state ${\bf{E}}_N = (N_\sharp, 0)$.

Let us now assume that $N = 0$, but $W\not = 0$. By replacing these values in (A.1b), we obtain the following equation.

Taking this expression, it is evident that the value of $W$ must be $W_\sharp = \dfrac{\rho_W-\alpha_W}{\beta_W}$, thus obtaining the third steady state ${\bf{E}}_W = (0, W_\sharp)$.

For the final steady state, we must assume both $N\not = 0$ and $W\not = 0$. By equation (A.1b), we obtain:

From this equation, we obtain that $W_c+N_c = \dfrac{\rho_W-\alpha_W}{\beta_W} = W_\sharp$. Replacing this value in (A.1b), we get:

Hence,

Finally, from this last expression, and using that $N_c+W_c = W_\sharp$ we obtain:

Thus, the last steady state is ${\bf{E}}_c = (N_c, W_c),$ and its coordinates match the formula (3.6).