Wave blocking in a bistable system by local introduction of a population: application to sterile insect techniques on mosquito populations

,

Technique (IIT), where released mosquitoes are not sterilized but rather infected with the bacterium Wolbachia, which shortens their lifespan and reduces their vector capacity (for instance for transmitting dengue). Indeed, when only Wolbachia-bearing males are released, the IIT is equivalent to the classical SIT; when both males and females are released, the expected result is that the existing population will be replaced by the Wolbachia infected population. In particular, for mosquito control, the SIT/IIT have been employed in multiple locations, such as Italy [4], Reunion Island [18], Polynesia [8] and China [25]. Aedes mosquitoes and the arboviroses for which they are transmission vectors have been a major problem in tropical regions since long ago. However, they are also becoming a major concern for public health authorities in temperate regions (including Europe) where, due to climate change and global exchanges, Aedes mosquitoes have become a successful invasive species..
Classical SIT and IIT have been modeled and studied in a large number of papers, from diverse mathematical viewpoints, such as discrete difference equations [16], as ODE systems by themselves [24] or linked with SIR systems of ODEs [10]. The spatial dynamics has been also studied thanks to PDE models to describe the spatial invasion of mosquitoes : for instance the recent work [21] proposes a study of the invasion of Aedes Albopictus in France. In [13], the authors propose a simple scalar reaction-diffusion equation to analyze the influence of releases of sterilized males on the traveling waves of invasion of mosquitoes.
In this paper, we are interested in the use of SIT/IIT locally as a way to block invasion or infestation of mosquitoes in a mosquito-free area. More precisely, we study the possibility of blocking a propagating front of mosquitoes by performing a local release of sterile mosquitoes in a band of width L. Such technique may be used as a sanitary cordon to avoid infestation of a mosquito-free area or re-infestation of a previously treated area. A numerical study of such strategy has been performed in [23] and in [1]. It is expected that when releasing enough sterilized mosquitoes in a wide enough area, it should be possible to block invading waves of mosquitoes. The aim of this paper is to justify rigorously these numerical results. From a mathematical point of view, it boils down to show the existence of stationary solutions for a system of reaction-diffusion equations.
More generally, let v(x, t) denote the density of a species, e.g. mosquitoes, at position x ∈ R and time t > 0. Let w(x, t) denote the density of an introduced species, e.g. sterilized males, which is released on a domain [0, L] at a constant density M 0 ; we denote µ the death rate of this species. Then, the population dynamics is governed by the following system of reaction-diffusion equations: (1.1b) In this system, the reaction term g describes the interaction between both species which is considered competitive in this study. It is assumed to be quite general (the complete statement of the assumptions is provided in the statement of Thm. 2.3). In the absence of the species w, i.e. when w = 0, we denote f (v) := g(v, 0) the reaction term for the first species. Since we want to study the elimination of the population v, we will assume that the steady state 0 is stable. Then we consider that f is a C 1 (R) bistable function, that is ∃ θ ∈ (0, 1), f (0) = f (1) = f (θ) = 0, f < 0 on (0, θ), f > 0 on (θ, 1).
Moreover, we assume that It is well-known (see e.g. [3,14]) that under this condition, there exist invading traveling waves for species v in absence of species w. The aim of this paper is to prove that, assuming appropriate conditions on the interaction function g, for L > 0 and M 0 large enough it is possible to block the invasion. In other words, we study the possibility to block the invasion by acting on a release function which is chosen here to be the piece-wise constant function M 0 1 [0,L] .
Wave-blocking in reaction-diffusion equations has been studied by several authors. For instance, the influence of the geometry on the propagation of a potential along a nerve axon has been considered in [19]. In the seminal paper [15], the authors show the existence of a wave-block due to heterogeneities in the media of propagation. To obtain their result, they introduce a geometrical technique by reasoning on the phase plan. This technique will be recalled below. In the present paper, we adapt and extend their technique to the problem at hand. In [7], the authors propose to study the absorption by the white matter of propagating waves in the brain. In [5], wave-block in a model for criminality is proved. Wave-block in bistable reaction-diffusion systems with a drift term, usually called gene-flow models, has been investigated in [11,17]. Finally, we mention that control strategies on reaction-diffusion equations with bistable nonlinearity has been recently investigated by several authors e.g. [20,22].
The outline of the paper is as follows. In Section 2, we state our main results and we describe our mathematical model. Section 3 is devoted to the proof of our main result for general bistable systems, i.e. the existence of a barrier when acting on a wide enough region with large enough intensity. We use a geometric method based on a phase plane analysis to construct such barrier. In Section 4 we apply this technique and result to our practical problem: the use of the sterile insect technique to avoid invasion of a species like mosquitoes. Then, we illustrate our results with some numerical applications in Section 5. We end this article with some conclusion in Section 6.

