<figure><img src="image_1.png" /></figure>

HAL
open science

# Efficacy of the Sterile Insect Technique in the presence of inaccessible areas: A study using two-patch models

Pierre-Alexandre Bliman, Nga Nguyen, Nicolas Vauchelet

► To cite this version:

Pierre-Alexandre Bliman, Nga Nguyen, Nicolas Vauchelet. Efficacy of the Sterile Insect Technique in the presence of inaccessible areas: A study using two-patch models. 2024. hal-04525237

HAL Id: hal-04525237

https://hal.science/hal-04525237

Preprint submitted on 28 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Efficacy of the Sterile Insect Technique in the presence of inaccessible areas: A study using two-patch models

Pierre-Alexandre Bliman, Nga Nguyen, Nicolas Vauchelet

► To cite this version:

Pierre-Alexandre Bliman, Nga Nguyen, Nicolas Vauchelet. Efficacy of the Sterile Insect Technique in the presence of inaccessible areas: A study using two-patch models. Mathematical Biosciences, inPress, pp.109290. 10.1016/j.mbs.2024.109290. hal-04525237v2

HAL Id: hal-04525237

https://hal.science/hal-04525237v2

Submitted on 10 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Efficacy of the Sterile Insect Technique in the presence of inaccessible areas: A study using two-patch models

Pierre-Alexandre Bliman<fnref n="1" />, Nga Nguyen<fnref n="2" />, Nicolas Vauchelet<fnref n="3" />

## Abstract

The Sterile Insect Technique (SIT) is one of the sustainable strategies for the control of disease vectors, which consists of releasing sterilized males that will mate with the wild females, resulting in a reduction and, eventually a local elimination, of the wild population. The implementation of the SIT in the field can become problematic when there are inaccessible areas where the release of sterile insects cannot be carried out directly, and the migration of wild insects from these areas to the treated zone may influence the efficacy of this technique. However, we can also take advantage of the movement of sterile individuals to control the wild population in these unreachable places. In this paper, we derive a two-patch model for *Aedes* mosquitoes where we consider the discrete diffusion between the treated area and the inaccessible zone. We investigate two different release strategies (constant and impulsive periodic releases), and by using the monotonicity of the model, we show that if the number of released sterile males exceeds some threshold, the technique succeeds in driving the whole population in both areas to extinction. This threshold depends on not only the biological parameters of the population but also the diffusion between the two patches.

**Keywords:** sterile insect technique, metapopulation model, monotone dynamical systems

## 1 Introduction

Mosquitoes of genus *Aedes aegypti* and *Aedes albopictus* play a crucial role in transmitting various arboviruses to humans including dengue, chikungunya, and Zika virus. Unfortunately, there are no specific vaccines or drugs available for these diseases. Consequently, the primary prevention lies in controlling the mosquito population [37]. However, traditional insecticide-based methods have limitations, prompting the need for innovative and sustainable strategies [2], [10]. Biological controls involve releasing large numbers of mosquitoes that are either sterile or incapable of transmitting diseases, which recently gained much attention. The Sterile Insect Technique is among these sustainable alternative methods which consist of the release of sterilized male mosquitoes that will mate with wild females [21], [19]. These wild females, unable to lay viable eggs, will gradually drive the wild population to decline. The efficacy of SIT relies on a comprehensive understanding of the vector behavior, as well as accurate modeling of its dispersal, to optimize the release strategies.

<footnote n="1">¹MAMBA, Inria Paris; LJLL, Sorbonne University, CNRS, 5 Place Jussieu, 75005 Paris, France (email : pierre-alexandre.bliman@inria.fr).</footnote>

<footnote n="2">²LAGA, CNRS UMR 7539, Institut Galilée, Université Sorbonne Paris Nord, 99 avenue Jean-Baptiste Clément, 93430 Villetaneuse; MAMBA, Inria Paris, Laboratoire Jacques Louis-Lions, Sorbonne Université, 5 place Jussieu, 75005 Paris, France (email : thiquynhnga.nguyen@math.univ-paris13.fr).</footnote>

<footnote n="3">³LAGA, CNRS UMR 7539, Institut Galilée, Université Sorbonne Paris Nord, 99 avenue Jean-Baptiste Clément, 93430 Villetaneuse, France (email : vauchelet@math.univ-paris13.fr).</footnote>
# Efficacy of the Sterile Insect Technique in the presence of inaccessible areas: A study using two-patch models

P.A. Bliman<fnref n="1" />, N. Nguyen<fnref n="2" />, <fnref n="1" />, and N. Vauchelet<fnref n="2" />

September 10, 2024

## Abstract

The Sterile Insect Technique (SIT) is one of the sustainable strategies for the control of disease vectors, which consists of releasing sterilized males that will mate with the wild females, resulting in a reduction and, eventually a local elimination, of the wild population. The implementation of the SIT in the field can become problematic when there are inaccessible areas where the release of sterile insects cannot be carried out directly, and the migration of wild insects from these areas to the treated zone may influence the efficacy of this technique. However, we can also take advantage of the movement of sterile individuals to control the wild population in these unreachable places. In this paper, we derive a two-patch model for *Aedes* mosquitoes where we consider the discrete diffusion between the treated area and the inaccessible zone. We investigate two different release strategies (constant and impulsive periodic releases), and by using the monotonicity of the model, we show that if the number of released sterile males exceeds some threshold, the technique succeeds in driving the whole population in both areas to extinction. This threshold depends on not only the biological parameters of the population but also the diffusion between the two patches.

**Keywords** – sterile insect technique, metapopulation model, monotone dynamical systems

## 1 Introduction

Mosquitoes of genus *Aedes aegypti* and *Aedes albopictus* play a crucial role in transmitting various arboviruses to humans including dengue, chikungunya, and Zika virus. Existing treatments are only symptomatic, and available vaccines (e.g. Dengvaxia for dengue) have a lot of constraints [1]. Consequently, the primary prevention lies in controlling the mosquito population [2]. However, traditional insecticide-based methods have limitations, prompting the need for innovative and sustainable strategies [3], [4]. Biological controls involve releasing large numbers of mosquitoes that are either sterile or incapable of transmitting diseases, which recently gained much attention. The Sterile Insect Technique is among these sustainable alternative methods which consist of the release of sterilized male mosquitoes that will mate with wild females [5], [6]. These wild females, unable to lay viable eggs, will gradually drive the wild population to decline. The efficacy of SIT relies on a comprehensive understanding of the vector behavior, as well as accurate modeling of its dispersal, to optimize the release strategies.

Spatial heterogeneity in mosquito populations and mosquito-borne diseases occurs due to differences in the quality and quantity of their habitats, as well as variations in host density, temperature, and rainfall [7], [8], [9]. Especially, the number and accessibility of sites where mosquitoes lay their eggs play a significant role in determining the size of adult mosquito populations by increasing the carrying capacity of the environment [10]. Developing models that capture mosquito behavior in response to environmental heterogeneity is crucial for designing effective control strategies, especially in the face of rapid global land-use changes. Models using monotone dynamical systems were introduced (see e.g. [11], [12], [13]) and applied efficiently (see e.g. [14], [15], [16]) to study the SIT. Not many mosquito modeling studies have incorporated migration or dispersal effects due to insufficient information on individual movement in the field as well as the complex analysis of models. Most of them used the diffusion approach, which considers space as a continuous variable. They were first developed in one-dimensional space using scalar reaction-diffusion equations [17], [18], then extended to sex-structured compartmental systems to consider the different behaviors of aquatic phases, wild females, males, and sterile

<footnote n="1">¹MAMBA, Inria Paris; LJLL, Sorbonne University, CNRS, 5 Place Jussieu, 75005 Paris, France</footnote>

<footnote n="2">²LAGA, CNRS UMR 7539, Institut Galilée, University Sorbonne Paris Nord, 99 avenue Jean-Baptiste Clément, 93430 Villetaneuse, France</footnote>
Spatial heterogeneity in mosquito populations and mosquito-borne diseases occurs due to differ-
ences in the quality and quantity of their habitats, as well as variations in host density, temperature, and
rainfall [15], [34], [31]. Especially, the number and accessibility of sites where mosquitoes lay their
eggs play a significant role in determining the size of adult mosquito populations by increasing the car-
rying capacity of the environment [1]. Developing models that capture mosquito behavior in response
to environmental heterogeneity is crucial for designing effective control strategies, especially in the face
of rapid global land-use changes. Models using monotone dynamical systems were introduced (see e.g.
[7], [18], [35]) and applied efficiently (see e.g. [11], [8], [12]) to study the SIT. Not many mosquito
modeling studies have incorporated migration or dispersal effects due to insufficient information on
individual movement in the field as well as the complex analysis of models. Most of them used the
diffusion approach, which considers space as a continuous variable. They were first developed in one-
dimensional space using scalar reaction-diffusion equations [30], [25], then extended to sex-structured
compartmental systems to consider the different behaviors of aquatic phases, wild females, males, and
sterile males (see e.g. [4], [6], [24]) and in higher dimension (see e.g. [17], [5]). However, it remains
challenging to explicitly incorporate the factors that affect the movement of sterile males. For instance,
when resources are concentrated in patches or distinct locations, a metapopulation approach in which
we treat space as a discrete set of patches and describe how the population on each patch varies with
time is more suitable for modeling mosquito dispersal [9], [28], [29].

The application of the SIT in the field encounters a difficulty of the limitation in space when there are
some inaccessible areas where people can not release sterile insects directly. For example, mosquitoes
of the genus *Aedes polynesiensis* primarily exploit land crab burrows for oviposition in certain French
Polynesian atolls [13], [22], [23], [20]. The larvae in the crab burrows emerge into adult mosquitoes
that can fly out to search for food and human blood for fertility. However, one standout advantage of the
SIT is that it relies on the natural ability of the male mosquitoes to move, locate, and mate with females.
This behavior will take place in those areas that cannot be reached with conventional control techniques
(i.e. insecticides). Therefore, we are interested in the mosquito population dynamics in the presence of
such reservoirs and the elimination of the whole population while considering that the released sterile
males can fly into unreachable sites. The patchy models with discrete diffusion mentioned above are a
useful approach to describe the mosquito dynamic taking into account the inaccessibility to the burrows.
We develop a two-patch model and in each patch, we consider a monotone dynamical system inspired
by the models in [35] where the population is divided into different compartments characterizing the
aquatic phase, wild females, wild males, and sterile males. Except for the aquatic phase, individuals in
other states move between patches at specific rates. The SIT is only carried out in the first patch and only
affects the second one through these natural movements. Two-patch models were used to study the same
problem in [38], where they considered a simple scalar equation to describe the population dynamics in
each patch. Our model provides a better understanding of how the dynamics of each stage influence the
result of the control method. However, the complexity of our system does not allow us to obtain the full
analysis of the model like what has been done in [38].

In the present work, we are interested in how to guarantee the successful elimination of the SIT
in both areas and how the diffusion rates as well as other biological parameters influence the efficacy.
To tackle this problem, we focus on studying the global stability of the extinction equilibrium in our
system. Results of global asymptotic behavior for the single-species model depending on the discrete
diffusion were provided in the literature[3], [36], [27]. Lyapunov's second method was used in [26]
to investigate the multi-species system with discrete diffusion. Many works have been done to design
robust strategies for releasing sterile males to drive a population to elimination [11], [12]. We extend
these control strategies to our two-patch system and prove the sufficient conditions for both constant
continuous and periodic impulsive releases to drive the whole system to extinction. We obtain that when
the number of released sterile males exceeds some threshold, the populations in both the treated and the
inaccessible zone reach elimination. We also show in the present work how the diffusion rates between
males (see e.g. [19], [20], [21]) and in higher dimension (see e.g. [22], [23]). However, it remains challenging to explicitly incorporate the factors that affect the movement of sterile males. For instance, when resources are concentrated in patches or distinct locations, a metapopulation approach in which we treat space as a discrete set of patches and describe how the population on each patch varies with time is more suitable for modeling mosquito dispersal [24], [25], [26].

The application of the SIT in the field encounters a difficulty of the limitation in space when there are some inaccessible areas where people can not release sterile insects directly. For example, mosquitoes of the genus *Aedes polynesiensis* primarily exploit land crab burrows for oviposition in certain French Polynesian atolls [27], [28], [29], [30]. The larvae in the crab burrows emerge into adult mosquitoes that can fly out to search for food and human blood for fertility. Another example is the inaccessibility in heavily forested or narrow mountainous terrain, where sterile insects must be released via helicopters owing to their better maneuverability [31]. However, one standout advantage of the SIT is that it relies on the natural ability of the male mosquitoes to move, locate, and mate with females. This behavior will take place in those areas that cannot be reached with conventional control techniques (i.e. insecticides). Therefore, we are interested in the mosquito population dynamics in the presence of such reservoirs and the elimination of the whole population while considering that the released sterile males can fly into unreachable sites. The patchy models with discrete diffusion mentioned above are a useful approach to describe the mosquito dynamic taking into account the inaccessibility to the burrows. We develop a two-patch model and in each patch, we consider a monotone dynamical system inspired by the models in [13] where the population is divided into different compartments characterizing the aquatic phase, wild females, wild males, and sterile males. Except for the aquatic phase, individuals in other states move between patches at specific rates. The SIT is only carried out in the first patch and only affects the second one through these natural movements. Two-patch models were used to study the same problem in [32], where they considered a simple scalar equation to describe the population dynamics in each patch. Our model provides a better understanding of how the dynamics of each stage influence the result of the control method. However, the complexity of our system does not allow us to obtain the full analysis of the model like what has been done in [32].

In the present work, we are interested in how to guarantee the successful elimination of the SIT in both areas and how the diffusion rates as well as other biological parameters influence the efficacy. To tackle this problem, we focus on studying the global stability of the extinction equilibrium in our system. Results of global asymptotic behavior for the single-species model depending on the discrete diffusion were provided in the literature [33], [34], [35]. Lyapunov's second method was used in [36] to investigate the multi-species system with discrete diffusion. Many works have been done to design robust strategies for releasing sterile males to drive a population to elimination [14], [16]. We extend these control strategies to our two-patch system and prove the sufficient conditions for both constant continuous and periodic impulsive releases to drive the whole system to extinction. We obtain that when the number of released sterile males exceeds some threshold, the populations in both the treated and the inaccessible zone reach elimination, and we show how this critical value depends on the diffusion rates between two areas and other biological intrinsic parameters. In the original mathematical model provided by Knipling [5], the population elimination depends on an overflooding ratio of the number of sterile males released to the initial wild male population size. In our result, the threshold number of sterile males does not depend on the initial level of the wild population due to the global stability of the zero steady state. Our results may help estimate the possibility and the amount of sterile mosquitoes necessary to complete elimination in the presence of hidden, inaccessible reservoirs.

The organization of the paper is as follows. In section 2, we present the formulation of the two-patch model and prove the monotonicity of the systems and some other preliminary results that will be applied in our proofs. Section 3 is devoted to the study of the system without sterile insects. In Theorem 3.1, we provide conditions for the persistence and extinction of the wild population on each patch. In section 4, we study the dynamics of mosquito population in the presence of the SIT with two release strategies: constant and impulsive releases. Theorem 4.2 presents sufficient conditions on the average number of sterile males released per time unit to drive the population to elimination. In section 4.1, we present the principle of our method and then apply it to prove Theorem 4.2. Section 5 is focused on the dependence of the critical number of sterile males on parameters. The results in 5.1.1 show that when the diffusion rates are large, the dynamics of the whole system are the same as in the case when there is no separation between the two sub-populations. Then, Theorem 5.3 shows that the critical number of released sterile males depends monotonically on the biological parameters. Finally, some numerical illustrations are provided in Section 6.

## 2 Model

In this section, we present the formulation of the model used to study the population dynamics in 2.1. Then, in 2.2, we provide some preliminary results that will be used later in the present work.
two areas and other biological parameters influence these conditions. These results may help estimate the possibility and the surplus of sterile mosquitoes necessary to complete elimination in the presence of hidden, inaccessible, reservoirs.

The organization of the paper is as follows. In section 2, we present the formulation of the two-patch model and prove the monotonicity of the systems and some other preliminary results that be applied in our proofs. Section 3 is devoted to the study of the system without sterile insects. In Theorem 3.1, we provide conditions for the persistence and extinction of the wild population on each patch. In section 4, we study the dynamics of mosquito population in the presence of the SIT with two release strategies: constant and impulsive releases. Theorem 4.2 presents sufficient conditions on the average number of sterile males released per time unit to drive the population to elimination. We provide the principle of the method used to treat the system in 4.1 and then apply this principle to prove Theorem 4.2. Section 5 is focused on the dependence of the critical number of sterile males on parameters. The results in 5.1.1 show that when the diffusion rates are large, the dynamics of the whole system are the same as in the case when there is no separation between the two sub-populations. Then, Theorem 5.3 shows that the critical number of released sterile males depends monotonically on the biological parameters. Finally, some numerical illustrations are provided in Section 6.

# 2 Model

In this section, we present the formulation of the model used to study the population dynamics in 2.1. Then, in 2.2, we provide some preliminary results that will be used later in the present work.

## 2.1 Formulation of the model

Consider two patches and denote $E_i, F_i, M_i$, and $M_i^s$ respectively the density of aquatic phase (eggs, larvae, pupae), fertile females, males, and sterile males on the patch $i$ depending on time $t$. We consider a two-patch model coupled by the diffusion terms as follows where the dynamic in each patch is inspired by the model in [35]

$$ \dot{E}_1 = bF_1 \left(1 - \frac{E_1}{K_1}\right) - (v_E + \mu_E)E_1, \quad (1a) $$

$$ \dot{F}_1 = r v_E E_1 \frac{M_1}{M_1 + \gamma M_1^s} - \mu_F F_1 - d_{12} F_1 + d_{21} F_2, \quad (1b) $$

$$ \dot{M}_1 = (1-r)v_E E_1 - \mu_M M_1 - \beta d_{12} M_1 + \beta d_{21} M_2, \quad (1c) $$

$$ \dot{M}_1^s = \Lambda - \mu_s M_1^s - \alpha d_{12} M_1^s + \alpha d_{21} M_2^s, \quad (1d) $$

$$ \dot{E}_2 = bF_2 \left(1 - \frac{E_2}{K_2}\right) - (v_E + \mu_E)E_2, \quad (1e) $$

$$ \dot{F}_2 = r v_E E_2 \frac{M_2}{M_2 + \gamma M_2^s} - \mu_F F_2 - d_{21} F_2 + d_{12} F_1, \quad (1f) $$

$$ \dot{M}_2 = (1-r)v_E E_2 - \mu_M M_2 - \beta d_{21} M_2 + \beta d_{12} M_1, \quad (1g) $$

$$ \dot{M}_2^s = -\mu_s M_2^s - \alpha d_{21} M_2^s + \alpha d_{12} M_1^s. \quad (1h) $$

The interpretation of the parameters used in the model, with $i, j \in \{1, 2\}$, is as below

* $\Lambda(t)$ is the number per time unit of sterile mosquitoes that are released at time $t$ on the first patch;

* the fraction $\frac{M_i}{M_i + \gamma M_i^s}$ corresponds to the probability that a female mates with a fertile male;
## 2.1 Formulation of the model

Consider two patches and denote $E_i$, $F_i$, $M_i$, and $M_i^s$ respectively the density of aquatic phase (eggs, larvae, pupae), fertilized females, wild males, and sterile males on the patch $i$ depending on time $t$. We consider a two-patch model coupled by the diffusion terms as follows where the dynamic in each patch is inspired by the model in [13, 15]

$$ \dot{E}_1 = bF_1 \left(1 - \frac{E_1}{K_1}\right) - (\nu_E + \mu_E)E_1, \quad (1a) $$

$$ \dot{F}_1 = r\nu_E E_1 \frac{M_1}{M_1 + \gamma M_1^s} - \mu_F F_1 - d_{12}F_1 + d_{21}F_2, \quad (1b) $$

$$ \dot{M}_1 = (1 - r)\nu_E E_1 - \mu_M M_1 - \beta d_{12} M_1 + \beta d_{21} M_2, \quad (1c) $$

$$ \dot{M}_1^s = \Lambda - \mu_s M_1^s - \alpha d_{12} M_1^s + \alpha d_{21} M_2^s, \quad (1d) $$

$$ \dot{E}_2 = bF_2 \left(1 - \frac{E_2}{K_2}\right) - (\nu_E + \mu_E)E_2, \quad (1e) $$

$$ \dot{F}_2 = r\nu_E E_2 \frac{M_2}{M_2 + \gamma M_2^s} - \mu_F F_2 - d_{21}F_2 + d_{12}F_1, \quad (1f) $$

$$ \dot{M}_2 = (1 - r)\nu_E E_2 - \mu_M M_2 - \beta d_{21} M_2 + \beta d_{12} M_1, \quad (1g) $$

$$ \dot{M}_2^s = -\mu_s M_2^s - \alpha d_{21} M_2^s + \alpha d_{12} M_1^s. \quad (1h) $$

In *Aedes* mosquitoes, fertilization takes place very quickly after the females have hatched (the males wait for the females to hatch near the breeding sites). So an important assumption we made in system (1) is that all female mosquitoes can mate with either wild or sterile males. The variable $F_i$ only characterizes the density of females fertilized by wild males and $E_i$ denotes the density of viable offspring in the aquatic phase. On patch $i$, the total amount of offspring emerging per time unit is $\nu_E E_i$. These correspond to the birth of males and females, with the respective quantities $(1 - r)\nu_E E_i$ and $r\nu_E E_i$. Among the $r\nu_E E_i$ females that hatch at each unit of time, a quantity $r\nu_E E_i \frac{M_i}{(M_i + \gamma M_i^s)}$ mate with a wild male and produce viable offspring (which corresponds to the first term of equations (1a) and (1e)); and the females that couple with a sterile male are no longer involved in the reproduction.

The interpretation of the parameters used in the model, with $i, j \in \{1, 2\}$, is as below

* $\Lambda(t)$ is the number per time unit of sterile mosquitoes that are released at time $t$ on the first patch;

* $\gamma \in [0, 1]$ characterizes the competitiveness of sterile males;

* $b > 0$ is the birth rate; $\mu_E > 0$, $\mu_M > 0$, and $\mu_F > 0$ denote the death rates for the mosquitoes in the aquatic phase, for adult males, and for adult females, respectively;

* $K_i$ is an environmental capacity for the aquatic phase on patch $i$, accounting also for the intraspecific competition;

* $\nu_E > 0$ is the rate of emergence;

* $r \in (0, 1)$ is the probability that a female emerges, then $(1 - r)$ is the probability that a male emerges.

* $d_{ij} > 0$ is the moving rate of female mosquitoes from patch $i$ to patch $j$; the fertile males and sterile males move slower but with proportional rates respectively $\beta d_{ij}$, $\alpha d_{ij}$ where typically $0 < \alpha < \beta < 1$ in practice.

Mosquito life cycle characteristics in both zones are linked to the zone sizes, comparable to the distance a mosquito can travel. In this way, the overall homogeneity of the species can be preserved, despite the possible effects of evolution, due to permanent mixing, so it is relevant to consider the same biological parameters $b$, $r$, $\mu_E$, $\mu_F$, $\mu_M$ in the two patches.

We recall the basic offspring number of the sub-population in one patch as introduced in [13]

$$ \mathcal{N} = \frac{br\nu_E}{\mu_F(\mu_E + \nu_E)}. \quad (2) $$

The persistence and extinction of the population in the patch depend strongly on the value of this number. In Section 3, we will show that $\mathcal{N}$ is also the basic offspring number of the whole two-patch system.

## 2.2 Preliminary results

First, we provide some definitions and denotations of the order used in the present work.
* $b > 0$ is the birth rate; $\mu_E > 0$, $\mu_M > 0$, and $\mu_F > 0$ denote the death rates for the mosquitoes in the aquatic phase, for adult males, and for adult females, respectively;

* $K_i$ is an environmental capacity for the aquatic phase on patch $i$th, accounting also for the in-traspecific competition;

* $v_E > 0$ is the rate of emergence;

* $r \in (0,1)$ is the probability that a female emerges, then $(1-r)$ is the probability that a male emerges.

* $d_{ij} > 0$ is the moving rate of female mosquitoes from patch $i$th to patch $j$th; the fertile males and sterile males move slower but with proportional rates respectively $\beta d_{ij}$, $\alpha d_{ij}$ where typically $0 < \alpha < \beta < 1$ in practice.

We recall the basic offspring number of the sub-population in one patch as introduced in [35]

$$ \mathcal{N} = \frac{brv_E}{\mu_F(\mu_E + v_E)}. \qquad (2) $$

The persistence and extinction of the population in the patch depend strongly on the value of this number. In Section, 3, we will show that $\mathcal{N}$ is also the basic offspring number of the whole two-patch system.

## 2.2 Preliminary results

First, we provide some definitions and denotations of the order used in the present work.

**Definition 2.1.** A matrix $A \in \mathcal{M}^{m \times n}$ is called non-negative, denote $A \ge 0$, if all of its entries are non-negative.

It is called positive, denote $A > 0$, if it is non-negative and there is at least one positive entry.

It is called strictly positive, denote $A \gg 0$, if all of its entries are strictly positive.

In the present work, we also use the above definition of order for vectors in $\mathbb{R}^n$. Next, we present a property of a Metzler matrix that will be used in this paper.

**Lemma 2.1.** Assume that a square matrix $A$ is Metzler and irreducible, then $e^A$ is strictly positive.

*Proof.* Since $A$ is Metzler, then there exists a constant $\delta > 0$ large enough such that $A + \delta I$ is a non-negative matrix with a positive element on the main diagonal. Moreover, $A$ is irreducible so $A + \delta I$ is also irreducible. Thus, $A + \delta I$ is primitive, that is, there exists an integer $n > 0$ such that $(A + \delta I)^n \gg 0$. Hence, we have $e^{A+\delta I} \gg 0$, and since $\delta I$ commutes with all matrices, one has $e^A = e^{A+\delta I}e^{-\delta I} \gg 0$. $\square$

We present in this section the so-called Kamke [16] or Chaplygin [14] lemma for a cooperative system (Lemma 2.2). Then, we apply this lemma to show the monotonicity of system (1) in Lemma 2.3.

**Lemma 2.2.** For any $n \in \mathbb{N}^*$, consider a smooth function $f: \mathbb{R}^n \to \mathbb{R}^n$, and a vector function $u(t)$ satisfying a differential equation

$$ \dot{u} = f(u). $$

Moreover, we assume that the above system is cooperative, that is,

$$ \frac{\partial f_i}{\partial u_j}(t) \ge 0, \quad \text{for } i \ne j, t > 0. \qquad (3) $$

