Global dynamics of a simple model for wild and sterile mosquitoes

1.   Introduction

Vector-borne diseases including dengue fever and malaria that mosquitoes transmit have been threats to human beings. Although the use of chemical insecticides is effective, from the standpoint of environmental friendliness, the sterile insect technique (SIT) is known. This technique has emerged as an innovative pest control technology to manage and reduce the wild populations. SIT was originally invented by American entomologist Edward F. Knipling and his collaborators to control populations by releasing sterile insects into the wild ^[1].

The genetically engineered sterile mosquitoes do not produce offspring, thus reducing and controlling the number of the next generation's population. However, because the SIT method involves genetic engineering, it is necessary to evaluate the effects of releasing sterile mosquitoes and to estimate the optimal number of sterile mosquitoes to release. Therefore, it is worthwhile to build a simple mathematical model to describe the population dynamics associated with the generic practical implementation of SIT methods, and to obtain qualitative information about the optimal sterile mosquito release strategy and the consequences of population dynamics through simulation and mathematical analysis ^{[2,3,4,5,6,7,8,9,10,11,12,13]}. Each of the previous studies is briefly presented below for the convenience of the reader. See also Chapter 2.5 in ^[1] for a comprehensive overview of deterministic mathematical models of SIT. Since it is crucial to assess the impact of the actual release of sterile mosquitoes, it is worthwhile to build a simple mathematical model to describe the population dynamics associated with the implementation of SIT methods.

In ^[6], the authors formulate continuous-time mathematical models describing the interaction of the wild and sterile mosquitoes with different strategies in releasing sterile mosquitoes. The authors study the following system of differential equations:

with $N = x+y$. $C(N)$ is the number of mates per individual per unit time and $B(\cdot)$ is the release rate of sterile mosquitoes. The authors investigate the existence and stability of equilibria for the following three cases. The first is the case where $C(N)c = a$ and $B = b$ and the constant release rate is considered. The second case is when $cC(N)x/(x+y)$ is replaced by $ax/(1+x+y)$ and $B = bx$, taking into account the Allee effect and the proportional release rate. Finally, the third case is the case where the Allee effect is continued and the Holling-Ⅱ type release rate, $B = bx/(1+x)$, is considered. For mathematical analysis, $\mu_{1} = \mu_{2}$ and $\xi_{1} = \xi_{2}$ are assumed. It is proven that, in each case, the model can have two positive equilibria and one positive equilibrium is stable while another equilibrium is a saddle point. Periodic solutions appear for the case of proportional release, but not for the case of constant release. Numerical simulations are also given. Population dynamics of the wild and sterile mosquitoes is revisited in ^[9].

In ^[5], a sex-structured mosquito population model is formulated. It is assumed that all females are equally able to mate and a different mortality rate for male and female mosquito is considered. Sterile male mosquitoes, which interfere with the reproduction of the wild mosquitoes, are incorporated into the sex-structured mosquito population model. The authors study the existence and stability of equilibria, assuming that the number of sterile male mosquitoes is at an equilibrium, thus the model is reduced to a two-dimensional system. The authors compare constant continuous release to periodic impulsive releases, and obtain conditions for successful elimination of the mosquito population. The engineered sterile mosquitoes are not necessarily completely sterilized. In ^[4], residual fertility is incorporated in the model developed in ^[5]. See also ^[7,8] for models including residual fertility. Stage-structured models are also formulated in ^[2,3,7,8,13], to take into account mosquitoes' metamorphic stage structure, see also ^[10].

Mathematical analysis generally becomes more difficult as models become more complex, and this is also true for mathematical models in the SIT method, a problem that must be overcome in the pursuit of reproducibility of biological situations. In this paper, we aim to establish a mathematical analysis that works with a simple model for the consideration of biological situations. We hope that the approach in this paper will be the foundation for further reproducing of the phenomenon and contributing to solving problems in the implementation of the SIT method.

This paper considers the following system of two-dimensional ordinary differential equations:

System (1.2), studied in ^[12], is a simplification of a model studied in ^[6], in the particular case where $C(N)c = a$, $B(\cdot) = bx$, $\mu_{1} = \mu_{2}$ and $\xi_{1} = \xi_{2}$, ignoring interspecific competition. Our aim in this paper is to complement the study in ^[12] elucidating the dynamics. According to ^[6,12], $x = x(t)$ is the density of wild mosquitoes at any time $t$ and $y = y(t)$ is the density of genetically engineered sterile male mosquitoes at any time $t$. The parameter $\mu$ is the natural death rate of both populations, $\xi$ is the intra-species competitive death rate of both populations, and $a$ is the total number of offspring produced by a female wild mosquito per unit time. In addition, the term $x/(x+y)$ is the probability of mating between two wild mosquitoes, that decreases as $y$ increases. The parameter $b$ is the release rate of genetically engineered sterile mosquitoes. The initial conditions are $(x(0), y(0)) = (x_{0}, y_{0})$, where $x_{0}\ge0$ and $y_{0}\ge 0$. For (1.2), we impose

Then the right hand side is continuous and local Lipschitz in $\mathbb{R}_{+}^{2}$. This assumption (1.3) is also imposed in ^[12].