Modelling and main results
In this section we state the main results of this paper for system (1.1), that is the blocking of the propagation under appropriate conditions. Then, we explain how to use such results to obtain in the application we have in mind, that is the local use of the sterile insect technique to control invasion by mosquito populations.

Main results for general bistable systems
We first start by defining a barrier for system (1.1): Definition 2.1. Let L ≥ 0 and M 0 ≥ 0. We call a barrier (of width L and height M 0 ) for system (1.1), any couple (v,w) ∈ (C 2 (R)) 2 which is a nonnegative solution of the following differential inequalities: We say that there exists a barrier to system (1.1) for L ≥ 0 and M 0 ≥ 0 if system (2.1) admits a nonnegative solution.
The following Lemma explains why we call solutions to (2.1) a barrier : indeed, the existence of such solutions may block the propagation. Lemma 2.2. Let g be Lipschitz-continuous and such that ∂ w g(v, w) ≤ 0 for all v, w. Let us assume that there exists a barrier for system (1.1) (as in Def. 2.1). If the initial data (v 0 , w 0 ) is such that v 0 ≤v and w 0 ≥w, then the solution to (1.1) is such that for all t > 0, v(·, t) ≤v.
Proof. It is straightforward to obtain this result using a comparison principle. Indeed, from (1.1) and (2.1), we deduce that We multiply the first equation by (v −v) + and the second equation by (w − w) + and integrate over R, where (·) + denotes the positive part. Summing the resulting equalities, we get From the assumptions on g, we deduce where C g is a Lipschitz constant for g. We conclude as usual by using a Gronwall type inequality and the fact that As a consequence of Lemma 2.2, if the species w is present in the environment in such a way that w(0, ·) ≥w, then any invasion front is blocked by the barrier, provided such a barrier exists. Under assumption (1.
Hence, when L = 0, there is no barrier. The following Theorem gives conditions for existence of barriers: Theorem 2.3. Let us assume that g satisfies the following properties: Remark 2.4. Let us mention that, from an experimental point of view, it is not difficult to guarantee the condition w 0 ≥w to build a barrier. Indeed, it is enough to make the release with M ≥ M 0 or on a domain large enough.

Modelling for the sterile insect technique
The situation to be modeled is the effect that a release of sterilized males in a limited region in space will have on a mosquito population that is already present in part of the domain. For the sake of simplicity, we will also consider only one dimension in space and consider all biological parameters to be constant over time, thus disregarding the effects of field heterogeneity and seasonal variations.
Following ideas in e.g. [6,24], we model the population of mosquitoes with four compartements : eggs, adult males, adult females, and sterilized mosquitoes. The dynamics is governed by the following system for t > 0 and x ∈ R: The variables and parameters included in our model are listed below: -E, M and F denote respectively the spatial density of mosquitoes in aquatic phase, adult males, and adult females; -M s (x, t) is the density of sterile mosquitoes which are released at point x and time t, the release function is denoted u; from now on, we will consider that the release function is u(x) = M 0 · 1 [0,L] (x); -the term (1 − e −β(M +γMs) ) is an Allee effect to model the fact that it may be difficult for a female to find a partner when the density of mosquitoes is small, similar terms have been introduced in [24]; β is a parameter to quantify this Allee effect : the smaller this parameter, the stronger the Allee effect, whereas it disappears when β tends to +∞. -the fraction M M +γMs corresponds to the probability that a female mates with a noninfected mosquito; -the parameter γ ∈ (0, 1] quantifies the mating preference of females between sterile and fertile males. If γ = 1 there is no preference between sterile and fertile male, whereas if γ = 0 females only consider fertile males for mating; -b > 0 is a birth rate; µ E > 0, µ M > 0, µ F > 0, and µ s > 0 denote the death rates for the mosquitoes in the aquatic phase, for adult males, for adult females, and for sterile males respectively; -K is an environmental capacity for the aquatic phase, accounting also for the intraspecific competition; -ν E > 0 is the rate of emergence; -r ∈ (0, 1) is the probability that a female emerges, then (1 − r) is the probability that a male emerges.
Since mathematical analysis for systems like (2.2) is quite complicated to handle, we will make some modelling assumptions to reduce (2.2) to a system such as (1.1) for which we will perform a full analysis. To do so, we assume that r = 1 2 , µ F < µ M , and M (x, 0) ≤ F (x, 0), so that the system rewrites Notice that the above assumptions are reasonable since the sex ratio is indeed expected to be close to 1 2 and the male mortality is higher than the female mortality (see e.g. [24]). Our second assumption is to consider that the dynamics of the egg compartment is fast, which boils down to assume that the first equation is at equilibrium: From here we can determine E: We may verify easily that ∂ ∂M E(F, M ) ≥ 0 and ∂ ∂F E(F, M ) ≥ 0. The second and third equations of (2.3) simplify into We mention that a similar minimalist model for the sterile insect technique, without the space dependence, has been introduced and studied in [6]. Notice that the main difference is that the system in [6] is not cooperative, whereas E being increasing with respect to F and M , system Hence we deduce that, if initially we have M (0, x) ≤ F (0, x), then by the comparison principle for reactiondiffusion equations, we have M ≤ F . Then, from (2.4b) and the fact that Finally, we will consider the following system This latter system has the same form as system (1.1), and any solution to (2.5)-(2.6) is a super-solution to (2.4)-(2.6). Hence, it is enough to prove the existence of a barrier for (2.5)-(2.6) to show that a blocking exists. We will show that we may apply Theorem 2.3 and obtain the existence of a barrier for L and M 0 large enough. Moreover, with the particular form of the reaction term g in this setting, we may obtain more information. More precisely, the main results are summarized in the following Then, for all L > 0, there exists M * (L) such that for all M 0 ≥M * (L), there exists a barrier for system (2.5)-(2.6).
And the function L → M * (L) has the properties (i), (ii), and (iii) in Theorem 2.3.
Remark 2.6. We will see in Lemma 4.2 that the condition on the coefficients in the statement of this Theorem guarantees that the function F → g(F, 0) is bistable and verifies (1.3). If this condition is not satisfied then there is no invasion by the population of mosquitoes and thus there is no need to perform a sterile insect technique.

A geometric construction of steady states
In order to build a barrier we will use the geometrical construction method which has been proposed by Lewis and Keener in [15]. Let us consider the system, for x ∈ R, The idea is to work on the phase plane (u, u x ). We know (see e.g. [14]) that the equation u xx + f (u) = 0 has a family of translated solutions any ξ ∈ R. This family traces out a homoclinic orbit associated with (u, u x ) = (0, 0), which we call curve A. In u xx + f (u) = 0, the steady state u ≡ 1 is stable, so there also exists a curve B which is its stable manifold, with points on it tending to (u, u x ) = (1, 0) when x → +∞.
We can derive an analytical expression for these two curves: Multiplying u xx + f (u) = 0 by u x and integrating we get where C is a suitable constant. Taking u = 1, u x = 0 gives us C = 1 0 f (s)ds and thus the expression for curve B is Likewise, taking u = u x = 0 gives us the expression for curve A: To summarize, we have where θ c is such that θc 0 f (s) ds = 0 and from (1.3) we have 0 < θ c < 1. We take a point on curve A which lies on the first quadrant and consider it as an initial condition (u(0), u x (0)), with u(0) > 0 and u x (0) > 0 (we can choose x = 0 as the starting point because ξ is arbitrary). We solve the second-order differential equation u (x) = −k(x, u) on (0, L) with this initial condition (u(0), u x (0)) to get (u(L), u x (L)) (assuming u does not blow up before x = L). If (u(L), u x (L)) falls on curve B, then for x > L we follow curve B. By doing so, we have constructed a stationary solution of (3.1) such that u(−∞) = 0 and u(+∞) = 1 (see Fig. 1). In order to determine when this happens, for a given L > 0, we define the following mapping ψ L by Proof. The idea of the proof is to show that for L large enough the curves ψ L (A) and B will intersect. We first notice that for any L > 0, we have ψ L (0, 0) = (0, 0). On [0, L], we have k(x, u) ≤ −αu. Then, solving u (x) − αu(x) = ν(x) where ν(x) := −αu(x) − k(x, u(x)) ≥ 0, with nonnegative initial data (u(0), u (0)), we get on [0, L], We have c 1 > |c 2 | ≥ 0. Then, since ν ≥ 0, the last term of the right hand side is clearly nonnegative. We deduce Then, since u (0) > 0, we have We also have Moreover, as above, we have is increasing and tends to +∞ when √ αL → +∞, tnot only does the mapping ψ L increase both coordinates, but u(L), u (L) → +∞ when L → +∞ for any u(0), u (0) > 0. This guarantees that there exists L * large enough such that, for any L ≥ L * , the curves ψ L (A) and B intersect, and we have existence of a solution of (3.1) with limit 0 at −∞ and with limit 1 at +∞.

Proof of Theorem 2.3
Before proving our first main result, we provide some useful informations.
After straightforward computations, we see that If L < L * , we obtain Thus, in both situations, we get that for M 0 large enough, We are now in position to prove our first main result. Proof of Theorem 2.3. We split the proof into several steps.
Step 1. Existence of a barrier.
Let L > 0, M 0 > 0, and w be a solution of the equation Let L * > 0 which will be fixed later. As a consequence of Lemma 3.2, we have that for M 0 ≥ M * large enough, , we clearly have w ≥ 0, and thus, g(v, w) ≤ g(v, 0) since, from the assumptions on g in the statement of Theorem 2.3, the function g is decreasing with respect to its second variable. Let us introduce v 2 solution of the problem (3.5) We clearly have from above remarks that v 2 + g( v 2 , w) ≤ 0 on R. Since f := g(·, 0) is bistable, we may apply Proposition 3.1 for system (3.5): There exists L * > 0 such that system (3.5) admits a solution with v 2 (−∞) = 0 and v 2 (+∞) = 1. Then, the couple ( v 2 , w) verifies (2.1). Hence, we have constructed a barrier for any L > 0 and for M 0 large enough. Let us recall that for L = 0 or M 0 = 0, we have w = 0 and we know that, due to Step 3. Proof of (ii).
As a consequence of (i), the limits lim From the first two assumptions on g in the statement of Theorem 2.3, there exists η > 0 small enough such that v → g(v, η) is bistable, admits 3 roots 0 < v 1,η < v 2,η and verifies v2,η 0 g(v, η) dv > 0. Hence there exists a traveling wave solution to the parabolic equation However, from (3.7), we deduce that for ε small enough, we have g(v, w ε ) ≥ g(v, η). Hence v ε is a super-solution for the parabolic equation (3.8) (provided the initial conditions are chosen correctly). This is not possible since it is a stationary solution and (3.8) admits traveling waves solutions with positive velocity. Thus, we obtain a contradiction. We deduce that lim L→L * M * (L) = +∞.
Thus, since g is decreasing with respect to its second variable, for any v, g v, M∞ It is well-known that there exists an invading traveling wave solution, denoted v ∞ , of the equation Since, by the comparison principle, for any M ≤ M ∞ and any L > 0, the solution of (1.1), with the same initial data, is such that ). Then, on the interval By conservation of energy, the solutions of this system are such that g(s, 0) ds = 0, for x ∈ ( 3L 4 , +∞), Hence, there exists a solution to (3.9) if there exists ν 1 < ν 2 such that 0 < ν 1 < θ c < ν 2 < 1, and g is decreasing with respect to its second variable. Finally, ( v, W ) is a barrier and we have constructed a barrier for any M > M ∞ .