If a vector function $v(t)$ satisfies a differential inequality $\dot{v} \le f(v)$ then, for initial data $v(0) \le u(0)$, we have $v(t) \le u(t)$ for all $t > 0$.
**Definition 2.1.** A matrix $A \in \mathcal{M}^{m \times n}$ is called non-negative, denoted $A \ge 0$, if all of its entries are non-negative. It is called positive, denoted $A > 0$, if it is non-negative and there is at least one positive entry. It is called strictly positive, denoted $A \gg 0$, if all of its entries are strictly positive.

In the present work, we also use the above definition of order for vectors in $\mathbb{R}^n$. Next, we recall the definition of a Metzler matrix and present a property of a Metzler matrix that will be used in this paper.

**Definition 2.2.** A matrix $A = (a_{ij})$ is called Metzler of all the off-diagonal components are nonnegative, i.e. $a_{ij} \ge 0$ for $i \neq j$.

A matrix $A \in \mathcal{M}^{m \times n}$ is irreducible if it is not similar via a permutation to a block upper triangular matrix.

**Lemma 2.1.** Assume that a square matrix $A$ is Metzler and irreducible, then $e^A$ is strictly positive.

*Proof.* Since $A$ is Metzler, then there exists a constant $\delta > 0$ large enough such that $A + \delta I$ is a non-negative matrix with a positive element on the main diagonal. Moreover, $A$ is irreducible so $A + \delta I$ is also irreducible. Thus, $A + \delta I$ is primitive, that is, there exists an integer $n > 0$ such that $(A + \delta I)^n \gg 0$. Hence, we have $e^{A+\delta I} \gg 0$, and since $\delta I$ commutes with all matrices, one has $e^A = e^{A+\delta I}e^{-\delta I} \gg 0$. $\square$

We present in this section the so-called Kamke [37] or Chaplygin [38] lemma for a cooperative system (Lemma 2.2). Then, we apply this lemma to show the monotonicity of system (1) in Lemma 2.3.

**Lemma 2.2.** For any $n \in \mathbb{N}^*$, consider a smooth function $f: \mathbb{R}^n \to \mathbb{R}^n$, and a vector function $u(t)$ satisfying a differential equation

$$ \dot{\mathbf{u}} = \mathbf{f}(\mathbf{u}). $$

Moreover, we assume that the above system is cooperative, that is,

$$ \frac{\partial f_i}{\partial u_j}(t) \ge 0, \quad \text{for } i \neq j, t > 0. \qquad (3) $$

If a vector function $v(t)$ satisfies a differential inequality $\dot{v} \le \mathbf{f}(v)$ then, for initial data $v(0) \le \mathbf{u}(0)$, we have $v(t) \le \mathbf{u}(t)$ for all $t > 0$.

To apply this Lemma to system (1), we first define the following order in $\mathbb{R}^8$ as follows

**Definition 2.3.** For any vectors $\mathbf{u}, \mathbf{v} \in \mathbb{R}^8$, we define an order $\leq$ such that $\mathbf{u} \leq \mathbf{v}$ if and only if

$$ \begin{cases} u_i \le v_i & \text{for } i \in \{1, 2, 3, 5, 6, 7\}, \\ u_i \ge v_i & \text{for } i \in \{4, 8\}. \end{cases} $$

Moreover, we write $\mathbf{u} \prec \mathbf{v}$ if $\mathbf{u} \leq \mathbf{v}$ and $\mathbf{u} \neq \mathbf{v}$.

The following result shows the comparison principle of system (1)

**Lemma 2.3.** By denoting $\mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s) \in \mathbb{R}^8$, we can write system (1) as the form $\dot{\mathbf{u}} = \mathbf{f}(\mathbf{u})$ with $\mathbf{f}$ is $C^1$ in $\mathbb{R}^8$. In the subset $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of $\mathbb{R}_+^8$, system (1) is monotone in the sense that if a vector function $v(t)$ satisfies a differential inequality $\dot{v} \le \mathbf{f}(v)$ then, for initial data $v(0) \le \mathbf{u}(0)$, we have $v(t) \le \mathbf{u}(t)$ for all $t > 0$.

*Proof.* By changing the variable to $\tilde{\mathbf{u}} = (E_1, F_1, M_1, -M_1^s, E_2, F_2, M_2, -M_2^s)$, we can write system (1) as

$$ \tilde{\mathbf{u}} = \tilde{\mathbf{f}}(\tilde{\mathbf{u}}). $$

This system is cooperative since in $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of $\mathbb{R}_+^8$, we have

$$ \frac{\partial \tilde{f}_1}{\partial \tilde{u}_2} = b \left(1 - \frac{E_1}{K_1}\right) \ge 0, \quad \frac{\partial \tilde{f}_1}{\partial \tilde{u}_j} = 0 \quad \text{for any } j > 2, $$

$$ \frac{\partial \tilde{f}_2}{\partial \tilde{u}_1} = r \nu_E \frac{M_1}{M_1 + \gamma M_1^s} \ge 0, \quad \frac{\partial \tilde{f}_2}{\partial \tilde{u}_3} = r \nu_E \frac{\gamma M_1^s}{(M_1 + \gamma M_1^s)^2} \ge 0, $$

$$ \frac{\partial \tilde{f}_2}{\partial \tilde{u}_4} = r \nu_E E_1 \frac{\gamma M_1}{(M_1 + \gamma M_1^s)^2} \ge 0, \quad \frac{\partial \tilde{f}_2}{\partial \tilde{u}_6} = d_{21} > 0, \quad \frac{\partial \tilde{f}_2}{\partial \tilde{u}_j} = 0 \quad \text{for } j \in \{5, 7, 8\}. $$

$$ \frac{\partial \tilde{f}_3}{\partial \tilde{u}_1} = (1 - r) \nu_E \ge 0, \quad \frac{\partial \tilde{f}_3}{\partial \tilde{u}_7} = \beta d_{21} > 0, \quad \frac{\partial \tilde{f}_1}{\partial \tilde{u}_j} = 0 \quad \text{for } j \in \{2, 4, 5, 6, 8\}, $$
To apply this Lemma to system (1), we first define the following order in $\mathbb{R}^8$ as follows

**Definition 2.2.** For any vectors $\mathbf{u}, \mathbf{v} \in \mathbb{R}^8$, we define an order $\preceq$ such that $\mathbf{u} \preceq \mathbf{v}$ if and only if

$$ \begin{cases} u_i \le v_i & \text{for } i \in \{1, 2, 3, 5, 6, 7\}, \\ u_i \ge v_i & \text{for } i \in \{4, 8\}. \end{cases} $$

Moreover, we write $\mathbf{u} \prec \mathbf{v}$ if $\mathbf{u} \preceq \mathbf{v}$ and $\mathbf{u} \neq \mathbf{v}$.

The monotonicity of system (1) is shown in the following result

**Lemma 2.3.** By denoting $\mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s) \in \mathbb{R}^8$, we can write system (1) as the form $\dot{\mathbf{u}} = \mathbf{f}(\mathbf{u})$ with $\mathbf{f}$ is $\mathcal{C}^1$ in $\mathbb{R}^8$. In the invariant subset $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of $\mathbb{R}_+^8$, system (1) is monotone in the sense that if a vector function $\mathbf{v}(t)$ satisfies a differential inequality $\dot{\mathbf{v}} \preceq \mathbf{f}(\mathbf{v})$ then, for initial data $\mathbf{v}(0) \preceq \mathbf{u}(0)$, we have $\mathbf{v}(t) \preceq \mathbf{u}(t)$ for all $t > 0$.

*Proof.* By changing the variable to $\tilde{\mathbf{u}} = (E_1, F_1, M_1, -M_1^s, E_2, F_2, M_2, -M_2^s)$, we can write system (1) as

$$ \tilde{\mathbf{u}} = \tilde{\mathbf{f}}(\tilde{\mathbf{u}}). $$

This system is cooperative since in $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of $\mathbb{R}_+^8$, we have

$$ \frac{\partial \tilde{f}_1}{\partial \tilde{u}_2} = b \left( 1 - \frac{E_1}{K_1} \right) \ge 0, \quad \frac{\partial \tilde{f}_1}{\partial \tilde{u}_j} = 0 \text{ for any } j > 2, $$

$$ \frac{\partial \tilde{f}_2}{\partial \tilde{u}_1} = r v_E \frac{M_1}{M_1 + \gamma M_1^s} \ge 0, \quad \frac{\partial \tilde{f}_2}{\partial \tilde{u}_3} = r v_E \frac{\gamma M_1^s}{(M_1 + \gamma M_1^s)^2} \ge 0, $$

$$ \frac{\partial \tilde{f}_2}{\partial \tilde{u}_4} = r v_E E_1 \frac{\gamma M_1}{(M_1 + \gamma M_1^s)^2} \ge 0, \quad \frac{\partial \tilde{f}_2}{\partial \tilde{u}_6} = d_{21} > 0, \quad \frac{\partial \tilde{f}_2}{\partial \tilde{u}_j} = 0 \text{ for } j \in \{5, 7, 8\}. $$

$$ \frac{\partial \tilde{f}_3}{\partial \tilde{u}_1} = (1 - r) v_E \ge 0, \quad \frac{\partial \tilde{f}_3}{\partial \tilde{u}_7} = \beta d_{21} > 0, \quad \frac{\partial \tilde{f}_1}{\partial \tilde{u}_j} = 0 \text{ for } j \in \{2, 4, 5, 6, 8\}, $$

$$ \frac{\partial \tilde{f}_4}{\partial \tilde{u}_8} = \alpha d_{21} > 0, \quad \frac{\partial \tilde{f}_4}{\partial \tilde{u}_j} = 0 \text{ for } j \in \{1, 2, 3, 5, 6, 7\}. $$

Similarly for $\tilde{f}_i$ with $i > 4$, so $\tilde{\mathbf{f}}$ is cooperative.

For any vector function $\mathbf{v}$ such that $\mathbf{v} \preceq \mathbf{f}(\mathbf{v})$, by the same variable change, one has $\tilde{\mathbf{v}} \preceq \tilde{\mathbf{f}}(\tilde{\mathbf{v}})$. The initial data $\mathbf{v}(0) \preceq \mathbf{u}(0)$ implies that $\tilde{\mathbf{v}}(0) \preceq \tilde{\mathbf{u}}(0)$. Therefore, by applying Lemma 2.2, one has $\tilde{\mathbf{v}}(t) \preceq \tilde{\mathbf{u}}(t)$ for any $t > 0$ which is equivalent to $\mathbf{v}(t) \preceq \mathbf{u}(t)$. $\square$

In order to define the solution of (1), we make some assumptions for the release function $\Lambda(t)$

**Assumption 2.1.** Assume that function $\Lambda(t)$ satisfies

$$ \Lambda(t) = \Lambda_1(t) + \Lambda_2(t), \quad (4) $$

where $\Lambda_1 \in L_{\text{loc}}^1(0, +\infty)$, $\Lambda_1(t) \ge 0$ for almost every $t$, and $\Lambda_2$ is a sum of Dirac masses with positive weights. Assume moreover that there exists a time $T > 0$ such that the average value of $\Lambda$ over any $T$-time interval is finite, that is,

$$ C_{\Lambda} := \frac{1}{T} \sup_{t \ge 0} \int_{t}^{t+T} \Lambda(s) ds < +\infty. \quad (5) $$
$$ \frac{\partial \tilde{f}_4}{\partial \tilde{u}_8} = \alpha d_{21} > 0, \frac{\partial \tilde{f}_4}{\partial \tilde{u}_j} = 0 \text{ for } j \in \{1, 2, 3, 5, 6, 7\}. $$

Similarly for $\tilde{f}_i$ with $i > 4$, so $\tilde{f}$ is cooperative.

For any vector function $\mathbf{v}$ such that $\dot{\mathbf{v}} \preceq \mathbf{f}(\mathbf{v})$, by the same variable change, one has $\dot{\tilde{\mathbf{v}}} \le \tilde{\mathbf{f}}(\tilde{\mathbf{v}})$. The initial data $\mathbf{v}(0) \preceq \mathbf{u}(0)$ implies that $\tilde{\mathbf{v}}(0) \le \tilde{\mathbf{u}}(0)$. Therefore, by applying Lemma 2.2, one has $\tilde{\mathbf{v}}(t) \le \tilde{\mathbf{u}}(t)$ for any $t > 0$ which is equivalent to $\mathbf{v}(t) \preceq \mathbf{u}(t)$. $\square$

In order to define the solution of (1), we make some assumptions for the release function $\Lambda(t)$

**Assumption 2.1.** Assume that function $\Lambda(t)$ satisfies

$$ \Lambda(t) = \Lambda_1(t) + \Lambda_2(t), \quad (4) $$

where $\Lambda_1 \in L_{\text{loc}}^1(0, +\infty)$, $\Lambda_1(t) \ge 0$ for almost every $t$, and $\Lambda_2$ is a sum of Dirac masses with positive weights. More precisely, for $0 < t_1 < t_2 < \dots$, one has

$$ \Lambda_2(t) = \sum_{i=0}^{+\infty} w_k \delta_{t_k}(t) \quad \text{with } w_k > 0. \quad (5) $$

Assume moreover that there exists a time $T > 0$ such that the average value of $\Lambda$ over any $T$-time interval is finite, that is,

$$ C_{\Lambda} := \frac{1}{T} \sup_{t \ge 0} \int_{t}^{t+T} \Lambda(s) ds < +\infty. \quad (6) $$

Assumption 2.1 is natural since in practice, the total amount of the sterile males released in a finite time interval is finite. The term $\Lambda_2$ corresponds to impulsive releases.

The next result shows that any trajectory of system (1) resulting from any non-negative initial data is bounded.

**Lemma 2.4.** Let $\Lambda$ satisfy Assumption 2.1. For any non-negative initial data $(E_1^0, F_1^0, M_1^0, M_1^{s,0}, E_2^0, F_2^0, M_2^0, M_2^{s,0})$, there exists a unique solution $(E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s)$ of system (1) that is smooth on each interval $(t_k, t_{k+1})$ of the impulsive function $\Lambda_2$ defined in (5) with $k = 0, 1, \dots$, and it is non-negative. If $E_i^0 < K_i$ with $i = 1, 2$, then $E_i(t) \le K_i$ for any $t > 0$. Moreover, for all $t \ge 0$, we have the uniform bounds

$$ F_1 + F_2 \le \max \{F_1^0 + F_2^0, C_F\}, \quad M_1 + M_2 \le \max \{M_1^0 + M_2^0, C_M\}, $$

where

$$ C_F := \frac{r\nu_E(K_1 + K_2)}{\mu_F}, \quad C_M := \frac{(1-r)\nu_E(K_1 + K_2)}{\mu_M}, $$

and

$$ M_1^s(t) + M_2^s(t) \le \max \left\{ M_1^{s,0} + M_2^{s,0}, \frac{TC_\Lambda}{1 - e^{-\mu_s T}} \right\} + TC_\Lambda, $$

with $T$ and $C_\Lambda$ defined in Assumption 2.1. One also has

$$ \limsup_{t \to +\infty} (F_1 + F_2)(t) \le C_F, \quad \limsup_{t \to +\infty} (M_1 + M_2)(t) \le C_M, \quad , $$

and

$$ \limsup_{t \to +\infty} (M_1^s + M_2^s)(t) \le \frac{TC_\Lambda}{1 - e^{-\mu_s T}} + TC_\Lambda =: C_{M^s}. $$

**Remark 2.1.** In the case $\Lambda \in L^\infty(0, +\infty)$, one can let $T$ tend to zero and obtain that $C_\Lambda = \sup_{t>0} \Lambda(t)$ and $\limsup_{t \to +\infty} (M_1^s + M_2^s) \le \frac{C_\Lambda}{\mu_s}$. The condition of $\Lambda$ that we made in Assumption 2.1 is weaker than the $L^\infty$ assumption since we also include impulsive releases, represented by a sum of Dirac masses.

**Proof of Lemma 2.4.** By denoting $\mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s) \in \mathbb{R}^8$, we can write system (1) as the form $\dot{\mathbf{u}} = \mathbf{f}(\mathbf{u})$ with $\mathbf{f} = (f_1, f_2, \dots, f_s)$ is $C^1$ in $\mathbb{R}^8$. We easily observe that in the subset $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of $\mathbb{R}_+^8$, one has $f_i(u_i = 0) \ge 0$ for $i = 1, \dots, 8$, and $f_1(E_1 = K_1) \le 0$, $f_5(E_2 = K_2) \le 0$. Therefore, the vector field $\mathbf{f}$ on the boundary of $\mathbb{R}_+^8 \cap \{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ is inward or tangential, so this set is positively invariant.

From equations (1b) and (1f), since for $i = 1, 2$, we have $\frac{M_i}{M_i + \gamma M_i^s} \le 1$, and

$$ \dot{F}_1 + \dot{F}_2 \le r\nu_E(E_1 + E_2) - \mu_F(F_1 + F_2). $$
Assumption 2.1 is natural since in practice, the total amount of the sterile males released in a finite time interval is finite. The term $\Lambda_2$ corresponds to impulsive releases.

The next result shows that any trajectory of system (1) resulting from any non-negative initial data is bounded.

**Lemma 2.4.** Let $\Lambda$ satisfy Assumption 2.1. For any non-negative initial data $(E_1^0, F_1^0, M_1^0, M_1^{s,0}, E_2, F_2, M_2, M_2^s)$ of system (1), and it is non-negative. If $E_i^0 < K_i$ with $i = 1, 2$, then $E_i(t) \le K_i$ for any $t > 0$. Moreover, for all $t \ge 0$, we have the uniform bounds

$$F_1 + F_2 \le \max \{F_1^0 + F_2^0, C_F\}, \quad M_1 + M_2 \le \max \{M_1^0 + M_2^0, C_M\},$$

where

$$C_F := \frac{r v_E (K_1 + K_2)}{\mu_F}, \quad C_M := \frac{(1-r) v_E (K_1 + K_2)}{\mu_M},$$

and

$$M_1^s(t) + M_2^s(t) \le \max \left\{ M_1^{s,0} + M_2^{s,0}, \frac{TC_\Lambda}{1 - e^{-\mu_s T}} \right\} + TC_\Lambda,$$

with $T$ and $C_\Lambda$ defined in Assumption 2.1. One also has

$$\limsup_{t \to +\infty} (F_1 + F_2)(t) \le C_F, \quad \limsup_{t \to +\infty} (M_1 + M_2)(t) \le C_M, \quad ,$$

and

$$\limsup_{t \to +\infty} (M_1^s + M_2^s)(t) \le \frac{TC_\Lambda}{1 - e^{-\mu_s T}} + TC_\Lambda := C_M^s.$$

**Remark 2.1.** In the case $\Lambda \in L^\infty(0, +\infty)$, one can let $T$ tend to zero and obtain that $C_\Lambda = \sup_{t>0} \Lambda(t)$ and $\limsup_{t \to +\infty} (M_1^s + M_2^s) \le \frac{C_\Lambda}{\mu_s}$. The condition of $\Lambda$ that we made in Assumption 2.1 is weaker than the $L^\infty$ assumption since we also include impulsive releases, represented by the Dirac masses.

**Proof of Lemma 2.4.** By applying the Lemma 2.3, one deduces that system (1) preserves the positivity. For $i = 1, 2$, we have $E_i(t = 0) = E_i^0 < K_i$ and assume that there exists a value $t_0 < \infty$ such that

$$t_0 = \inf\{t > 0 : E_i(t) = K_i\}$$

then $\dot{E}_i(t_0) > 0$ but from (1a) and (1e), one has $\dot{E}_i(t_0) = -(v_E + \mu_E)K_i < 0$ (contradictory). Then we deduce that $E_i(t) \le K_i$ for any $t > 0$.

From equations (1b) and (1f), since for $i = 1, 2$, $\frac{M_i}{M_i + \gamma M_i^s} \le 1$ one has

$$\dot{F}_1 + \dot{F}_2 \le r v_E (E_1 + E_2) - \mu_F (F_1 + F_2).$$

Since $E_1, E_2$ are bounded then we deduce that

$$(F_1 + F_2)(t) \le (F_1^0 + F_2^0)e^{-\mu_F t} + \frac{r v_E (K_1 + K_2)}{\mu_F} (1 - e^{-\mu_F t}) \le \max \left\{ F_1^0 + F_2^0, \frac{r v_E (K_1 + K_2)}{\mu_F} \right\},$$

for any $t \ge 0$. For $i = 1, 2$, one has $F_i \ge 0$, thus $F_i(t) \le \max \{F_1^0 + F_2^0, C_F\}$ for any $t \ge 0$. Let $t$ goes to infinity we get $\limsup_{t \to +\infty} (F_1 + F_2)(t) \le C_F$. One obtains similarly the inequalities for $M_1, M_2$.

For $M_1^s$ and $M_2^s$, by denoting $X_s(t) = M_1^s(t) + M_2^s(t)$, then from equations (1d) and (1h), one has

$$\dot{X}_s(t) = -\mu_s X_s(t) + \Lambda(t).$$
Since $E_1, E_2$ are bounded then we deduce that

$$ (F_1 + F_2)(t) \le (F_1^0 + F_2^0)e^{-\mu_F t} + \frac{r\nu_E(K_1 + K_2)}{\mu_F}(1 - e^{-\mu_F t}) \le \max \left\{ F_1^0 + F_2^0, \frac{r\nu_E(K_1 + K_2)}{\mu_F} \right\}, $$

for any $t \ge 0$. For $i = 1, 2$, one has $F_i \ge 0$, thus $F_i(t) \le \max\{F_1^0 + F_2^0, C_F\}$ for any $t \ge 0$. Let $t$ go to infinity we get $\limsup(F_1 + F_2)(t) \le C_F$. One obtains similarly the inequalities for $M_1, M_2$.

For $M_1^s$ and $M_2^s$, by denoting $X_s(t) = M_1^s(t) + M_2^s(t)$, then from equations (1d) and (1h), one has

$$ \dot{X}_s(t) = -\mu_s X_s(t) + \Lambda(t). $$

For any integer $k$, by integrating both sides of this equality in $((k-1)T, kT)$ with $T$ defined in Assumption 2.1, one gets

$$ X_s(kT) = e^{-\mu_s T} X_s((k-1)T) + \int_{(k-1)T}^{kT} e^{-\mu_s(t-(k-1)T)} \Lambda(t) dt. $$

Since $e^{-\mu_s(t-(k-1)T)} < 1$ for any $t \in ((k-1)T, kT)$ and by Assumption 2.1, we deduce that

$$ X_s(kT) \le e^{-\mu_s T} X_s((k-1)T) + TC_\Lambda, $$

with $C_\Lambda$ defined in (6). Using the iteration with respect to $k$, we deduce that

$$ X_s(kT) \le e^{-\mu_s kT} X_s^0 + TC_\Lambda (1 + e^{-\mu_s T} + \dots + e^{-\mu_s (k-1)T}) = e^{-\mu_s kT} X_s^0 + TC_\Lambda \frac{1 - e^{-\mu_s kT}}{1 - e^{-\mu_s T}}. $$

Now for any time $t > 0$, there exists an integer $k$ such that $t \in [kT, (k+1)T)$. Then, we obtain that

$$ \begin{aligned} X_s(t) &= e^{-\mu_s(t-kT)} X_s(kT) + \int_{kT}^{t} e^{-\mu_s(t-s)} \Lambda(s) ds \\ &\le e^{-\mu_s t} X_s^0 + TC_\Lambda \frac{e^{-\mu_s(t-kT)} - e^{-\mu_s t}}{1 - e^{-\mu_s T}} + TC_\Lambda. \\ &\le e^{-\mu_s t} X_s^0 + \frac{TC_\Lambda}{1 - e^{-\mu_s T}} (1 - e^{-\mu_s t}) + TC_\Lambda \end{aligned} $$

since $e^{-\mu_s(t-kT)} < 1$. The inequality of $M_1^s + M_2^s$ follows. $\square$

### 3 Mosquito dynamics without sterile males

First, we describe the dynamics of wild mosquitoes in the two areas by considering the following system which is re-obtained from system (1) in the absence of sterile males

$$ \dot{E}_1 = bF_1 \left(1 - \frac{E_1}{K_1}\right) - (\nu_E + \mu_E)E_1, \quad (7a) $$

$$ \dot{F}_1 = r\nu_E E_1 - \mu_F F_1 - d_{12} F_1 + d_{21} F_2, \quad (7b) $$

$$ \dot{M}_1 = (1-r)\nu_E E_1 - \mu_M M_1 - \beta d_{12} M_1 + \beta d_{21} M_2, \quad (7c) $$

$$ \dot{E}_2 = bF_2 \left(1 - \frac{E_2}{K_2}\right) - (\nu_E + \mu_E)E_2, \quad (7d) $$

$$ \dot{F}_2 = r\nu_E E_2 - \mu_F F_2 - d_{21} F_2 + d_{12} F_1, \quad (7e) $$

$$ \dot{M}_2 = (1-r)\nu_E E_2 - \mu_M M_2 - \beta d_{21} M_2 + \beta d_{12} M_1, \quad (7f) $$

It is clear that the subset $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of the positive cone of $\mathbb{R}^6$ is positively invariant over time. The following result shows the nature of the equilibrium points of system (7).

**Theorem 3.1.** For $\mathcal{N} \le 1$, zero is the unique equilibrium of system (7), and all trajectories of (7) resulting from non-negative initial data converge to zero as time evolves.

For $\mathcal{N} > 1$, system (7) has two equilibrium points: zero and $\mathbf{u}^+ = (E_1^+, F_1^+, M_1^+, E_2^+, F_2^+, M_2^+)$ strictly positive. Moreover, the zero equilibrium is unstable. All trajectories of (7) resulting from any positive initial data $(E_1^0, F_1^0, M_1^0, E_2^0, F_2^0, M_2^0)$ such that $(E_1^0, F_1^0, E_2^0, F_2^0) > 0$ converge to $\mathbf{u}^+$ when $t \to +\infty$.
For any integer $k$, by integrating both sides of this equality in $((k-1)T, kT)$ with $T$ defined in Assumption 2.1, one gets

$$X_s(kT) = e^{-\mu_s T} X_s((k-1)T) + \int_{(k-1)T}^{kT} e^{-\mu_s (t-(k-1)T)} \Lambda(t) dt.$$

Since $e^{-\mu_s(t-(k-1)T)} < 1$ for any $t \in ((k-1)T, kT)$ and by Assumption 2.1, we deduce that

$$X_s(kT) \le e^{-\mu_s T} X_s((k-1)T) + TC_\Lambda,$$

with $C_\Lambda$ defined in (5). Using the iteration with respect to $k$, we deduce that

$$X_s(kT) \le e^{-\mu_s kT} X_s^0 + TC_\Lambda (1 + e^{-\mu_s T} + \dots + e^{-\mu_s (k-1)T}) = e^{-\mu_s kT} X_s^0 + TC_\Lambda \frac{1 - e^{-\mu_s kT}}{1 - e^{-\mu_s T}}.$$

Now for any time $t > 0$, there exists an integer $k$ such that $t \in [kT, (k+1)T)$. Then, we obtain that

$$\begin{aligned} X_s(t) &= e^{-\mu_s(t-kT)}X_s(kT) + \int_{kT}^{t} e^{-\mu_s(t-s)}\Lambda(s)ds \\ &\le e^{-\mu_s t}X_s^0 + TC_\Lambda \frac{e^{-\mu_s(t-kT)} - e^{-\mu_s t}}{1 - e^{-\mu_s T}} + TC_\Lambda. \\ &\le e^{-\mu_s t}X_s^0 + \frac{TC_\Lambda}{1 - e^{-\mu_s T}}(1 - e^{-\mu_s t}) + TC_\Lambda \end{aligned}$$

since $e^{-\mu_s(t-kT)} < 1$. The inequality of $M_1^s + M_2^s$ follows. $\square$

### 3 Mosquito dynamics without sterile males

First, we describe the dynamics of wild mosquitoes in the two areas by considering the following system which is re-obtained from system (1) in the absence of sterile males

$$\dot{E}_1 = bF_1 \left(1 - \frac{E_1}{K_1}\right) - (v_E + \mu_E)E_1, \quad (6a)$$

$$\dot{F}_1 = r v_E E_1 - \mu_F F_1 - d_{12} F_1 + d_{21} F_2, \quad (6b)$$

$$\dot{M}_1 = (1-r)v_E E_1 - \mu_M M_1 - \beta d_{12} M_1 + \beta d_{21} M_2, \quad (6c)$$

$$\dot{E}_2 = bF_2 \left(1 - \frac{E_2}{K_2}\right) - (v_E + \mu_E)E_2, \quad (6d)$$

$$\dot{F}_2 = r v_E E_2 - \mu_F F_2 - d_{21} F_2 + d_{12} F_1, \quad (6e)$$

$$\dot{M}_2 = (1-r)v_E E_2 - \mu_M M_2 - \beta d_{21} M_2 + \beta d_{12} M_1, \quad (6f)$$

It is clear that the subset $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of the positive cone of $\mathbb{R}^6$ is positively invariant over time. The following result shows the nature of the equilibrium points of system (6).

**Theorem 3.1.** For $\mathcal{N} < 1$, zero is the unique equilibrium of system (6), and all trajectories of (6) resulting from non-negative initial data converge to zero as time evolves.

For $\mathcal{N} > 1$, system (6) has two equilibrium points: zero and $\mathbf{u}^+ = (E_1^+, F_1^+, M_1^+, E_2^+, F_2^+, M_2^+)$ strictly positive. Moreover, the zero equilibrium is unstable. All trajectories of (6) resulting from any positive initial data $(E_1^0, F_1^0, M_1^0, E_2^0, F_2^0, M_2^0)$ such that $(E_1^0, F_1^0, E_2^0, F_2^0) > 0$ converge to $\mathbf{u}^+$ when $t \to +\infty$.
Theorem 3.1 shows that the constant $\mathcal{N}$ defined in (2) is the basic offspring number of the whole two-patch system (1). When $\mathcal{N} > 1$, the populations in both areas remain persistent for any diffusion rates as time evolves. In the rest of the paper, we only consider the case $\mathcal{N} > 1$.

To prove this theorem, we first consider the sub-system of $E_1, E_2, F_1, F_2$. From equations (7b) and (7e), the positive equilibrium satisfies

$$r\nu_E \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix} = \begin{pmatrix} \mu_F + d_{12} & -d_{21} \\ -d_{12} & \mu_F + d_{21} \end{pmatrix} \begin{pmatrix} F_1^+ \\ F_2^+ \end{pmatrix},$$