The boundedness of the solution in the case of (1.2), the existence of positive equilibria and their local stability are then investigated in ^[12] with numerical simulation results, based on the case where the proportional release rate in ^[6] is considered. For (1.2), the boundedness of the solution, the existence and stability of positive equilibria are investigated in ^[12] with numerical simulation results. Sasmal et. al. ^[12] provided a mathematical insight into the sterile mosquito release strategy in the simplest situation possible, without considering the constant rate and Holling-Ⅱ type release rate of sterile mosquitoes and the interspecific competition between wild and sterile mosquitoes. The existence of equilibria in the $\mu b$-parameter plane is also discussed. See Section 2 and ^[12] for details.

For (1.2), the origin plays the role of a singularity due to the ratio dependence, and even if (1.3) is imposed, only numerical simulation results are given for the dynamical system near the origin, and no mathematical analysis of the system has been performed. As shown in the numerical simulation results, when the initial condition is chosen from $(x(0), y(0)) = (x_{0}, y_{0})\in \mathbb{R}_{+}^{2}$, convergence to the origin is confirmed by the values of $\mu$ and $b$. Considering the original problem, the essential question is how the population dynamics of mosquitoes in the wild change with the chosen parameters, and it is important to know when convergence to the origin occurs. However, in (1.2), studying the behavior of solutions near the origin is a challenging problem, due to the difficulty of linear stability analysis, as the right-hand side is not differentiable at the origin.

Analysis on the neighborhood of the origin of differential equations where the derivative is not defined at the origin, such as the system (1.2), also appears in the predator-prey system with ratio-dependence. For instance, the following system of equations is considered in ^[14]:

The time-scale transformation $dt = (my+x)d\tau$ is introduced and the dynamical system around the equilibrium $(x, y) = (0, 0)$ is investigated using the blow-up technique with the polar coordinates and the Briot-Bouquet transformation.

In this paper, we apply the idea of analyzing dynamical systems near the singular point by means of a time-scale transformation and the blow-up technique in ^[14] to (1.2).

The blow-up technique is a method for studying dynamical systems near non-hyperbolic equilibria with double-zero eigenvalues. The details of the blow-up technique are described in Section 4. Also, see, for instance, ^[15,16,17] and references therein. In this paper, instead of using the polar coordinate, we use a directional blow-up in which the information on the singularity is divided into two local coordinates. The resulting dynamical system near the equilibria in the blow-up vector field is investigated using the center manifold theory, and the dynamical system near the origin of (1.2) can be fully understood.

In this paper, we analyze the dynamics near the origin and combine the informations, the boundedness of the solution, the existence and local stability of the positive equilibria, and the non-existence of periodic solutions, to obtain the global dynamics. This allows us to obtain a classification of the behavior in $(x(t), y(t))\in \mathbb{R}_{+}^{2}$ depending on the choice of the sterile mosquito release rate $b$ and the natural death rate $\mu$ for both populations. The qualitative dynamics is illustrated in a parameter plane. By providing the global behavior of the solution, we can provide qualitative information on the optimal sterile mosquito release strategy and the consequences of population dynamics, as well as considerations of the next generation population dynamics associated with the sterile mosquito release rate $b$ and validation of the mathematical model. This also improves the discussion of the existence and stability of the internal and positive equilibria with sharper results than ^[12]. The results are obtained not by solving the algebraic equations that appear when finding the positive equilibria, but by considering it as an extreme value problem, which gives a clearer view of the discussion of stability as well as the existence classification of the positive equilibria. See Appendix A for details.

The model used in this paper to evaluate the effectiveness of the SIT method and to predict population dynamics makes a number of assumptions regarding model simplification. Although the natural death and intraspecific competition rates of wild and sterile mosquitoes are assumed to be the same, the authors hope that the mathematical analysis methods in this paper will be effective even with more complex models and will provide insight into the release strategies of sterile mosquitoes in the SIT method.

This paper is organized as follows. In the next section, we summarize the results obtained in ^[12]. Those results are used in this paper. Section 3 describes the main results. In Section 4, we study the behavior of the system (1.2) near the origin, applying the blow-up technique. The proof of the main result is given in Subsection 4.3. Finally, we offer a discussion of the population dynamics of wild and genetically engineered mosquitoes and the effect of the SIT method based on the results obtained in this paper in Section 5.

2.   Preliminaries

In this section, we discuss the known results in (1.2) in ^[12] and some preparations for the mathematical analysis in this paper.

2.1.   Boundedness of the solutions and nonexistence of periodic solutions

In (1.2), the time-scale transformation $d\tilde{t}/dt = a$ allows us to set $a = 1$ by setting $\tilde{t}$ to $t$ again. Also, with $\tilde{x} = \xi x$ and $\tilde{y} = \xi y$, we can set $\xi = 1$ in (1.2). Since we obtain the differential equations

for $\tilde{x}$ and $\tilde{y}$, we can rewrite $\tilde{x}$ and $\tilde{y}$, as $x$ and $y$ respectively. This simplification is not made in ^[12] and simplifies the mathematical analysis in this paper. Therefore, the following systems are considered:

Also, for $\mu \ge 1$,

holds, as shown in Lemma 3.1 of ^[12], so

can be obtained. Note that $0 < \mu < 1$ is the case we are interested in. In ^[12], the authors obtain the following results on the boundedness of solutions for (2.2) and (1.3) in $\mathbb{R}_{+}^{2}$. Since they also play an important role in this paper, we mention them for the convenience of the reader.

Proposition 2.1 (^[12], Theorem 3.1). System (2.2) with (1.3) is positively invariant in $\mathbb{R}_{+}^{2}$. Let $\mu < 1$ hold. Then the solution of the system (2.2) with (1.3) is ultimately bounded in $\mathbb{R}_{+}^{2}$, satisfying the following properties:

(i) $\limsup\limits_{t\to \infty} x(t) \le 1-\mu$

(ii) $\limsup\limits_{t\to \infty} y(t) \le b(1-\mu)\mu^{-1}$

(iii) $\limsup\limits_{t\to \infty} [x(t)+y(t)] \le (1-\mu)(\mu+b)\mu^{-1}$

It has also been shown in ^[12] that there are no periodic solutions for the system (2.2). This result is important to discuss the global behavior of the solutions in this paper and is presented for the convenience of the reader.

Proposition 2.2 (^[12], Theorem 3.3). System (2.2) does not have any periodic/closed orbit in $\mathbb{R}_{+}^{2}$.

2.2.   Existence of equilibria and stability of positive equilibria

The internal and positive equilibria of (2.2) are also discussed in ^[12] and are important information for understanding the global behavior of the solutions in this paper. First, $E_{0}:(x, y) = (0, 0)$ always exists from (1.3) and is called an extinction equilibrium. Next, we review the existence of positive equilibria and their local stability in $\mu \ge 1$, $2^{-1} \le \mu < 1$, and $\mu < 2^{-1}$ based on ^[12].

Proposition 2.3 (^[12], Lemma 3.1, 3.2). When $\mu \ge 1$, system (2.2) has no positive equilibrium and the extinction equilibrium is globally asymptotically stable. When $2^{-1}\le \mu < 1$, system (2.2) has at most one positive equilibrium, moreover,

(i) if $b\ge 1-\mu$, then system (2.2) has no positive equilibrium,

(ii) if $b < 1-\mu$, then system (2.2) has unique positive equilibrium, $E_{1}: (x, y) = (x_{1}^{*}, y_{1}^{*})$.

In ^[12], conditions for the existence of internal and positive equilibria in (2.2) are given. In this paper, we obtain explicit conditions for the existence and local stability of the positive equilibria. The positive equilibria of (2.2) are found as the solutions $(x, y)$ of the following equations:

with $0 < \mu < 1$, $0 < x < 1-\mu$, and $0 < y$. Substituting the second equation into the first one leads to the problem of considering the roots of the following equation:

In this paper, the variable transformation:

is used to reduce to the problem of examining the roots of the following equation derived from (2.4) in the interval $(0, \mu^{-1}-1)$:

with

As described in Appendix A, $g^\prime (z)$ has only one zero point at $z > 0$, which is denoted by $z^{*}(\mu)$. The $g(z)$ takes its maximum value at $z = z^{*}(\mu)$. In Appendix A, the following results for positive equilibria for $0 < \mu < 2^{-1}$ are shown.

Proposition 2.4. Let us assume that $0 < \mu < 2^{-1}$ holds.

(iii) Suppose that $0 < \mu < 2^{-1}(3-\sqrt{5})$ holds.

(iii-1) If $b\leq1-\mu$, then system (2.2) has a unique positive equilibrium in $\mathit{\text{int}}\, \mathbb{R}_{+}^{2}$.

(iii-2) If $1-\mu < b < g(z^{*}(\mu))$, then system (2.2) has exactly two positive equilibria in $\mathit{\text{int}}\, \mathbb{R}_{+}^{2}$.

(iii-3) If $b = g(z^{*}(\mu))$, then system (2.2) has a unique positive equilibrium in $\mathit{\text{int}}\, \mathbb{R}_{+}^{2}$.

(iii-4) If $b > g(z^{*}(\mu))$, then system (2.2) has no positive equilibria in $\mathit{\text{int}}\, \mathbb{R}_{+}^{2}$.

(iv) Suppose that $\mu\geq 2^{-1}(3-\sqrt{5})$ holds.

(iv-1) If $b < 1-\mu$, then system (2.2) has a unique positive equilibrium in $\mathit{\text{int}}\, \mathbb{R}_{+}^{2}$.

(iv-2) If $1-\mu\leq b$, then system (2.2) has no positive equilibria in $\mathit{\text{int}}\, \mathbb{R}_{+}^{2}$.

In Theorem 3.2 in ^[12], the Jacobian matrices for the asymptotic stability of positive equilibria of (2.2) are obtained. The conditions given in Theorem 3.2, a local stability result, is implicit in terms of the parameters. In this paper, therefore, we obtain a simplified representation of the Jacobian matrices from the conditions satisfied by the equilibria, and obtain the following results on the asymptotic stability of the positive equilibria.

Proposition 2.5. Let us assume that $0 < \mu < 1$ holds.

(v) Suppose that $0 < \mu < 2^{-1}(3-\sqrt{5})$ holds.

(v-1) If $b\leq 1-\mu$, then system (2.2) has a unique positive equilibrium and it is asymptotically stable.

(v-2) If $1-\mu < b < g(z^{*}(\mu))$, then system (2.2) has exactly two positive equilibria. We denote the two positive equilibria by $(\overline{x}_{1}, \overline{y}_{1})$ and $(\overline{x}_{2}, \overline{y}_{2})$ satisfying

The positive equilibria $(\overline{x}_{1}, \overline{y}_{1})$ and $(\overline{x}_{2}, \overline{y}_{2})$ are unstable and asymptotically stable, respectively.

(vi) Suppose that $\mu\geq 2^{-1}(3-\sqrt{5})$ holds. If $b < 1-\mu$, then system (2.2) has a unique positive equilibrium and it is asymptotically stable.