Application to the sterile insect technique
In this section, we will apply previous results to system (2.5)-(2.6) which has been introduced to model the sterile insect technique. We first show that the reaction term fulfills the hypotheses of Theorem 2.3. Then, we describe the barriers.

Bistability
Let us denote g the right hand side of (2.5). We have We display in Figure 2 the function F → g(F, M s ) for some values 0 < M 1 s < M 2 s . In this case, it is bistable for 0 and M 1 s 1 and mono-stable for M 2 s . The parameter values we will use in this work are given in Table 1 below and are extracted from [24]. Clearly, with these parameters values, the function g(·, 0) is bistable.

We have
1. Differentiating (4.1) with respect to M s and disregarding all positive factors, we have that the sign of ∂g ∂Ms depends on the sign of (β(F + γM s ) + 1)e −β(F +γMs) − 1. Putting x = β(F + γM s ), h(x) = (x + 1)e −x − 1, we have that h(x) < 0 for all x > 0, which proves that ∂g ∂Ms (F, M s ) < 0 for all F > 0, M s ≥ 0. 2. The sign of g(F, M s ) depends on the numerator We compute also Since for all M s ≥ 0, lim Since ∂g ∂F is continuous and has finite limit when F → 0 + and F → +∞, then ∂g ∂F is bounded for every F ≥ 0 and M s ≥ 0.