then

$$\begin{pmatrix} F_1^+ \\ F_2^+ \end{pmatrix} = \frac{r\nu_E}{\mu_F(\mu_F + d_{12} + d_{21})} \begin{pmatrix} \mu_F + d_{21} & d_{21} \\ d_{12} & \mu_F + d_{12} \end{pmatrix} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}. \quad (8)$$

On the other hand, from equation (7a) and (7d), we also have

$$F_1^+ = \frac{\nu_E + \mu_E}{b} \frac{E_1^+}{1 - \frac{E_1^+}{K_1}}, \quad F_2^+ = \frac{\nu_E + \mu_E}{b} \frac{E_2^+}{1 - \frac{E_2^+}{K_2}}. \quad (9)$$

From (8) and (9), we deduce that

$$E_2^+ = \frac{\mu_F + d_{12} + d_{21}}{d_{21}} \frac{1}{\mathcal{N}} \frac{E_1^+}{1 - \frac{E_1^+}{K_1}} - \frac{\mu_F + d_{21}}{d_{21}} E_1^+ =: f_{21}(E_1^+), \quad (10)$$

$$E_1^+ = \frac{\mu_F + d_{12} + d_{21}}{d_{12}} \frac{1}{\mathcal{N}} \frac{E_2^+}{1 - \frac{E_2^+}{K_2}} - \frac{\mu_F + d_{12}}{d_{12}} E_2^+ =: f_{12}(E_2^+). \quad (11)$$

The following lemma provides information for these functions.

**Lemma 3.2.** For $\{i, j\} = \{1, 2\}$, function $f_{ij}(x)$ is defined and convex on $(0, K_j) \subset \mathbb{R}$.
If $\mathcal{N} \le \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$, then $f_{ij}$ has no positive root and it is increasing on $(0, K_j)$.
Otherwise, it has a unique positive root

$$K_j^+ := K_j \left( 1 - \frac{\mu_F + d_{ji} + d_{ij}}{\mu_F + d_{ij}} \frac{1}{\mathcal{N}} \right) < K_j \quad (12)$$

Moreover, $f_{ij} < 0$ on $(0, K_j^+)$, $f_{ij} > 0$ and increasing on $(K_j^+, K_j)$.

**Proof of Lemma 3.2.** We recall function $f_{ij}(x) := \frac{\mu_F + d_{ij} + d_{ji}}{d_{ij}} \frac{1}{\mathcal{N}} \frac{x}{1 - \frac{x}{K_j}} - \frac{\mu_F + d_{ij}}{d_{ij}} x$. One has $f_{ij} = 0$ if and only if

$$\frac{\mu_F + d_{ij}}{K_j} x^2 + \left( \frac{\mu_F + d_{ij} + d_{ji}}{\mathcal{N}} - \mu_F - d_{ij} \right) x = 0.$$

We deduce that $f_{ij} = 0$ at 0 and $K_j^+$ as in (12), and $K_j^+ > 0$ if and only if $\mathcal{N} > \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$. Moreover, $f_{ij} < 0$ on $(0, K_j^+)$, $f_{ij} > 0$ and is increasing on $(K_j^+, K_j)$. It is defined and convex on $(0, K_j) \subset \mathbb{R}$ since

$$f_{ij}''(x) = 2 \frac{\mu_F + d_{ij} + d_{ji}}{d_{ij}} \frac{1}{K_j \mathcal{N}} \left(1 - \frac{x}{K_j}\right)^3 > 0,$$

for any $x \in (0, K_j)$. We also have $f_{ij}(0) = 0$, $\lim_{x \to K_j} f_{ij}(x) = +\infty$. $\square$

**Proof of Theorem 3.1. Existence and uniqueness of the positive equilibrium.** System (7) has a positive equilibrium iff system (10)-(11) has a solution $(E_1^+, E_2^+)$ in $(0, K_1) \times (0, K_2)$.

**Case 1:** $0 < \mathcal{N} \le \min_{\substack{i,j \in \{1,2\} \\ i \ne j}} \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$.

According to Lemma 3.2, we have $f_{12} : [0, K_2) \to [0, +\infty)$ is positive and increasing, so this function is invertible (see Figure 1). We denote $f_{12}^{-1} : [0, K_1] \to [0, K_2)$ the restriction of the inverse function of $f_{12}$ on $[0, K_1]$, then

$$E_2^+ = f_{21}(E_1^+) = f_{12}^{-1}(E_1^+).$$
Theorem 3.1 shows that the constant $\mathcal{N}$ defined in (2) is the basic offspring number of the whole two-patch system (1). When $\mathcal{N} > 1$, the populations in both areas remain persistent for any diffusion rates as time evolves. In the rest of the paper, we only consider the case $\mathcal{N} > 1$.

To prove this theorem, we first consider the sub-system of $E_1, E_2, F_1, F_2$. From equations (6b) and (6e), the positive equilibrium satisfies

$$r v_E \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix} = \begin{pmatrix} \mu_F + d_{12} & -d_{21} \\ -d_{12} & \mu_F + d_{21} \end{pmatrix} \begin{pmatrix} F_1^+ \\ F_2^+ \end{pmatrix},$$

then

$$\begin{pmatrix} F_1^+ \\ F_2^+ \end{pmatrix} = \frac{r v_E}{\mu_F(\mu_F + d_{12} + d_{21})} \begin{pmatrix} \mu_F + d_{21} & d_{21} \\ d_{12} & \mu_F + d_{12} \end{pmatrix} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}. \quad (7)$$

On the other hand, from equation (6a) and (6d), we also have

$$F_1^+ = \frac{v_E + \mu_E}{b} \frac{E_1^+}{1 - \frac{E_1^+}{K_1}}, \quad F_2^+ = \frac{v_E + \mu_E}{b} \frac{E_2^+}{1 - \frac{E_2^+}{K_2}}. \quad (8)$$

From (7) and (8), we deduce that

$$E_2^+ = \frac{\mu_F + d_{12} + d_{21}}{d_{21}} \frac{1}{\mathcal{N}} \frac{E_1^+}{1 - \frac{E_1^+}{K_1}} - \frac{\mu_F + d_{21}}{d_{21}} E_1^+ := f_{21}(E_1^+), \quad (9)$$

$$E_1^+ = \frac{\mu_F + d_{12} + d_{21}}{d_{12}} \frac{1}{\mathcal{N}} \frac{E_2^+}{1 - \frac{E_2^+}{K_2}} - \frac{\mu_F + d_{12}}{d_{12}} E_2^+ := f_{12}(E_2^+). \quad (10)$$

The following lemma provides information for these functions.

**Lemma 3.2.** For $\{i, j\} = \{1, 2\}$, function $f_{ij}(x)$ is defined and convex on $(0, K_j) \subset \mathbb{R}$.

If $\mathcal{N} \le \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$, then $f_{ij}$ has no positive root and it is increasing on $(0, K_j)$.

Otherwise, it has a unique positive root

$$K_j^+ := K_j \left( 1 - \frac{\mu_F + d_{ji} + d_{ij}}{\mu_F + d_{ij}} \frac{1}{\mathcal{N}} \right) < K_j \quad (11)$$

Moreover, $f_{ij} < 0$ on $(0, K_j^+)$, $f_{ij} > 0$ and increasing on $(K_j^+, K_j)$.

**Proof of Lemma 3.2.** We recall function $f_{ij}(x) := \frac{\mu_F + d_{ij} + d_{ji}}{d_{ij}} \frac{x}{\mathcal{N} \left(1 - \frac{x}{K_j}\right)} - \frac{\mu_F + d_{ij}}{d_{ij}} x$. One has $f_{ij} = 0$ if and only if

$$\frac{\mu_F + d_{ij}}{K_j} x^2 + \left( \frac{\mu_F + d_{ij} + d_{ji}}{d_{ij}} - \mu_F - d_{ij} \right) x = 0.$$

We deduce that $f_{ij} = 0$ at 0 and $K_j^+$ as in (11), and $K_j^+ > 0$ if and only if $\mathcal{N} \le \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$. Moreover, $f_{ij} < 0$ on $(0, K_j^+)$, $f_{ij} > 0$ and is increasing on $(K_j^+, K_j)$. It is defined and convex on $(0, K_j) \subset \mathbb{R}$ since

$$f_{ij}''(x) = 2 \frac{\mu_F + d_{ij} + d_{ji}}{d_{ij}} \frac{1}{K_j \mathcal{N}} \frac{1}{\left(1 - \frac{x}{K_j}\right)^3} > 0,$$

for any $x \in (0, K_j)$. We also have $f_{ij}(0) = 0$, $\lim_{x \to K_j} f_{ij}(x) = +\infty$. $\square$
<figure><img src="image_2.png" /><figcaption>Figure 1: Relations $E_2^+ = f_{21}(E_1^+)$ (in green) and $E_1^+ = f_{12}(E_2^+)$ (in red) presented in the $E_1 - E_2$ plane when $\mathcal{N} > 1$. In case (a), function $f_{12}$ and $f_{21}$ are increasing and have a unique root zero, so $f_{12}$ is invertible on $[0, K_2)$ and the inverse function $f_{12}^{-1} : [0, K_1] \to [0, K_2)$ is plotted in the solid red line. In case (b), function $f_{12}$ has another positive root $K_2^+$ and is invertible on $[K_2^+, K_2)$. The inverse function $f_{12}^{-1} : [0, K_1] \to [K_2^+, K_2)$ is plotted in the solid red line.</figcaption></figure>

Thus, $E_1^+$ is a positive root of function $f_{21} - f_{12}^{-1}$. For any $x \in (0, K_1)$, one has

$$ (f_{21} - f_{12}^{-1})'(x) = f'_{21}(x) - \frac{1}{f'_{12}(f_{12}^{-1}(x))}, $$

then

$$ (f_{21} - f_{12}^{-1})''(x) = f''_{21}(x) + \frac{f''_{12}(f_{12}^{-1}(x))}{(f'_{12}(f_{12}^{-1}(x)))^3} > 0 $$

since $f_{ij}$ is convex on $(0, K_j)$. Hence, $f_{21} - f_{12}^{-1}$ is convex on $(0, K_1)$. Moreover, we have $(f_{21} - f_{12}^{-1})(0) = 0$, and $\lim_{x \to K_1} (f_{21} - f_{12}^{-1})(x) = +\infty$. Therefore, this function has a unique positive root if and only if the derivative at zero is negative. We have

$$ (f_{21} - f_{12}^{-1})'(0) = \frac{\frac{1}{\mathcal{N}}(\mu_F + d_{12} + d_{21}) - \mu_F - d_{21}}{d_{21}} - \frac{d_{12}}{\frac{1}{\mathcal{N}}(\mu_F + d_{12} + d_{21}) - \mu_F - d_{12}}. $$

Since $\mathcal{N} \le \min_{\substack{i,j \in \{1,2\} \\ i \ne j}} \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}} < \frac{\mu_F + d_{12} + d_{21}}{\mu_F}$, then $(f_{21} - f_{12}^{-1})'(0) < 0$ if and only if $\mathcal{N} > 1$.

**Case 2:** $\mathcal{N} > \min_{\substack{i,j \in \{1,2\} \\ i \ne j}} \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$.

Without loss of generality, we assume that $d_{12} > d_{21}$, then $\frac{\mu_F + d_{12} + d_{21}}{\mu_F + d_{12}} = \min_{\substack{i,j \in \{1,2\} \\ i \ne j}} \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$.

If $\mathcal{N} > \frac{\mu_F + d_{12} + d_{21}}{\mu_F + d_{12}} > 1$, again according to Lemma 3.2, function $f_{12}$ has a unique positive root $K_2^+$ and is invertible on $[K_2^+, K_2)$ (see Figure 1(b)). We denote again $f_{12}^{-1} : [0, K_1] \to [K_2^+, K_2)$ the restriction of the inverse function of $f_{12}$ on $[0, K_1]$, then we also have $f_{21} - f_{12}^{-1}$ convex on $(0, K_1)$, and $(f_{21} - f_{12}^{-1})(0) = -f_{12}^{-1}(0) = -K_2^+ < 0$, $\lim_{x \to K_1} (f_{21} - f_{12}^{-1})(x) = +\infty > 0$. We can deduce that $f_{21} - f_{12}^{-1}$ has a unique positive root on $(0, K_1)$.

**Instability of the zero equilibrium.** At the origin $0_6 = (0, 0, 0, 0, 0, 0)$ of $\mathbb{R}^6$, the Jacobian matrix of system (7) is

$$ J(0) = \begin{pmatrix} -\nu_E - \mu_E & b & 0 & 0 & 0 & 0 \\ r\nu_E & -\mu_F - d_{12} & 0 & 0 & d_{21} & 0 \\ (1-r)\nu_E & 0 & -\mu_M - \beta d_{12} & 0 & 0 & \beta d_{21} \\ 0 & 0 & 0 & -\nu_E - \mu_E & b & 0 \\ 0 & d_{12} & 0 & r\nu_E & -\mu_F - d_{21} & 0 \\ 0 & 0 & \beta d_{12} & (1-r)\nu_E & 0 & -\mu_F - \beta d_{21} \end{pmatrix}, $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Behaviors of $f_{21}$ and $f_{12}^{-1}$</figcaption></figure>

*Proof of Theorem 3.1. Existence and uniqueness of the positive equilibrium.* System (6) has a positive equilibrium iff system (9)-(10) has a solution $(E_1^+, E_2^+)$ in $(0, K_1) \times (0, K_2)$.

First, we study the case where $0 < \mathcal{N} \le \min_{\substack{i,j \in \{1,2\} \\ i \ne j}} \frac{\mu_F + d_{ij} + d_{ji}}{\mu_F + d_{ij}}$, then according to Lemma 3.2, we

have $f_{12} : [0, K_2) \to [0, +\infty)$ is positive and increasing, so this function is invertible (see Figure 1). We denote $f_{12}^{-1} : [0, K_1] \to [0, K_2)$ the restriction of the invert function of $f_{12}$ on $[0, K_1]$, then

$$E_2^+ = f_{21}(E_1^+) = f_{12}^{-1}(E_1^+).$$

Thus, $E_1^+$ is a positive root of function $f_{21} - f_{12}^{-1}$. For any $x \in (0, K_1)$, one has

$$(f_{21} - f_{12}^{-1})'(x) = f'_{21}(x) - \frac{1}{f'_{12}(f_{12}^{-1}(x))},$$

then

$$(f_{21} - f_{12}^{-1})''(x) = f''_{21}(x) + \frac{f''_{12}(f_{12}^{-1}(x))}{(f'_{12}(f_{12}^{-1}(x)))^3} > 0$$

since $f_{ij}$ is convex on $(0, K_j)$. Hence, $f_{21} - f_{12}^{-1}$ is convex on $(0, K_1)$. Moreover, we have $(f_{21} - f_{12}^{-1})(0) = 0$, and $\lim_{x \to K_1} (f_{21} - f_{12}^{-1})(x) = +\infty$. Therefore, this function has a unique positive root if and only if the derivative at zero is negative. We have

$$(f_{21} - f_{12}^{-1})'(0) = \frac{\frac{1}{\mathcal{N}}(\mu_F + d_{12} + d_{21}) - \mu_F - d_{21}}{d_{21}} - \frac{d_{12}}{\frac{1}{\mathcal{N}}(\mu_F + d_{12} + d_{21}) - \mu_F - d_{12}}.$$

Then, $(f_{21} - f_{12}^{-1})'(0) < 0$ if and only if $1 < \mathcal{N} < \frac{\mu_F + d_{12} + d_{21}}{\mu_F}$.

Now, without loss of generality, we assume that $d_{12} > d_{21}$, then $\frac{\mu_F + d_{12} + d_{21}}{\mu_F + d_{12}} < \frac{\mu_F + d_{12} + d_{21}}{\mu_F + d_{21}}$.

If $\mathcal{N} > \frac{\mu_F + d_{12} + d_{21}}{\mu_F + d_{12}} > 1$, again according to Lemma 3.2, function $f_{12}$ has a unique positive root $K_2^+$ and is invertible on $[K_2^+, K_2]$ (see Figure 1). We denote again $f_{12}^{-1} : [0, K_1] \to [K_2^+, K_2)$ the restriction of
with the characteristic polynomial

$$
\det(J(0) - \lambda I) = [(\lambda + \mu_F)(\lambda + \mu_M) + \beta(d_{12} + d_{21} + d_{12}\mu_F + d_{21}\mu_M)][(\lambda + \nu_E + \mu_E)(\lambda + \mu_F) - \mathrm{br}\nu_E] \\
\times [(\lambda + \nu_E + \mu_E)(\lambda + \mu_F + d_{12} + d_{21}) - \mathrm{br}\nu_E].
$$

Since $\mathcal{N} > 1$, we have $\mu_F(\nu_E + \mu_E) - \mathrm{br}\nu_E < 0$. Thus, we can deduce that the factor $(\lambda + \nu_E + \mu_E)(\lambda + \mu_F) - \mathrm{br}\nu_E = \lambda^2 + \lambda(\nu_E + \mu_E + \mu_F) + \mu_F(\nu_E + \mu_E) - \mathrm{br}\nu_E$ has one positive root $\lambda > 0$. Hence, the Jacobian at zero has at least one positive eigenvalue so the zero equilibrium is unstable.

**Stability of the positive equilibrium.** First, we can see that the system (7a)-(7b), (7d)-(7e) of ($E_1, F_1, E_2, F_2$) does not depend on $M_1, M_2$, and it is cooperative and irreducible. By applying Theorem 1.1 in Chapter 4 of [39], one deduces that this system is strongly monotone. When $\mathcal{N} > 1$, this system admits exactly two equilibria: $\mathbf{0}_4 = (0,0,0,0)$, and $(E_1^+, F_1^+, E_2^+, F_2^+)$. But the zero equilibrium is unstable, so by Theorem 2.2 in Chapter 2 of [39], if the initial data satisfies that $\mathbf{0}_4 < (E_1^0, F_1^0, E_2^0, F_2^0) \le (E_1^+, F_1^+, E_2^+, F_2^+)$, the solution $(E_1, F_1, E_2, F_2)$ converges to $(E_1^+, F_1^+, E_2^+, F_2^+)$ when $t \to +\infty$.

Now if the initial data satisfies that $(E_1^0, F_1^0, E_2^0, F_2^0) > (E_1^+, F_1^+, E_2^+, F_2^+)$, then there exists a constant $\lambda > 1$ large enough such that $\lambda(E_1^+, F_1^+, E_2^+, F_2^+) \ge (E_1^0, F_1^0, E_2^0, F_2^0)$. Since $1 - \frac{\lambda E_i^+}{K_i} < 1 - \frac{E_i^+}{K_i}$, one has

$$
b\lambda F_i^+ \left( 1 - \frac{\lambda E_i^+}{K_i} \right) - (\nu_E + \mu_E)\lambda E_i^+ < \lambda \left[ bF_i^+ \left( 1 - \frac{E_i^+}{K_i} \right) - (\nu_E + \mu_E)E_i^+ \right] = 0,
$$

and the right-hand side of system (7a)-(7b), (7d)-(7e) at $\lambda(E_1^+, F_1^+, E_2^+, F_2^+)$ is non positive. Thus, the trajectory resulting from the initial data $\lambda(E_1^+, F_1^+, E_2^+, F_2^+)$ is non-increasing, and therefore converges to $(E_1^+, F_1^+, E_2^+, F_2^+)$. By applying the Lemma 2.2 to system (7a)-(7b), (7d)-(7e), we deduce that the trajectory resulting from the initial data $(E_1^0, F_1^0, E_2^0, F_2^0)$ lies between $(E_1^+, F_1^+, E_2^+, F_2^+)$ and the trajectories resulting from $\lambda(E_1^+, F_1^+, E_2^+, F_2^+)$. Hence, it also converges to $(E_1^+, F_1^+, E_2^+, F_2^+)$ when time $t$ goes to infinity.

Moreover, since the trajectories issued from the initial data above and below $(E_1^+, F_1^+, E_2^+, F_2^+)$ all converge to the same limit, then by the comparison principle, we deduce that the trajectory resulting from any positive initial data with values between these initial values converges to this equilibrium.

Secondly, if we denote matrix $A = \begin{pmatrix} -\mu_M - \beta d_{12} & \beta d_{21} \\ \beta d_{12} & -\mu_M - \beta d_{21} \end{pmatrix}$, this matrix is Hurwitz. Functions $M_1, M_2$ satisfy $\begin{pmatrix} \dot{M}_1 \\ \dot{M}_2 \end{pmatrix} = A \begin{pmatrix} M_1 \\ M_2 \end{pmatrix} + (1-r)\nu_E \begin{pmatrix} E_1 \\ E_2 \end{pmatrix}$. Thus, for any $t > 0$,

$$
\begin{pmatrix} M_1(t) \\ M_2(t) \end{pmatrix} = e^{tA} \begin{pmatrix} M_1^0 \\ M_2^0 \end{pmatrix} + (1-r)\nu_E \int_0^t e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds. \quad (13)
$$

Moreover, the equilibrium satisfies

$$
\begin{pmatrix} M_1^+ \\ M_2^+ \end{pmatrix} = -(1-r)\nu_E A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}. \qquad (14)
$$

Hence, from (14) and (13), we deduce that

$$
\begin{pmatrix} M_1(t) \\ M_2(t) \end{pmatrix} - \begin{pmatrix} M_1^+ \\ M_2^+ \end{pmatrix} = e^{tA} \begin{pmatrix} M_1^0 \\ M_2^0 \end{pmatrix} + (1-r)\nu_E \left[ \int_0^t e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds + A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix} \right]. \quad (15)
$$

Moreover, when $t \to +\infty$, we have that $\begin{pmatrix} E_1(t) \\ E_2(t) \end{pmatrix}$ converges to $\begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}$ and $e^{tA} \to 0$ since A is Hurwitz. Thus, for any $\varepsilon > 0$, there exists a time $T_\varepsilon > 0$ large enough such that for any $t > T_\varepsilon$,

$$
E_i^+ - \varepsilon < E_i(t) < E_i^+ + \varepsilon, \quad i = 1, 2, \qquad (16)
$$

and $\|e^{tA}\| < \|e^{T_\varepsilon A}\| \le \varepsilon$.

Since matrix $A$ is Metzler and irreducible, then by applying Lemma 2.1, one has that $e^{At}$ is strictly positive for any $t > 0$. Moreover, one has $E_i \in (0, K_i)$ in $(0, +\infty)$, then for any $t > 2T_\varepsilon$,

$$
0 < \int_{0}^{T_{\varepsilon}} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds < \int_{0}^{T_{\varepsilon}} e^{(t-s)A} ds \begin{pmatrix} K_1 \\ K_2 \end{pmatrix},
$$

then

$$
\left\| \int_{0}^{T_{\varepsilon}} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds \right\| < \left\| \int_{0}^{T_{\varepsilon}} e^{(t-s)A} ds \right\| \left\| \begin{pmatrix} K_1 \\ K_2 \end{pmatrix} \right\| = \left\| A^{-1} \left( e^{tA} - e^{(t-T_{\varepsilon})A} \right) \right\| \left\| \begin{pmatrix} K_1 \\ K_2 \end{pmatrix} \right\| < \varepsilon C_1,
$$
the invert function of $f_{12}$ on $[0, K_1]$, then we also have $f_{21} - f_{12}^{-1}$ convex on $(0, K_1)$, and $(f_{21} - f_{12}^{-1})(0) = -f_{12}^{-1}(0) = -K_2^+ < 0$, $\lim_{x \to K_1} (f_{21} - f_{12}^{-1})(x) = +\infty > 0$. We can deduce that $f_{21} - f_{12}^{-1}$ has a unique positive root on $(0, K_1)$.

**Instability of the zero equilibrium.** At the origin $\mathbf{0} = (0,0,0,0,0,0)$ of $\mathbb{R}^6$, the Jacobian matrix of system (6) is

$$ J(\mathbf{0}) = \begin{pmatrix} -v_E - \mu_E & b & 0 & 0 & 0 & 0 \\ r v_E & -\mu_F - d_{12} & 0 & 0 & d_{21} & 0 \\ (1-r)v_E & 0 & -\mu_M - \beta d_{12} & 0 & 0 & \beta d_{21} \\ 0 & 0 & 0 & -v_E - \mu_E & b & 0 \\ 0 & d_{12} & 0 & r v_E & -\mu_F - d_{21} & 0 \\ 0 & 0 & \beta d_{12} & (1-r)v_E & 0 & -\mu_F - \beta d_{21} \end{pmatrix}, $$

with the characteristic polynomial

$$ \det(J(0) - \lambda I) = [(\lambda + \mu_F)(\lambda + \mu_M) + \beta(d_{12} + d_{21} + d_{12}\mu_F + d_{21}\mu_M)][(\lambda + v_E + \mu_E)(\lambda + \mu_F) - brv_E] \\ \times [(\lambda + v_E + \mu_E)(\lambda + \mu_F + d_{12} + d_{21}) - brv_E]. $$

Since $\mathcal{N} > 1$, we have $\mu_F(v_E + \mu_E) - brv_E < 0$. Thus, we can deduce that the factor $(\lambda + v_E + \mu_E)(\lambda + \mu_F) - brv_E = \lambda^2 + \lambda(v_E + \mu_E + \mu_F) + \mu_F(v_E + \mu_E) - brv_E$ has one positive root $\lambda > 0$. Hence, the Jacobian at zero has at least one positive eigenvalue so the zero equilibrium is unstable.

**Stability of the positive equilibrium.** First, we can see that the system (6a)-(6b), (6d)-(6e) of $(E_1, F_1, E_2, F_2)$ does not depend on $M_1, M_2$, and it is cooperative and irreducible. By applying Theorem 1.1 in Chapter 4 of [33], one deduces that this system is strongly monotone. When $\mathcal{N} > 1$, this system admits exactly two equilibria: $(0,0,0,0)$, and $(E_1^+, F_1^+, E_2^+, F_2^+)$. But the zero equilibrium is unstable, so by Theorem 2.2 in Chapter 2 of [33], if the initial data satisfies that $0 < (E_1^0, F_1^0, E_2^0, F_2^0) \le (E_1^+, F_1^+, E_2^+, F_2^+)$, the solution $(E_1, F_1, E_2, F_2)$ converges to $(E_1^+, F_1^+, E_2^+, F_2^+)$ when $t \to +\infty$.

Now if the initial data satisfies that $(E_1^0, F_1^0, E_2^0, F_2^0) > (E_1^+, F_1^+, E_2^+, F_2^+)$, then there exists a constant $\lambda > 1$ large enough such that $\lambda(E_1^+, F_1^+, E_2^+, F_2^+) \ge (E_1^0, F_1^0, E_2^0, F_2^0)$. Since $1 - \frac{\lambda E_i^+}{K_i} < 1 - \frac{E_i^+}{K_i}$, one has

$$ b\lambda F_i^+ \left(1 - \frac{\lambda E_i^+}{K_i}\right) - (v_E + \mu_E)\lambda E_i^+ < \lambda \left[bF_i^+ \left(1 - \frac{E_i^+}{K_i}\right) - (v_E + \mu_E)E_i^+\right] = 0, $$

and the right-hand side of system (6a)-(6b), (6d)-(6e) at $\lambda(E_1^+, F_1^+, E_2^+, F_2^+) = (0,0,0,0)$ is non-positive. Thus, the trajectory resulting from the initial data $\lambda(E_1^+, F_1^+, E_2^+, F_2^+)$ is non-increasing, and therefore converges to $(E_1^+, F_1^+, E_2^+, F_2^+)$. By applying the Lemma 2.2 to system (6a)-(6b), (6d)-(6e), we deduce that the trajectory resulting from the initial data $(E_1^0, F_1^0, E_2^0, F_2^0)$ lies between $(E_1^+, F_1^+, E_2^+, F_2^+)$ and the trajectories resulting from $\lambda(E_1^+, F_1^+, E_2^+, F_2^+)$. Hence, it also converges to $(E_1^+, F_1^+, E_2^+, F_2^+)$ when time $t$ goes to infinity.

Moreover, since the trajectories issued from the initial data above and below $(E_1^+, F_1^+, E_2^+, F_2^+)$ all converge to the same limit, then by the comparison principle, we deduce that the trajectory resulting from any positive initial data with values between these initial values converges to this equilibrium.

Secondly, if we denote matrix $A = \begin{pmatrix} -\mu_M - \beta d_{12} & \beta d_{21} \\ \beta d_{12} & -\mu_M - \beta d_{21} \end{pmatrix}$, this matrix is Hurwitz. Functions $M_1, M_2$ satisfy $\begin{pmatrix} M_1 \\ M_2 \end{pmatrix} = A \begin{pmatrix} M_1 \\ M_2 \end{pmatrix} + (1-r)v_E \begin{pmatrix} E_1 \\ E_2 \end{pmatrix}$. Thus, for any $t > 0$,

$$ \begin{pmatrix} M_1(t) \\ M_2(t) \end{pmatrix} = e^{tA} \begin{pmatrix} M_1^0 \\ M_2^0 \end{pmatrix} + (1-r)v_E \int_0^t e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds. \quad (12) $$
with some positive constant $C_1$ not depending on $\varepsilon$. Using the second inequality in (16), one has

$$ \int_{T_{\varepsilon}}^{t} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds \leq A^{-1}(e^{(t-T_{\varepsilon})A} - I) \begin{pmatrix} E_1^+ + \varepsilon \\ E_2^+ + \varepsilon \end{pmatrix} = A^{-1}e^{(t-T_{\varepsilon})A} \begin{pmatrix} E_1^+ + \varepsilon \\ E_2^+ + \varepsilon \end{pmatrix} - \varepsilon A^{-1} - A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}. $$