Appendix A outlines the proofs of Propositions 2.4 and 2.5. In ^[12], the existence or nonexistence of positive equilibria and their local stability are organized in the $\mu b$-parameter plane. Since the classification of the solution behavior in this parameter plane plays an important role in the discussion of the global behavior of solutions and the global stability of equilibria including the origin, which is the main objective of this paper, we reproduce Figure 1 of ^[12] as shown in Figure 1 below.

3.   Main results

In this section, we present the main results on the global behavior of the solutions of (2.2). These main results are obtained by combining the information on the existence and stability of the positive equilibria, classified in the $\mu b$-parameter plane given in ^[12] (see also Figure 1), with the dynamics near the origin of (1.2), which we clarify in this paper. Finally, the global behavior of the solution is classified in the $\mu b$-plane, see Figure 2.

The curve $b = f(\mu)$ is given as the unique positive root of the following equation:

for $0 < \mu < \mu^*$, where $\mu^{*} = 2^{-1}(3-\sqrt{5})$. One can see that $f(\mu)$ coincides with $g(z^*(\mu))$.

In the $\mu b$-plane, we define the region (A) as

The region (B) is the region denoted by $\{(b, \mu) \mid b = f(\mu), 0 < \mu < \mu^* \}$ and the region (C) is the region denoted by $\{(b, \mu) \mid b \leq f(\mu), \ 0 < \mu < \mu^* \}$. The region (D) is $\{(b, \mu) \mid 1-b-\mu = 0, b^{*} < b < 1\}$ with $b^{*} = 2^{-1}(\sqrt{5}-1)$ and the region (E) is $\{(b, \mu) \mid 0 < \mu < 1, 0 < b < 1\}$. See also Figure 2.

Theorem 3.1. Assume that $b > 0$ and $\mu > 0$. The extinction equilibrium $E_{0}: (x, y) = (0, 0)$ of system (1.2) is global asymptotically stable if $b$ and $\mu$ are in region (A) in the $\mu b$-plane.

Theorem 3.2. Assume that $b > 0$ and $\mu > 0$. The equilibrium that is not the extinction equilibrium $E_{0}: (x, y) = (0, 0)$ of system (2.2) is global asymptotically stable if $b$ and $\mu$ are in region (D). The unique positive equilibrium $E_{1}^{*} = (x_{1}^{*}, y_{1}^{*})$ of system (2.2) is global asymptotically stable if $b$ and $\mu$ are in region (E).

When $b$ and $\mu$ belong to region (B) in Figure 2, we know that solutions converge to the origin or a positive equilibria depending on the initial values. The same is true for region (C). In this paper, we do not obtain the basin of the attraction of the origin and positive equilibrium. This issue would be an important indicator for evaluating the effectiveness of the SIT method. Providing an answer to this question is a subject for future work.

4.   Behavior near the origin and proof of theorems

First, for (2.2), we introduce the following time-scale desingularization:

By (4.1), we are led to the problem of considering the following differential equation:

The time-scale desingularization (4.1) simply multiplies the vector field of (2.2) by $(x+y)$. Except for the singularity $\{x = y = 0\}$, the solution curve of the system (vector field) remains the same, although the parameterization is different. See Section 7.7 of ^[18] and references therein for an analytical treatment of desingularization with time scaling. The basic idea for (4.1) is the same as the transformation applied to (4.3) in ^[14]. The system (4.2) has an equilibrium $(x, y) = (0, 0)$. Therefore, the problem of understanding the dynamical system near the origin, which is the singular point of (2.2), that was left unsolved in ^[12], has been reduced to the problem of investigating the dynamical system near the origin, which is the equilibrium of (4.2). The Jacobian matrix (linearized matrix) at the equilibrium $(x, y) = (0, 0)$ of (4.2) is

Thus, the equilibrium $(x, y) = (0, 0)$ is not hyperbolic. To obtain the dynamical system near the origin, the following transformation is used to perform the directional blow-up technique:

For the blow-up technique, see, for instance, ^[15,16,17] and references therein. Although the polar blow-up technique using the polar coordinates is used in ^[14], the directional blow-up technique is used in this paper for the sake of computational simplicity.

Since we are interested in the dynamical system near the origin in $\{(x, y) \mid x\ge 0, y\ge 0\}$, we consider the dynamical system in $\{\overline{x} = 1\}$ and $\{\overline{y} = 1\}$ combine them. See Figure 3 for a schematic diagram of the directional blow-up technique.

The directional blow-up technique is a geometric approach that covers all of the information at $\{r = 0\}$ by virtually inflating the origin (corresponding to the dotted line $\{r = \sqrt{x^{2}+y^{2}} = 0\}$ in Figure 3) and projecting it onto the local coordinates $\{\bar{x} = 1\}$ and $\{\bar{y} = 1\}$. It corresponds to (4.3) as a transformation of the differential equation. In the next subsection, we examine the dynamical system at the local coordinate $\{\bar{y} = 1\}$ and give the dynamical system at $\{r = 0\}\, \backslash\, \{y = 0\}$. Then, to examine the dynamical system in the excluded neighborhood, the dynamical system at $\{\bar{x} = 1\}$ is studied in Subsection 4.2.

4.1.   Dynamics on the chart $\{\bar{y}=1\}$

By the change of coordinates $x = \varepsilon \bar{x}$ and $y = \varepsilon$,

holds. The time-rescaling $ds/d\tau = \varepsilon$ yields

where $\varepsilon_{s} = d\varepsilon/ds$ and $\bar{x}_{s} = d\bar{x}/ds$. The equilibria of (4.4) on $\{\varepsilon = 0\}$ are