Wave-blocking region
This part is devoted to the proof of Theorem 2.5.
Proof of Theorem 2.5. Under assumption (4.2) the function F → g(F, 0) is bistable (see Lem. 4.2) and thus, it admits three nonnegative roots 0 < F − < F + . If Let us assume now that (4.2) holds and that . This means that no matter how many sterile males we introduce in this interval, we cannot make F tend to zero asymptotically faster than e −µ F t . This result can be seen as a consequence of the fact that introducing the sterile males only perturbs the reproduction and does not impact the mortality of the females. Therefore, the females cannot decrease with a mortality rate larger than µ F which is, in this simple model, the mortality rate they would have in isolation. However, when we introduce a large number of sterile males, we can choose our linear bound α ≈ µ F , which reflects that the probability of females mating with fertile males, and laying viable eggs, decreases, and therefore the female population decreases at a rate close to µ F .  [24].

Numerical simulations
To illustrate our theoretical results, we provide some numerical simulations for the mathematical model used for the sterile insect technique (2.5)-(2.6). The function g defining the right hand side is given in (4.1). The values of the numerical parameters are taken from [24] and are given in Table 1. The carrying capacity is computed as in [24] by taking the value of the number of males at equilibrium (without performing the SIT). We may assume that when no control technique is performed, the system (2.2) is at equilibrium with value (E, M , F ) which verifies (see system (2.2)) Once the value of the equilibrium is determined, it is easy to compute the carrying capacity since it should be such that g(F , 0) = 0. Therefore, with the expression of g in (4.1) we get For instance, for a value M = 5106 km −2 (see [24]) we find K = 20793 km −2 .
With these parameter values, we first verify easily that the conditions on the parameters in Theorem 2.5 are satisfied. Then, we perform numerical simulations of system (2.5)-(2.6) in one space dimension. More precisely, we consider a domain of width 50 km discretized by 1000 points. We discretize this system with a finite difference scheme where the right hand side is treated explicitly. Since, for numerical purposes we are in a bounded domain, we impose Neumann boundary conditions at the boundary of the numerical domain.
We show in Figure 3 the density of the female mosquitoes F computed by this discetization of system (2.5)-(2.6) for L = 10 km, and for the two values M 0 = 10 000 km −2 (Fig. 3-left) and M 0 = 30 000 km −2 (Fig. 3-right). We observe that there is an invading wave of mosquitoes coming from the left of the domain. When performing the SIT with the value M 0 = 10 000 km −2 in the domain [0, L], the wave is slowing down but can cross this region and continue to invade the whole domain. Increasing M 0 sufficiently, we notice that the invasion seems to be blocked, as illustrated in Figure 3-right for M 0 = 30 000 km −2 .
In order to illustrate Theorems 2.3 and 2.5, we display in Figure 4 the curve L → M * (L) obtained for the function g used to model the sterile insect technique (see (4.1)) with the values of Table 1. To obtain this curve, we test numerically, for a given value of the width L of the band, for which value of M 0 the invading wave is blocked at least numerically. We consider that the invading wave is blocked when the difference between two consecutive time iterations is smaller than a given threshold (here we choose 10 −6 ). We observe that this curve L → M * (L) verifies the properties stated in Theorem 2.3.
Moreover, we also compute the limit value M ∞ , defined in Theorem 2.3, by a dichotomy method. More precisely, for a given value of M 0 , we compute the maximal root F + of the function g(·, M0 µs ) thanks to a Newton method. Then we compute the integral F+ 0 g(F, M0 µs ) dF by a trapezoidal rule. By dichotomy, we obtain an approximation of M ∞ such that this latter integral vanishes. With the numerical values at hand, we find M ∞ = 25083, 58 km −2 . We verify in Figure 4 that the curve seems to converge to this asymptotic value M ∞ .

Conclusion
We provide in this paper a rigorous study of the feasability to block an invading species by performing a sterile insect technique in a band, which may be seen as a saniraty cordon to protect an area. The main result of this paper may be summarized as follows : for any width of the band, there exists a minimal value of the number of sterile insects to release continuously in this band to block the invasion of the species. We provide some numerical simulations for this strategy for the case of mosquito invasion.
This study raises several questions. One question is directly related to our observation from numerical experiments. Indeed it will be interesting to study the slowing down of the wave since for practical situations it may be interesting to only slow down the invasion to protect an area for several years. Another extension of this work is the so-called "rolling carpet" strategy which consists in moving the band where the blocking occurs in the opposite direction of the wave to eradicate the insect species in a whole area. Such a study has been performed in [2]. Moreover, we have assumed in this study that all biological parameters are constant in space and time. Obviously, this is a rough approximation and a more precise study should take into account the variation due to the spatial heterogeneity and due to the seasonality.
This journal is currently published in open access under a Subscribe-to-Open model (S2O). S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

Please help to maintain this journal in open access!
Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting subscribers@edpsciences.org More information, including a list of sponsors and a financial transparency report, available at: https://www.edpsciences.org/en/maths-s2o-programme