Proving similarly for the other inequality, we can deduce that

$$ \left\| \int_{T_{\varepsilon}}^{t} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds + A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix} \right\| \leq \varepsilon C_2, $$

with some positive constant $C_2$ not depending on $\varepsilon$. Hence, from (15), we deduce that for any $t > 2T_{\varepsilon}$,

$$ \left\| \begin{pmatrix} M_1(t) \\ M_2(t) \end{pmatrix} - \begin{pmatrix} M_1^+ \\ M_2^+ \end{pmatrix} \right\| < \varepsilon \left[ \left\| \begin{pmatrix} M_1^0 \\ M_2^0 \end{pmatrix} \right\| + C_1 + C_2 \right]. $$

Therefore, $(M_1(t), M_2(t))$ converges to $(M_1^+, M_2^+)$ when $t$ tends to $+\infty$. $\square$

In the following section, by considering the releases of sterile males, we look for a condition of release functions $\Lambda$ such that the positive equilibrium disappears.

# 4 Elimination with releases of sterile males

In this section, we consider $\Lambda(t)$ in system (1) the number of sterile males released per time unit, and our goal is to adjust its values such that the wild population reaches elimination. We consider two release strategies as follows

**Constant release:** Let the release function $\Lambda(t) \equiv \Lambda > 0$. As time goes to infinity, the density of sterile males $(M_1^s, M_2^s)$ converges to $(M_1^{s*}, M_2^{s*})$ that is the solution of system

$$ \begin{cases} \Lambda - \mu_s M_1^{s*} - \alpha d_{12} M_1^{s*} + \alpha d_{21} M_2^{s*} = 0 \\ -\mu_s M_2^{s*} - \alpha d_{21} M_2^{s*} + \alpha d_{12} M_1^{s*} = 0 \end{cases} \quad (17) $$

By denoting

$$ \tau_1 := \frac{(\mu_s + \alpha d_{21})}{\mu_s(\mu_s + \alpha d_{12} + \alpha d_{21})}, \quad \tau_2 := \frac{\alpha d_{12}}{\mu_s(\mu_s + \alpha d_{12} + \alpha d_{21})}, \quad (18) $$

we have $M_1^{s*} = \tau_1 \Lambda$, and $M_2^{s*} = \tau_2 \Lambda$ and $M_1^{s*} + M_2^{s*} = \frac{\Lambda}{\mu_s}$.

**Impulsive periodic releases:** Consider the release function

$$ \Lambda(t) = \sum_{k=0}^{+\infty} \tau \Lambda_k^{\text{per}} \delta_{k\tau}, \quad (19) $$

with period $\tau > 0$ and $\Lambda_k^{\text{per}}$ is the average number of sterile males released per time unit during the time interval $(k\tau, (k+1)\tau)$ for $k = 0, 1, \dots$. We choose in this work $\Lambda_k^{\text{per}}$ constant and drop consequently the sub-index $k$. The release function $\Lambda(t)$ in (19) means that we release a total amount of $\tau \Lambda^{\text{per}}$ mosquitoes at the beginning of each time period $(t = k\tau)$.

Denote vector $X(t) = \begin{pmatrix} M_1^s(t) \\ M_2^s(t) \end{pmatrix}$, then with $k = 0, 1, \dots$, the density of sterile males satisfies the following system

$$ X'(t) = A_s X(t) \quad \text{for any } t \in \bigcup_{k=0}^{\infty} (k\tau, (k+1)\tau), \quad (20a) $$

$$ X(k\tau^+) = X(k\tau^-) + \begin{pmatrix} \tau\Lambda_k^{\text{per}} \\ 0 \end{pmatrix}, \quad (20b) $$

with matrix $A_s = \begin{pmatrix} -\alpha d_{12} - \mu_s & \alpha d_{21} \\ \alpha d_{12} & -\alpha d_{21} - \mu_s \end{pmatrix}$, and $X(k\tau^\pm)$ denote the right and left limits of $X(t)$ at time $k\tau$ and by convention, we set $X(0^-) = 0$.

The following lemma shows that $X(t)$ converges to a periodic solution as time evolves and this periodic solution is discontinuous and bounded from below by a positive value.
Moreover, the equilibrium satisfies

$$ \begin{pmatrix} M_1^+ \\ M_2^+ \end{pmatrix} = -(1-r)\nu_E A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}. \quad (13) $$

Hence, from (13) and (12), we deduce that

$$ \begin{pmatrix} M_1(t) \\ M_2(t) \end{pmatrix} - \begin{pmatrix} M_1^+ \\ M_2^+ \end{pmatrix} = e^{tA} \begin{pmatrix} M_1^0 \\ M_2^0 \end{pmatrix} + (1-r)\nu_E \left[ \int_0^t e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds + A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix} \right]. \quad (14) $$

Moreover, when $t \to +\infty$, we have that $\begin{pmatrix} E_1(t) \\ E_2(t) \end{pmatrix}$ converges to $\begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}$ and $e^{tA} \to 0$ since A is Hurwitz. Thus, for any $\varepsilon > 0$, there exists a time $T_\varepsilon > 0$ large enough such that for any $t > T_\varepsilon$,

$$ E_i^+ - \varepsilon < E_i(t) < E_i^+ + \varepsilon, \quad i = 1, 2, \quad (15) $$

and $\|e^{tA}\| < \|e^{T_\varepsilon A}\| \le \varepsilon$.

Since matrix A is Metzler and irreducible, then by applying Lemma 2.1, one has that $e^{At}$ is strictly positive for any $t > 0$. Moreover, one has $E_i \in (0, K_i)$ in $(0, +\infty)$, then for any $t > 2T_\varepsilon$,

$$ 0 < \int_0^{T_\varepsilon} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds < \int_0^{T_\varepsilon} e^{(t-s)A} ds \begin{pmatrix} K_1 \\ K_2 \end{pmatrix}, $$

then

$$ \left\| \int_{0}^{T_{\varepsilon}} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds \right\| < \left\| \int_{0}^{T_{\varepsilon}} e^{(t-s)A} ds \right\| \left\| \begin{pmatrix} K_1 \\ K_2 \end{pmatrix} \right\| = \left\| A^{-1} \left( e^{tA} - e^{(t-T_{\varepsilon})A} \right) \right\| \left\| \begin{pmatrix} K_1 \\ K_2 \end{pmatrix} \right\| < \varepsilon C_1, $$

with some positive constant $C_1$ not depending on $\varepsilon$. Using the second inequality in (15), one has

$$ \int_{T_{\varepsilon}}^{t} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds \le A^{-1} (e^{(t-T_{\varepsilon})A} - I) \begin{pmatrix} E_1^+ + \varepsilon \\ E_2^+ + \varepsilon \end{pmatrix} = A^{-1} e^{(t-T_{\varepsilon})A} \begin{pmatrix} E_1^+ + \varepsilon \\ E_2^+ + \varepsilon \end{pmatrix} - \varepsilon A^{-1} - A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix}. $$

Proving similarly for the other inequality, we can deduce that

$$ \left\| \int_{T_{\varepsilon}}^{t} e^{(t-s)A} \begin{pmatrix} E_1(s) \\ E_2(s) \end{pmatrix} ds + A^{-1} \begin{pmatrix} E_1^+ \\ E_2^+ \end{pmatrix} \right\| \le \varepsilon C_2, $$

with some positive constant $C_2$ not depending on $\varepsilon$. Hence, from (14), we deduce that for any $t > 2T_\varepsilon$,

$$ \left\| \begin{pmatrix} M_1(t) \\ M_2(t) \end{pmatrix} - \begin{pmatrix} M_1^+ \\ M_2^+ \end{pmatrix} \right\| < \varepsilon \left[ \left\| \begin{pmatrix} M_1^0 \\ M_2^0 \end{pmatrix} \right\| + C_1 + C_2 \right]. $$

Therefore, $(M_1(t), M_2(t))$ converges to $(M_1^+, M_2^+)$ when $t$ tends to $+\infty$. $\square$

In the following section, by considering the releases of sterile males, we look for a condition of release functions $\Lambda$ such that the positive equilibrium disappears.

# 4 Elimination with releases of sterile males

In this section, we consider $\Lambda(t)$ in system (1) the number of sterile males released per time unit and our goal is to adjust its values such that the wild population reaches elimination. We consider two release strategies as follows
**Lemma 4.1.** The solution $X(t)$ of system (20) converges to a periodic solution $X^{\text{per}}(t)$ that satisfies

$$X^{\text{per}}(k\tau^+) = (I - e^{A_s\tau})^{-1} \begin{pmatrix} \tau\Lambda^{\text{per}} \\ 0 \end{pmatrix}, \quad X^{\text{per}}(t) = e^{A_s t} X^{\text{per}}(k\tau^+) \quad \text{for any } t \in \bigcup_{k=0}^{\infty} (k\tau, (k+1)\tau). \quad (21)$$

Moreover, there exist positive constants $\tau_1^{\text{per}}$, $\tau_2^{\text{per}}$ which depend on $A_s$ and period $\tau$ such that for any $t > 0$, one has $X^{\text{per}}(t) \ge \Lambda^{\text{per}} \begin{pmatrix} \tau_1^{\text{per}} \\ \tau_2^{\text{per}} \end{pmatrix}$.

*Proof.* The densities of the sterile males evolve according to (20a) on the union of open intervals $(k\tau, (k+1)\tau)$ while $X$ is submitted to jump at each point $k\tau$ as in (20b). For such a release schedule, the solution of system (20) satisfies

$$X(k\tau^+) = \sum_{i=0}^{k} e^{iA_s\tau} \begin{pmatrix} \tau\Lambda^{\text{per}} \\ 0 \end{pmatrix}, \quad X(t) = e^{A_s t} X(k\tau^+) \quad \text{for any } t \in \bigcup_{k=0}^{\infty} (k\tau, (k+1)\tau).$$

Since matrix $A_s$ is Hurwitz, when $t \to +\infty$, we have that $X$ converges to the periodic solution $X^{\text{per}} = \begin{pmatrix} M_1^{s,\text{per}} \\ M_2^{s,\text{per}} \end{pmatrix}$ that satisfies that (21) for $k = 0, 1, \dots$.

We have matrix $A_s$ is Metzler and irreducible, then by applying Lemma 2.1, we deduce that $e^{A_s} \gg 0$.

On the other hand, matrix $A_s$ is Hurwitz so $(I - e^{A_s\tau})^{-1} = \sum_{i=0}^{+\infty} e^{iA_s\tau} \gg 0$ for any $\tau > 0$. Moreover, we have $e^{A_s t}$ is also strictly positive for any $t > 0$, so

$$\inf_{t>0} X^{\text{per}}(t) = \min_{t \in [0, \tau]} X^{\text{per}}(t) = \min_{t \in [0, \tau]} e^{A_s t} (I - e^{A_s \tau})^{-1} \begin{pmatrix} \tau \Lambda^{\text{per}} \\ 0 \end{pmatrix}.$$

By taking $\begin{pmatrix} \tau_1^{\text{per}} \\ \tau_2^{\text{per}} \end{pmatrix} = \min_{t \in [0, \tau]} e^{A_s t} (I - e^{A_s \tau})^{-1} \begin{pmatrix} \tau \\ 0 \end{pmatrix}$, the result of Lemma 4.1 follows. $\square$

**Remark 4.1.** The parameters $\tau_i$ defined in (18) and $\tau_i^{\text{per}}$ play a similar role to each other: they define a relationship between an average release rate $\Lambda$ per time unit and a (minimum) level of the sterile mosquito density. These parameters depend on the diffusion rates, the death rate of sterile males, and the release period $\tau$ in the periodic case.

We provide in the following result a condition on $\Lambda(t)$ for the wild population to reach elimination.

**Theorem 4.2.** Consider system (1) with the release function $\Lambda(t)$. Then

* In the constant release case, for $\Lambda(t) \equiv \Lambda$, there exists a positive number $\bar{\Lambda}$ satisfying

$$\bar{\Lambda} \le \max_{i=1,2} \frac{1}{\gamma \tau_i} (\mathcal{N} - 1) C_M \quad i = 1, 2,$$

with $\tau_i$ defined in (18), $C_M$ defined in Lemma 2.4 such that if $\Lambda > \bar{\Lambda}$, system (1) has a unique equilibrium

$$\mathbf{u}_0^* = (0, 0, 0, M_1^{s*}, 0, 0, 0, M_2^{s*}).$$

Moreover, in this case, for any non-negative initial data, the solution of (1) converges to this equilibrium when $t \to +\infty$.

* In the periodic release case, for $\Lambda(t)$ defined in (20), There exists a positive constant $\bar{\Lambda}^{\text{per}}$ satisfying

$$\bar{\Lambda}^{\text{per}} \le \max_{i=1,2} \frac{1}{\gamma \tau_i^{\text{per}}} (\mathcal{N} - 1) C_M,$$

with $\tau_i^{\text{per}}$ defined in Lemma 4.1, $C_M$ defined in Lemma 2.4 such that if $\Lambda^{\text{per}} > \bar{\Lambda}^{\text{per}}$, then for any non-negative initial data, the solution of the initial value problem of system (1) converges to the unique steady state

$$\mathbf{u}_0^{\text{per}} = (0, 0, 0, M_1^{s,\text{per}}, 0, 0, 0, M_2^{s,\text{per}}),$$

as time $t \to +\infty$.

This result shows that with a sufficiently large number of sterile males released in the first zone, we can succeed in driving the wild population in both areas to elimination. In the following, we describe the principle idea to prove this result.
**Constant release:** Let the release function $\Lambda(t) \equiv \Lambda > 0$. As time goes to infinity, the density of sterile males $(M_1^s, M_2^s)$ converges to $(M_1^{s*}, M_2^{s*})$ that is the solution of system

$$ \begin{cases} \Lambda - \mu_s M_1^{s*} - \alpha d_{12} M_1^{s*} + \alpha d_{21} M_2^{s*} = 0 \\ -\mu_s M_2^{s*} - \alpha d_{21} M_2^{s*} + \alpha d_{12} M_1^{s*} = 0 \end{cases} \quad (16) $$

By denoting

$$ \tau_1 := \frac{(\mu_s + \alpha d_{21})}{\mu_s(\mu_s + \alpha d_{12} + \alpha d_{21})}, \quad \tau_2 := \frac{\alpha d_{12}}{\mu_s(\mu_s + \alpha d_{12} + \alpha d_{21})}, \quad (17) $$

we have $M_1^{s*} = \tau_1 \Lambda$, and $M_2^{s*} = \tau_2 \Lambda$ and $M_1^{s*} + M_2^{s*} = \frac{\Lambda}{\mu_s}$.

**Impulsive periodic releases:** Consider the release function

$$ \Lambda(t) = \sum_{k=0}^{+\infty} \tau \Lambda_k^{\text{per}} \delta_{k\tau}, \quad (18) $$

with period $\tau > 0$ and $\Lambda_k^{\text{per}}$ is the average number of sterile males released per time unit during the time interval $(k\tau, (k+1)\tau)$ for $k = 0, 1, \dots$. We choose in this work $\Lambda_k^{\text{per}}$ constant and drop consequently the sub-index $k$. The release function $\Lambda(t)$ in (18) means that we release a total amount of $\tau\Lambda^{\text{per}}$ mosquitoes at the beginning of each time period $(t = k\tau)$.

Denote vector $X(t) = \begin{pmatrix} M_1^s(t) \\ M_2^s(t) \end{pmatrix}$, then with $k = 0, 1, \dots$, the density of sterile males satisfies the following system

$$ X'(t) = A_s X(t) \quad \text{for any } t \in \bigcup_{k=0}^{\infty} (k\tau, (k+1)\tau), \quad (19a) $$

$$ X(k\tau^+) = X(k\tau^-) + \begin{pmatrix} \tau\Lambda^{\text{per}} \\ 0 \end{pmatrix}, \quad (19b) $$

with matrix $A_s = \begin{pmatrix} -\alpha d_{12} - \mu_s & \alpha d_{21} \\ \alpha d_{12} & -\alpha d_{21} - \mu_s \end{pmatrix}$, and $X(k\tau^\pm)$ denote the right and left limits of $X(t)$ at time $k\tau$ and by convention, we set $X(0^-) = 0$. The densities of the sterile males evolve according to (19a) on the union of open intervals $(k\tau, (k+1)\tau)$ while $X$ is submitted to jump at each point $k\tau$ as in (19b). For such a release schedule, the solution of system (19) satisfies

$$ X(k\tau^+) = \sum_{i=0}^{k} e^{iA_s\tau} \begin{pmatrix} \tau\Lambda^{\text{per}} \\ 0 \end{pmatrix}, \quad X(t) = e^{A_s t} X(k\tau^+) \quad \text{for any } t \in \bigcup_{k=0}^{\infty} (k\tau, (k+1)\tau). $$

Since matrix $A_s$ is Hurwitz, when $t \to +\infty$, we have that $X$ converges to the periodic solution $X^{\text{per}} = \begin{pmatrix} M_1^{s,\text{per}} \\ M_2^{s,\text{per}} \end{pmatrix}$ that satisfies that, for $k = 0, 1, \dots$

$$ X^{\text{per}}(k\tau^+) = (I - e^{A_s \tau})^{-1} \begin{pmatrix} \tau \Lambda^{\text{per}} \\ 0 \end{pmatrix}, \quad X^{\text{per}}(t) = e^{A_s t} X^{\text{per}}(k\tau^+) \quad \text{for any } t \in \bigcup_{k=0}^{\infty} (k\tau, (k+1)\tau). \quad (20) $$

The following lemma shows that the periodic solution $X^{\text{per}}$ is strictly positive at any time $t$.

**Lemma 4.1.** There exists positive constant $\tau_1^{\text{per}}, \tau_2^{\text{per}}$ which depends on $A_s$ and period $\tau$ such that for any $t > 0$, one has $X^{\text{per}}(t) \ge \Lambda^{\text{per}} \begin{pmatrix} \tau_1^{\text{per}} \\ \tau_2^{\text{per}} \end{pmatrix}$.
## 4.1 Principle of the method

To provide conditions for the release $\Lambda$ to stabilize the zero equilibrium, our strategy is as follows:

Step 1: We consider $\rho_i = \sup_{t>0} \frac{M_i(t)}{M_i(t) + \gamma M_i^s(t)}$, for $i = 1, 2$, in system (1) to be smaller than some level, then we study the system with the fractions replaced by some constant.

Step 2: We show how to realize, through an adequate choice of $\Lambda$, the above behavior of $M_i^s$.

### 4.1.1 Step 1: Setting the sterile population level directly

Theorem 3.1 shows us that when the basic offspring number is smaller than 1, the zero equilibrium is globally asymptotically stable. For the controlled system, the basic offspring numbers are smaller than $\rho_i N$. It suggests that for stabilizing the origin of system (1), it is sufficient to ensure $\rho_i N \le 1$.

**Proposition 4.3.** If the trajectory resulting from any positive initial data of system (1) satisfies that for $N$ defined in (2),

$$ \frac{M_i(t)}{M_i(t) + \gamma M_i^s(t)} \le \frac{1}{N}, \quad t \ge 0, i = 1, 2. \quad (22) $$

then $\mathbf{u}' = (E_1, F_1, M_1, E_2, F_2, M_2)$ converges to $\mathbf{0}_6$ as time $t \to +\infty$.

*Proof.* Assume that we can set $M_i^s$ to be large enough such that (22) holds, and we consider the following system

$$ \dot{E}_1 = bF_1 \left(1 - \frac{E_1}{K_1}\right) - (\nu_E + \mu_E)E_1, \quad (23a) $$

$$ \dot{F}_1 = \frac{1}{N} r \nu_E E_1 - \mu_F F_1 - d_{12} F_1 + d_{21} F_2, \quad (23b) $$