● $E_{1}: (\varepsilon, \bar{x}) = (0, 0)$ for $b\ge 1$.

● $E_{1}: (\varepsilon, \bar{x}) = (0, 0)$ and $E_{2}: (\varepsilon, \bar{x}) = (0, b^{-1}(1-b))$ for $0 < b < 1$.

Let $J_{1}$ (resp. $J_{2}$) be the Jacobian matrix of the vector fields (4.4) at $E_{1}$ (resp. $E_{2}$). So, we have

In the neighborhood of the equilibrium $E_{1}: (\varepsilon, \bar{x}) = (0, 0)$ in (4.4), $\{\varepsilon = 0\}$ ($\bar{x}$-axis) is a center manifold. Also, since

holds, we see that in (4.4), if we have an initial value on $\{(\varepsilon, \bar{x}) \mid \varepsilon > 0, \bar{x} > 0\}$, the solution can not enter $\{(\varepsilon, \bar{x}) \mid \varepsilon > 0, \bar{x} < 0\}$. When $b\ge 1$,

holds. Therefore, the solution with initial values on $\{(\varepsilon, \bar{x}) \mid \varepsilon > 0, \bar{x} > 0\}$ can only go to $E_{1}$. For $0 < b < 1$, the case is divided by $b$ and $\mu$, which is determined from the dynamical system of $E_{2}$ described next.

The eigenvalues of $J_{2}$ are $b^{-1}(1-b-\mu)$ and $-b^{-1}(1-b)^{2} < 0$. The corresponding eigenvectors are

with $T$ as the transpose matrix. In the following, we examine the dynamics near the equilibrium point $E_{2}$ for $0 < b < 1$. The dynamics are inherently different depending on the sign of $1-b-\mu$, because $b^{-1}(1-b-\mu)$ is the eigenvalue of $J_2$.

(ⅰ) If $1-b-\mu < 0$ holds, then $E_{2}$ is asymptotically stable.

(ⅱ) If $1-b-\mu > 0$ holds, then $E_{2}$ is a saddle.

(ⅲ) If $1-b-\mu = 0$ holds, then $E_{2}$ is not hyperbolic.

Let us consider the case $1-b-\mu = 0$. The dynamical system near the equilibrium $E_{2}$ can be understood by the center manifold theorem. For details on the center manifold theorem, see ^[19,20].

For the application of the center manifold theorem, we first introduce the following transformation that shifts the equilibrium $E_{2}$ to $(\varepsilon, \bar{x}) = (0, 0)$:

We then have

Therefore, the problem reduces to considering a dynamical system near the equilibrium $(\tilde{\varepsilon}, \tilde{x}) = (0, 0)$ in (4.5). The linearized matrix at this equilibrium is

the eigenvalues are $0$ and $-b^{-1}\mu^{2}$, and the corresponding eigenvectors are ${\mathbf{v}}_{1} = (b^{-1}\mu^{2}, b^{-3}\mu(2b-1))^{T}$ and ${\mathbf{v}}_{2} = (0, 1)^{T}$, respectively. Let $P = \left(\begin{array}{cc} {\mathbf{v}}_{1} & {\mathbf{v}}_{2} \end{array}\right)$, then

The matrix $P$ is invertible and

Next, the center manifold theorem is put into a form that can be applied to understand dynamical systems near the origin, which is the equilibrium of (4.5). We then have

by (4.5). By multiplying both sides by $P^{-1}$ from the left,

holds. Setting $(\hat{\varepsilon}, \hat{x})^{T} = P^{-1}(\tilde{\varepsilon}, \tilde{x})^{T}$, we obtain

and

with

The center manifold theorem is applied to understand the dynamical system near the origin of equation (4.6). The approximation of the center manifold of (4.6) is

The dynamics on the center manifold is topologically equivalent to the dynamics of the following equation:

Thus, by following back the previous transformation above, the approximation of the center manifold near $E_{2}$ of (4.4) at $1-b-\mu = 0$ is

The dynamics on the center manifold is topologically equivalent to the dynamics of the following equation:

Let $h(b) = b^{2}+b-1 = 0$. The coefficient of $\varepsilon^{2}$ is $h(b)/[b^{2}(1-b)]$. One can see that $h(b) = 0$ for $b = b^{*} = 2^{-1}(\sqrt{5}-1)\in (2^{-1}, 1)$. One has $h(b) < 0$ when $0 < b < b^{*}$, $h(b) > 0$ when $b^{*} < b < 1$, and $h(b) = 0$ when $b = b^{*}$. Note that

holds when $b = b^{*}$.

From the above, the schematic diagram of the dynamical system of (4.4) is shown in Figure 4. In ^[14], the degenerate case is not studied. The present paper uses a directional blow-up technique and precisely computes the approximation of the center manifold to give the dynamical systems of (4.4), including the degenerate case $b = b^{*}$.

4.2.   Dynamics on the chart $\{\bar{x}=1\}$

By the transformation $x = \varepsilon$, $y = \varepsilon \bar{y}$, and the time-scale transformation $ds/d\tau = \varepsilon$, we obtain

where $\varepsilon_{s} = d\varepsilon/ds$ and $\bar{x}_{s} = d\bar{x}/ds$. The equilibria of (4.9) on $\{\varepsilon = 0\}$ are:

● When $b\ge 1$, there is no equilibrium in the range $\bar{y}\ge 0$. In addition, we obtain the following:

● When $0 < b < 1$, there exists the equilibrium $E_{3}: (\varepsilon, \bar{y}) = (0, b(1-b)^{-1})$. Let $J_{3}$ be the Jacobian matrix of the vector fields (4.9) at $E_{3}$. Then, we have

The eigenvalues of $J_{3}$ are $(1-b-\mu)(1-b)^{-1}$ and $-(1-b) < 0$, and the dynamical system in the neighborhood of $E_{3}$ varies with the sign of $1-b-\mu$. When $1-b-\mu < 0$, $E_{3}$ is asymptotically stable, and when $1-b-\mu > 0$, $E_{3}$ is a saddle. When $1-b-\mu = 0$, the dynamical system near $E_{3}$ can be understood by the same argument as in Subsection 4.1.

From the above, combining the dynamical system in the local coordinates $\{\bar{x} = 1\}$ and $\{\bar{y} = 1\}$, we obtain the dynamical system on $\{r = 0\}$. Finally, we obtain the dynamical system near the origin of (2.2) and (4.2). See also Figure 5.

4.3.   Proof of theorems

A combination of the results obtained in Subsections 4.1 and 4.2 gives the dynamical system near the origin of (2.2) and (4.2). See also Figure 5. The previous study ^[12] shows the boundedness of the solution (Proposition 2.1), the non-existence of periodic orbits (Proposition 2.2), the classification of the existence of positive equilibria (Proposition 2.4), and their local stability (Proposition 2.5). The classification of the existence and stability of positive equilibria is described in 2. See Appendix A for details.

With this information and the fact that (2.2) and (4.2) are two-dimensional systems, we can use Poincaré-Bendixson's theorem (see, for instance, ^[20] and references therein). Application of this theorem yields a complete understanding of the global behavior of the solution (convergence to the origin and global stability of the positive equilibria). Consequently, Theorem 3.1 and Theorem 3.2 are derived. This completes the proof of the main result in this paper.

5.   Discussion

In this paper, the dynamics of the system (1.2) near the origin $(x, y) = (0, 0)$ is investigated using a time-scale transformation, directional blow-up technique, and the center manifold theorem. System (1.2) is a mathematical model that describes the population dynamics of wild and sterile mosquitoes and estimates the proportional release of sterile mosquitoes in the infanticide technique.

Proportional releasing has been studied in the models of SIT ^[6,9,10,11]. The proportional release would be a possible control strategy based on the surveillance of the wild mosquito population. In the paper ^[6], it is shown that periodic oscillation of the population is possible for the case of proportional release, but not for the case of constant release. In the system (1.2), due to the proportional release and mating function, the trivial equilibrium is a singular point, where the standard linearized stability analysis does not determine the stability of the extinction equilibrium. In this paper, we treat this problem by using the idea of time-scale transformation in ^[14] and the method of a directional blow-up technique. System (1.2) is a simple mathematical model for population dynamics of the wild and sterile mosquitoes. The modeling of the sterile mosquito population dynamics could be improved. In ^[4,5], no intraspecific competition is assumed in sterile mosquitoes and the mortality rate of sterile mosquitoes is larger than that of the wild mosquitoes. There are also several controlling strategies of the sterile male mosquitoes. Periodically impulsive releases and feedback control (closed-loop control) are considered in ^[4,5].

The dynamical system near the origin is derived from the directional blow-up technique and is precisely analyzed using the standard form transformation and the center manifold theorem. The boundedness of solutions, the non-existence of periodic solutions, the classification of positive equilibria excluding the origin, and their local stability are discussed in ^[12]. The combination of these arguments with the Poincaré-Bendixson theorem gives a classification of the global behavior of the solutions of (1.2) in the $\mu b$-parameter plane. See Theorem 3.1, Theorem 3.2, and Figure 2 for details.

From the results of ^[12], one can conjecture that in the parameter region (Ⅰ) in Figure 1, the solution converges to the origin irrespective of the initial values. However, the convergence has not been rigorously proven. In this paper, by combining a number of nontrivial methods from dynamical systems theory as described above, we prove that the solution converges to the origin regardless of the initial condition. We also obtain rigorous results for the parameter regions (Ⅱ, Ⅴ, Ⅵ) in Figure 1, where we prove the global asymptotic stability of the positive equilibrium. The classification of the global behavior of the solution provides a qualitative understanding of the population dynamics of wild and sterile mosquitoes. Although the mathematical model is simplified to demonstrate the effectiveness of the SIT method in actual ecology, this paper is significant in that it provides results of mathematical analysis that are sufficiently robust to future model complexity.

From the global behavior of the solution in the simple model, we can obtain insight into the release strategy of infertile mosquitoes by $b$ and $\mu$. That is, since $\mu$ is not a parameter we can control, we can show the release strategy for the release rate $b$ with initial values $x(0)$ and $y(0)$ given $\mu$.

(ⅰ) For $b$ and $\mu$ in the region (A), the population density of wild and sterile mosquitoes converge to $0$. Thus pest control is successful, irrespective of the initial conditions.

(ⅱ) The regions (B) and (C) correspond to the case where the natural death rate of mosquitoes $\mu$ is small (strong mosquito survivability), and a positive equilibrium is stable in addition to the origin, so, the pest control may not be successful. This may give us a direction to improve the mathematical model so that the SIT method can be successfully used for pest control. On the other hand, if the initial number of wild mosquitoes is small, the release rate $b$ can be increased to achieve extinction. The mathematical analysis suggests that such a bistable scenario is possible in the actual ecology, and the implementation of the SIT method would be an issue for such a case. The mathematical analysis in this paper does not provide a criterion for the critical amount of initial wild mosquitoes.