$$ \dot{M}_1 = (1-r)\nu_E E_1 - \mu_M M_1 - \beta d_{12} M_1 + \beta d_{21} M_2, \quad (23c) $$

$$ \dot{E}_2 = bF_2 \left(1 - \frac{E_2}{K_2}\right) - (\nu_E + \mu_E)E_2, \quad (23d) $$

$$ \dot{F}_2 = \frac{1}{N} r \nu_E E_2 - \mu_F F_2 - d_{21} F_2 + d_{12} F_1, \quad (23e) $$

$$ \dot{M}_2 = (1 - r)\nu_E E_2 - \mu_M M_2 - \beta d_{12} M_2 + \beta d_{21} M_1, \quad (23f) $$

Denote $\tilde{\mathbf{u}} = (\tilde{E}_1, \tilde{F}_1, \tilde{M}_1, \tilde{E}_2, \tilde{F}_2, \tilde{M}_2)$ solution of system (23). Since system (23) is cooperative and the inequality (22) holds, one obtains that $\tilde{\mathbf{u}}$ is a super-solution of the system (1a)-(1c), (1e)-(1g), and by applying Lemma 2.2, we have $\tilde{\mathbf{u}} \ge \mathbf{u}'$.

Denote $\mathbf{u}^* = (E_1^*, F_1^*, M_1^*, E_2^*, F_2^*, M_2^*)$ a positive equilibrium of system (23) if exists. Similar to the previous section, we have

$$ E_2^* = \frac{\mu_F + d_{12} + d_{21}}{d_{21}} \frac{E_1^*}{1 - \frac{E_1^*}{K_1}} - \frac{\mu_F + d_{21}}{d_{21}} E_1^* =: g_{21}(E_1^*), \quad (24) $$

$$ E_1^* = \frac{\mu_F + d_{12} + d_{21}}{d_{12}} \frac{E_2^*}{1 - \frac{E_2^*}{K_2}} - \frac{\mu_F + d_{12}}{d_{12}} E_2^* =: g_{12}(E_2^*). \quad (25) $$

The analysis of $g_{ij}$ is analogous to $f_{ij}$ in Lemma 3.2. It is easy to check that $g_{12}$ is increasing on $(0, K_2)$, so it is invertible. Then, $E_1^*$ satisfies $(g_{21} - g_{12}^{-1})(E_1^*) = 0$. Function $g_{21}$ is convex and $g_{12}^{-1}$ is concave in $(0, K_1)$, and

$$ g'_{21}(0) = \frac{\mu_F + d_{12} + d_{21}}{d_{21}} - \frac{(\mu_F + d_{21})}{d_{21}} = \frac{d_{12}}{d_{21}}, \\ (g_{12}^{-1})'(0) = \frac{1}{g'_{12}(0)} = \frac{1}{\frac{\mu_F + d_{12} + d_{21}}{d_{12}} - \frac{\mu_F + d_{12}}{d_{12}}} = \frac{d_{12}}{d_{21}}. $$

We obtain that $g'_{21}(0) = (g_{12}^{-1})'(0)$ (see Figure 2), so zero is the unique equilibrium of system (23). By applying Theorem 3.1 in Chapter 2 of [39], we deduce that when $t \to +\infty$, the solution $\tilde{\mathbf{u}}(t)$ converges to the equilibrium zero. Since $\mathbf{u}'(t) \le \tilde{\mathbf{u}}(t)$ for all $t > 0$, we deduce that the $\mathbf{u}'$ also converges to zero when $t$ large. $\square$

### 4.1.2 Step 2: Shaping the release function

We now want to choose $\Lambda$ such that the condition (22) holds, which means

$$ \gamma M_i^s(t) \ge (N - 1) M_i(t), \quad t \ge 0, i = 1, 2 $$

The upper bound of $M_i$ can be obtained from Lemma 2.4. For time $t > 0$ large enough, it is sufficient to choose $\Lambda$ such that $\gamma M_i^s(t) \ge (N - 1) C_M$.
*Proof.* We have matrix $A_s$ is Metzler and irreducible, then by applying Lemma 2.1, we deduce that $e^{A_s} \gg 0$.

On the other hand, matrix $A_s$ is Hurwitz so $(I - e^{A_s} \tau)^{-1} = \sum_{i=0}^{+\infty} e^{iA_s \tau} \gg 0$ for any $\tau > 0$. Moreover, we have $e^{A_s t}$ is also strictly positive for any $t > 0$, thus there exist positive constants $\tau_1^{\text{per}}, \tau_2^{\text{per}}$ depending on $\tau$ and $A_s$ such that

$$ \inf_{t>0} X^{\text{per}}(t) = \min_{t \in [0, \tau]} X^{\text{per}}(t) = \min_{t \in [0, \tau]} e^{A_s t} (I - e^{A_s} \tau)^{-1} \begin{pmatrix} \tau \Lambda^{\text{per}} \\ 0 \end{pmatrix} \geq \Lambda^{\text{per}} \begin{pmatrix} \tau_1^{\text{per}} \\ \tau_2^{\text{per}} \end{pmatrix}. $$

The result of Lemma 4.1 follows.

□

**Remark 4.1.** The parameters $\tau_i$ defined in (17) and $\tau_i^{\text{per}}$ play a similar role to each other: they define a relationship between an average release rate $\Lambda$ per time unit and a (minimum) level of the sterile mosquito density.

We provide in the following result a condition on $\Lambda(t)$ for the wild population to reach elimination.

**Theorem 4.2.** Consider system (1) with the release function $\Lambda(t)$. Then

* In the constant release case, for $\Lambda(t) \equiv \Lambda$, there exists a positive number $\bar{\Lambda}$ satisfying

$$ \bar{\Lambda} \leq \max_{i=1,2} \frac{1}{\gamma \tau_i} (\mathcal{N} - 1) C_M \quad i = 1,2, $$

with $\tau_i$ defined in (17), $C_M$ defined in Lemma 2.4 such that if $\Lambda > \bar{\Lambda}$, system (1) has a unique equilibrium

$$ \mathbf{u}_0^* = (0, 0, 0, M_1^{s*}, 0, 0, 0, M_2^{s*}). $$

Moreover, in this case, for any non-negative initial data, the solution of (1) converges to this equilibrium when $t \to +\infty$.

* In the periodic release case, for $\Lambda(t)$ defined in (19), There exists a positive constant $\bar{\Lambda}^{\text{per}}$ satisfying

$$ \bar{\Lambda}^{\text{per}} \leq \max_{i=1,2} \frac{1}{\gamma \tau_i^{\text{per}}} (\mathcal{N} - 1) C_M, $$

with $\tau_i^{\text{per}}$ defined in Lemma 4.1, $C_M$ defined in Lemma 2.4 such that if $\Lambda^{\text{per}} > \bar{\Lambda}^{\text{per}}$, then for any non-negative initial data, the solution of the initial value problem of system (1) converges to the unique steady state

$$ \mathbf{u}_0^{\text{per}} = (0, 0, 0, M_1^{s, \text{per}}, 0, 0, 0, M_2^{s, \text{per}}), $$

as time $t \to +\infty$.

This result shows that with a sufficiently large number of sterile males released in the first zone, we can succeed in driving the wild population in both areas to elimination. In the following, we describe the principle idea to prove this result.

## 4.1 Principle of the method

To provide conditions for the release $\Lambda$ to stabilize the zero equilibrium, our strategy is as follows:

Step 1: We consider $\rho_i = \sup_{t>0} \frac{M_i(t)}{M_i(t) + \gamma M_i^s(t)}$, for $i = 1, 2$, in system (1) to be smaller than some level, then we study the system with the fractions replaced by some constant.

Step 2: We show how to realize, through an adequate choice of $\Lambda$, the above behavior of $M_i^s$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Relations $E_2^* = g_{21}(E_1^*)$ (in green) and $E_1^* = g_{12}(E_2^*)$ (in red) presented in the $E_1 - E_2$ plane when $\mathcal{N} > 1$. The two curves only intersect at the origin.</figcaption></figure>

## 4.2 Proof of Theorem 4.2

*Proof of Theorem 4.2. Constant release:* For $\Lambda(t) \equiv \Lambda$, let us recall $A_s = \begin{pmatrix} -\mu_s - \alpha d_{12} & \alpha d_{21} \\ \alpha d_{12} & -\mu_s - \alpha d_{21} \end{pmatrix}$, and it is a Hurwitz matrix. And we have

$$ \begin{pmatrix} M_1^s(t) \\ M_2^s(t) \end{pmatrix} = e^{tA_s} \begin{pmatrix} M_1^{s0} \\ M_2^{s0} \end{pmatrix} + (1-r)\nu_E \int_0^t e^{(t-s)A_s} ds \begin{pmatrix} \Lambda \\ 0 \end{pmatrix}. $$

Thus, when $t \to +\infty$, we can deduce that $(M_1^s(t) / M_2^s(t))$ converges to $(M_1^{s*} / M_2^{s*})$ for any initial data since $e^{tA_s} \to 0$ when $t \to +\infty$. Hence, for any $\varepsilon > 0$, there exists a value $T_\varepsilon > 0$ such that for any $t > T_\varepsilon$, $i = 1, 2$,

$$ M_i^s(t) \geq M_i^{s*} - \varepsilon. $$

If we take $\Lambda$ such that $\gamma(M_i^{s*} - \varepsilon) \geq (\mathcal{N} - 1)C_M$ with $C_M$ defined in Lemma 2.4, then condition (22) holds. By applying Proposition 4.3, we deduce that for $i = 1, 2$, if $\Lambda > \max_{i=1,2} \frac{1}{\gamma\tau_i}(\mathcal{N} - 1)C_M$, system (1) has a unique equilibrium point

$$ \mathbf{u}_0^* = (0, 0, 0, M_1^{s*}, 0, 0, 0, M_2^{s*}), $$

and every trajectory converges to this equilibrium when $t \to +\infty$. The dynamics of system (1) depend continuously and monotonically on $\Lambda$, then we deduce that there exists a positive critical value $\bar{\Lambda} \leq \max_{i=1,2} \frac{1}{\gamma\tau_i}(\mathcal{N} - 1)C_M$ such that for any $\Lambda > \bar{\Lambda}$, and for any non-negative initial data, solution $\mathbf{u}' = (E_1, F_1, M_1, E_2, F_2, M_2)(t)$ converges to $\mathbf{0}_6$ when $t \to +\infty$.

**Impulsive periodic releases:** Consider $\Lambda(t)$ defined in (20), denote $(E_1^{\text{per}}, F_1^{\text{per}}, M_1^{\text{per}}, E_2^{\text{per}}, F_2^{\text{per}}, M_2^{\text{per}})$ a solution of (1a)-(1c), (1e)-(1g) with $M_i^s \equiv M_i^{s,\text{per}}$ defined in (21). From Lemma 4.1, one has $M_i^{s,\text{per}}(t) \geq \Lambda^{\text{per}}\tau_i^{\text{per}}$ for all $t > 0$. Therefore, if we take $\Lambda^{\text{per}}$ such that $\gamma\Lambda^{\text{per}}\tau_i^{\text{per}} \geq (\mathcal{N} - 1)C_M$, then condition (22) holds. By applying Proposition 4.3, we deduce that $\mathbf{u}^{\text{per}} = (E_1^{\text{per}}, F_1^{\text{per}}, M_1^{\text{per}}, E_2^{\text{per}}, F_2^{\text{per}}, M_2^{\text{per}})$ converges to zero as $t$ grows. Since $(M_1^{s,\text{per}}, M_2^{s,\text{per}})$ converges to $(M_1^{s,\text{per}}, M_2^{s,\text{per}})$ as $t \to +\infty$, we have $\mathbf{u}' = (E_1, F_1, M_1, E_2, F_2, M_2)$ approaches $\mathbf{u}^{\text{per}}$ and thus converges to $\mathbf{0}_6$ as time $t$ goes to infinity.

Since the dynamics of system (1) depends continuously and monotonically on $\Lambda$, there exists a positive critical value $\bar{\Lambda}^{\text{per}} \leq \max_{i=1,2} \frac{1}{\gamma\tau_i}(\mathcal{N} - 1)C_M$ such that if $\Lambda^{\text{per}} > \bar{\Lambda}^{\text{per}}$ the equilibrium $\mathbf{u}_0^{\text{per}}$ of (1) with $\Lambda(t)$ defined in (19) is globally asymptotically stable. $\square$

# 5 Parameter dependence of the critical values of the release rate

In this section, we consider the constant release case and examine how the critical value $\bar{\Lambda}$ depends on the parameters of system (1). In this model, the elimination of the population depends not only on the diffusion rate between the inaccessible area and the treated area, but also on the biological intrinsic values like the birth/death rates, and the carrying capacity.
### 4.1.1 Step 1: Setting the sterile population level directly

Theorem 3.1 shows us that when the basic offspring number is smaller than 1, the zero equilibrium is globally asymptotically stable. For the controlled system, the basic offspring numbers are smaller than $\rho_i \mathcal{N}$. It suggests that for stabilizing the origin of system (1), it is sufficient to ensure $\rho_i \mathcal{N} \le 1$.

**Proposition 4.3.** If the trajectory resulting from any positive initial data of system (1) satisfies that for $\mathcal{N}$ defined in (2),

$$ \frac{M_i(t)}{M_i(t) + \gamma M_i^s(t)} \le \frac{1}{\mathcal{N}}, \quad t \ge 0, i = 1, 2. \qquad (21) $$

then $\mathbf{u}' = (E_1, F_1, M_1, E_2, F_2, M_2)$ converges to $\mathbf{0}_6$ as time $t \to +\infty$.

*Proof.* Assume that we can set $M_i^s$ to be large enough such that (21) holds, and we consider the following system

$$ \dot{E}_1 = bF_1 \left(1 - \frac{E_1}{K_1}\right) - (v_E + \mu_E)E_1, \qquad (22a) $$

$$ \dot{F}_1 = \frac{1}{\mathcal{N}} r v_E E_1 - \mu_F F_1 - d_{12} F_1 + d_{21} F_2, \qquad (22b) $$

$$ \dot{M}_1 = (1 - r)v_E E_1 - \mu_M M_1 - \beta d_{12} M_1 + \beta d_{21} M_2, \qquad (22c) $$

$$ \dot{E}_2 = bF_2 \left(1 - \frac{E_2}{K_2}\right) - (v_E + \mu_E)E_2, \qquad (22d) $$

$$ \dot{F}_2 = \frac{1}{\mathcal{N}} r v_E E_2 - \mu_F F_2 - d_{21} F_2 + d_{12} F_1, \qquad (22e) $$

$$ \dot{M}_2 = (1 - r)v_E E_2 - \mu_M M_2 - \beta d_{21} M_2 + \beta d_{12} M_1, \qquad (22f) $$

Denote $\tilde{\mathbf{u}} = (\tilde{E}_1, \tilde{F}_1, \tilde{M}_1, \tilde{E}_2, \tilde{F}_2, \tilde{M}_2)$ solution of system (22). Since system (22) is cooperative and the inequality (21) holds, one obtains that $\tilde{\mathbf{u}}$ is a super-solution of the system (1a)-(1c), (1e)-(1g), and by applying Lemma 2.2, we have $\tilde{\mathbf{u}} \ge \mathbf{u}'$.

Denote $\mathbf{u}^* = (E_1^*, F_1^*, M_1^*, E_2^*, F_2^*, M_2^*)$ a positive equilibrium of system (22) if exists. Similar to the previous section, we have

$$ E_2^* = \frac{\mu_F + d_{12} + d_{21}}{d_{21}} \frac{E_1^*}{1 - \frac{E_1^*}{K_1}} - \frac{\mu_F + d_{21}}{d_{21}} E_1^* =: g_{21}(E_1^*), \qquad (23) $$

$$ E_1^* = \frac{\mu_F + d_{12} + d_{21}}{d_{12}} \frac{E_2^*}{1 - \frac{E_2^*}{K_2}} - \frac{\mu_F + d_{12}}{d_{12}} E_2^* =: g_{12}(E_2^*). \qquad (24) $$

The analysis of $g_{ij}$ is analogous to $f_{ij}$ in Lemma 3.2. It is easy to check that $g_{12}$ is increasing on $(0, K_2)$, so it is invertible. Then, $E_1^*$ satisfies $(g_{21} - g_{12}^{-1})(E_1^*) = 0$. Function $g_{21}$ is convex and $g_{12}^{-1}$ is concave in $(0, K_1)$, and

$$ g'_{21}(0) = \frac{\mu_F + d_{12} + d_{21}}{d_{21}} - \frac{(\mu_F + d_{21})}{d_{21}} = \frac{d_{12}}{d_{21}}, \\ (g_{12}^{-1})'(0) = \frac{1}{g'_{12}(0)} = \frac{1}{\frac{\mu_F + d_{12} + d_{21}}{d_{12}} - \frac{\mu_F + d_{12}}{d_{12}}} = \frac{d_{12}}{d_{21}}. $$

We obtain that $g'_{21}(0) = (g_{12}^{-1})'(0)$ (see Figure 2), so zero is the unique equilibrium of system (22). By applying Theorem 3.1 in Chapter 2 of [33], we deduce that when $t \to +\infty$, the solution $\tilde{\mathbf{u}}(t)$ converges to the equilibrium zero. Since $\mathbf{u}'(t) \le \tilde{\mathbf{u}}(t)$ for all $t > 0$, we deduce that the $\mathbf{u}'$ also converges to zero when $t$ large. $\square$
## 5.1 Diffusion rates

In this part, we want to compare the critical values of $\Lambda$ corresponding to different values of $d_{12}$, $d_{21}$. We show that when the diffusion rates are large enough, the critical number of sterile males released is the same as in the case when there is no separation between the two sub-populations.

### 5.1.1 The case $d_{12}, d_{21}$ large

First, we present a result of uniform convergence of system (1) when $d_{12}$, $d_{21}$ go to $+\infty$ and $d_{12}$ is proportional to $d_{21}$.

**Proposition 5.1.** For $\varepsilon > 0$, consider the diffusion rates $d_{12} = \frac{1}{\varepsilon}$, $d_{21} = \frac{\eta}{\varepsilon}$ with $\eta = \frac{d_{21}}{d_{12}} > 0$. Denote $\mathbf{u}^{\varepsilon} = (E_1^{\varepsilon}, F_1^{\varepsilon}, M_1^{\varepsilon}, M_1^{s,\varepsilon}, E_2^{\varepsilon}, F_2^{\varepsilon}, M_2^{\varepsilon}, M_2^{s,\varepsilon})$ the solution of system (1) with the initial data $\mathbf{u}^{\varepsilon,0}$ satisfying that
$$ \{E_i^{\varepsilon,0}\}, \{F_i^{\varepsilon,0}\}, \{M_i^{\varepsilon,0}\} \text{ converge to } E_i^{0,0}, F_i^{0,0}, M_i^{0,0} \text{ respectively as } \varepsilon \to 0, \text{ with } i = 1, 2, $$

and

$$ E_1^{\varepsilon,0} - \eta E_2^{\varepsilon,0} = \mathcal{O}(\varepsilon), \quad F_1^{\varepsilon,0} - \eta F_2^{\varepsilon,0} = \mathcal{O}(\varepsilon), \quad M_1^{\varepsilon,0} - \eta M_2^{\varepsilon,0} = \mathcal{O}(\varepsilon), \quad M_1^{s,\varepsilon,0} = M_2^{s,\varepsilon,0} = 0. \quad (26) $$

Then, when $\varepsilon \to 0$, the sequence $\{\mathbf{u}^{\varepsilon}\}$ converges uniformly to a limit $(E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s)$ on $[0, +\infty)$. Moreover, we have

$$ F_1 = \eta F_2, \quad M_1 = \eta M_2, \quad M_1^s = \eta M_2^s. \quad (27) $$

If we denote $F = F_1 + F_2$, $M = M_1 + M_2$, $M^s = M_1^s + M_2^s$, then $(E_1, E_2, F, M, M^s)$ solves the following system

$$ \dot{E}_1 = \frac{\eta}{\eta+1} bF \left(1 - \frac{E_1}{K_1}\right) - (\nu_E + \mu_E)E_1, \quad (28a) $$

$$ \dot{E}_2 = \frac{1}{\eta+1} bF \left(1 - \frac{E_2}{K_2}\right) - (\nu_E + \mu_E)E_2, \quad (28b) $$

$$ \dot{F} = r\nu_E(E_1 + E_2) \frac{M}{M + \gamma M^s} - \mu_F F, \quad (28c) $$

$$ \dot{M} = (1 - r)\nu_E(E_1 + E_2) - \mu_M M, \quad (28d) $$

$$ \dot{M}^s = \Lambda_\infty - \mu_s M^s, \quad (28e) $$

with the corresponding initial data $E_1^{0,0}$, $E_2^{0,0}$, $F^{0,0} = F_1^{0,0} + F_2^{0,0}$, $M^{0,0} = M_1^{0,0} + M_2^{0,0}$, $M^{s,0,0} = 0$.

It is straightforward to see that the previous result implies that $F_1, F_2, M_1, M_2$ satisfy

$$ F_1 = \frac{\eta}{1 + \eta} F, \quad F_2 = \frac{1}{1 + \eta} F, \quad M_1 = \frac{\eta}{1 + \eta} M, \quad M_2 = \frac{1}{1 + \eta} M. $$

*Proof of Proposition 5.1.* To prove this result, we first apply the Arzela-Ascoli theorem for the sequence of smooth solution $\{\mathbf{u}^{\varepsilon}\}_{\varepsilon \in \text{a close interval } [0, T]}$ with any $T > 0$. Then, we extend the convergence at infinity.

* Uniform convergence on $[0, T]$:* First, we check the uniform boundedness of this sequence. For $i = 1, 2$, from Lemma 2.4, one has $E_i^{\varepsilon}(t) \le K_i$ for all $t > 0$ and $\varepsilon > 0$. Again by this Lemma, for any $t > 0$, one has

$$ F_i^{\varepsilon}(t) \le \max \left\{ F_1^{\varepsilon,0} + F_2^{\varepsilon,0}, \frac{r\nu_E(K_1 + K_2)}{\mu_F} \right\} \le C_F^0 $$

where $C_F^0$ does not depend on $\varepsilon$ since the initial data converge as $\varepsilon$ goes to zero. Similarly, we can apply Lemma 2.4 to show that there are positive constants $C_M^0$, $C_M^s$ not depending on $\varepsilon$ such that for any $t > 0$, one has $M_i^{\varepsilon}(t) < C_M^0$, $M_i^{s,\varepsilon}(t) < C_M^s$.

Next, we prove that the sequence of derivative $\{\dot{\mathbf{u}}^{\varepsilon}\}_{\varepsilon \in \text{a uniform bounded. For any } t > 0}$,

$$ \dot{F}_1^{\varepsilon} < r \nu_E K_1 - \frac{1}{\varepsilon} (F_1^{\varepsilon} - \eta F_2^{\varepsilon}). $$

We show that $\frac{F_1^{\varepsilon} - \eta F_2^{\varepsilon}}{\varepsilon}$ is uniformly bounded on $[0, T]$. Indeed, we have for all $t > 0$ and $\varepsilon > 0$,

$$ \begin{aligned} \dot{F}_1^{\varepsilon} - \eta \dot{F}_2^{\varepsilon} &= r \nu_E \left( E_1^{\varepsilon} \frac{M_1^{\varepsilon}}{M_1^{\varepsilon} + \gamma M_1^{s,\varepsilon}} - \eta E_2^{\varepsilon} \frac{M_2^{\varepsilon}}{M_2^{\varepsilon} + \gamma M_2^{s,\varepsilon}} \right) - \left( \mu_F + \frac{\eta+1}{\varepsilon} \right) (F_1^{\varepsilon} - \eta F_2^{\varepsilon}) \\ &= A^{\varepsilon} - \left( \frac{\eta+1}{\varepsilon} \right) (F_1^{\varepsilon} - \eta F_2^{\varepsilon}), \end{aligned} $$
<figure><img src="image_3.png" /><figcaption>Figure 2: The case when the positive root disappears</figcaption></figure>

### 4.1.2 Step 2: Shaping the release function

We now want to choose $\Lambda$ such that the condition (21) holds, which means

$$ \gamma M_i^s(t) \geq (\mathcal{N} - 1) M_i(t), \quad t \geq 0, i = 1, 2 $$

The upper bound of $M_i$ can be obtained from Lemma 2.4. For time $t > 0$ large enough, it is sufficient to choose $\Lambda$ such that $\gamma M_i^s(t) \geq (\mathcal{N} - 1) C_M$.

## 4.2 Proof of Theorem 4.2

*Proof of Theorem 4.2. Constant release:* For $\Lambda(t) \equiv \Lambda$, let us recall $A_s = \begin{pmatrix} -\mu_s - \alpha d_{12} & \alpha d_{21} \\ \alpha d_{12} & -\mu_s - \alpha d_{21} \end{pmatrix}$, and it is a Hurwitz matrix. And we have

$$ \begin{pmatrix} M_1^s(t) \\ M_2^s(t) \end{pmatrix} = e^{tA_s} \begin{pmatrix} M_1^{s0} \\ M_2^{s0} \end{pmatrix} + (1-r)v_E \int_0^t e^{(t-s)A_s} ds \begin{pmatrix} \Lambda \\ 0 \end{pmatrix}. $$

Thus, when $t \to +\infty$, we can deduce that $\begin{pmatrix} M_1^s(t) \\ M_2^s(t) \end{pmatrix}$ converges to $\begin{pmatrix} M_1^{s*} \\ M_2^{s*} \end{pmatrix}$ for any initial data since $e^{tA_s} \to 0$ when $t \to +\infty$. Hence, for any $\varepsilon > 0$, there exists a value $T_\varepsilon > 0$ such that for any $t > T_\varepsilon$, $i = 1, 2$,

$$ M_i^s(t) \geq M_i^{s*} - \varepsilon. $$

If we take $\Lambda$ such that $\gamma(M_i^{s*} - \varepsilon) \geq (\mathcal{N} - 1) C_M$ with $C_M$ defined in Lemma 2.4, then condition (21) holds. By applying Proposition 4.3, we deduce that for $i = 1, 2$, if $\Lambda > \max_{i=1,2} \frac{1}{\gamma\tau_i}(\mathcal{N} - 1)C_M$, system (1) has $\mathbf{u}_0^*$ as a unique equilibrium point, and every trajectory converges to this equilibrium when $t \to +\infty$. The dynamics of system (1) depend continuously and monotonically on $\Lambda$, then we deduce that there exists a positive critical value $\bar{\Lambda} \leq \max_{i=1,2} \frac{1}{\gamma\tau_i}(\mathcal{N} - 1)C_M$ such that for any $\Lambda > \bar{\Lambda}$, and for any non-negative initial data, solution $\mathbf{u}' = (E_1, F_1, M_1, E_2, F_2, M_2)(t)$ converges to $\mathbf{0}_6$ when $t \to +\infty$.

**Impulsive periodic releases:** Consider $\Lambda(t)$ defined in (19), denote $(E_1^{\text{per}}, F_1^{\text{per}}, M_1^{\text{per}}, E_2^{\text{per}}, F_2^{\text{per}}, M_2^{\text{per}})$ a solution of (1a)-(1c), (1e)-(1g) with $M_i^s \equiv M_i^{s,\text{per}}$ defined in (20). From Lemma 4.1, one has $M_i^{s,\text{per}}(t) \geq \Lambda^{\text{per}}\tau_i^{\text{per}}$ for all $t > 0$. Therefore, if we take $\Lambda^{\text{per}}$ such that $\gamma\Lambda^{\text{per}}\tau_i^{\text{per}} \geq (\mathcal{N} - 1)C_M$, then condition
where $A^e := r\nu_E \left( E_1^e \frac{M_1^e}{M_1^e + \gamma M_1^{s,e}} - \eta E_2^e \frac{M_2^e}{M_2^e + \gamma M_2^{s,e}} \right) - \mu_F (F_1^e - \eta F_2^e)$ is uniformly bounded since we already proved that $\mathbf{u}^e$ is uniformly bounded. Then, for any $\varepsilon > 0$, we have $|A^e(t)| < C$ for any $t > 0$ and some constant $C > 0$. By the Duhamel formula, we obtain

$$ (F_1^e - \eta F_2^e)(t) = (F_1^{e,0} - \eta F_2^{e,0})e^{-\frac{\eta+1}{\varepsilon}t} + \int_0^t A^e(s)e^{-\frac{\eta+1}{\varepsilon}(t-s)}ds. $$

So for all $t \ge 0$, one has

$$ \frac{|F_1^e - \eta F_2^e|(t)}{\varepsilon} \le \frac{|F_1^{e,0} - \eta F_2^{e,0}|}{\varepsilon} e^{-\frac{\eta+1}{\varepsilon}t} + \frac{C}{\eta+1} \left(1 - e^{-\frac{\eta+1}{\varepsilon}t}\right). $$

For any $t \in [0, +\infty)$ and $\varepsilon > 0$, one has $0 < e^{-\frac{\eta+1}{\varepsilon}t} < 1$. And due to the Assumption (26) for the initial data, the right-hand side is uniformly bounded with respect to $\varepsilon$. Hence, we deduce that $\dot{F}_i^e$ is uniformly bounded on $[0, T]$. We obtain analogously the uniform boundedness of $\dot{F}_i^e$, $M_i^e$, and $M_i^{s,e}$. Due to the positivity of system (1), one has $\dot{E}_i^e(t) < bC_{F_i}^0$ for all $t > 0$ and $\varepsilon > 0$.

Since the sequence of derivatives $\{\mathbf{u}^e\}_\varepsilon$ is uniformly bounded on $[0, T]$, we deduce the equicontinuity of the sequence $\{\mathbf{u}^e\}_\varepsilon$. Hence, by the Arzela-Ascoli theorem, this sequence has a uniformly convergent subsequence. We denote its limit $\mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s)$. If we multiply system (1) with $\varepsilon$ and let it go to zero, we obtain the equalities (27) and system (28).

With the initial data satisfying the assumptions in Proposition 5.1, the solution of system (28) on $(0, +\infty)$ is unique. Since all the subsequence of $\{\mathbf{u}^e\}_\varepsilon$ converge to the same limit, we deduce that the whole sequence converges uniformly to this limit on $[0, T]$.

• Extension to $+\infty$: For all $t \ge 0$, we prove that for all $\delta > 0$, there exists $\varepsilon_0 > 0$ such that for all $\varepsilon \in (0, \varepsilon_0)$, we have $\|\mathbf{u}^e(t) - \mathbf{u}(t)\| < \delta$.

Indeed, the solution of both (1) and (28) converges to a constant as time $t$ goes to infinity, then there exists a time $T > 0$ large enough and $\varepsilon_1 > 0$ such that for any $\varepsilon \in (0, \varepsilon_1)$ and all $t > T$, one has

$$ \|\mathbf{u}^e(t) - \mathbf{u}^e(T)\| < \frac{\delta}{3}, \quad \|\mathbf{u}(t) - \mathbf{u}(T)\| < \frac{\delta}{3}. $$

Moreover, we have that the sequence $\{\mathbf{u}^e\}_\varepsilon$ converges uniformly to $\mathbf{u}$ in the closed interval $[0, T]$. Thus, there exists a positive value $\varepsilon_0 < \varepsilon_1$ such that for all $\varepsilon \in (0, \varepsilon_0)$, $\sup_{[0,T]} \|\mathbf{u}^e - \mathbf{u}\| < \frac{\delta}{3}$. Hence, we have $\|\mathbf{u}^e(T) - \mathbf{u}(T)\| < \frac{\delta}{3}$ and we deduce that

$$ \|\mathbf{u}^e(t) - \mathbf{u}(t)\| \le \|\mathbf{u}^e(t) - \mathbf{u}^e(T)\| + \|\mathbf{u}^e(T) - \mathbf{u}(T)\| + \|\mathbf{u}(t) - \mathbf{u}(T)\| < \delta. $$

It is clear that for $t \le T$, one has $\|\mathbf{u}^e(t) - \mathbf{u}(t)\| < \frac{\delta}{3} < \delta$. So we obtain the convergence on $[0, +\infty)$. $\square$

In the next result, we study the limit system (28).

**Theorem 5.2.** Consider system (28) with the release function given by

(i) constant release $\Lambda_\infty(t) \equiv \Lambda_\infty$.

Then there exists $\bar{\Lambda}_\infty > 0$ such that for any $\Lambda_\infty > \bar{\Lambda}_\infty$, system (28) has a unique equilibrium $u_\infty^0 = (0, 0, 0, 0, \frac{\Lambda_\infty}{\mu_\infty})$ and it is globally asymptotically stable.

(ii) impulsive periodic release $\Lambda_\infty(t) = \sum_{k=0}^{+\infty} \tau \Lambda_\infty^{\text{per}} \delta_{k\tau}$ with period $\tau$.

Then there exists $\bar{\Lambda}_\infty^{\text{per}} > 0$ such that for any $\Lambda_\infty^{\text{per}} > \bar{\Lambda}_\infty^{\text{per}}$, all trajectories of (28) resulting from any non-negative initial data satisfy that $(E_1, E_2, F, M)$ converges to the equilibrium $\mathbf{0}_4 \in \mathbb{R}^4$.

*Proof.* Firstly, we show that for $\Lambda_\infty(t)$ large enough such that $\frac{M}{M + \gamma M^s} \le \frac{1}{N}$, then all trajectories of (28) resulting from any non-negative initial data satisfy that $(E_1, E_2, F, M)$ converges to the equilibrium $\mathbf{0}_4 \in \mathbb{R}^4$. Indeed, consider the first four equations of system (28) with $\frac{M}{M + \gamma M^s}$ replaced by $\frac{1}{N}$, and we denote the equilibrium $(E_1^*, E_2^*, F^*, M^*)$ of this system satisfy

$$ F^* = \frac{r\nu_E(E_1^* + E_2^*)}{N\mu_F}, \quad M^* = \frac{(1-r)\nu_E(E_1^* + E_2^*)}{\mu_M}, $$
(21) holds. By applying Proposition 4.3, we deduce that $\mathbf{u}^{\text{per}} = (E_1^{\text{per}}, F_1^{\text{per}}, M_1^{\text{per}}, E_2^{\text{per}}, F_2^{\text{per}}, M_2^{\text{per}})$ converges to zero as $t$ grows. Since $(M_1^s, M_2^s)$ converges to $(M_1^{s,\text{per}}, M_2^{s,\text{per}})$ as $t \to +\infty$, we have $\mathbf{u}' = (E_1, F_1, M_1, E_2, F_2, M_2)$ approaches $\mathbf{u}^{\text{per}}$ and thus converges to $\mathbf{0}_6$ as time $t$ goes to infinity.

Since the dynamics of system (1) depends continuously and monotonically on $\Lambda$, there exists a positive critical value $\bar{\Lambda}^{\text{per}} \le \max_{i=1,2} \frac{1}{\gamma\tau_i^{\text{per}}}(\mathcal{N}-1)C_M$ such that if $\Lambda^{\text{per}} > \bar{\Lambda}^{\text{per}}$ the equilibrium $\mathbf{u}_0^{\text{per}}$ of (1) with $\Lambda(t)$ defined in (18) is globally asymptotically stable. $\square$

# 5 Parameter dependence of the critical values of the release rate

In this section, we consider the constant release case and examine how the critical value $\bar{\Lambda}$ depends on the parameters of system (1). In this model, the elimination of the population depends not only on the diffusion rate between the inaccessible area and the treated area, but also on the biological intrinsic values like the birth/death rates, and the carrying capacity.

## 5.1 Diffusion rates

In this part, we want to compare the critical values of $\Lambda$ corresponding to different values of $d_{12}$, $d_{21}$. We show that when the diffusion rates are large enough, the critical number of sterile males released is the same as in the case when there is no separation between the two sub-populations.

### 5.1.1 The case $d_{12}, d_{21}$ large

First, we present a result of uniform convergence of system (1) when $d_{12}, d_{21}$ go to $+\infty$ and $d_{12}$ is proportional to $d_{21}$.

**Proposition 5.1.** For $\varepsilon > 0$, consider the diffusion rates $d_{12} = \frac{1}{\varepsilon}$, $d_{21} = \frac{\eta}{\varepsilon}$ with $\eta = \frac{d_{21}}{d_{12}} > 0$. Denote $\mathbf{u}^{\varepsilon} = (E_1^{\varepsilon}, F_1^{\varepsilon}, M_1^{\varepsilon}, M_1^{s,\varepsilon}, E_2^{\varepsilon}, F_2^{\varepsilon}, M_2^{\varepsilon}, M_2^{s,\varepsilon})$ the solution of system (1) with the initial date $\mathbf{u}^{0,0}$ satisfying that

$\{E_i^{\varepsilon,0}\}_\varepsilon, \{F_i^{\varepsilon,0}\}_\varepsilon, \{M_i^{\varepsilon,0}\}_\varepsilon$ converge to $E_i^{0,0}, F_i^{0,0}, M_i^{0,0}$ respectively as $\varepsilon \to 0$, with $i = 1, 2$,

and

$$E_1^{\varepsilon,0} - \eta E_2^{\varepsilon,0} = \mathcal{O}(\varepsilon), \quad F_1^{\varepsilon,0} - \eta F_2^{\varepsilon,0} = \mathcal{O}(\varepsilon), \quad M_1^{\varepsilon,0} - \eta M_2^{\varepsilon,0} = \mathcal{O}(\varepsilon), \quad M_1^{s,\varepsilon,0} = M_2^{s,\varepsilon,0} = 0. \quad (25)$$

Then, when $\varepsilon \to 0$, the sequence $\{\mathbf{u}^\varepsilon\}_\varepsilon$ converges uniformly to a limit $(E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s)$ on $[0, +\infty)$. Moreover, we have

$$F_1 = \eta F_2, \quad M_1 = \eta M_2, \quad M_1^s = \eta M_2^s. \qquad (26)$$

If we denote $F = F_1 + F_2$, $M = M_1 + M_2$, $M^s = M_1^s + M_2^s$, then $(E_1, E_2, F, M, M^s)$ solves the following system

$$\dot{E}_1 = \frac{\eta}{\eta+1} bF \left(1 - \frac{E_1}{K_1}\right) - (v_E + \mu_E)E_1, \qquad (27a)$$

$$\dot{E}_2 = \frac{1}{\eta+1} bF \left(1 - \frac{E_2}{K_2}\right) - (v_E + \mu_E)E_2, \qquad (27b)$$
and

$$ (E_1^* + E_2^*) = \frac{E_1^* + E_2^*}{\frac{\eta}{\eta+1} \left(1 - \frac{E_1^*}{K_1}\right) + \frac{1}{\eta+1} \left(1 - \frac{E_2^*}{K_2}\right)}. $$

This is equivalent to either $E_1^* + E_2^* = 0$ or

$$ \frac{\eta}{\eta+1} \left(1 - \frac{E_1^*}{K_1}\right) + \frac{1}{\eta+1} \left(1 - \frac{E_2^*}{K_2}\right) = 1. $$

The left-hand side of this equality is smaller than 1 since $E_i^* < K_i$ with $i = 1, 2$, thus we deduce that $E_1^* = E_2^* = F^* = M^* = 0$. Hence, this system has exactly one equilibrium $\mathbf{0}_4$ and all trajectories converge to this steady state by using Theorem 3.1 in Chapter 2 of [39]. Then, by applying the comparison Lemma 2.2, we deduce the convergence of system (28).

Analogously to system (1), we have the boundedness for the solution of (28) and the monotonicity of the system with respect to $\Lambda_\infty$. Therefore, we can deduce the existence of the critical values for both the constant and periodic cases. $\square$

Next, we make a comparison between the previous case and the case where there is no separation between the two sub-populations.

### 5.1.2 The non-separation case

When there is no separation between the two sub-populations of mosquitoes, we consider one population $(E, F, M, M^s)$ in a habitat with aquatic carrying capacity $K = K_1 + K_2$. Then $(E, F, M, M^s)$ satisfies the following system

$$ \dot{E} = bF \left(1 - \frac{E}{K}\right) - (\nu_E + \mu_E)E, \qquad (29a) $$

$$ \dot{F} = r\nu_E E \frac{M}{M + \gamma M^s} - \mu_F F, \qquad (29b) $$

$$ \dot{M} = (1 - r)\nu_E E - \mu_M M, \qquad (29c) $$

$$ \dot{M}^s = \Lambda - \mu_s M^s. \qquad (29d) $$

For the constant release, the positive equilibrium $(E^*, F^*, M^*, M^{s*})$ satisfies

$$ M^* = \frac{(1-r)\nu_E}{\mu_M} E^*, \quad M^{s*} = \frac{\Lambda}{\mu_s}, \quad F^* = \frac{r\nu_E}{\mu_F} \frac{E^*}{1 + \frac{\mu_M \gamma \Lambda}{(1-r)\nu_E \mu_s E^*}}; $$

and from (29a), we deduce that

$$ \frac{br\nu_E}{\mu_F} \frac{E^*}{1 + \frac{\mu_M \gamma \Lambda}{(1-r)\nu_E \mu_s E^*}} \left(1 - \frac{E^*}{K}\right) - (\nu_E + \mu_E)E^* = 0. $$

This equation has no positive solution if and only if $\Lambda > \bar{\Lambda}_0 = \frac{(1-r)\nu_E K \mu_s (1-\mathcal{N})^2}{4\mathcal{N}\mu_M \gamma}$.

**Remark 5.1.** We can see that in the special case where $K_1 = \eta K_2$, by taking $E = E_1 + E_2$, we can write system (28) as system (29) for $(E, F, M, M^s)$ with carrying capacity $K = K_1 + K_2$. Hence, we deduce that $\bar{\Lambda}_\infty = \bar{\Lambda}_0$. This suggests that the critical number of sterile males released in the case with very large diffusion rate is the same as in the non-separation case in 5.1.2.

## 5.2 Biological intrinsic values

In this section, we compare the critical value of $\Lambda$ corresponding to different values of the parameters namely the birth rate $b$, the death rates $\mu_E, \mu_F, \mu_M, \mu_s$, and the carrying capacities $K_1, K_2$. In this section, we show that the critical value $\bar{\Lambda}$ is monotone with respect to these parameters. To prove this claim, we first define in $\mathbb{R}_+^7$ an order such that $(\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \triangleq (\mu'_E, \mu'_F, \mu'_M, \mu'_s, b', K'_1, K'_2)$ if and only if

$$ \mu_E \le \mu'_E, \mu_F \le \mu'_F, \mu_M \le \mu'_M, \mu_s \ge \mu'_s, b \ge b', K_1 \ge K'_1, K_2 \ge K'_2. $$

Moreover, we write $(\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \triangleq (\mu'_E, \mu'_F, \mu'_M, \mu'_s, b', K'_1, K'_2)$ if the two vectors are not identical. With this order relation, we have the following result
$$ \dot{F} = r v_E (E_1 + E_2) \frac{M}{M + \gamma M^s} - \mu_F F, \qquad (27c) $$

$$ \dot{M} = (1 - r) v_E (E_1 + E_2) - \mu_M M, \qquad (27d) $$

$$ \dot{M}^s = \Lambda_{\infty} - \mu_s M^s, \qquad (27e) $$

with the corresponding initial data $E_1^{0,0}$, $E_2^{0,0}$, $F^{0,0} = F_1^{0,0} + F_2^{0,0}$, $M^{0,0} = M_1^{0,0} + M_2^{0,0}$, $M^{s,0,0} = 0$.

It is straightforward to see that the previous result implies that the functions $F_1, F_2, M_1, M_2$ fulfill the following identities:

$$ F_1 = \frac{\eta}{1 + \eta} F, \quad F_2 = \frac{1}{1 + \eta} F, \quad M_1 = \frac{\eta}{1 + \eta} M, \quad M_2 = \frac{1}{1 + \eta} M. $$

*Proof of Proposition 5.1.* To prove this result, we first apply the Arzela-Ascoli theorem for the sequence of smooth solution $\{\mathbf{u}^\epsilon\}_\epsilon$ on a close interval $[0, T]$ with any $T > 0$. Then, we extend the convergence at infinity.

• Uniform convergence on $[0, T]$ : First, we check the uniform boundedness of this sequence. For $i = 1, 2$, from Lemma 2.4, one has $E_i^\epsilon(t) \le K_i$ for all $t > 0$ and $\epsilon > 0$. Again by this Lemma, for any $t > 0$, one has

$$ F_i^\epsilon(t) \le \max \left\{ F_1^{\epsilon,0} + F_2^{\epsilon,0}, \frac{r v_E (K_1 + K_2)}{\mu_F} \right\} \le C_F^0 $$

where $C_F^0$ does not depend on $\epsilon$ since the initial data converge as $\epsilon$ goes to zero. Similarly, we can apply Lemma 2.4 to show that there are positive constants $C_M^0, C_M^0$ not depending on $\epsilon$ such that for any $t > 0$, one has $M_i^\epsilon(t) < C_M^0, M_i^{s,\epsilon}(t) < C_M^0$.

Next, we prove that the sequence of derivative $\{\mathbf{u}^\epsilon\}_\epsilon$ is also uniformly bounded. For any $t > 0$,

$$ \dot{F}_1^\epsilon < r v_E K_1 - \frac{1}{\epsilon} (F_1^\epsilon - \eta F_2^\epsilon). $$

We show that $\frac{F_1^\epsilon - \eta F_2^\epsilon}{\epsilon}$ is uniformly bounded on $[0, T]$. Indeed, we have for all $t > 0$ and $\epsilon > 0$,

$$ \begin{aligned} \dot{F}_1^\epsilon - \eta \dot{F}_2^\epsilon &= r v_E \left( E_1^\epsilon \frac{M_1^\epsilon}{M_1^\epsilon + \gamma M_1^{s,\epsilon}} - \eta E_2^\epsilon \frac{M_2^\epsilon}{M_2^\epsilon + \gamma M_2^{s,\epsilon}} \right) - \left( \mu_F + \frac{\eta+1}{\epsilon} \right) (F_1^\epsilon - \eta F_2^\epsilon) \\ &= A^\epsilon - \left( \frac{\eta+1}{\epsilon} \right) (F_1^\epsilon - \eta F_2^\epsilon), \end{aligned} $$

where $A^\epsilon := r v_E \left( E_1^\epsilon \frac{M_1^\epsilon}{M_1^\epsilon + \gamma M_1^{s,\epsilon}} - \eta E_2^\epsilon \frac{M_2^\epsilon}{M_2^\epsilon + \gamma M_2^{s,\epsilon}} \right) - \mu_F (F_1^\epsilon - \eta F_2^\epsilon)$ is uniformly bounded since we already proved that $\mathbf{u}^\epsilon$ is uniformly bounded. Then, for any $\epsilon > 0$, we have $|A^\epsilon(t)| < C$ for any $t > 0$ and some constant $C > 0$. By the Duhamel formula, we obtain

$$ (F_1^\epsilon - \eta F_2^\epsilon)(t) = (F_1^{\epsilon,0} - \eta F_2^{\epsilon,0})e^{-\frac{\eta+1}{\epsilon}t} + \int_0^t A^\epsilon(s)e^{-\frac{\eta+1}{\epsilon}(t-s)}ds. $$

So for all $t \ge 0$, one has

$$ \frac{|F_1^\epsilon - \eta F_2^\epsilon|(t)}{\epsilon} \le \frac{|F_1^{\epsilon,0} - \eta F_2^{\epsilon,0}|}{\epsilon} e^{-\frac{\eta+1}{\epsilon}t} + \frac{C}{\eta+1} \left(1 - e^{-\frac{\eta+1}{\epsilon}t}\right). $$

For any $t \in [0, +\infty)$ and $\epsilon > 0$, one has $0 < e^{-\frac{\eta+1}{\epsilon}t} < 1$. And due to the Assumption (25) for the initial data, the right-hand side is uniformly bounded with respect to $\epsilon$. Hence, we deduce that $\dot{F}_1^\epsilon$ is uniformly
<table><caption>Table 1: Parameter values of *Aedes albopictus* mosquitoes used for the numerical simulation</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>b</td><td>Birth rate of fertile females</td><td>10</td><td>day<sup>-1</sup></td></tr><tr><td>ν<sub>E</sub></td><td>Emerging rate of viable eggs</td><td>0.08</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>E</sub></td><td>Death rate of aquatic phase</td><td>0.05</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>F</sub></td><td>Female death rate</td><td>0.1</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>M</sub></td><td>Wild male death rate</td><td>0.14</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>s</sub></td><td>Sterile male death rate</td><td>0.14</td><td>day<sup>-1</sup></td></tr><tr><td>K<sub>1</sub></td><td>Carrying capacity of aquatic phase in patch 1</td><td>200</td><td>-</td></tr><tr><td>K<sub>2</sub></td><td>Carrying capacity of aquatic phase in patch 2</td><td>180</td><td>-</td></tr><tr><td>γ</td><td>Mating competitiveness of sterile male</td><td>1</td><td>-</td></tr><tr><td>r</td><td>Ratio of female hatch</td><td>0.5</td><td>-</td></tr><tr><td>α</td><td>Ratio between diffusion rates of sterile males and female</td><td>0.5</td><td>-</td></tr><tr><td>β</td><td>Ratio between diffusion rates of sterile males and female</td><td>0.8</td><td>-</td></tr></tbody></table>

**Theorem 5.3.** Consider system (1) and the basic offspring number $N > 1$, consider the critical values $\bar{\Lambda}$ and $\bar{\Lambda}^{\text{per}}$ as defined in Theorem 4.2, then we have the mappings from $\mathbb{R}_+^7$ to $\mathbb{R}_+$

$$ (\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \mapsto \bar{\Lambda}, \quad (\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \mapsto \bar{\Lambda}^{\text{per}}, $$

are non-increasing with respect to the order $\leq$.

*Proof.* First, we consider system (1) with two sets of parameters

$$ \Theta = (\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2), \quad \Theta' = (\mu'_E, \mu'_F, \mu'_M, \mu'_s, b', K'_1, K'_2), $$

where $\Theta \leq \Theta'$. We fix the same value of $\Lambda$ in both cases and consider

$$ \mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s), \quad \mathbf{v} = (\tilde{E}_1, \tilde{F}_1, \tilde{M}_1, \tilde{M}_1^s, \tilde{E}_2, \tilde{F}_2, \tilde{M}_2, \tilde{M}_2^s) $$

where $\mathbf{u}, \mathbf{v}$ are the solutions of (1) using the same initial data with the parameters $\Theta, \Theta'$, respectively. We have $\dot{\mathbf{u}} = \mathbf{f}_{\Theta}(\mathbf{u})$, and $\dot{\mathbf{v}} = \mathbf{f}_{\Theta'}(\mathbf{v}) \leq \mathbf{f}_{\Theta}(\mathbf{v})$ in the subset $\{0 \leq E_1 \leq K_1\} \cap \{0 \leq E_2 \leq K_2\}$ of $\mathbb{R}_+^8$. Moreover, functions $\mathbf{f}_{\Theta}$ and $\mathbf{f}_{\Theta'}$ satisfy the assumptions in Lemma 2.3, then by applying this lemma, we obtain that $\mathbf{v} \leq \mathbf{u}$ for the same initial data, so

$$ E_i(t) \geq \tilde{E}_i(t), \quad F_i(t) \geq \tilde{F}_i(t), \quad M_i(t) \geq \tilde{M}_i(t) \quad \text{for all } t > 0, i = 1, 2. $$

On the other hand, for any $\Lambda > \bar{\Lambda}_{\Theta}$, by Theorem 4.2 we have that $E_i(t), F_i(t), M_i(t)$ converge to zero as $t$ goes to infinity. As a consequence of the above inequalities, we deduce that $\tilde{E}_i(t), \tilde{F}_i(t), \tilde{M}_i(t)$ also converge to zero for all initial data. So $\Lambda > \bar{\Lambda}_{\Theta'}$, and we can deduce that $\bar{\Lambda}_{\Theta} \geq \bar{\Lambda}_{\Theta'}$. $\square$

# 6 Numerical simulations

Following [13, 22], we consider the parameters as in Table 1.

## 6.1 Trajectories and Equilibria

We fix the moving rate $d_{12} = 0.06$, $d_{21} = 0.04$ (day$^{-1}$), and plot the numerical solutions of system (1) with different release functions $\Lambda(t)$. In the case $\Lambda = 0$, system (7) has a unique equilibrium

$$ (E_1^+, F_1^+, M_1^+, E_2^+, F_2^+, M_2^+) = (192.63, 67.92, 50.29, 174.83, 79.06, 54.69). $$

To highlight the global stability of equilibria, we numerically solve the system with initial data in three levels: close to zero, intermediate, and close to the positive equilibrium when $\Lambda = 0$. For the level close to zero, the initial data is taken between 1% to 10% of the value of the above equilibrium, between 40% and 50% for the intermediate level, and between 90% to 100% for the third case. More precisely, we take

$$ (E_1^0, F_1^0, M_1^0, M_1^{s,0}, E_2^0, F_2^0, M_2^0, M_2^{s,0}) \in \{(2, 5, 4, 0, 3, 5, 3, 0), (80, 30, 20, 0, 70, 30, 30, 0), (160, 60, 50, 0, 155, 70, 50, 0)\}. $$

When the largest wild mosquito population was less than $10^{-2}$, we considered the wild population to be eliminated. The following section presents several numerical simulations showing the trajectories and approximated equilibria according to different release strategies.
bounded on $[0, T]$. We obtain analogously the uniform boundedness of $\dot{E}_i^\epsilon$, $\dot{M}_i^\epsilon$, and $\dot{M}_i^{\epsilon, \epsilon}$. Due to the positivity of system (1), one has $\dot{E}_i^\epsilon(t) < bC_{F_i}^0$ for all $t > 0$ and $\epsilon > 0$.

Since the sequence of derivatives $\{\dot{\mathbf{u}}^\epsilon\}_\epsilon$ is uniformly bounded on $[0, T]$, we deduce the equicontinuity of the sequence $\{\mathbf{u}^\epsilon\}_\epsilon$. Hence, by the Arzela-Ascoli theorem, this sequence has a uniformly convergent subsequence. We denote its limit $\mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s)$. If we multiply system (1) with $\epsilon$ and let it go to zero, we obtain the equalities (26) and system (27).

With the initial data satisfying the assumptions in Proposition 5.1, the solution of system (27) on $(0, +\infty)$ is unique. Since all the subsequence of $\{\mathbf{u}^\epsilon\}_\epsilon$ converge to the same limit, we deduce that the whole sequence converges uniformly to this limit on $[0, T]$.

• Extension to $+\infty$: For all $t \ge 0$, we prove that for all $\delta > 0$, there exists $\epsilon_0 > 0$ such that for all $\epsilon \in (0, \epsilon_0)$, we have $\|\mathbf{u}^\epsilon(t) - \mathbf{u}(t)\| < \delta$.

Indeed, the solution of both (1) and (27) converges to a constant as time $t$ goes to infinity, then there exists a time $T > 0$ large enough and $\epsilon_1 > 0$ such that for any $\epsilon \in (0, \epsilon_1)$ and all $t > T$, one has

$$ \|\mathbf{u}^\epsilon(t) - \mathbf{u}^\epsilon(T)\| < \frac{\delta}{3}, \quad \|\mathbf{u}(t) - \mathbf{u}(T)\| < \frac{\delta}{3}. $$

Moreover, we have that the sequence $\{\mathbf{u}^\epsilon\}_\epsilon$ converges uniformly to $\mathbf{u}$ in the closed interval $[0, T]$. Thus, there exists a positive value $\epsilon_0 < \epsilon_1$ such that for all $\epsilon \in (0, \epsilon_0)$,

$$ \sup_{[0,T]} \|\mathbf{u}^\epsilon - \mathbf{u}\| < \frac{\delta}{3}. $$

Hence, we have $\|\mathbf{u}^\epsilon(T) - \mathbf{u}(T)\| < \frac{\delta}{3}$ and we deduce that

$$ \|\mathbf{u}^\epsilon(t) - \mathbf{u}(t)\| \le \|\mathbf{u}^\epsilon(t) - \mathbf{u}^\epsilon(T)\| + \|\mathbf{u}^\epsilon(T) - \mathbf{u}(T)\| + \|\mathbf{u}(t) - \mathbf{u}(T)\| < \delta. $$

It is clear that for $t \le T$, one has $\|\mathbf{u}^\epsilon(t) - \mathbf{u}(t)\| < \frac{\delta}{3} < \delta$. So we obtain the convergence on $[0, +\infty)$. $\square$

In the next result, we study the limit system (27).

**Theorem 5.2.** Consider system (27) with the release function given by

(i) constant release $\Lambda_\infty(t) \equiv \Lambda_\infty$.

Then there exists $\bar{\Lambda}_\infty > 0$ such that for any $\Lambda_\infty > \bar{\Lambda}_\infty$, system (27) has a unique equilibrium $u_0^0 = (0, 0, 0, 0, \frac{\Lambda_\infty}{\mu_s})$ and it is globally asymptotically stable.

(ii) impulsive periodic release $\Lambda_\infty(t) = \sum_{k=0}^{+\infty} \tau \Lambda_\infty^{\text{per}} \delta_{k\tau}$ with period $\tau$.

Then there exists $\bar{\Lambda}_\infty^{\text{per}} > 0$ such that for any $\Lambda_\infty^{\text{per}} > \bar{\Lambda}_\infty^{\text{per}}$, all trajectories of (27) resulting from any non-negative initial data satisfy that $(E_1, E_2, F, M)$ converges to the equilibrium $\mathbf{0}_4 \in \mathbb{R}^4$.

*Proof.* Firstly, we show that for $\Lambda_\infty(t)$ large enough such that $\frac{M}{M + \gamma M^s} \le \frac{1}{\mathcal{N}}$, then all trajectories of (27) resulting from any non-negative initial data satisfy that $(E_1, E_2, F, M)$ converges to the equilibrium $\mathbf{0}_4 \in \mathbb{R}^4$. Indeed, consider the first four equations of system (27) with $\frac{M}{M + \gamma M^s}$ replaced by $\frac{1}{\mathcal{N}}$, and we denote the equilibrium $(E_1^*, E_2^*, F^*, M^*)$ of this system satisfy

$$ F^* = \frac{r v_E (E_1^* + E_2^*)}{\mathcal{N} \mu_F}, \quad M^* = \frac{(1 - r) v_E (E_1^* + E_2^*)}{\mu_M}, $$
### 6.1.1 Constant continuous releases

We take three different constant values of $\Lambda \in \{0, 200, 500\}$ (day$^{-1}$). The initial density of sterile males is equal to zero. We approximate the positive equilibria in each case and plot the trajectories of $E_1$ and $E_2$ in Figures 3 according to different values of $\Lambda$. We observe the following:

* When $\Lambda = 0$, there is one positive equilibrium

$$ (E_1^*, E_2^*) = (192.62, 174.82). $$

All positive trajectories converge to the positive steady state $(E_1^*, E_2^*)$.

* When $\Lambda = 200$ (day$^{-1}$), there are two positive equilibria

$$ (e_1^*, e_2^*) = (17.29, 49.98), \quad (E_1^*, E_2^*) = (85.79, 130.02). $$

All positive trajectories also converge to the larger positive steady state $(E_1^*, E_2^*)$.

* When $\Lambda = 500$ (day$^{-1}$), there is no positive equilibrium. All the trajectories converge to the zero equilibrium.

This validates the result in Theorem 4.2 that when $\Lambda$ exceeds some critical value, zero is the unique equilibrium of system (1). The observation for $\Lambda = 0$ illustrates the result in Theorem 3.1 that there is one positive equilibrium and it is globally asymptotically stable. The introduction of sterile males ($\Lambda = 200 > 0$) reduces the value of the positive steady state (see Figure 3b), and when $\Lambda = 500$ (day$^{-1}$) exceeds some critical value (at most equal to 500), all trajectories converge to the zero equilibrium (see Figure 3c). This illustrates the first point of Theorem 4.2. To approximate the critical value of $\Lambda$, we provide some numerical bifurcation diagrams in Section 6.2.

### 6.1.2 Periodic impulsive releases

In this part, we consider the periodic impulsive releases with $\Lambda(t)$ defined in (19), with $\Lambda^{\text{per}}$ equal to 200 and 300 (day$^{-1}$), the period $\tau = 10$ (days). The trajectories of $E_1, E_2$ shown in Figure 4 converge to the periodic solution when $\Lambda^{\text{per}} = 200$ (day$^{-1}$) and go to zero when $\Lambda^{\text{per}} = 300$ (day$^{-1}$). This illustrates the second point of Theorem 4.2 that when the number of sterile males released exceeds a critical value $\bar{\Lambda}^{\text{per}}$, the wild populations of mosquitoes in both areas reach elimination.

## 6.2 Critical values and bifurcation

Our aim in this section is to approximate the critical value of $\Lambda$ where the bifurcation occurs.

### 6.2.1 Bifurcation diagram in the constant release case

We solve a system of nonlinear stationary problems $\mathcal{F}(u; \Lambda) = 0$ for all values of the parameter $\Lambda$, knowing that the solutions are continuous with respect to $\Lambda$. Solving by numerical approximations can be done using numerical continuation methods (see [40]).

Here we present the simplest method called *Natural Parameter Continuation* (incremental methods, see [40]): Iteratively find approximate roots of $\mathcal{F}(u; \Lambda) = 0$ for several values of $\Lambda_i$ with index $i \in \mathbb{N}^*$. The root of step $i$ is used as an initial guess for the numerical solver at step $i + 1$. The first initial guess is the root for the smallest $\Lambda$. To approximate the critical value $\bar{\Lambda}$ in the constant case and examine what happens when $0 < \Lambda \le \bar{\Lambda}$, we draw the bifurcation diagram for $\Lambda \in [0.1, 500]$. The initial positions of the numerical continuation are taken at the approximated equilibria when $\Lambda = 0.1$.

We obtain the bifurcation diagrams in Figure 5 for two scenarios. We observed that the critical value of $\Lambda$ decreases when the diffusion rates increase.

* For $d_{12} = 1$, $d_{21} = 2$, the critical value $\bar{\Lambda} = 106.45$ (day$^{-1}$).

* For $d_{12} = 0.06$, $d_{21} = 0.04$, the critical value $\bar{\Lambda} = 250.88$ (day$^{-1}$).

Values of $\bar{\Lambda}$ varying with respect to different values of $d_{12}$, $d_{21} \in [0.05, 2]$ were depicted in Figure 6a. We observe that as the diffusion rates get larger, the critical value $\bar{\Lambda}$ is decreasing and converges to some value $\bar{\Lambda}_{\infty}$. This validates the result provided by Proposition 5.2 where $\bar{\Lambda}_{\infty}$ is the critical value of $\Lambda$ corresponding to system (28).

By considering the case $d_{21} = \eta d_{12}$, we plot values of $\bar{\Lambda}$ corresponding to different ratios $\eta$ while fixing $d_{12} = 0.6$ in Figure 6b. We can see that $\bar{\Lambda}$ reach a local maximum as $\eta$ passes through the ratio $\frac{K_1}{K_2}$. In practice, since the mosquitoes will likely move to areas with more breeding sites, the ratio $\eta = \frac{d_{21}}{d_{12}}$ is close to the carrying capacity $\frac{K_1}{K_2}$. This indicates that intervening on the breeding sites to increase the difference between these two ratios can help decrease the number of sterile males needing release.

Moreover, as $d_{12}$ get larger, this value of $\bar{\Lambda}$ gets closer to the critical value $\bar{\Lambda}_0$ of the system when there is no separation between the two sub-populations defined in 5.1.2.
and

$$ (E_1^* + E_2^*) = \frac{E_1^* + E_2^*}{\frac{\eta}{\eta+1} \left(1 - \frac{E_1^*}{K_1}\right) + \frac{1}{\eta+1} \left(1 - \frac{E_2^*}{K_2}\right)}. $$

This is equivalent to either $E_1^* + E_2^* = 0$ or

$$ \frac{\eta}{\eta+1} \left(1 - \frac{E_1^*}{K_1}\right) + \frac{1}{\eta+1} \left(1 - \frac{E_2^*}{K_2}\right) = 1. $$

The left-hand side of this equality is smaller than 1 since $E_i^* < K_i$ with $i = 1, 2$, thus we deduce that $E_1^* = E_2^* = F^* = M^* = 0$. Hence, this system has exactly one equilibrium $0_4$ and all trajectories converge to this steady state by using Theorem 3.1 in Chapter 2 of [33]. Then, by applying the comparison Lemma 2.2, we deduce the convergence of system (27).

Analogously to system (1), we have the boundedness for the solution of (27) and the monotonicity of the system with respect to $\Lambda_\infty$. Therefore, we can deduce the existence of the critical values for both the constant and periodic cases. $\square$

Next, we make a comparison between the previous case and the case where there is no separation between the two sub-populations.

### 5.1.2 The non-separation case

When there is no separation between the two sub-populations of mosquitoes, we consider one population $(E, F, M, M^s)$ in a habitat with aquatic carrying capacity $K = K_1 + K_2$. Then $(E, F, M, M^s)$ satisfies the following system

$$ \dot{E} = bF \left(1 - \frac{E}{K}\right) - (v_E + \mu_E)E, \qquad (28a) $$

$$ \dot{F} = r v_E E \frac{M}{M + \gamma M^s} - \mu_F F, \qquad (28b) $$

$$ \dot{M} = (1 - r) v_E E - \mu_M M, \qquad (28c) $$

$$ \dot{M}^s = \Lambda - \mu_s M^s. \qquad (28d) $$

For the constant release, the positive equilibrium $(E^*, F^*, M^*, M^{s*})$ satisfies

$$ M^* = \frac{(1-r)v_E}{\mu_M} E^*, \quad M^{s*} = \frac{\Lambda}{\mu_s}, \quad F^* = \frac{r v_E}{\mu_F} \frac{E^*}{1 + \frac{\mu_M \gamma \Lambda}{(1-r) v_E \mu_s E^*}}; $$

and from (28a), we deduce that

$$ \frac{brv_E}{\mu_F} \frac{E^*}{1 + \frac{\mu_M \gamma \Lambda}{(1-r) v_E \mu_s E^*}} \left(1 - \frac{E^*}{K}\right) - (v_E + \mu_E)E^* = 0. $$

This equation has no positive solution if and only if $\Lambda > \bar{\Lambda}_0 = \frac{(1-r)v_E K \mu_s (1 - \mathcal{N})^2}{4 \mathcal{N} \mu_M \gamma}$.

**Remark 5.1.** We can see that in the special case where $K_1 = \eta K_2$, by taking $E = E_1 + E_2$, we can write system (27) as system (28) for $(E, F, M, M^s)$ with carrying capacity $K = K_1 + K_2$. Hence, we deduce that $\bar{\Lambda}_\infty = \bar{\Lambda}_0$. This suggests that the critical number of sterile males released in the case with very large diffusion rate is the same as in the non-separation case in 5.1.2.
<figure><img src="image_4.png" /><figcaption>Figure 3: Trajectories of E₁ and E₂ in the constant release case with diffusion rates $d_{12} = 0.06$, $d_{21} = 0.04$ (day⁻¹).</figcaption></figure>
## 5.2 Biological intrinsic values

In this section, we compare the critical value of $\Lambda$ corresponding to different values of the parameters namely the birth rate $b$, the death rate $\mu_E$, $\mu_F$, $\mu_M$, $\mu_s$, and the carrying capacities $K_1, K_2$. In this section, we show that the critical value $\bar{\Lambda}$ is monotone with respect to these parameters. To prove this claim, we first define in $\mathbb{R}_+^7$ an order such that $(\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \triangleleft (\mu'_E, \mu'_F, \mu'_M, \mu'_s, b', K'_1, K'_2)$ if and only if

$$ \mu_E \le \mu'_E, \mu_F \le \mu'_F, \mu_M \le \mu'_M, \mu_s \ge \mu'_s, b \ge b', K_1 \ge K'_1, K_2 \ge K'_2. $$

Moreover, we write $(\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \triangleleft (\mu'_E, \mu'_F, \mu'_M, \mu'_s, b', K'_1, K'_2)$ if the two vectors are not identical. With this order relation, we have the following result

**Theorem 5.3.** Consider system (1) and the basic offspring number $\mathcal{N} > 1$, consider the critical values $\bar{\Lambda}$ and $\bar{\Lambda}^{\text{per}}$ as defined in Theorem 4.2, then we have the mappings from $\mathbb{R}_+^7$ to $\mathbb{R}_+$

$$ (\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \mapsto \bar{\Lambda}, \quad (\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2) \mapsto \bar{\Lambda}^{\text{per}}, $$

are non-increasing with respect to the order $\triangleleft$.

*Proof.* First, we consider system (1) with two sets of parameters

$$ \Theta = (\mu_E, \mu_F, \mu_M, \mu_s, b, K_1, K_2), \quad \Theta' = (\mu'_E, \mu'_F, \mu'_M, \mu'_s, b', K'_1, K'_2), $$

where $\Theta \triangleleft \Theta'$. We fix the same value of $\Lambda$ in both cases and consider

$$ \mathbf{u} = (E_1, F_1, M_1, M_1^s, E_2, F_2, M_2, M_2^s), \quad \mathbf{v} = (\tilde{E}_1, \tilde{F}_1, \tilde{M}_1, \tilde{M}_1^s, \tilde{E}_2, \tilde{F}_2, \tilde{M}_2, \tilde{M}_2^s) $$

where $\mathbf{u}, \mathbf{v}$ are the solutions of (1) with the parameters $\Theta, \Theta'$, respectively. We have $\dot{\mathbf{u}} = \mathbf{f}_{\Theta}(\mathbf{u})$, and $\dot{\mathbf{v}} = \mathbf{f}_{\Theta'}(\mathbf{v}) \preceq \mathbf{f}_{\Theta}(\mathbf{v})$ in the subset $\{0 \le E_1 \le K_1\} \cap \{0 \le E_2 \le K_2\}$ of $\mathbb{R}_+^8$. Moreover, functions $\mathbf{f}_{\Theta}$ and $\mathbf{f}_{\Theta'}$ satisfy the assumptions in Lemma 2.3, then by applying this lemma, we obtain that $\mathbf{v} \preceq \mathbf{u}$ for the same initial data, so

$$ E_i(t) \ge \tilde{E}_i(t), F_i(t) \ge \tilde{F}_i(t), M_i(t) \ge \tilde{M}_i(t) \quad \text{for all } t > 0, i = 1, 2. $$

On the other hand, for any $\Lambda > \bar{\Lambda}_{\Theta}$, by Theorem 4.2 we have that $E_i(t), F_i(t), M_i(t)$ converge to zero as $t$ goes to infinity. As a consequence of the above inequalities, we deduce that $\tilde{E}_i(t), \tilde{F}_i(t), \tilde{M}_i(t)$ also converge to zero for all initial data. So $\Lambda > \bar{\Lambda}_{\Theta'}$, and we can deduce that $\bar{\Lambda}_{\Theta} \ge \bar{\Lambda}_{\Theta'}$. $\square$

# 6 Numerical simulations

Following [17, 35], we consider the parameters as in Table 1.

## 6.1 Trajectories and Equilibria

We fix the moving rate $d_{12} = 0.06$, $d_{21} = 0.04$ (day$^{-1}$), and plot the numerical solutions of system (1) with different release functions $\Lambda(t)$. In each case, we numerically solve the system with different initial data $(E_1^0, F_1^0, M_1^0, E_2^0, F_2^0, M_2^0) : \{(2, 5, 6, 3, 5, 6), (10, 20, 60, 25, 40, 60), (100, 50, 60, 120, 80, 60)\}$. In the following section, we present several numerical simulations showing the trajectories and approximated equilibria according to different release strategies.
<figure><img src="image_5.png" /><figcaption>Figure 4: Trajectories of $E_1$ and $E_2$ in the periodic release case with period $\tau = 10$ (days), diffusion rates $d_{12} = 0.06$, $d_{21} = 0.04$ ($\text{day}^{-1}$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Bifurcation diagrams of $E_1^*$ with parameter $\Lambda$ in the constant continuous release case.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Trajectories of E₁ and E₂ in the constant release case with diffusion rates d₁₂ = 0.06, d₂₁ = 0.04 (day⁻¹).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Values of $\bar{\Lambda}$ varies with respect to diffusion rates $d_{12}$, $d_{21} \in [0.05, 2]$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Densities of $M_1^s$ and $M_2^s$ in both cases. Left: constant continuous releases with $\Lambda = \bar{\Lambda} = 250.88$ (day$^{-1}$), Right: periodic impulsive releases with $\Lambda^{\text{per}} = \bar{\Lambda}^{\text{per}} = 255.15$ (day$^{-1}$).</figcaption></figure>
<table><caption>Table 1: Parameter values of *Aedes albopictus* mosquitoes used for the numerical simulation</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>b</td><td>Birth rate of fertile females</td><td>10</td><td>day<sup>-1</sup></td></tr><tr><td>ν<sub>E</sub></td><td>Emerging rate of viable eggs</td><td>0.08</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>E</sub></td><td>Death rate of aquatic phase</td><td>0.05</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>F</sub></td><td>Female death rate</td><td>0.1</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>M</sub></td><td>Wild male death rate</td><td>0.14</td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>s</sub></td><td>Sterile male death rate</td><td>0.14</td><td>day<sup>-1</sup></td></tr><tr><td>K<sub>1</sub></td><td>Carrying capacity of aquatic phase in patch 1</td><td>200</td><td>—</td></tr><tr><td>K<sub>2</sub></td><td>Carrying capacity of aquatic phase in patch 2</td><td>180</td><td>—</td></tr><tr><td>γ</td><td>Mating competitiveness of sterile male</td><td>1</td><td>—</td></tr><tr><td>r</td><td>Ratio of female hatch</td><td>0.5</td><td>—</td></tr><tr><td>α</td><td>Ratio between diffusion rates of sterile males and female</td><td>0.5</td><td>—</td></tr><tr><td>β</td><td>Ratio between diffusion rates of sterile males and female</td><td>0.8</td><td>—</td></tr></tbody></table>

### 6.1.1 Constant continuous releases

We take three different constant values of $\Lambda \in \{0, 200, 500\}$ (day$^{-1}$). The initial density of sterile males is equal to zero. We approximate the positive equilibria in each case and plot the trajectories of $E_1$ and $E_2$ in Figures 3 according to different values of $\Lambda$. We observe the following:

* When $\Lambda = 0$, there is one positive equilibrium

$$ (E_1^*, E_2^*) = (192.62, 174.82). $$

All positive trajectories converge to the positive steady state $(E_1^*, E_2^*)$.

* When $\Lambda = 200$ (day$^{-1}$), there are two positive equilibria

$$ (E_1^+, E_2^+) = (17.29, 49.98), \quad (E_1^*, E_2^*) = (85.79, 130.02). $$

All positive trajectories also converge to the larger positive steady state $(E_1^*, E_2^*)$.

* When $\Lambda = 500$ (day$^{-1}$), there is no positive equilibrium. All the trajectories converge to the zero equilibrium.

This validates the result in Theorem 4.2 that when $\Lambda$ exceeds some critical value, zero is the unique equilibrium of system (1). The observation for $\Lambda = 0$ illustrates the result in Theorem 3.1 that there is one positive equilibrium and it is globally asymptotically stable. The introduction of sterile males ($\Lambda = 200 > 0$) reduces the value of the positive steady state (see Figure 3b), and when $\Lambda = 500$ (day$^{-1}$) exceeds some critical value (at most equal to 500), all trajectories converge to the zero equilibrium (see Figure 3c). This illustrates the first point of Theorem 4.2. To approximate the critical value of $\Lambda$, we provide some numerical bifurcation diagrams in Section 6.2.

### 6.1.2 Periodic impulsive releases

In this part, we consider the periodic impulsive releases with $\Lambda(t)$ defined in (18), with $\Lambda^{\text{per}}$ equal to 200 and 300 (day$^{-1}$), the period $\tau = 10$ (days). The trajectories of $E_1, E_2$ shown in Figure 4 converge to the periodic solution when $\Lambda^{\text{per}} = 200$ (day$^{-1}$) and go to zero when $\Lambda^{\text{per}} = 300$ (day$^{-1}$). This illustrates the second point of Theorem 4.2 that when the number of sterile males released exceeds a critical value $\bar{\Lambda}^{\text{per}}$, the wild populations of mosquitoes in both areas reach elimination.
<figure><img src="image_9.png" /><figcaption>Figure 8: Densities of wild mosquitoes in two patches in both cases. Left: constant continuous release with $\Lambda = 300$ mosquitoes released per day, Right: periodic impulsive releases with $\Lambda^{\text{per}} = 300$ mosquitoes released at the beginning of each time period.</figcaption></figure>

### 6.2.2 Comparison of release strategies

In practice, the strategy using impulsive releases is more realistic than the constant strategy. In this section, we make a comparison between these two strategies.

For the fixed diffusion rates $d_{12} = 0.06$, $d_{21} = 0.04$, we approximated the critical number of sterile males released in both cases using the method in 6.2

* When $\Lambda(t) \equiv \Lambda$ constant, the critical value $\bar{\Lambda} \approx 250.88 (\text{day}^{-1})$;

* When $\Lambda(t) = \sum_{k=0}^{+\infty} \tau \Lambda^{\text{per}} \delta_{k\tau}$ with period $\tau = 10$, the critical value of $\Lambda^{\text{per}}$ is $\bar{\Lambda}^{\text{per}} \approx 255.15 (\text{day}^{-1})$.

We can see that $\bar{\Lambda}$ and $\bar{\Lambda}^{\text{per}}$ are consistent. We also present numerical simulations in both cases with the same total amount of sterile males released where $\Lambda^{\text{per}} = \Lambda = 300$. The densities of sterile males in both cases are shown in Figure 7. In the constant release case, one observes that the density of sterile males in both zones converges to an equilibrium $(M_1^{s,*}, M_2^{s,*})$. In the periodic case, Figure 7 illustrates the results provided in Lemma 4.1 in which the density of sterile males converges to a periodic solution that is bounded from below.

We obtained in Figure 8 that in both cases, the wild mosquito population reaches elimination at time $t \approx 300$. Again we can see that the two strategies provide the same performance.

## 7 Discussion and conclusion

The existence of some hidden areas (e.g. crab burrows) that can not be accessed by the SIT hinders the population from reaching elimination. Without the implementation of this technique, Theorem 3.1 showed that the wild populations in both areas are persistent and converge towards the unique positive equilibrium (see Figure 3a) and are independent of the diffusion rates between them. The main results obtained in the present work indicated that with a sufficient number of sterile males released, the SIT succeeds in driving both sub-populations to extinction. We investigated both continuous constant releases and impulsive periodic releases in Theorem 4.2. The two strategies provided almost similar performance but the periodic release is more realistic in practice. The idea in our proof can also be used to design a feedback release strategy and this could be studied in future works.
<figure><img src="image_5.png" /><figcaption>Figure 4: Trajectories of E₁ and E₂ in the periodic release case with period τ = 10 (days), diffusion rates d₁₂ = 0.06, d₂₁ = 0.04 (day⁻¹).</figcaption></figure>
The results in Theorem 5.3 also pointed out that the critical numbers of released sterile males are monotone with respect to the biological parameters of the population (see Section 5.2). A population with a larger birth rate and a bigger environment carrying capacity requires more sterile males to reach elimination. A larger death rate in any compartment of the wild mosquitoes reduces this critical value, and on the contrary, a larger death rate for the sterile males increases this value. From the control measure point of view, one can lower the threshold value of $\Lambda$ by killing mosquitoes to increase mosquito mortality or removing breeding sites, even only in the accessible zone to reduce the carrying capacity ($K_1$). This indicates that other conventional control measures can be combined with the SIT control to make the SIT elimination threshold easier to attain.

Moreover, the critical number of sterile males also depends on the diffusion rates between the treated area and the inaccessible zone. More precisely, if the diffusion rates are large, this system approaches the case when there is no separation between two sub-populations (see Theorem 5.2). Numerically, we showed that the larger the values of diffusion rates, the smaller the threshold we need to exceed to obtain elimination (see Figure 6). This also showed that when the movement is at a low level, the leak of wild mosquitoes from the inaccessible area impedes the eradication in the treated zone and it requires a larger number of sterile males to break through this obstacle. In practice, this could be an unrealistic amount of sterile mosquitoes. It is not surprising that the scenario with larger diffusion between two areas is better since more sterile males can arrive at the unreachable zone. Further study on dispersal rates of mosquitoes is necessary to estimate the release rates for SIT elimination in the presence of inaccessible zones.

## Acknowledgments

This study is supported by the STIC AmSud project BIO-CIVIP 23-STIC-02.

<figure><img src="image_10.png" /></figure>

N. Nguyen has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 945332.

## References

[1] CDC, "Dengue Vaccination: What Everyone Should Know | CDC," Jan. 2022.

[2] "Global vector control response 2017-2030.," 2017. World Health Organization.

[3] N. L. Achee, J. P. Grieco, H. Vatandoost, G. Seixas, J. Pinto, L. Ching-NG, A. J. Martins, W. Juntarajum-nong, V. Corbel, C. Gouagna, J.-P. David, J. G. Logan, J. Orsborne, E. Marois, G. J. Devine, and J. Vontas, "Alternative strategies for mosquito-borne arbovirus control," *PLOS Neglected Tropical Diseases*, vol. 13, p. e0006822, Jan. 2019. Publisher: Public Library of Science.

[4] N. Becker, D. Petrić, M. Zgomba, C. Boase, M. B. Madon, C. Dahl, and A. Kaiser, *Mosquitoes: Identification, Ecology and Control*. Fascinating Life Sciences, Cham: Springer International Publishing, 2020.

[5] E. F. Knipling, "Sterile-Male Method of Population Control," *Science*, vol. 130, pp. 902–904, Oct. 1959. Publisher: American Association for the Advancement of Science.

[6] V. A. Dyck, J. Hendrichs, and A. S. Robinson, eds., *Sterile Insect Technique: Principles and Practice in Area-Wide Integrated Pest Management*. Boca Raton: CRC Press, 2 ed., Jan. 2021.

[7] L. F. Chaves, G. L. Hamer, E. D. Walker, W. M. Brown, M. O. Ruiz, and U. D. Kitron, "Climatic variability and landscape heterogeneity impact urban mosquito diversity and vector abundance and infection," *Eco-sphere*, vol. 2, no. 6, p. art70, 2011. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1890/ES11-00088.1.

[8] S. T. Stoddard, B. M. Forshey, A. C. Morrison, V. A. Paz-Soldan, G. M. Vazquez-Prokopec, H. Astete, R. C. Reiner, S. Vilcarromero, J. P. Elder, E. S. Halsey, T. J. Kochel, U. Kitron, and T. W. Scott, "House-to-house human movement drives dengue virus transmission," *Proceedings of the National Academy of Sciences*, vol. 110, pp. 994–999, Jan. 2013. Publisher: Proceedings of the National Academy of Sciences.

[9] P. Martens and L. Hall, "Malaria on the move: human population movement and malaria transmission.," *Emerging Infectious Diseases*, vol. 6, no. 2, pp. 103–109, 2000.

[10] F. V. S. d. Abreu, M. M. Morais, S. P. Ribeiro, and A. E. Eiras, "Influence of breeding site availability on the oviposition behaviour of *Aedes aegypti*," *Memorias do Instituto Oswaldo Cruz*, vol. 110, pp. 669–676, July 2015. Publisher: Instituto Oswaldo Cruz, Ministério da Saúde.

[11] R. Anguelov, Y. Dumont, and J. Lubuma, "Mathematical modeling of sterile insect technology for control of anopheles mosquito," *Computers & Mathematics with Applications*, vol. 64, pp. 374–389, Aug. 2012.

[12] Y. Dumont and J. M. Tchuenche, "Mathematical studies on the sterile insect technique for the Chikungunya disease and Aedes albopictus," *Journal of Mathematical Biology*, vol. 65, pp. 809–854, Nov. 2012.
<figure><img src="image_6.png" /><figcaption>Figure 5: Bifurcation diagrams of $E_1^*$ with parameter $\Lambda$ in the constant continuous release case.</figcaption></figure>

## 6.2 Critical values and bifurcation

Our aim in this section is to approximate the critical value of $\Lambda$ where the bifurcation occurs.

### 6.2.1 Bifurcation diagram in the constant release case

We solve a system of nonlinear stationary problem $\mathcal{F}(u;\Lambda) = 0$ for all values of the parameter $\Lambda$, knowing that the solutions are continuous with respect to $\Lambda$. Solving by numerical approximations can be done using numerical continuation methods (see [32]).

Here we present the simplest method called *Natural Parameter Continuation* (incremental methods, see [32]): Iteratively find approximate roots of $\mathcal{F}(u, \Lambda) = 0$ for several values of $\Lambda_i$ with index $i \in \mathbb{N}^*$. The root of step $i$ is used as an initial guess for the numerical solver at step $i+1$. The first initial guess is the root for the smallest $\Lambda$. To approximate the critical value $\bar{\Lambda}$ in the constant case and examine what happens when $0 < \Lambda \le \bar{\Lambda}$, we draw the bifurcation diagram for $\Lambda \in [0.1, 500]$. The initial positions of the numerical continuation are taken at the approximated equilibria when $\Lambda = 0.1$.

We obtain the bifurcation diagrams in Figure 5 for two scenarios. We observed that the critical value of $\Lambda$ decreases when the diffusion rates increase.

* For $d_{12} = 1$, $d_{21} = 2$, the critical value $\bar{\Lambda} = 106.45$ (day$^{-1}$).

* For $d_{12} = 0.06$, $d_{21} = 0.04$, the critical value $\bar{\Lambda} = 250.88$ (day$^{-1}$).

Taking $d_{12} = d_{21}$, we plot the critical value $\bar{\Lambda}$ corresponding to the moving rates $d_{12}$ (see Figure 6). This shows that the value of $\bar{\Lambda}$ decreases when the diffusion rate gets larger, and converges to a value $\bar{\Lambda}_\infty \approx 109.45$ (day$^{-1}$) as $d_{12}$ goes to infinity. This validates the result provided by Proposition 5.2 where $\bar{\Lambda}_\infty$ is the critical value of $\Lambda$ corresponding to system (27). We also found that $\bar{\Lambda}_\infty = \bar{\Lambda}_0$ where $\bar{\Lambda}_0$ is the critical value of the system when there is no separation between the two sub-populations defined in 5.1.2.

### 6.2.2 Comparison of release strategies

In practice, the strategy using impulsive releases is more realistic than the constant strategy. In this section, we make a comparison between these two strategies.

For the fixed diffusion rates $d_{12} = 0.06$, $d_{21} = 0.04$, we approximated the critical number of sterile males released in both cases using the method in 6.2
[13] M. Strugarek, H. Bossin, and Y. Dumont, “On the use of the sterile insect release technique to reduce or eliminate mosquito populations,” *Applied Mathematical Modelling*, vol. 68, pp. 443–470, Apr. 2019.

[14] P.-A. Bliman, D. Cardona-Salgado, Y. Dumont, and O. Vasilieva, “Implementation of control strategies for sterile insect techniques,” *Mathematical Biosciences*, vol. 314, pp. 43–60, Aug. 2019.

[15] R. Anguelov, Y. Dumont, and I. V. Yatat Djeumen, “Sustainable vector/pest control using the permanent sterile insect technique,” *Mathematical Methods in the Applied Sciences*, vol. 43, no. 18, pp. 10391–10412, 2020. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/mma.6385.

[16] P.-A. Bliman and Y. Dumont, “Robust control strategy by the Sterile Insect Technique for reducing epidemiological risk in presence of vector migration,” *Mathematical Biosciences*, vol. 350, p. 108856, Aug. 2022.

[17] V. S. Manoranjan and P. Van Den Driessche, “On a diffusion model for sterile insect release,” *Mathematical Biosciences*, vol. 79, pp. 199–208, June 1986.

[18] M. A. Lewis and P. Van Den Driessche, “Waves of extinction from sterile insect release,” *Mathematical Biosciences*, vol. 116, pp. 221–247, Aug. 1993.

[19] L. Almeida, J. Estrada, and N. Vauchelet, “Wave blocking in a bistable system by local introduction of a population: application to sterile insect techniques on mosquito populations,” *Mathematical Modelling of Natural Phenomena*, vol. 17, p. 22, 2022. Publisher: EDP Sciences.

[20] L. Almeida, A. Léculier, and N. Vauchelet, “Analysis of the Rolling Carpet Strategy to Eradicate an Invasive Species,” *SIAM Journal on Mathematical Analysis*, vol. 55, pp. 275–309, Feb. 2023. Publisher: Society for Industrial and Applied Mathematics.

[21] A. Leculier and N. Nguyen, “A control strategy for the sterile insect technique using exponentially decreasing releases to avoid the hair-trigger effect,” *Mathematical Modelling of Natural Phenomena*, vol. 18, pp. 1–25, 2023. Publisher: EDP Sciences.

[22] C. Dufourd and Y. Dumont, “Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control,” *Computers & Mathematics with Applications*, vol. 66, pp. 1695–1715, Nov. 2013.

[23] L. Almeida, A. Leculier, N. Nguyen, and N. Vauchelet, “Rolling carpet strategies to eliminate mosquitoes: two-dimensional case.” (paper in preparation), 2023.

[24] P. Auger, E. Kouokam, G. Sallet, M. Tchuente, and B. Tsanou, “The Ross–Macdonald model in a patchy environment,” *Mathematical Biosciences*, vol. 216, pp. 123–131, Dec. 2008.

[25] A. M. Lutambi, M. A. Penny, T. Smith, and N. Chitnis, “Mathematical modelling of mosquito dispersal in a heterogeneous environment,” *Mathematical Biosciences*, vol. 241, pp. 198–216, Feb. 2013.

[26] M. L. Mann Manyombe, B. Tsanou, J. Mbang, and S. Bowong, “A metapopulation model for the population dynamics of anopheles mosquito,” *Applied Mathematics and Computation*, vol. 307, pp. 71–91, Aug. 2017.

[27] D. D. Bonnet and H. Chapman, “The larval habitats of Aedes polynesiensis Marks in Tahiti and methods of control.,” *American Journal of Tropical Medicine and Hygiene*, vol. 7, no. 5, 1958. Publisher: Baltimore, Md.

[28] F. Lardeux, F. Riviere, Y. Sechan, and B. H. Kay, “Release of Mesocyclops aspericornis (Copepoda) for Control of Larval Aedes polynesiensis (Diptera: Culicidae) in Land Crab Burrows on an Atoll of French Polynesia,” *Journal of Medical Entomology*, vol. 29, pp. 571–576, July 1992.

[29] F. Lardeux, Y. Sechan, and M. Faaruia, “Evaluation of Insecticide Impregnated Baits for Control of Mosquito Larvae in Land Crab Burrows on French Polynesian Atolls,” *Journal of Medical Entomology*, vol. 39, pp. 658–661, July 2002.

[30] K. Heath, M. B. Bonsall, J. Marie, and H. C. Bossin, “Mathematical modelling of the mosquito Aedes polynesiensis in a heterogeneous environment,” *Mathematical Biosciences*, vol. 348, p. 108811, June 2022.

[31] R. I. Vargas, L. Whitehand, W. A. Walsh, J. P. Spencer, and C.-L. Hsu, “Aerial Releases of Sterile Mediterranean Fruit Fly (Diptera: Tepbritidae) by Helicopter: Dispersal, Recovery, and Population Suppression,” *Journal of Economic Entomology*, vol. 88, pp. 1279–1287, Oct. 1995.

[32] C. Yang, X. Zhang, and J. Li, “Dynamics of two-patch mosquito population models with sterile mosquitoes,” *Journal of Mathematical Analysis and Applications*, vol. 483, p. 123660, Mar. 2020.

[33] L. J. S. Allen, “Persistence, extinction, and critical patch number for island populations,” *Journal of Mathematical Biology*, vol. 24, pp. 617–625, Feb. 1987.

[34] Y. Takeuchi, “Cooperative Systems Theory and Global Stability of Diffusion Models,” in *Evolution and Control in Biological Systems* (A. B. Kurzhanski and K. Sigmund, eds.), (Dordrecht), pp. 49–57, Springer Netherlands, 1989.
<figure><img src="image_7.png" /><figcaption>Figure 6: Dependence of $\bar{\Lambda}$ on the diffusion rate $d_{12}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Densities of $M_1^s$ and $M_2^s$ in both cases. Left: constant continuous releases with $\Lambda = \bar{\Lambda} = 250.88$ (day$^{-1}$), Right: periodic impulsive releases with $\Lambda^{\text{per}} = \bar{\Lambda}^{\text{per}} = 255.15$ (day$^{-1}$).</figcaption></figure>
[35] Z. Lu and Y. Takeuchi, "Global asymptotic behavior in single-species discrete diffusion systems," *J. Math. Biol.*, vol. 32, pp. 67-77, Nov. 1993.

[36] M. Y. Li and Z. Shuai, "Global-stability problem for coupled systems of differential equations on networks," *Journal of Differential Equations*, vol. 248, pp. 1-20, Jan. 2010.

[37] W. A. Coppel, *Stability and Asymptotic Behavior of Differential Equations*. Heath, 1965.

[38] S. Chaplygin, "Selected works on mechanics and mathematics," *State Publ. House, Technical-Theoretical Literature, Moscow*, 1954.

[39] H. Smith, *Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems*, vol. 41 of *Mathematical Surveys and Monographs*. Providence, Rhode Island: American Mathematical Society, Mar. 2008.

[40] W. C. Rheinboldt, "Numerical continuation methods: a perspective," *Journal of Computational and Applied Mathematics*, vol. 124, pp. 229-244, Dec. 2000.
<figure><img src="image_9.png" /><figcaption>Figure 8: Densities of wild mosquitoes in two patches in both cases. Left: constant continuous release with $\Lambda = 300$ mosquitoes released per day, Right: periodic impulsive releases with $\tau\Lambda^{\text{per}} = 3000$ mosquitoes released at the beginning of each time period.</figcaption></figure>

* When $\Lambda(t) \equiv \Lambda$ constant, the critical value $\bar{\Lambda} \approx 250.88 (\text{day}^{-1})$;

* When $\Lambda(t) = \sum_{k=0}^{+\infty} \tau\Lambda^{\text{per}}\delta_k\tau$ with period $\tau = 10$, the critical value of $\Lambda^{\text{per}}$ is $\bar{\Lambda}^{\text{per}} \approx 255.15 (\text{day}^{-1})$.

We can see that $\bar{\Lambda}$ and $\bar{\Lambda}^{\text{per}}$ are consistent. We also present numerical simulations in both cases with the same total amount of sterile males released where $\Lambda^{\text{per}} = \Lambda = 300$. The densities of sterile males in both cases are shown in Figure 7. We obtained in Figure 8 that in both cases, the wild mosquito population reaches elimination at time $t \approx 300$. Again we can see that the two strategies provide the same performance.

# 7 Discussion and conclusion

The existence of some hidden areas (e.g. crab burrows) that can not be accessed by the SIT hinders the population from reaching elimination. Without the implementation of this technique, Theorem 3.1 showed that the wild populations in both areas are persistent and converge towards the unique positive equilibrium (see Figure 3a) and are independent of the diffusion rates between them. The main results obtained in the present work indicated that with a sufficient number of sterile males released, the SIT succeeds in driving both sub-populations to extinction. We investigated both continuous constant releases and impulsive periodic releases in Theorem 4.2. The two strategies provided almost similar performance but the periodic release is more realistic in practice. The idea in our proof can also be used to design a feedback release strategy and this could be studied in future works.
The results also pointed out that the critical numbers of released sterile males are monotone with respect to the biological parameters of the population (see Section 5.2). A population with a larger birth rate of wild mosquitoes and a bigger carrying capacity of the environment requires more sterile males to reach elimination. A larger death rate in any compartment of the wild mosquitoes reduces this critical value, and on the contrary larger death rate for the sterile males increases this value.

Moreover, the critical number of sterile males also depends on the diffusion rates between the treated area and the inaccessible zone. More precisely, if the diffusion rates are large, this system approaches the case when there is no separation between two sub-populations (see Theorem 5.2). Numerically, we showed that the larger the values of diffusion rates, the smaller the threshold we need to exceed to obtain elimination (see Figure 6). This also showed that when the movement is at a low level, the leak of wild mosquitoes from the inaccessible area impedes the eradication in the treated zone and it requires a large number of sterile males to break through this obstacle. In practice, this could be an unrealistic amount of sterile mosquitoes. It is not surprising that the scenario with larger diffusion between two areas is better since more sterile males can arrive at the unreachable zone.

## Acknowledgments

This study is supported by the STIC AmSud project BIO-CIVIP 23-STIC-02.

<figure><img src="image_10.png" /></figure>

N. Nguyen has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 945332.

## References

[1] ABREU, F. V. S. D., MORAIS, M. M., RIBEIRO, S. P., AND EIRAS, A. E. Influence of breeding site availability on the oviposition behaviour of *Aedes aegypti*. *Memorias do Instituto Oswaldo Cruz* 110 (July 2015), 669–676. Publisher: Instituto Oswaldo Cruz, Ministério da Saúde.

[2] ACHEE, N. L., GRIECO, J. P., VATANDOOST, H., SEIXAS, G., PINTO, J., CHING-NG, L., MARTINS, A. J., JUNTARAJUMNONG, W., CORBEL, V., GOUAGNA, C., DAVID, J.-P., LOGAN, J. G., ORSBORNE, J., MAROIS, E., DEVINE, G. J., AND VONTAS, J. Alternative strategies for mosquito-borne arbovirus control. *PLOS Neglected Tropical Diseases* 13, 1 (Jan. 2019), e0006822. Publisher: Public Library of Science.

[3] ALLEN, L. J. S. Persistence, extinction, and critical patch number for island populations. *Journal of Mathematical Biology* 24, 6 (Feb. 1987), 617–625.

[4] ALMEIDA, L., ESTRADA, J., AND VAUCHELET, N. Wave blocking in a bistable system by local introduction of a population: application to sterile insect techniques on mosquito populations. *Mathematical Modelling of Natural Phenomena* 17 (2022), 22. Publisher: EDP Sciences.

[5] ALMEIDA, L., LECULIER, A., NGUYEN, N., AND VAUCHELET, N. Rolling carpet strategies to eliminate mosquitoes: two-dimensional case. (paper in preparation), 2023.

[6] ALMEIDA, L., LÉCULIER, A., AND VAUCHELET, N. Analysis of the Rolling Carpet Strategy to Eradicate an Invasive Species. *SIAM Journal on Mathematical Analysis* 55, 1 (Feb. 2023), 275–309. Publisher: Society for Industrial and Applied Mathematics.

[7] ANGUELOV, R., DUMONT, Y., AND LUBUMA, J. Mathematical modeling of sterile insect technology for control of anopheles mosquito. *Computers & Mathematics with Applications* 64, 3 (Aug. 2012), 374–389.
[8] ANGUELOV, R., DUMONT, Y., AND YATAT DJEUMEN, I. V. Sustainable vector/pest control using the permanent sterile insect technique. *Mathematical Methods in the Applied Sciences* 43, 18 (2020), 10391–10412. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/mma.6385.

[9] AUGER, P., KOUOKAM, E., SALLET, G., TCHUENTE, M., AND TSANOU, B. The Ross–Macdonald model in a patchy environment. *Mathematical Biosciences* 216, 2 (Dec. 2008), 123–131.

[10] BECKER, N., PETRIĆ, D., ZGOMBA, M., BOASE, C., MADON, M. B., DAHL, C., AND KAISER, A. *Mosquitoes: Identification, Ecology and Control*. Fascinating Life Sciences. Springer International Publishing, Cham, 2020.

[11] BLIMAN, P.-A., CARDONA-SALGADO, D., DUMONT, Y., AND VASILIEVA, O. Implementation of control strategies for sterile insect techniques. *Mathematical Biosciences* 314 (Aug. 2019), 43–60.

[12] BLIMAN, P.-A., AND DUMONT, Y. Robust control strategy by the Sterile Insect Technique for reducing epidemiological risk in presence of vector migration. *Mathematical Biosciences* 350 (Aug. 2022), 108856.

[13] BONNET, D. D., AND CHAPMAN, H. The larval habitats of Aedes polynesiensis Marks in Tahiti and methods of control. *American Journal of Tropical Medicine and Hygiene* 7, 5 (1958). Publisher: Baltimore, Md.

[14] CHAPLYGIN, S. Selected works on mechanics and mathematics. *State Publ. House, Technical-Theoretical Literature, Moscow* (1954).

[15] CHAVES, L. F., HAMER, G. L., WALKER, E. D., BROWN, W. M., RUIZ, M. O., AND KITRON, U. D. Climatic variability and landscape heterogeneity impact urban mosquito diversity and vector abundance and infection. *Ecosphere* 2, 6 (2011), art70. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1890/ES11-00088.1.

[16] COPPEL, W. A. *Stability and Asymptotic Behavior of Differential Equations*. Heath, 1965.

[17] DUFOURD, C., AND DUMONT, Y. Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control. *Computers & Mathematics with Applications* 66, 9 (Nov. 2013), 1695–1715.

[18] DUMONT, Y., AND TCHUENCHE, J. M. Mathematical studies on the sterile insect technique for the Chikungunya disease and Aedes albopictus. *Journal of Mathematical Biology* 65, 5 (Nov. 2012), 809–854.

[19] DYCK, V. A., HENDRICHS, J., AND ROBINSON, A. S., Eds. *Sterile Insect Technique: Principles and Practice in Area-Wide Integrated Pest Management*, 2 ed. CRC Press, Boca Raton, Jan. 2021.

[20] HEATH, K., BONSALL, M. B., MARIE, J., AND BOSSIN, H. C. Mathematical modelling of the mosquito Aedes polynesiensis in a heterogeneous environment. *Mathematical Biosciences* 348 (June 2022), 108811.

[21] KNIPLING, E. F. Sterile-Male Method of Population Control. *Science* 130, 3380 (Oct. 1959), 902–904. Publisher: American Association for the Advancement of Science.

[22] LARDEUX, F., RIVIERE, F., SECHAN, Y., AND KAY, B. H. Release of Mesocyclops aspericornis (Copepoda) for Control of Larval Aedes polynesiensis (Diptera: Culicidae) in Land Crab Burrows on an Atoll of French Polynesia. *Journal of Medical Entomology* 29, 4 (July 1992), 571–576.
[23] LARDEUX, F., SECHAN, Y., AND FAARUIA, M. Evaluation of Insecticide Impregnated Baits for Control of Mosquito Larvae in Land Crab Burrows on French Polynesian Atolls. *Journal of Medical Entomology* 39, 4 (July 2002), 658–661.

[24] LECULIER, A., AND NGUYEN, N. A control strategy for the sterile insect technique using exponentially decreasing releases to avoid the hair-trigger effect. *Mathematical Modelling of Natural Phenomena* 18 (2023), 1–25. Publisher: EDP Sciences.

[25] LEWIS, M. A., AND VAN DEN DRIESSCHE, P. Waves of extinction from sterile insect release. *Mathematical Biosciences* 116, 2 (Aug. 1993), 221–247.

[26] LI, M. Y., AND SHUAI, Z. Global-stability problem for coupled systems of differential equations on networks. *Journal of Differential Equations* 248, 1 (Jan. 2010), 1–20.

[27] LU, Z., AND TAKEUCHI, Y. Global asymptotic behavior in single-species discrete diffusion systems. *J. Math. Biol.* 32, 1 (Nov. 1993), 67–77.

[28] LUTAMBI, A. M., PENNY, M. A., SMITH, T., AND CHITNIS, N. Mathematical modelling of mosquito dispersal in a heterogeneous environment. *Mathematical Biosciences* 241, 2 (Feb. 2013), 198–216.

[29] MANN MANYOMBE, M. L., TSANOU, B., MBANG, J., AND BOWONG, S. A metapopulation model for the population dynamics of anopheles mosquito. *Applied Mathematics and Computation* 307 (Aug. 2017), 71–91.

[30] MANORANJAN, V. S., AND VAN DEN DRIESSCHE, P. On a diffusion model for sterile insect release. *Mathematical Biosciences* 79, 2 (June 1986), 199–208.

[31] MARTENS, P., AND HALL, L. Malaria on the move: human population movement and malaria transmission. *Emerging Infectious Diseases* 6, 2 (2000), 103–109.

[32] RHEINBOLDT, W. C. Numerical continuation methods: a perspective. *Journal of Computational and Applied Mathematics* 124, 1 (Dec. 2000), 229–244.

[33] SMITH, H. *Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems*, vol. 41 of *Mathematical Surveys and Monographs*. American Mathematical Society, Providence, Rhode Island, Mar. 2008.

[34] STODDARD, S. T., FORSHEY, B. M., MORRISON, A. C., PAZ-SOLDAN, V. A., VAZQUEZ-PROKOPEC, G. M., ASTETE, H., REINER, R. C., VILCARROMERO, S., ELDER, J. P., HALSEY, E. S., KOCHEL, T. J., KITRON, U., AND SCOTT, T. W. House-to-house human movement drives dengue virus transmission. *Proceedings of the National Academy of Sciences* 110, 3 (Jan. 2013), 994–999. Publisher: Proceedings of the National Academy of Sciences.

[35] STRUGAREK, M., BOSSIN, H., AND DUMONT, Y. On the use of the sterile insect release technique to reduce or eliminate mosquito populations. *Applied Mathematical Modelling* 68 (Apr. 2019), 443–470.

[36] TAKEUCHI, Y. Cooperative Systems Theory and Global Stability of Diffusion Models. In *Evolution and Control in Biological Systems* (Dordrecht, 1989), A. B. Kurzhanski and K. Sigmund, Eds., Springer Netherlands, pp. 49–57.

[37] Global vector control response 2017-2030., 2017. World Health Organization.

[38] YANG, C., ZHANG, X., AND LI, J. Dynamics of two-patch mosquito population models with sterile mosquitoes. *Journal of Mathematical Analysis and Applications* 483, 2 (Mar. 2020), 123660.