(ⅲ) The regions (D) and (E) correspond to the case where the natural death rate of mosquitoes $\mu$ is small (mosquito viability is strong). In these regions, the positive equilibrium is asymptotically stable and the origin is unstable indicating that wild and sterile mosquitoes cannot be exterminated.

If the natural death rate of mosquitoes $\mu$ is fixed, varying the release rate $b$ makes the following considerations.

(ⅳ) When $\mu\ge \mu^{*} = 2^{-1}(3-\sqrt{5})$, mosquitoes will not become extinct if the release rate $b$ is small, but mosquitoes will become extinct for any initial condition if $b$ is $b\ge 1-\mu$. This means that the larger the mosquito's natural mortality rate $\mu$ (i.e., the weaker the mosquito's ability to survive), the more it can be exterminated by simply releasing a small amount of sterile mosquitoes.

(ⅴ) When $\mu < \mu^{*}$, no matter how much the release rate $b$ is increased, it is not always possible to exterminate wild mosquitoes. If we increase the release rate $b$ and enter the region (A), we can successfully exterminate the pests no matter what the initial value is. However, if we enter the region (B) or (C), the release strategy depends on the initial population density of wild and infertile mosquitoes. This is as described in (ⅱ) above.

The model (1.2) is a simplified mathematical model of population dynamics for wild and sterile mosquitoes. The question of whether the mathematical model can be improved so that the natural death rate $\mu$ and the intraspecific competitive death rate $\xi$ for wild and sterile mosquitoes can be varied, and the global behavior of the solution in such cases, is an important one, and the approach taken in this paper is expected to be effective. The introduction of models that reproduce the biological situation, the presence or absence of intraspecific competition, the method of introducing release rates, and the discussion of models in previous studies ^[4,5] are issues for future work.

This paper provides a simple and meaningful analysis of the behavior of the solution of a system of differential equations with a ratio-dependent mathematical model for pest control related to the SIT method. The estimation of the effect of the SIT method on the growth process of mosquitoes and the implementation of simulations based on the model in this paper are expected to be further developed. We believe that the results of this paper will be meaningful as a mathematical analysis when the model becomes more realistic through a biological viewpoint.

Use of AI tools declaration

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

Acknowledgments

The authors would like to thank the anonymous referees who provided thoughtful comments, which improved the original manuscript. YI was partially supported by JSPS KAKENHI Grant Number JP22KJ2844. YN was partially supported by JSPS KAKENHI Grant-in-Aid for Scientific Research (C) 20K03734. This work is partially supported by the Meiji Institute for Advanced Study of Mathematical Sciences (MIMS) Center for Mathematical Modeling and Applications (CMMA) Interdisciplinary Research Support Program for Life Sciences and Mathematical Sciences.

Conflict of interest

The authors declare there are no conflicts of interest.

Appendix

A.   The existence of positive equilibrium

To find equilibria of the system (2.2), we consider (2.3) for $x > 0$ and $y > 0$. We call the solution of (2.3) in $\text{int}\, \mathbb{R}_+^2$ the positive equilibrium of (2.2). We state the proof of Proposition 2.4.

Lemma A.1. If $\mu\geq1$, then there is no positive equilibria.

Consider the case that $0 < \mu < 1$ holds. Note that the positive equilibrium satisfies the following inequality:

Since (2.3) holds, we obtain (2.4):

noting that $x \neq 0$. In the paper ^[12], the authors study the root of (2.4).

We here consider the transformation of a variable (2.5):

Then, from (2.5), we obtain (2.6) and (2.7):

We consider the solution of (2.6) for $z\in (0, \mu^{-1}-1)$, corresponding to $x\in (0, 1-\mu)$.

Lemma A.2. It holds that

Thus $g^{\prime}(0) = \mu$ holds. If $0 < \mu < 1$ holds, then there exists a unique zero of $g^{\prime}(z)$ for $z > 0$.

Proof. Direct computation from (2.7) shows the formula (A.1). From (A.1), we have

Since $0 < \mu < 1$ holds, $g^{\prime\prime}(0) = 2-2\mu > 0$. It is easy to see that $g^{\prime\prime}(z)$ is monotone decreasing and $\lim_{z\to\infty}g^{\prime\prime}(z) = -\infty$. Thus $g^{\prime\prime}(z)$ has a unique zero for $z\geq0$, which we denote by $\overline{z}$. One can see that $g^{\prime}(0) = \mu > 0$ and $g^{\prime}(z)$ attains a unique local maximum at $z = \overline{z} > 0$. Therefore, $g^{\prime}(z)$ monotonically decreases for $z\geq\overline{z}$. Since $\lim_{z\to\infty}g^{\prime\prime}(z) = -\infty$, there exists a unique zero of $g^{\prime}(z)$ for $z > \overline{z} > 0$. □

Let us define the interval $I$ as

We consider the solution of (2.7) on the interval $I$. Note that

We now obtain the following result.

Lemma A.3. Let $0 < \mu < 1$ hold. If $0 < b < 1-\mu$, then (2.7) has a unique root on $I$.

Proof. From Lemma A.2, $g(z)$ either monotonically increases on $I$ or attains a local maximum on $I$. From (A.2), in any case, if $0 < b < 1-\mu$, then (2.6) has a unique root on $I$. □

We denote by $z^{*}(\mu)$ the unique zero of $g^{\prime}(z) = 0$, i.e.,

holds. Note that, from Lemma A.3, for $0 < \mu < 1$, $z^{*}(\mu)$ exists and $z^*(\mu) > 0$.

From (A.1), one can compute

and that $-\mu^{2}+3\mu-1 = 0$ has a root

for $0 < \mu < 1$. Note that $2^{-1}(3-\sqrt{5}) = 0.381966\dots$. The sign of $g^{\prime}\left(\mu^{-1}-1\right)$ determines the order of $z^{*}(\mu)$ and $\mu^{-1}-1$ as follows.

Lemma A.4. Let $0 < \mu < 1$ hold. Then it holds that

We now have the following propositions concerning the number of roots of (2.6). Since the proof is straightforward, we omit the proof.

Lemma A.5. Let us assume that $0 < \mu < 1$ and $b = 1-\mu$ hold. Then (2.6) has a root $z = \mu^{-1}-1$.

(a) Suppose that $0 < \mu < 2^{-1}(3-\sqrt{5})$ holds. Then (2.6) has one root on $I$.

$(b)$ Suppose that $2^{-1}(3-\sqrt{5}) \leq \mu < 1$ holds. Then (2.6) has no roots on $I$.

Lemma A.6. Let us assume that $0 < \mu < 1$ and $b > 1-\mu$ hold.

(c) Suppose that $0 < \mu < 2^{-1}(3-\sqrt{5})$ holds.

(c-1) If $1-\mu < b < g(z^{*}(\mu))$, then (2.6) has exactly two positive roots on $I$.

(c-2) If $b = g(z^{*}(\mu))$, then (2.6) has a unique root $z^{*}(\mu)\in I$.

(c-3) If $b > g(z^{*}(\mu))$, then (2.6) has no roots on $I$.

(d) Suppose that $2^{-1}(3-\sqrt{5}) \leq\mu < 1$ holds. Then (2.6) has no roots on $I$.

Finally, combining Lemmas A.3, A.5, and A.6, we obtain the following result for the existence of the positive equilibrium for (2.2).

In the $\mu b$-parameter plane, the region enclosed by the line $b = 1-\mu$ and the curve $b = g(z^{*}(\mu))$ for $0 < \mu < 2^{-1}(3-\sqrt{5})$ is the parameter region where the system has two positive equilibria. The curve $b = g(z^{*}(\mu))$ for $0 < \mu < 2^{-1}(3-\sqrt{5})$ can be parametrized by $z$ as follows. Since $z = z^{*}(\mu)$ is a unique root of $g^{\prime}(z) = 0$ for $0 < \mu < 1$, from (A.1), its inverse function is given as

Furthermore, substituting (A.3) into $g(z)$ in (2.6), we obtain

From Lemma A.4, the range of $z$ is

So we complete the proof of Proposition 2.4.

B.   Stability of the positive equilibrium

In this section, we consider the locally asymptotic stability of the positive equilibrium for (2.2). The Jacobian matrix for the positive equilibrium is computed in ^[12]. Using the condition (2.3) that the equilibrium satisfies, we obtain a simpler expression for the Jacobian matrix. Our objective in this section is to give the proof of Proposition 2.5.

Lemma B.1. The Jacobian matrix evaluated at the positive equilibrium of (2.2) is given as

where $(x, y)$ denotes the positive equilibrium of (2.2), where $x > 0$ and $y > 0$.

Proof. The Jacobian matrix evaluated at the positive equilibrium $(x, y)$ is given as

From (2.3), it holds that $x+y = x(\mu+x)^{-1}$ and $y = x[(\mu+ x)^{-1}-1]$. Thus

Therefore, we obtain the Jacobian matrix as in (B.1). □

We denote by $J$ the Jacobian matrix given as in (B.1). We obtain the following result.

Lemma B.2. For the Jacobian matrix $J$ given as in (B.1), it holds that

and

Proof. It is straightforward to obtain (B.2) from (B.1). We compute $D: = \text{det}\, J$.

From the conditions (2.3) and (2.4), we have

and

Therefore, we see

Hence, we obtain (B.3). □

Since $\text{trace}\, J < 0$, the positive equilibrium is asymptotically stable if $\text{det}\, J > 0$ and it is unstable if $\text{det}\, J < 0$. We recall that $z$ is related to $x$ via the transformation (2.6), where $x$ denotes the component of the positive equilibrium of (2.2). We now obtain the following result.

Lemma B.3. Let $z$ be given by (2.6), where $x$ denotes the component of the positive equilibrium of (2.2). If $g^\prime(z) > 0$, then the positive equilibrium is asymptotically stable, and if $g^\prime(z) < 0$, then it is unstable.

Proof. We express $g^{\prime}(z)$ in terms of $x$, using (2.6), as follows:

Therefore, from (B.3), we see the relation that $g^\prime (z) > 0\, (<0)$ if and only if $\text {det}\, J > 0\, (<0)$. □

Consider the case that (2.2) has only one positive equilibrium $(x, y)$. We can see that, for the case, $g^\prime(z) > 0$, thus the positive equilibrium is asymptotically stable. Suppose that (2.2) has two positive equilibria, $({x}_{1}, {y}_{1})$ and $({x}_{2}, {y}_{2})$. Let

It is easy to see that $z_{2} < z_{1}$. Therefore, we have

which implies the consequence that $({x}_{1}, {y}_{1})$ is unstable while $({x}_{2}, {y}_{2})$ is asymptotically stable. Therefore, we complete the proof of Proposition 2.5.