Martin Strugarek

Hervé Bossin

Yves Dumont

On the use of the sterile insect technique or the incompatible insect technique to reduce or eliminate mosquito populations

Keywords: Vector control, elimination, sterile insect technique, monotone dynamical system, basin of attraction, numerical simulation, Aedes spp MSC Classification: 34A12, 34C12, 34C60, 34K45, 92D25

Vector control is critical to limit the circulation of vector-borne diseases like chikungunya, dengue or zika which have become important issues around the world. Among them the Sterile Insect Technique (SIT) and the Incompatible Insect Technique (IIT) recently aroused a renewed interest. In this paper we derive and study a minimalistic mathematical model designed for Aedes mosquito population elimination by SIT/IIT. Contrary to most of the previous models, it is bistable in general, allowing simultaneously for elimination of the population and for its survival. We consider different types of releases (constant, periodic or impulsive) and show necessary conditions to reach elimination in each case. We also estimate both sufficient and minimal treatment times. Biological parameters are estimated from a case study of an Aedes polynesiensis population, for which extensive numerical investigations illustrate the analytical results. The applications of this work are two-fold: to help identifying some key parameters that may need further field investigations, and to help designing release protocols.

Introduction

Sterile insect technique (SIT) is a promising technique that has been first studied by E. Knipling and collaborators and first experimented successfully in the early 50's by nearly eradicating screwworm population in Florida. Since then, SIT has been applied on different pest and disease vectors, like fruit flies or mosquitoes (see [START_REF] Dyck | The Sterile Insect Technique, Principles and Practice in Area-Wide Integrated Pest Management[END_REF] for an overall presentation of SIT and its applications). The classical SIT relies on the mass releases of males sterilized by ionizing radiations. The released sterile males transfer their sterile sperms to wild females, which results in a progressive reduction of the target population. For mosquito control in particular, new approaches stemming from SIT have emerged, namely the RIDL technique, and the Wolbachia technique. Wolbachia is a bacterium that infects many Arthropods, and among them some mosquito species in nature. It was discovered in 1924 [START_REF] Hertig | Studies on rickettsia-like micro-organisms in insects[END_REF]. Since then, particular properties of these bacteria have been unveiled. One of of these properties is particularly useful for vector control: the cytoplasmic incompatibility (CI) property [START_REF] Sinkins | Wolbachia and cytoplasmic incompatibility in mosquitoes[END_REF][START_REF] Bourtzis | Wolbachia-based technologies for insect pest population control[END_REF]. CI can serve two different control strategies:

• Incompatible Insect Technique (IIT): the sperm of W-males (males infected with CI-inducing Wolbachia) is altered so that it can no longer successfully fertilize uninfected eggs. This can result in a progressive reduction of the target population. Thus, when only W-males are released the IIT can be seen as classical SIT. This also supposes that releases are made regularly until extinction is achieved (when possible) or until a certain threshold is reached (in order to reduce exposure to mosquito bites and the epidemiological risk).

• Population replacement: when males and W-females are released in a susceptible (uninfected) population, due to CI, W-females will typically produce more offspring than uninfected females. Because Wolbachia is maternally inherited this will result in a population replacement by Wolbachia infected mosquitoes (such replacements or invasions have been observed in natural population, see [START_REF] Rasgon | Wolbachia and cytoplasmic incompatibility in the california culex pipiens mosquito species complex: parameter estimates and infection dynamics in natural populations[END_REF] for the example of Californian Culex pipiens). It has been showed that this technique may be very interesting with Aedes aegypti, shortening their lifespan (see for instance [START_REF] Schraiber | Constraints on the use of lifespan-shortening wolbachia to control dengue fever[END_REF]), or more importantly, cutting down their competence for dengue virus transmission [START_REF] Moreira | A wolbachia symbiont in aedes aegypti limits infection with dengue, chikungunya, and plasmodium[END_REF]. However, it is also acknowledged that Wolbachia infection can have fitness costs, so that the introgression of Wolbachia into the field can fail [START_REF] Schraiber | Constraints on the use of lifespan-shortening wolbachia to control dengue fever[END_REF].

Based on these biological properties, classical SIT and IIT (see [START_REF] Dufourd | Modeling and simulations of mosquito dispersal. the case of aedes albopictus[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and Aedes albopictus[END_REF][START_REF]Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control[END_REF][START_REF] Li | Modelling releases of sterile mosquitoes with different strategies[END_REF][START_REF] Huang | Modelling and analysis of impulsive releases of sterile mosquitoes[END_REF] and references therein) or population replacement (see [START_REF] Farkas | Structured and unstructured continuous models for wolbachia infections[END_REF][START_REF] Fenton | Solving the Wolbachia paradox: modeling the tripartite interaction between host, Wolbachia, and a natural enemy[END_REF][START_REF] Schraiber | Constraints on the use of lifespan-shortening wolbachia to control dengue fever[END_REF][START_REF] Hughes | Modeling the Use of Wolbachia to Control Dengue Fever Transmission[END_REF][START_REF] Farkas | Modelling wolbachia infection in a sex-structured mosquito population carrying west nile virus[END_REF][START_REF] Nadin | Hindrances to bistable front propagation: application to wolbachia invasion[END_REF][START_REF] Strugarek | Quantifying the survival uncertainty of Wolbachia-infected mosquitoes in a spatial model[END_REF] and referencs therein) have been modeled and studied theoretically in a large number of papers in order to derive results to explain the success or not of these strategies using discrete, continuous or hybrid modeling approaches, temporal and spatio-temporal models. Recently, the theory of monotone dynamical systems has been applied efficiently to study SIT [START_REF] Anguelov | Mathematical modeling of sterile insect technology for control of anopheles mosquito[END_REF] or population replacement [START_REF] Sallet | Monotone dynamical systems and some models of wolbachia in aedes aegypti populations[END_REF][START_REF] Bliman | Ensuring successful introduction of wolbachia in natural populations of aedes aegypti by means of feedback control[END_REF] systems.

Here, we derive a monotone dynamical system to model the release and elimination process for SIT/IIT. The analytical study of this model is complemented by a detailed parametrization to describe real-life settings, and a thorough investigation of numerical scenarios.

The outline of the paper is as follows. First, we explain in Section 1 the biological situation we consider and the practical questions we want to answer, namely: how to quantify the release effort required to eliminate an Aedes population using SIT/IIT, with particular emphasis on the timing and size of the releases. We also justify our modeling choices and give value intervals deduced from experimental results for most biological parameters in Table 1. Then, we perform the theoretical analysis of a simple, compartimentalized population model featuring an Allee effect and a constant sterilizing male population in Section 2. Proposition 2.4 gives the bistable asymptotic behavior of the system, and introduces the crucial separatrix between extinction and survival of the population. We also provide analytical inequalities on the entrance time of a trajectory into the extinction set (Proposition 2.6), which is extremely useful to understand what parameters are really relevant and how they interact. We then analyze the model as a control system, after adding a release term. Finally, Section 3 exposes numerical investigations of the various models, and applies them to a specific case study (a pilot field trial led by one of the authors).

In general, all mathematical results are immediately interpreted biologically. To keep the exposition as readable as possible, we gather all technical developments of the proofs into Appendices.

1 Modeling and biological parameter estimation

Modeling context

Our modeling effort is oriented towards an understanding of large-scale time dynamics of a mosquito population in the Aedes genus exposed to artificial releases of sterilizing males. These males can be either sterilized by irradiation (Sterile Insect Technique approach) or simply have a sterile crossing with wild females due, for instance, to incompatible strains of Wolbachia bacteria (Incompatible Insect Technique approach). In either techniques (SIT or IIT), the released males are effectively sterilizing the wild females they mate with.

Eggs from mosquitoes of various species in the Aedes genus resist to dessication and can wait for months before hatching. Due to rainfall-dependency of natural breeding sites availability, this feature allows for maintaining a large quiescent egg stock through the dry season, which triggers a boom in mosquito abundance when the rainy season resumes. For the populations we model here, natural breeding sites are considered to be prominent, and therefore it is absolutely necessary that our models take the egg stock into account.

We use a system biology approach to model population dynamics. In the present work we neglect the seasonal variations and assume all biological parameters to be constant over time.

Our first compartmental model features egg, larva, adult male and adult female (fertile or sterile) populations. Most transitions between compartments are assumed to be linear. Only three non-linear effects are accounted for.

First, the population size is bounded due to an environmental carrying capacity for eggs, which we model by a logistic term. Secondly, the sterilizing effect creates two sub-populations among inseminated females. Some are inseminated by wild males and become fertile while the others are inseminated by sterilizing males and become sterile. Hence the relative abundance (or more precisely the relative mating power) of sterilizing males with respect to wild males must appear in the model, and is naturally a nonlinear ratio. Many other parameters may interfere with the mating process for Aedes mosquitoes, but this process is not currently totally understood in particular from the male point of view [START_REF] Lees | Review: Improving our knowledge of male mosquito biology in relation to genetic control programmes[END_REF][START_REF] Oliva | Male reproductive biology of aedes mosquitoes[END_REF], and we stick here to the simplest possible modeling. Thirdly, as a result of sterilizing matings, we expect that the male population can drop down to a very low level. We introduce an Allee effect which come into play in this near-elimination regime. This effect reduces the insemination rate at low male density, as a consequence of difficult mate-finding. It can also be interpreted as a quantification of the size of the mating area relative to the total size of the domain, and compensates in some ways the intrinsic limitations of a mean-field model for a small and dispersed population (cf. [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF] and see Remark 1.1). Indeed, we model here temporal dynamics by neglecting spatial variations and assuming homogeneous spatial distribution of the populations. In nature, the distribution of Aedes mosquitoes is mostly heterogenous, depending on environmental factors such as vegetation coverage, availability of breeding containers and blood hosts. The proposed simplified homogenous model will thus be exposed to potential criticism.

Models and their basic properties

We denote by E the eggs, L the larvae, M the fertile males, F the fertile females and F st the sterile females (either inseminated by sterilizing males or not inseminated at all, due to male scarcity). The time-varying sterilizing male population is denoted M i . We use Greek letters µ for mortality rates, ν for transition rates and denote fecundity by b (viable eggs laid per female and per unit of time) and egg carrying capacity by K. The full model reads:

                                 dE dt = bF (1 - E K ) -( ν E + µ E )E, dL dt = ν E E -(ν L + µ L )L, dM dt = (1 -r)ν L L -µ M M, dF dt = rν L L(1 -e -β(M +γ i M i ) ) M M + γ i M i -µ F F, dF st dt = rν L L e -β(M +γ i M i ) + γ i M i M + γ i M i (1 -e -β(M +γ i M i ) ) -µ F F st . (1)
Dynamics of the full system (1) is not different from that of the following simplified, threepopulations system. We only keep egg, fertile and sterilizing male, and fertile female populations. The value of the hatching parameter ν E must be updated to take into account survivorship and development time in the larval stage.

               dE dt = bF (1 - E K ) -(ν E + µ E )E, dM dt = (1 -r)ν E E -µ M M, dF dt = rν E E(1 -e -β(M +γ i M i ) ) M M + γ i M i -µ F F. (2)
The following straightforward lemma means that (1) and ( 2) are well-suited for population dynamics modeling since all populations, in these systems, remain positive and bounded.

Lemma 1 Let M i be a non-negative, piecewise continuous function on R + . The solution to the Cauchy problems associated with (1), (2) and non-negative initial data is unique, exists on R + , is continuous and piecewise continuously differentiable. This solution is also forward-bounded and remains non-negative. It is positive for all positive times if F (0) or E(0) (or also L(0) in the case of (1)) is positive.

In addition, these systems are monotone in the sense of the monotone systems theory (see [START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF]).

Lemma 2

The system (2) is monotone on the set E 3 := {E ≤ K} ⊂ R 3 + for the order induced by R 3 + and the restriction of system (1) to the four first coordinates (omitting F st , which does not appear in any other compartment) is monotone on the set

E 4 := {E ≤ K} ⊂ R 4
+ for the order induced by R 4 + . Moreover, E 3 (respectively E 4 ) is forward invariant for (2) (respectively for the restriction of (1) to the four first coordinates), and any trajectory enters it in finite time.

Proof. We compute the Jacobian matrix of the system (2):

J =   -bF K -(ν E + µ E ) 0 b(1 -E K ) (1 -r)ν E -µ M 0 rν E (1 -e -β(M +γ i M i ) ) M M +γ i M i rν E E M +γ i M i βM e -β(M +γ i M i ) + (1 -e -β(M +γ i M i ) ) γ i M i M +γ i M i -µ F   .
It has non-negative extra-diagonal coefficients on E 3 , which proves that the system is indeed monotone on this set. In addition, if

E(t 0 ) > K then let T [t 0 ] := {t ≥ t 0 , ∀t ′ ∈ [t 0 , t), E(t ′ ) > K} ⊂ R. Let T + [t 0 ] := sup T [t 0 ]. For any t ∈ T [t 0 ] we have Ė(t) ≤ -(ν E + µ E )E(t).
Hence by integration we find that T + [t 0 ] ≤ t 0 + 1 ν E +µ E log(K/E(t 0 )) < +∞, which proves Lemma 2 (the proof being similar for the claims on (1)).

Remark 1.1 The Allee effect term 1exp(-βM ) can also be interpreted in the light of [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF]. This is the probability that an emerging female finds a male to mate with in her neighborhood.

Using a "mean-field" model of ordinary differential equations here is certainly debatable, since in the case of population extinction the individuals may eventually be very dispersed, and heterogeneity would play a very important role. However, we think that getting a neat mathematical understanding of the simplest system we study here is a necessary first step before moving to more complex systems. The Allee term compensates, as far as the qualitative behavior is concerned, what the model structurally lacks. Here, we are able to perform proofs and analytical computations. This gives a starting point for benchmarking what to expect as an output of release programs using sterilizing males, according to the models. For numerical simulations, we use experimental (lab and field) values of the biological parameters in (1)- [START_REF] Anguelov | On nonstandard finite difference schemes in biosciences[END_REF]. We consider specifically a population of Aedes polynesiensis in French Polynesia which has been studied in [START_REF] Jachowski | Filariasis in american samoa. y. bionomics of the principal vector, aedes polynesiensis marks[END_REF][START_REF] Suzuki | Breeding habits of vector mosquitoes of filariasis and dengue fever in western samoa[END_REF][START_REF] Rivière | Ecologie de Aedes (Stegomyia) polynesiensis, Marks[END_REF], and more recently in [START_REF] Chambers | Male mating competitiveness of a wolbachia-introgressed aedes polynesiensis strain under semi-field conditions[END_REF][START_REF] Hapairai | Effect of temperature and larval density on aedes polynesiensis (diptera: Culicidae) laboratory rearing productivity and male characteristics[END_REF][START_REF] Hapairai | Population studies of the filarial vector aedes polynesiensis (diptera: Culicidae) in two island settings of french polynesia[END_REF][START_REF] Hapairai | Studies on Aedes polynesiensis introgression and ecology to facilitate lymphatic filariasis control[END_REF].

Parameter estimation from experimental data

Values of most parameters are given in Table 2, and are deduced from experimental data gathered in Table 1. Some data come from unpublished results obtained at Institut Louis Malardé during the rearing of Aedes polynesiensis for a pilot IIT program. They are labelled as "Production data (ILM)". Note that we do not give values for β and ν E because they are very hard to estimate. Ongoing experiments of one of the author may help approximating them in the future for this Aedes polynesiensis population. Finally when it exists, we use the knowledge about population size (male and female) granted by mark-release-recapture experiments to adjust the environmental carrying capacity K for population and season.

≤ i ≤ d, x i ≤ y i ,
• x < y if and only if x ≤ y and x = y,

• x ≪ y if and only if for all 1 ≤ i ≤ d, x i < y i .

Constant incompatible male density

First we study system [START_REF] Anguelov | On nonstandard finite difference schemes in biosciences[END_REF] with constant incompatible male density M i (t) ≡ M i .

We introduce the three scalars

N := brν E µ F (ν E + µ E ) , λ := µ M (1 -r)ν E K , ψ := λ β (3)
and define the function f : R 2 + → R, with the two parameters N and ψ:

f (x, y) := x(1 -ψx)(1 -e -(x+y) ) - 1 N (x + y). (4)
The two aggregated numbers, N and ψ essentially contain all the information about system (2): N is the classical basic offspring number, ψ is the ratio between the typical male population size at which the Allee effect comes into play and the male population size at wild equilibrium, as prescribed by the egg carrying capacity.

The ODE system (2) has simple dynamical properties because it is monotone and we can count its steady states and even know their local stability. Let M i ≥ 0. It is straightforward to show that system (2) always admits a trivial steady-state (0, 0, 0) and eventually one (at least) non-trivial steady state (E * , M * , F * ) ∈ R 3 + solution of

E = b ν E + µ E F (1 - E K ), E = µ M (1 -r)ν E M, F = rν E µ F E(1 -e -β(M +γ i M i ) ) M M + γ i M i .
Using the first two equation into the third one yields

µ F (ν E + µ E ) brν E (M + γ i M i ) = M (1 - µ M (1 -r)ν E K M )(1 -e -β(M +γ i M i ) ), from which we deduce        E * = KλM * , F * = K(ν E + µ E ) b λM * 1 -λM * , f (βM * , γ i βM i ) = 0.
Hence for a given value M i ≥ 0, the number of steady states of ( 2) is equal to the number of positive solutions M * to f (βM * , βγ i M i ) = 0, plus 1. The trivial steady state (0, 0, 0) is also locally asymptotically stable (LAS). The following lemma give us additional informations about the positive steady state(s):

Lemma 3 Assume N > 4ψ. Let θ 0 ∈ (0, 1) be the unique solution to 1 -θ 0 = -4ψ
N log(θ 0 ), and

M crit i := 1 γ i β max θ∈[θ 0 ,1] -log(θ) - 1 2ψ 1 -1 + 4ψ N log(θ) 1 -θ . If M crit i > 0 then (2) has: • 0 positive steady state if M i > M crit i ,
• 2 positive steady states

E -≪ E + if M i ∈ [0, M crit i ),
• 1 positive steady state

E if M i = M crit i .
In addition, E -is unstable and E + is locally asymptotically stable. If M crit i < 0 then (2) has no positive steady state, and if M crit i = 0 then there exists a unique positive steady state. In particular, if N ≤ 1 then M crit i < 0. On the contrary, if N ≤ 4ψ then there is no positive steady state.

Proof. Let us give a quick overview of the remainder of the proof, which is detailed in Appendix A, page 23. We are going to study in details the solutions (x, y) to f (x, y) = 0. First, we prove that x < 1/ψ. Then, we check that for any y > 0, x → f (x, y) is either concave or convex-concave. In addition, it is straightforward that f (0, y) < 0 and lim x→+∞ f (x, y) = -∞, so that for any y > 0, we conclude that there are either 0, 1 or 2 real numbers x > 0 such that f (x, y) = 0.

Then, we introduce ξ = 4ψ/N . In fact, in order to determine (x, y) ∈ R 2 + such that f (x, y) = 0 we can introduce θ = e -(x+y) and then check easily that y = h ± (θ), where

h ± (θ) = -log(θ) - 1 2ψ ± 1 2ψ 1 + ξ log(θ) 1 -θ . ( 5
)
Let θ 0 (ξ) be the unique solution in (0, 1) to 1θ 0 (ξ) = -ξ log(θ 0 (ξ)), and

α crit (ξ, N ) := max θ∈[θ 0 (ξ),1] -log(θ) - 1 2ψ 1 -1 + ξ log(θ) 1 -θ . (6)
Collecting the previous facts, and studying the function h ± (see Appendix A.2, page 24), we can prove that the next point of Lemma 3 holds with the threshold

M crit i = N 4ψβγ i α crit (ξ, N ).
We remark that if N ≤ 1 then it is easily checked that M crit i < 0, using the fact that if α ∈ (0, 1)

then √ 1 -α ≤ (1 -α)/2. If θ ∈ (θ 0 , 1) then 4ψ log(θ) N (1-θ) < 1, and therefore -log(θ) - 1 4ψ 1 -1 + 4ψ N log(θ) 1 -θ ≤ -log(θ) 1 - 1 N - 1 4ψ < 0.
In the final part of the proof, we show that 0 is always locally stable and then treat separately the cases M i = 0 and M i > 0, showing that, when they exist, the greater positive steady state is locally stable while the smaller one is unstable.

Remark 2.1 In Lemma 3, the condition to have at least one positive equilibrium, N > 4ψ, is very interesting and particularly makes sense when rewritten as

N λ > 4 β
. Indeed N λ can be seen as the theoretical male progeny at next generation, starting from wild equilibrium. If this amount is large enough (larger than some constant times the population size at which the Allee effect comes into play) then the population can maintain. In any case, if this condition is not satisfied, then the population collapses. For the population to maintain: either the fitness is good and thus N is very large, or the probability of one female to mate is high and thus 1/β is small. However, whatever the values taken by N and β, if, for any reason, the male population at equilibrium decays, the population can be controlled and possibly collapses.

Remark 2.2 If β is not too small, then the "wild" steady state is approximately given by

M * (M i = 0) ≃ 1 λ (1 -1 N ) and the critical sterilizing level is approximately M crit i ≃ y = N 4λγ i 1 -1 N 2 (see the definition in Appendix A, in particular we know that M crit i ≤ y).
As a consequence, the target minimal constant density of sterilizing males compared to wild males in order to get unconditional extinction (i.e. to make (0, 0, 0) globally asymptotically stable, see Proposition 2.3, page 9) is well approximated by the simple formula

ρ * := M crit i M * (M i = 0) ≃ N -1 4γ i .
With the values from Tables 1 and2, for γ i = 1 (this means that introduced male are as competitive as wild ones for mating with wild females), we find

ρ * ∈ 7.46 • 0.46 • ν E 4 • 0.046 • (ν E + 0.046) -0.25, 14.85 • 0.48 • ν E 4 • 0.033 • (ν E + 0.023) -0.25
For instance, if ν E = 0.01 then this interval is (3.5, 22, 7), if ν E = 0.05 then this interval is (10.6, 51.7) and if ν E = 0.1 then this interval is (14.1, 61.4). As ν E goes to +∞, the interval goes to (20.7, 75.7). This example agrees with standard SIT Protocol that indicates to release at least 10 times more sterile males than wild males, recalling that here we deal with a highly reproductive species (with the above values, the lowest estimated basic reproduction number is 14.9, obtained for ν E = 0.01).

Asymptotic dynamics are easily deduced from the characterization of steady states and local behavior of the system (Lemma 3), because of the monotonicity (see [START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF]).

Proposition 2.3 If (2) has only the steady state (0, 0, 0) then it is globally asymptotically stable.

If there are two other steady states E -≪ E + then almost every orbit converges to E + or (0, 0, 0). Let K + := [(0, 0, 0), E + ]. The compact set K + is globally attractive and positively invariant. The basin of attraction of (0, 0, 0) contains [0, E -) and the basin of attraction of E + contains (E -, ∞). Now that we have established that the system is typically bistable, the main object to investigate is the separatrix between the two basins of attraction. This is the aim of the next proposition.

Proposition 2.4 Assume M crit i > 0 and M i ∈ [0, M crit i ).
Then there exists a separatrix Σ ⊂ R 3 + , which is a sub-manifold of dimension 2, such that for all X = Y ∈ Σ, X ≤ Y and Y ≤ X, and for all X ∈ Σ, X 0 > X implies that X(t) converges to E + , and X 0 < X implies that X(t) converges to 0. In particular, E -∈ Σ.

Let

Σ + := X ∈ R 3 + , ∃ X ∈ Σ, X > X and Σ -:= X ∈ R 3 + , ∃ X ∈ Σ, X < X . Then R 3 + = Σ -∪ Σ ∪ Σ + , Σ +
is the basin of attraction of E + and Σ -is the basin of attraction of 0.

In addition, there exists E M , F M > 0 such that

Σ -⊂ X ∈ R 3 + , X 1 ≤ E M , X 3 ≤ F M .
Remark 2.5 In order to reach extinction, the last point of Proposition 2.4 states that both egg and fertile female populations must stand simultaneously below given thresholds. This obvious fact receives here a mathematical quantification. With simple words: no matter how low the fertile female population F has dropped, if there remains at least E M eggs then the wild population will recover.

Proof. [Proposition 2.4] We state a preliminary fact: For all

v 0 ∈ {v ∈ R 3 + , ∀i, v i > 0, i v i = 1} =: S 2
+ , there exists a unique ρ 0 (v 0 ) such that the solution to (2) with initial data ρv 0 converges to 0 if ρ < ρ 0 (v 0 ) and to

E + if ρ > ρ 0 (v 0 ).
This fact comes from the strict monotonicity of the system, and from the estimate

ρ 0 (v 0 ) ≤ max i v 0 i (E -) i < +∞, combined with Proposition 2.3. Then we claim that Σ = {ρ 0 (v 0 )v 0 , v 0 ∈ S 2 + }.
The direct inclusion is a corollary of the previous fact. The converse follows from the fact that Σ ± , being the basins of attraction of attracting points, are open sets.

The remainder of the proof consists of a simple computation showing that if F 0 or E 0 is large enough then for some t > 0 we have (E, M, F )(t) > E -. In details, we can prove that if F 0 is large enough then for any E 0 , M 0 and ǫ > 0, we can get

E(s) ≥ (1-ǫ)K for s ∈ (t 0 (ǫ, E 0 , F 0 ), t 1 (ǫ, E 0 , F 0 )),
where t 0 is decreasing in F 0 and t 1 is increasing in F 0 and unbounded as F 0 goes to +∞. Then, if E > (1ǫ)K for ǫ small enough on a large enough time-interval, we deduce

M (t) > (1 -ǫ) 2 (1 -r) ν E
µ M K for some t > 0. Upon choosing ǫ small enough and F 0 large enough we finally get (E, M, F )(t) > E -. The scheme is similar when taking E 0 large enough.

At this stage, we know that starting from the positive equilibrium, and assuming that the population of sterile males M i is greater than M crit i , the solution will reach the basin of attraction of the trivial equilibrium in a finite time, τ (M i ). We obtain now quantitative estimates on the duration of this transitory regime. Rigorously, we define

τ (M i ) := inf t ≥ 0, (E, M, F )(t) ∈ Σ -(M i = 0), where (E, M, F )(0) = E + (M i = 0) and (E, M, F ) satisfies (2) . (7)
We obtain simple upper and lower bounds for τ (M i ) in terms of various parameters:

Proposition 2.6 Let M i > M crit
i , and Z = Z(ψ) be the unique real number in (0, 1 2ψ ) such that

e -Z = ψ 1 + ψ -ψZ ,
and

Z 0 := 1 + ψ -ψZ. Then τ (M i ) ≥ 1 µ F log 1 + N 2 (1 -ψZ) 3 ψZZ 2 0 - N (1 -ψZ) ψZZ 0 . ( 8
)
Let σ = sgn(ν E + µ E -µ F ), σ E := µ M /(ν E + µ E ) and σ F := µ M /µ F . If ǫ := M * + M * + +M i < 1/N , let g(ε) := 1 + 4N σ E σ F ε (σ F -σ E ) 2 . Assume that σ F , σ E > 1, g(ε)σ(σ F -σ E ) < max (2N -1)σ F + σ E , (2σ E -1)σ F , (σ F -1)(σ E -1) > εN .
Then

τ (M i ) ≤ 2σ E µ F σ F + σ E -g(ǫ)σ(σ F -σ E ) log N -1 ψ (N -1)σ F + 1 -εN (σ F -1)(σ E -1) -εN + σ E σ F g(ε)σ(σ F -σ E ) + (2N -1)σ F + σ E 2σ E σ F -(σ E + σ F ) + σ(σ F -σ E )g(ε) g(ε)σ(σ F -σ E ) . (9)
Proof. The proof relies on explicit computation of sub-and super-solutions, detailed in Appendix B.

Remark 2.7

The dependency in ψ of Proposition 2.6's upper estimate on τ is approximately equal to 1 min(ν E +µ E ,µ F ) . One order of magnitude of ψ (the ratio between the wild population size and the Allee population size) therefore typically corresponds to the maximum of one adult female and one egg lifespan in terms of release duration needed to get extinction. Remark 2.8 At this stage, we obtain an analytic upper bound only in the case of massive releases (ǫ small enough). A more refined upper bound could theoretically be obtained, see the derivation in Appendix B, in particular Lemma 12.

Adding a control by means of releases

In a slightly more realistic model, the level of sterilizing male population should vary with time, depending on the releases t → u(t) ≥ 0 and on a fixed death rate µ i . This model reads

                         dE dt = bF (1 - E K ) -(ν E + µ E )E, dM dt = (1 -r)ν E E -µ M M, dM i dt = u(t) -µ i M i , dF dt = rν E E(1 -e -β(M +γ i M i ) ) M M + γ i M i -µ F F. (10)
In [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF], the number of sterilizing males released between times t 1 and t 2 > t 1 is simply equal to

t 2 t 1 u(t)dt. First, if the release is constant, say u(t) ≡ u 0 , then M i (t) = e -µ i t M 0 i + u 0 µ i (1-e -µ i t
). The special case M 0 i = u 0 µ i leads back to system (2), with M i ≡ M 0 i . For general M 0 i ≥ 0, we notice that M i (t) converges to u 0 µ i as t goes to +∞.

Proposition 2.9 Assume u(t) ≡ u 0 . If u 0 > µ i M crit i (defined in Lemma 3) then 0 is globally asymptotically stable. If u 0 < µ i M crit i , then there exists open sets Σ -(u 0 ), Σ + (u 0 ) ⊂ R 4 +
, respectively the basins of attraction of 0 and E + (defined for (2) with M i = u 0 µ i ), separated by a set Σ(u 0 ) which enjoys the same properties as those of Σ(0), listed in Proposition 2.4.

(We do not treat the case u 0 = µ i M crit i ). Proof. Since system [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF] is monotone with respect to the control u (with sign pattern (-, -, -, +)), we can use Lemma 3 and Proposition 2.4 with sub-and super-solution to get this result in a straightforward way.

From now on we will restrict ourselves to (possibly truncated) time-periodic controls, which means that we assume that there exists N r ∈ Z + ∪ {+∞} (the number of release periods), a period T > 0 and a function u

0 : [0, T ] → R + such that u(t) = u 0 (t -nT ) if nT ≤ t < (n + 1)T for some N r > n ∈ Z + , 0 otherwise. (11)
We use the notation u ≡ [T, u 0 , N r ] to describe this control u.

As before, we can compute in case ( 11)

M i (t) = e -µ i t M 0 i + t 0 u(t ′ )e -µ i (t-t ′ ) dt ′ = e -µ i t M 0 i + e µ i (⌊ t T ⌋∧Nr)T -1 e µ i T -1 T 0 u 0 (t ′ )e µ i t ′ dt ′ + t T (⌊ t T ⌋∧Nr) u(t ′ )e µ i t ′ dt ′ (Here, for a, b ∈ Z, we let a ∧ b = min(a, b)).
If N r = +∞, for any u 0 = 0 there exists a unique periodic solution M i , uniquely defined by its initial value

M 0,per i = 1 1 -e -µ i T T 0 u 0 (t ′ )e µ i t ′ dt ′ ,
and which we denote by

M per i [u 0 ].
Lemma 4 Solutions to [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF]

with u ≡ [T, u 0 , +∞] are such that M i converges to M per i [u 0 ],
and the other compartments converge to a solution of

                 dE dt = bF (1 - E K ) -(ν E + µ E )E, dM dt = (1 -r)ν E E -µ M M, dF dt = rν E E(1 -e -β(M +γ i M per i [u 0 ]) ) M M + γ i M per i [u 0 ] -µ F F. ( 12
)
Convergence takes place in the sense that the L ∞ norm on (t, +∞) of the difference converges to 0 as t goes to +∞.

Proof. Convergence of M i is direct from the previous formula. Then, as for Proposition 2.9 the monotonicity of the system implies the convergence.

Let M i [u 0 ] := max M per i [u 0 ] and M i [u 0 ] := min M per i [u 0 ]. Proposition 2.10 If M i [u 0 ] > M crit i
then 0 is globally asymptotically stable for [START_REF] Farkas | Modelling wolbachia infection in a sex-structured mosquito population carrying west nile virus[END_REF]. On the contrary, if M i [u 0 ] < M crit i then (12) has at least one positive periodic orbit. In this case the basin of attraction of 0 contains the interval 0, E -(M i = M i [u 0 ]) , and any initial data above

E + (M i = M i [u 0 ]) converges to X per [u 0 ].
Proof. System ( 12) is a periodic monotone dynamical system. It admits a unique non-negative solution X = (E, M, F ). In fact, we consider the constant sterile population model

                 dE m dt = bF m 1 - E m K -(ν E + µ E )E m , dM m dt = (1 -r)ν E E m -µ M M m , dF m dt = rν E M m M m + M i [u 0 ] 1 -e -β(Mm+γ i M i [u 0 ]) E m -µ F F m . (13)
such that, using a comparison principle, the solution

X m = (E m , M m , F m ) verifies X m ≥ X for all time t > 0.
Thus if X m converges to 0, so will X. The behavior of system (13) follows from the results obtained in the previous section. A sufficient condition to have 0 globally asymptotically stable in ( 12) is therefore given by M per i > M crit i . The remainder of the claim is better seen at the level of the discrete dynamical system defined by [START_REF] Farkas | Modelling wolbachia infection in a sex-structured mosquito population carrying west nile virus[END_REF]. Periodic orbits are in one-to-one correspondence with the fixed points of the monotone mapping Φ[u 0 ] : R 3 + → R 3 + defined as the Poincaré application of (12) (mapping an initial data to the solution at time t = T ). Now, if X * := (E * , M * , F * ) is the biggest (i.e. stable) steady state of (2) at level

M i = M i [u 0 ] < M crit
i , then for any (E, M, F ) ≫ (E * , M * , F * ) and M ′ i ≤ M i , writing the right-hand side as Ψ = (Ψ 1 , Ψ 2 , Ψ 3 ) we have

Ψ 1 (E * , M, F, M ′ i ) > 0, Ψ 2 (E, M * , F, M ′ i ) > 0, Ψ 3 (E, M, F * , M ′ i ) > 0.
In other words, the interval X * , +∞ is a positively invariant set. Therefore, Φ[u 0 ](X * ) > X * . Thus the sequence Φ[u 0 ] k (X * ) k is increasing and bounded in R 3 + : it must converge to some X * > X * . The same reasoning (with reversed inequalities) applies with the sequence starting at the stable equilibrium associated with M i = M i [u 0 ]: it must decrease, and thus converge to some X * ≥ X * .

By our proof we have shown that the open interval E + (M i = M i [u 0 ]), +∞ belongs to the basin of attraction of X per , and we can also assert that

E -(M i = M i [u 0 ]), E + (M i = M i [u 0 ])
belongs to the basin of attraction of X per , while as usual 0, E -

(M i = M i [u 0 ]
) is in the basin of attraction of 0.

By a direct application of the previous results

Lemma 5

If M i [u 0 ] > M crit i
then the control u ≡ [T, u 0 , n] (with n ∈ Z + ) leads to extinction (i.e. the solution with initial data E + goes to 0 as t goes to +∞) as soon as

n ≥ τ (M i [u 0 ]) T . (14)
A special case of ( 10)-( 12) is obtained by choosing u 0 = u ǫ 0 = Λ ǫ 1 [0,ǫ] for some Λ > 0 and letting ǫ go to 0. Then there exists a unique limit as ǫ goes to 0, which is given by the following impulsive differential system derived from [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF]:

                               dE dt = bF (1 - E K ) -(ν E + µ E )E, dM dt = (1 -r)ν E E -µ M M, dM i dt = -µ i M i , M i (nT + ) = M i (nT ) + Λ for n ∈ Z + with 0 ≤ n < N r , dF dt = rν E E(1 -e -β(M +γ i M i ) ) M M + γ i M i -µ F F. (15)
In [START_REF] Hapairai | Effect of temperature and larval density on aedes polynesiensis (diptera: Culicidae) laboratory rearing productivity and male characteristics[END_REF], M i converges to the periodic solution

M imp i (t) := lim ǫ→0 M per i [u ǫ 0 ] = Λe -µ i (t-⌊ t T ⌋T ) 1 -e -µ i T
We can compute explicitly

M imp i := Λe -µ i T 1-e -µ i T and M imp i := Λ 1-e -µ i T
, respectively the minimum and the maximum of M imp i . We also define the following periodic monotone system as a special case of [START_REF] Farkas | Modelling wolbachia infection in a sex-structured mosquito population carrying west nile virus[END_REF]:

                 dE dt = bF (1 - E K ) -(ν E + µ E )E, dM dt = (1 -r)ν E E -µ M M, dF dt = rν E E(1 -e -β(M +γ i M imp i ) ) M M + γ i M imp i -µ F F. (16)
The right-hand side of system ( 15) is locally Lipschitz continuous on R 3 . Thus, using a classic existence theorem (Theorem 1.1, p. 3 in [START_REF] Bainov | Impulsive Differential Equations: Periodic Solutions and Applications[END_REF] ), there exists T e > 0 and a unique solution defined from (0, T e ) → R 3 . Using standard arguments, it is straightforward to show that the positive orthant R 3 + is an invariant region for system [START_REF] Hapairai | Effect of temperature and larval density on aedes polynesiensis (diptera: Culicidae) laboratory rearing productivity and male characteristics[END_REF]. We estimate the (minimum) size of the releases Λ and periodicity T , such that the wild population goes to extinction.

Proposition 2.11 Let

S := 1 -r ν E N 4µ M γ i 1 - 1 N 2 K. If T ≤ 1 µ i log 1 + Λ S ( 17
)
then 0 is globally asymptotically stable in [START_REF] Hapairai | Population studies of the filarial vector aedes polynesiensis (diptera: Culicidae) in two island settings of french polynesia[END_REF].

Condition (17) is equivalent to Λ ≥ S e µ i T -1 .
Proof. We know (see Appendix A and Remark 2.2) that

M crit i ≤ N 4λγ i 1-1 N 2 .
Hence the following is a sufficient condition for global asymptotic stability of 0:

M imp i ≥ N 4λγ i 1 - 1 N 2 = (1 -r) ν E N 4µ M γ i 1 - 1 N 2 K.

That is

Λe -µ i τ 1 -e -µ i τ ≥ (1 -r) ν E N 4µ M γ i 1 - 1 N 2 K,
and the result is proved.

Remark 2.12 As a continuation of Remark 2.2, we note that Proposition 2.11 gives a very simple estimate for the target ratio of sterilizing males per release over initial wild male population as a function of the period between impulsive releases in the form

ρ(T ) := Λ M * (M i = 0) ≃ e µ i T -1 N -1 4γ i .
We can specify Lemma 5 for impulses and combine it with Proposition 2.6 to get a sufficient condition for extinction in the impulsive cases: Proposition 2. [START_REF] Farkas | Structured and unstructured continuous models for wolbachia infections[END_REF] The impulsive control of amplitude Λ > 0 and period

T > 0 satisfying Λ ≥ S e µ i T -1 leads to extinction in n impulses if n ≥ τ (M imp i ) T , where
M imp i = Λe -µ i T 1 -e -µ i T . ( 18
)
3 Numerical study

Numerical method and parametrization

In order to preserve positivity of solutions and comparison principle, we use a nonstandard finitedifferences (NSFD) scheme to integrate the differential systems (see for instance [START_REF] Anguelov | On nonstandard finite difference schemes in biosciences[END_REF] for an overview).

For system [START_REF] Durrett | The importance of being discrete (and spatial)[END_REF], it reads

                           E n+1 -E n Φ(∆t) = bF n S (1 - E n+1 K ) -(ν E + µ E )E n , M n+1 -M n Φ(∆t) = (1 -r)ν E E n -µ M M n , M n+1 i -M n i Φ(∆t) = -µ i M n i + u n , F n+1 -F n Φ(∆t) = rν E M n+1 M n+1 + M n+1 i 1 -e -β(M n+1 +M n+1 i ) E n -µ F F n , ( 19
)
where ∆t is the time discretization parameter,

Φ(∆t) = 1-e -Q∆t Q , Q = max{µ M , µ F , ν E + µ E , µ i } and X n (respectively u n ) is the approximation of X(n∆t) (respectively u(n∆t)) for n ∈ N. Parameter β b r µ E ν E µ F µ M γ i µ i ∆t Value
10 -4 -1 10 0.49 0.03 0.001 -0.25 0.04 0.1 1 0.12 0.1

Table 3: Numerical values fixed for the simulations.

We fix the value of some parameters using the values from Tables 1 and2 (see Table 3). Then, in order to get results relevant for an island of 74 ha with an estimated male population of about 69 ha -1 , we let ν E and β vary, and fix K such that

M * + = 69 • 74 = 5106, that is K = 5106 • µ M (1 -r)ν E (1 - 1 N (1-e -β•5106 ) )
.

Recall that for the choice from The other aggregated value of interest,

ψ = µ M (1-r)ν E βK = N -(1-e -βM * + ) N M * + β
, ranges from 1.4 • 10 -4 to 2, approximately.

All computations were performed using Python programming language (version 3.6.2). The most costly operation was the separatrix approximation, which needed to be done once for each set of parameter values. We first compute points close to the separatrix (see details in Section 3.3), starting from a regular triangular mesh with 40 points on each side, then we reduce the points if any comparable pairs appeared. From these (at most 861) scattered points we build recursively a comparison tree by selecting the point P which minimizes the distance to all other points, and distributing the remaining points into six subtrees, corresponding to each affine orthant whose vertex is P . Each tree was saved using pickle module, and loaded when necessary. This was done to reduce the number of operations for checking if a point is below the separatrix, as this needs to be done several times along each computed trajectory. Indeed, using the fact that two points on the separatrix cannot be related by the partial order, one only needs to investigate 3 of the 6 remaining orthants to determine if the candidate point is below any of the scattered points or not. For any given input of released sterilizing males, the computation of a trajectory ended either when the maximal number of iterations was reached (here, we fixed that value at 3 • 10 5 ) or when it was found below the separatrix, using the comparison tree. Trial CPU times (on a laptop computer with Intel R Core TM i5-2410M CPU @ 2.30GHz x 4 processor) for all these operations are given in Table 4.

Equilibria and effort ratio

We first compute the position of equilibria for a range of values of β and ν E . This enables us to compute the effort ratio ρ * , defined in Remark 2.2 as the ratio between the wild steady state male population M * (M i = 0) and the critical constant value of sterilizing males M crit i necessary in order to make 0 globally asymptotically stable. Values are shown in Table 5. We note that ρ * depends practically only on ν E , because the Allee (with parameter β) does not apply at high population levels. In fact the ratio (and thus the control effort) increases with increasing values of ν E , that favor the maintenance of the wild population (the larger the value of ν E , the larger the value of N and the shorter the period in the eggs compartment).

ν E 0.

Computation of the basin of attraction of 0 for (2)

We start from a regular triangular mesh of the triangle {(E, M, F ) ∈ R 3 + , E + M + F = 1}, with 40 points on each side. Given ǫ > 0, for each vertex V of this mesh we compute λ ∈ (0, +∞) such that λV ∈ Σ -and (1 + ǫ)λV ∈ Σ + . The points λV (which are numerically at distance at most ǫ of the separatrix Σ) are then plotted. Figure 1 is typically the kind of figure that we can draw for each set of parameters. Depending on the parameters values, the basin of attraction of 0 can be tiny, or not. Its shape emphasizes the important role of eggs and, even, males abundance in the maintenance of the wild population. In fact, even if almost all females have disappeared, the control must go on in order to further reduce the stock of eggs before eventually reaching the separatrix.

Constant releases and entrance time into basin

For the same set of parameters as before, we compute the entrance time into the basin of 0.

First, we use Proposition 2.6 to get in Table 6 an underestimation of the entrance time, whatever the releasing effort could be, these entrance times represent the minimal time under which the SIT control cannot be successful (in fact, this under-estimation corresponds to the situation where M i = +∞, that is an infinite releasing effort). Table 7: Entrance time into the basin of 0 (in days) for various values of (ν E , β), with

N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A 10 -
N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A
M i = 1.2M crit i (left), M i = 2M crit i (middle) and M i = 8M crit i (right).
We notice that the entrance times corresponding to the biggest effort ratio are of the same order of magnitude as the analytic under-estimation from formula [START_REF]Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control[END_REF].

Another interesting output of Table 7 is that the release effort ratio is not so important in terms of duration of the control: depending on the values taken by ν E and β, the lowest ratio needs between 4 to 7 more weeks to reach the basin, than the largest ratio. Contrary to what could have been expected, there is no linear relationship. This can be explained by the fact that a female mates only once. Thus if males are in abundance, all females have mated, and then many released males become useless with regards to sterilization. Of course, this has to be mitigated taking into account that our model implicitly assumes a homogeneous distribution, while in real, environmental parameters (like vegetation, climate, etc.) have to be taken into account [START_REF]Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control[END_REF]

i = 1.2M crit i (left), M i = 2M crit i (middle) and M i = 8M crit i (right)
. The total effort ratio in this case is defined as M i /M * + multiplied by µ i times the entrance time, and corresponds to the number of males that should be released at a constant level, divided by the initial male population. Table 8, page 18, clearly emphasizes that a large effort ratio, i.e. φ = 8, means the use (and then the production) of a large number of sterile males with a really small time-saving compared to the case φ = 2. For instance with ν E = 0.05 and β = 10 -2 , the total effort ratio for φ = 8 is approximately 6 times larger than for φ = 2 (24395 against 4059), with a time-saving of 37 days, that is approximately one tenth of the total protocol duration (336 days against 373).

In other words, releasing a large number of sterile males is not necessarily a good strategy, from the economical point of view, but also from the control point of view.

In the next subsection, We consider a more realistic scenario, where sterile males are released periodically and instantaneously (system ( 16)).

Periodic releases

In the case of periodic releases by pulses u = [T, Λδ 0 , ∞], for a given couple (ν E , β) we compute the first time t > 0 such that (E, M, F )(t) is below one point of the previously computed separatrix.

We performed the computations with T ∈ {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, choosing

Λ = K φ(1 -r)ν E N 4µ M 1 - 1 N 2 e µ i T -1 for φ ∈ {1.2, 1.4, 1.6, 1.8, 2, 4, 8}.
For all combinations of (ν E , β), we indicate in Table 9 the maximal and minimal (with respect to (T, φ)) total effort ratio ρ tot defined as the number of released mosquitoes at the time when the basin of 0 is reached, divided by the initial male population that is:

ρ tot := n tot Λ/M * + , n tot = min{⌊t/T ⌋, (E, M, F )(t) ∈ Σ -}.
These extremal values are obtained for a period T and with an entrance time t * that are shown in parentheses. We also indicate in Table 10 the maximal and minimal entrance times, obtained for a period T and an effort ratio ρ tot that are shown in parentheses. Note that consistently, the minimal entrance time is always obtained for φ = 8 and corresponds to the maximal effort ratio. Maximal entrance time is obtained for T = 1 (minimal tested period) and the minimal entrance time is obtained for T = 10 (maximal tested period). However, the minimal effort ratio is sometimes obtained with T = 2.

Comparing Tables 8 and9 shows that in general, a periodic control achieves the target of bringing the population into Σ -at a smaller cost than the constant control (in terms of total number of released mosquitoes, counted with respect to the wild population). Table 9: Minimal (left) and maximal (right) total effort ratio to get into the basin of 0 (in days) for various values of (ν E , β), the minimum and maximum being taken with respect to (T, φ), with a period and an entrance time shown in parentheses. The total effort ratio is defined as the total number of released male mosquitoes divided by the initial (wild) male mosquito population. Table 10: Minimal (left) and maximal (right) entrance time into the basin of 0 (in days) for various values of (ν E , β), the minimum and maximum being taken with respect to (T, φ), with a period and a total effort ratio shown in parentheses. 11: Entrance time into the basin of 0 (in days) for various values of (ν E , β) with constant weekly (T = 7 days) releases at p = 4 (left), p = 6 (center) or p = 8 (right).

Case study: Onetahi motu

We now parametrize explicitly our model to the case of Onetahi motu in Tetiaroa atoll (French Polynesia), where weekly (T = 7 days) releases have been performed over a year. Male population was estimated at 69 • 74 ≃ 5000 individuals, and the initial effort ratio p := Λ/M * + was estimated at 8.

For p ∈ {4, 6, 8}, entrance times (in days) are shown in Table 11 and final total female ratio in Table 12. This last quantity is important for practical purposes to help answering the question: when is it time to stop the releases? The trap counts during the experiment are to be compared with the initial trap counts (before the releases), and roughly, the process can be stopped once the ratio between the counts goes below the values in Table 12. Interestingly, β determines the order of magnitude of this final ratio.

Table 11 provides us interesting information on the entrance time versus the transition rate ν E and the mating parameter β. If the effort ratio p is not large enough, the SIT treatment can fail, and even if it is large enough (say p = 8) the time to reach the basin of 0 can be very large.

In the 3-dimensional state space (E, M, F ) we draw the full trajectory for the same sample value ν E \β 10 -4 10 -3 10 -2 10 -1 10 0 0.001 0.943252 0.147678 0.020134 0.002495 0.000283 0.002 0.567382 0.071552 0.009875 0.001247 0.000141 0.005 0.205116 0.031070 0.004439 0.000568 0.000069 0.008 0.133889 0.021388 0.003170 0.000425 0.000052 0.010 0.111803 0.018284 0.002779 0.000380 0.000047 0.015 N/A N/A N/A N/A N/A Table 12: Final total female ratio (F +Fst)(t)

F * +F * st at time t when the trajectory enter into the basin of 0 for various values of (ν E , β) with constant weekly (T = 7 days) releases at p = 8. (ν E = 0.008, β = 10 -3 , p = 8) along with a zoom in the last 30 days of treatment showing also the separatrix between the basins of E + and 0 as dots in Figure 2. According to Table 11, page 19, the entrance time is 541, which justifies that the control should last for more than one year. Our system being monotone, the trajectory is monotone decreasing (see Figure 2 (left), page 20). However, the rate of the decrease is relatively large at the beginning of the treatment, and then becomes small and, almost, constant. We also show time dynamics of four relevant normalized quantities, for the same sample value (ν E = 0.008, β = 10 -3 , p = 8) in Figure 3.

Conclusion

In this paper we have derived a minimalistic model to control mosquito population by Sterile Insect Technique, using either irradiation or the cytoplasmic incompatibility of Wolbachia to release sterilizing males. We particularly focus on the chance of collapsing the wild population, provided that the selected area allows elimination. Thus contrary to previous SIT and IIT models, the trivial equilibrium, 0 is always Locally Asymptotically Stable, at least. We consider different type of releases (constant, continuous, or periodic and instantaneous) and show necessary conditions to reach elimination, in each case. We also derived the minimal time under which elimination cannot occur, (i.e. entrance into the basin of attraction of 0 is impossible), whatever the control effort. Obviously, the knowledge on the mosquito parameters are very important, particularly the duration of the egg compartment,

1 µ E + ν E
and the mating parameter, β. Surprisingly, mosquito entomologists have not yet really focused their experiments on β or the probability of meeting/mating between one male and one female according to the size of the domains. Our model illustrates the importance of this parameter (and others) in the duration of the SIT control. In general, SIT entomologists recommend to release a minimum of ten times more sterile males than (estimated) wild males: this can be necessary if the competitiveness of the sterile male is weak compared to the wild ones (this can be the case with irradiation SIT approach). Our approach may help standardizing and quantifying this estimated ratio. Finally, we focus on a real case scenario, the Onetahi motu, where a experiment has been conducted by Bossin and collaborators, driving the local mosquito, Aedes polynesiensis, to nearly elimination. Our preliminary results show some good agreement with field observation (mainly trapping).

Our results also show the importance of eggs in the survival of the wild population. If the egg stock is sufficiently large, and depending on weather parameters, the wild population can re-emerge after the control has stopped. That is why, according to our model and numerical results, it is recommended to pursue the release of sterilizing males even after wild mosquito females are no longer collected in monitoring traps.

Last but not least, we hope that our theoretical results will be helpful to improve future SIT experiments and particularly to take into account the long term dynamics of eggs.

A Proof of Lemma 3

A.1 Study of f

We first study function f defined in [START_REF] Bliman | Ensuring successful introduction of wolbachia in natural populations of aedes aegypti by means of feedback control[END_REF]. For any y ≥ 0, if x ≥ 1 ψ then f (x, y) < -1 N (x + y) so in particular f (x, y) < 0. Therefore all steady states must satisfy βM * < 1 ψ . Likewise,

y ≥ 0, 0 ≤ x < 1 ψ =⇒ (1 -ψx)(1 -e -(x+y) ) < 1.
Hence for all x < 1 ψ we find f (x, y) < (1 -1 N )x -1 N y. As a consequence, if N ≤ 1 then f (x, y) < 0 for all (x, y) ∈ R 2 + \{0}, and system (2) has no positive steady state. From now on we assume that N > 1.

We also compute directly f (0, y) = -1 N y < 0 and lim x→+∞ f (x, y) = -∞.

Remark A.1 For all x ∈ (0, 1/ψ), we notice that

f (x, y) < Q y (x) = -ψx 2 + (1 - 1 N )x - y N .
The discriminant of the second-order polynomial Q y is

∆ y = (1 - 1 N ) 2 - 4yψ N . Let y := N 4ψ (1 -1 N ) 2 . If y ≥ y then ∆ y ≤ 0, hence f < 0.
At this stage we know that if βγ i M i ≥ y then there is no positive steady state.

The quantity y is used in Remark 2.2 to obtain a first-order approximation of the target release ration.

We now compute the derivatives of f :

∂ x f = (1 -2ψx)(1 -e -(x+y) ) - 1 N + x(1 -ψx)e -(x+y) , ∂ 2 xx f = -2ψ(1 -e -(x+y) ) + e -(x+y) 2 -(4ψ + 1)x + ψx 2 ∂ 3 xxx f = e -(x+y) -6ψ -3 + (6ψ + 1)x -ψx 2 =: e -(x+y) Q 3 (x) ∂ y f = x(1 -ψx)e -(x+y) - 1 N , ∂ 2 yy f = -x(1 -ψx)e -(x+y) < 0 for x ∈ (0, 1/ψ). Obviously, ∂ x f (x, y) < 0 if x ≥ 1 ψ and ∂ x f (0, y) = 1 -e -y -1 N , which is positive if and only if y > -log(1 -1 N ) = log(1 + 1 N -1
). In order to know the variations of ∂ 3 xxx f we study the second-order polynomial

Q 3 (x) = -6ψ -3β + x 6ψ + 1 -ψx 2 .
Its discriminant is

∆ 3 = (6ψ + 1) 2 -4ψ(6ψ + 3) = 1 + 12ψ 2 ,
which is positive. Therefore ∂ 3 xxx f is negative-positive-negative. More precisely, Q 3 is positive on

(w -, w + ) := 6ψ + 1 -1 + 12ψ 2 2ψ , 6ψ + 1 + 1 + 12ψ 2 2ψ .
To go one step further, we need to know the signs of ∂ This is always true:

Lemma 6 For all ψ > 0,

log 1 + 1 ψ -(2 + 1 2ψ )w + + 1 2 w 2 + -w + < 0.
Proof. To prove it, we introduce γ = 1 2ψ so that we are left with

log 7 + 3γ + γ 2 + (4 + γ) 3 + γ 2 -(3 + γ + 3 + γ 2 ) < 0.
To check this we introduce g(x) := log(7

+ 3x + x 2 + (4 + x) 3 + x 2 ) -(3 + x + 3 + x 2 ),
and we want to prove that g is negative. We compute that the sign of g ′ (x) is equal to that of

-(4 + x)(3 + x 2 ) -2x -3 + x 2 (8 + 2x + x 2 ) < 0. It remains to check that g(0) = log(7 + 4 √ 3) -(3 + √ 3) < 0, which is true since e 3+ √ 3 > e 4 > 2 4 > 7 + 8 > 7 + 4 √ 3,
where we used e > 2 and 1 < √ 3 < 2. Thus we obtain that x → ∂ 2 xx f (x, y) is either positive-negative (if y < log(1+ 1 ψ )) or always negative (otherwise).

The conclusion of all these computations is that in both cases (f is either convex-concave or simply concave), for any y, f (0, y) < 0, f (+∞, y) = -∞ so that all in all there are either 0, 1 or 2 solutions to f (x, y) = 0, depending merely on the sign of the maximum of x → f (x, y).

A.2 Study of functions h ±

We move on to the next step of the proof, studying the functions h ± defined in [START_REF] Bourtzis | Wolbachia-based technologies for insect pest population control[END_REF]. Recall that solving f (x, y) = 0 (for x, y > 0) is equivalent to picking θ = e -(x+y) ∈ (0, 1) and y = h ± (θ).

First, to check that h + and h -are well-defined we need to check that 1 + ξ log(θ) 1-θ > 0 for some θ ∈ (0, 1). It is easily checked that this is the case on (θ 0 (ξ), 1), and θ 0 (ξ) is well-defined as soon as ξ < 1.

Hence if ξ ≥ 1 then there is no nonzero steady state. Assume therefore that ξ < 1. Then there exists a unique θ 0 (ξ) ∈ (0, 1) such that 1θξ log(θ) has the same sign as θθ 0 on (0, 1), that is,

1 -θ 0 = 4ψ N log(θ 0 ). We can check that h -is decreasing, h -< h + on (θ 0 , 1], h ± (θ 0 ) = - 1 2ψ -log(θ 0 ), and
h -(1) < h + (1) = 1 2ψ -1 + 1 -ξ < 0. Indeed (recall that N ξ = 4ψ), h ′ -(θ) = - 1 θ - 1 N 1 θ(1-θ) + log(θ) (1-θ) 2 1 + ξ log(θ) 1-θ < 0, since - log(θ) 1 -θ < 1 θ .
Let y crit := max θ∈[θ 0 (ξ),1] h + (θ). If y = y crit then there is exactly one solution to f (x, y) = 0. For any y ∈ [0, y crit ), there are at least two solutions. By the previous computations we know that there are at most two solutions. So in this case there are exactly two solutions. To describe them one should consider

I 1 := [0, h -(θ 0 (ξ))], if h -(θ 0 (ξ)) > 0 (I 1 = ∅ otherwise)
, and I 2 = (max(θ 0 (ξ), 0), y crit ). If y ∈ I 1 then there is a solution of the form h -(θ -) and one of the form h + (θ -). If y ∈ I 2 then both solutions are of the form h + (θ), for two values of θ whose range contains the argument of y crit . And for y > y crit there is no solution.

At this stage we proved that if ξ ≥ 1 then there is no positive steady state; if ξ < 1 then if y crit > 0 then there are two positive steady states for βγ i M i ∈ [0, y crit ), 1 for βγ i M i = y crit and 0 for βγ i M i > y crit . If y crit = 0 then there is a unique positive steady state and if y crit < 0 then there is no positive steady state for any M i ≥ 0.

A.3 Stability

Finally, in order to compute the linearized stability of the steady states, we decompose J = M 0 +N 0 , where M 0 is non-negative and N 0 is diagonal non-positive. Then J Metzler, since E < K at steady states) is stable if and only if ρ(-N -1 0 M 0 ) < 1. We compute

N 0 =   -bF K -(ν E + µ E ) 0 0 0 -µ M 0 0 0 -µ F   and M 0 =   0 0 b(1 -E K ) (1 -r)ν E 0 0 rν E M M +γ i M i (1 -e -β(M +γ i M i ) ) rν E E M +γ i M i βM e -β(M +γ i M i ) + γ i M i M +γ i M i (1 -e -β(M +γ i M i ) ) 0  
so that for some X 1 , X 2 ∈ R (which we compute below at steady states) we have

-N -1 0 M 0 =     0 0 b(1-E K ) b F K +ν E +µ E (1-r)ν E µ M 0 0 X 1 µ F X 2 µ F 0     .
At the steady state (0, 0, 0), we have directly unconditional stability as

J =   -(ν E + µ E ) 0 b (1 -r)ν E -µ M 0 0 0 -µ F   ,
In particular, φ ′ (M + ) > 0 so

M + > 1 λ + 1 β (1 -e βM + ) = 1 λ - 1 β 1 (1 -λM + )N -1 .
Multiplying this inequality by λβ (1 -λM + )N -1 yields exactly the stability of M + , since (1 -λM + )N > 1. Indeed,

N (1 -λM ± ) = e βM ± e βM ± -1 > 1.
By a similar computation one can show that the smaller steady state M -is unstable.

We move now to the general case M i ≥ 0, assume M i < M crit i and write that ∂ x f < 0 (which was proved to hold at the bigger steady state) is equivalent to

(1 -2λM )(1 -e -β(M +γ i M i ) ) + βM (1 -λM )e -β(M +γ i M i ) < 1 N .
Using as before the fact that M is a steady state allows us to rewrite this last inequality as

(1 -2λM ) 1 N M + γ i M i M (1 -λM ) + βM (1 -λM ) 1 - M + γ i M i N (1 -λM )M < 1 N .

Multiplying this inequality by

N (1 -λM ) M M +γ i M i yields (1 -2λM ) + β(1 -λM ) N M 2 1 -λM M + γ i M i -M < (1 -λM ) M M + γ i M i , that is (1 -λM ) 2 - M M + γ i M i + βM -1 + N M 1 -λM M + γ i M i < 1,
whence the stability of the bigger steady state, since we recover [START_REF] Jachowski | Filariasis in american samoa. y. bionomics of the principal vector, aedes polynesiensis marks[END_REF]. Likewise, at the smaller steady state we have ∂ x f > 0, and the reverse inequality holds. This concludes the proof.

B Basin entrance time approximation

B.1 Bounds on the wild equilibria

For M i = 0, under the assumptions of Lemma 3 such that there are two positive steady states E -≪ E + for (2), we get explicit bounds on these states. In particular, we assume N > 4ψ.

We recall that the positive equilibria can be expressed as an increasing function of their second coordinate M ∈ (0, 1/λ):

E(M ) :=   KλM M ν E +µ E b λM (1-λM )   ,
and E(M ) is an equilibrium if and only if f (βM ) = 0, where

f (x) = (1 -ψx)(1 -e -x ) - 1 N . ( 21
)
Lemma 7 The function f (defined in [START_REF] Hapairai | Studies on Aedes polynesiensis introgression and ecology to facilitate lymphatic filariasis control[END_REF]) is concave on [0, 1/ψ]. It reaches its maximum value on this interval at Z(ψ) ∈ (0, 1 2ψ ), where we define

e -Z(ψ) = ψ 1 + ψ -ψZ(ψ) , F (ψ) := 1 + ψ -ψZ(ψ) (1 -ψZ(ψ)) 2 . ( 22
)
Then f on [0, 1/ψ] has no zero if N < F (ψ), exactly 1 zero if N = F (ψ) and exactly 2 zeros if N > F (ψ).

In addition, Z and F have the following asymptotics:

Z(ψ) ∼ ψ→+∞ 1 2ψ , Z(ψ) ∼ ψ→0 log 1 ψ , F (ψ) ∼ ψ→+∞ 4ψ, F ---→ ψ→0 1.
Proof. We compute

f ′ (x) = e -x 1 + ψ -ψx -ψ, f ′′ (x) = e -x ψx -1 -2ψ , hence f ′′ < 0 on [0, 1/ψ]. Since f (0) = f (1/ψ) = -1/N < 0, f reaches a unique maximum at the (necessarily unique) point Z(ψ) ∈ (0, 1/ψ) such that f ′ (Z(ψ)) = 0. The claim that Z(ψ) < 1/(2ψ)
follows from the inequality e x > 1 + x, which implies that

1 ψ f ′ ( 1 2ψ ) = e -1/2ψ 1 + 1 2ψ -1 < 0.
Moreover, the sign of f (Z(ψ)) is exactly that of N -F (ψ). The equivalents and limit follow from straightforward computations.

Remark B.1

We notice that Z is related to a well-known special function: let us introduce the (principal branch of the) special Lambert W function, that is: W (y) = z, z ≥ -1 ⇐⇒ ze z = y.

Since if y > 1 then z > 0, we obtain Z(ψ) = log W (e 1+1/ψ ) .

Assume N > F (ψ) (defined in ( 22)), and denote by x -< x + the two positive zeros of f .

Lemma 8

We have x -> 1/N .

1 N < x -< 1 ψ 1 - κ * N < Z(ψ) < 1 ψ 1 - κ * N < x + ,
where

κ * = 1 + ψ 1 -ψZ(ψ)
, κ * = N -ψZ(ψ)(1 + ψ -ψZ(ψ))

(1 -ψZ(ψ)) 2 .

If in addition N > 2 then x + < 1 ψ 1 -1 N .

Proof. The first inequality is obtained by using the inequalities 1e -x ≤ x and 1 -√ 1x > x/2 for x ∈ (0, 1). The first one implies that f (x) ≤ x(1ψx) -1/N , which is a second order polynomial equal to f at 0 and at 1/ψ, with roots located at 1 ± 1 -4ψ/N /(2ψ) (recall that we have N > 4ψ). Hence x -> 1 -1 -4ψ/N /(2ψ) > 1/N by the second inequality.

The upper bound on x + comes from the fact that if N > 2 then by Lemma 7

1 - 1 N 1 ψ > 1 2ψ > Z(ψ).
Finally to get the two other bounds, we introduce

H(κ) := f 1 ψ (1 - κ N ) = κ 1 -e -1 ψ (1-κ N ) -1.
By Lemma 7, it is concave on [0, N ], equal to -1 at 0 and N and reaches its maximum at κ := N (1 -ψZ(ψ)). To get κ * and κ * , we simply use the fact that the graph of H is above the segments from (0, -1) to ( κ, H( κ)) on the first hand, and from ( κ, H( κ)) to (N , 0) on the other hand, so that we define

-1 + H( κ) + 1 κ κ * = 0 = -1 -(κ * -N ) H( κ) + 1 N -κ ,
and the expressions of κ * < κ < κ * follow from a straightforward computation.

Back to the steady states of (2), we deduce from Lemma 8 the following bounds, assuming N > 2:

E -:=    λK N β 1 N β ν E +µ E b λK N β    ≤E -≤ 1 - κ * N   K 1 λ ν E +µ E b KN κ *   =: E - (23)
E + := 1 - κ * N   K 1 λ ν E +µ E b KN κ *   ≤E + ≤ 1 - 1 N   K 1 λ KN (ν E +µ E ) b   =: E + . (24)

B.2 Results

A lower bound. First, we give a lower bound on the entrance times. We consider the fact that for a solution to [START_REF] Anguelov | On nonstandard finite difference schemes in biosciences[END_REF] with initial data given by E + , thanks to the overestimation in [START_REF] Moreira | A wolbachia symbiont in aedes aegypti limits infection with dengue, chikungunya, and plasmodium[END_REF],

F (t) ≥ F + e -µ F t =: F ♭ (t).
This implies

E(t) ≥ e -(ν E +µ E )t-b F + K (1-e -µ F t ) E + + b F + t 0
e -µ F t ′ e -(ν E +µ E )(t-t ′ ) e -b F + K (e -µ F t ′ -e -µ F t ) dt ′ =: E ♭ (t), and M (t) ≥ e -µ M t M + + (1r)ν E t 0 e -µ M ((t-t ′ ) E ♭ (t ′ )dt ′ =: M ♭ (t).

Using the underestimation of E -from (23), we define t Z ♭ := min{t ≥ 0, Z ♭ (t) ≤ Z -} for Z ∈ {E, M, F }.

Then, we construct the explicit solution (E, M ) = ( E ♯ , M ♯ ) to

Ė = b F ♯ -ν E + µ E + F ♯ K E, E(0) = E + , Ṁ = (1 -r)ν E E -µ M M, M (0) = M + .
In details:

F ♯ (t) = E + rν E ǫ + e -µ F t F + -rν E ǫ E + , E ♯ (t) = e -(ν E +µ E + E + rν E ǫ K )t- F + -rν E ǫ E + Kµ F (1-e -µ F t ) E + + t 0 b E + rν E ǫ +be -µ F t ′ ( F + -rν E ǫ E + ) e (ν E +µ E + E + rν E ǫ K )t ′ - F + +rν E ǫ E + Kµ F (1-e -µ F t ′ ) dt ′ , M ♯ (t) = e -µ F t M + + (1 -r)ν E t 0 e µ F t ′ E ♯ (t ′ )dt ′ .
We use this super-solution on [0, t 0 ] (for some t 0 > 0 to be determined), and then glue the solution on [t 0 , +∞) of

             Ė = bF -(ν E + µ E )E, E(t 0 ) = E ♯ (t 0 ), Ṁ = (1 -r)ν E E -µ M M, M (t 0 ) = M ♯ (t 0 ), Ḟ = rν E ǫ 0 E -µ F F, F (t 0 ) = F ♯ (t 0 ),
with ǫ 0 = M ♯ (t 0 )/( M ♯ (t 0 ) + M i ) < ǫ.

For Z ∈ {E, M, F } we let t Z ♯ (t 0 ) := min{t ≥ t 0 , Z ♯ ≤ Z -}.

Then as before:

Lemma 12 For all t 0 > 0, τ (M i ) ≤ t ♯ (t 0 ) := max(t E ♯ (t 0 ), t M ♯ (t 0 ), t F ♯ (t 0 )). By using Lemma 12, we can theoretically obtain a finite upper bound for τ (M i ) (upon choosing a suitable t 0 ) as soon as ǫ 0 < 1/N for t 0 large enough, that is if and only if

ρ * (ρ * + 1) µ M (1 -r)ν E + N -1) > N -1. (27)
Condition ( 27) is weaker than (26) (and in general, much weaker). It holds if and only if

ρ * > -(N -1 + φ) + (N -1 + φ) 2 + 4φ(N -1) 2φ , φ := λK = µ M (1 -r)ν E ,
which is true for instance if ρ * > (N -1)/φ. However, we do not develop any further these analytic computations in the present paper.

B.3 Analytic computations

Applying Lemma 11, in order to express analytically the solution X e := (E e , M e , F e ) of ( 25), we only need to diagonalize the matrix

R e := -(ν E + µ E ) b rν E ǫ -µ F .
R e has negative trace, and positive determinant if and only if 1 N > ǫ. Hence if N ǫ(M i ) < 1 then 0 is globally asymptotically stable for [START_REF] Nadin | Hindrances to bistable front propagation: application to wolbachia invasion[END_REF].

In this case its eigenvalues are real, negative and equal to κ ± associated respectively with eigenvectors 1

x ± , where

κ ± := -(ν E + µ E + µ F ) ± (ν E + µ E -µ F ) 2 + 4brν E ǫ 2 , x ± := ν E + µ E -µ F ± (ν E + µ E -µ F ) 2 + 4brν E ǫ 2b .
Then we deduce that for some real numbers (r 0 ± , s 0 ± ) ∈ R 4 , E e (t) = r 0 + e κ + t + r 0 -e κ -t , F e (t) = s 0 + e κ + t + s 0 -e κ -t , M e (t) = e -µ M t M 0 e + (1r)ν E t 0 e -µ M (t-t ′ ) r 0 + e κ + t ′ + r 0 -e κ -t ′ dt ′ .

In details, we find

r 0 + = x - x --x + E 0 e - 1 x --x + F 0 e , r 0 -= -x + x --x + E 0 e + 1 x --x + F 0 e s 0 + = x + x - x --x + E 0 e - x + x --x +
F 0 e , s 0 -= -x + x - x -x + E 0 e + x - x -x + F 0 e .

Assuming κ + = -µ M and κ -= -µ M (which must hold generically since these are biological parameters), we get M e (t) = e -µ M t M 0 e + (1r)ν E r 0 + e κ + te -µ M t µ M + κ + + r 0 -e κ -te -µ M t µ M + κ - .

Assuming N > 2, we use the overestimation [START_REF] Moreira | A wolbachia symbiont in aedes aegypti limits infection with dengue, chikungunya, and plasmodium[END_REF] of E + as an initial data (E 0 e , M 0 e , F 0 e ), and with the notations

g(ǫ) = 1 + 4brν E ǫ (ν E + µ E -µ F ) 2 , σ = sgn(ν E + µ E -µ F ), we deduce r 0 ± = K 2 1 - 1 N 1 ± (2N -1)(ν E + µ E ) + µ F g(ǫ)|ν E + µ E -µ F | , s 0 ± = K|ν E + µ E -µ F | 4bg(ε) 1 - 1 N σ ± g(ε) g(ε) ± (2N -1)(ν E + µ E ) + µ F |ν E + µ E -µ F | .

Figure 1 :

 1 Figure 1: Two viewpoints of scattered points lying around the separatrix (ε = 10 -2 ) for ν E = 0.1 and β = 10 -4 . In this case, 5 females or 900 eggs are enough to prevent population elimination.

ν

  E \β 10 -4 10 -3 10 -2 10 -1 10 0 N/A N/A N/A N/A N/A 0.010 N/A N/A N/A N/A N/A 0.015 N/A N/A N/A N/A N/A ν E \β 10 -4 10 -3 10 -2 10 -1 10 0 N/A N/A N/A N/A N/A 0.015 N/A N/A N/A N/A N/A ν E \β 10 -4 10 -3 10 -2 10 -1 10 0 N/A N/A N/A N/A N/A Table

Figure 2 :

 2 Figure 2: Trajectory t → (E(t), M (t), F (t)) for ν E = 0.008 and β = 10 -3 (left) and a zoom in the last 30 days of treatment displaying also the separatrix as dots (right).

  Relative total female count F (t)+Fst(t)

Figure 3 :

 3 Figure 3: Time dynamics of different ratio for ν E = 0.008 and β = 10 -3 .

Table 1 :

 1 Parameter values for some populations of Aedes polynesiensis in French Polynesia at a temperature of 27 • C.

	Symbol Name	Value interval	Source
	r viable N eggs τ gono τ E τ L r L r	Proportion of viable eggs Number of eggs laid per laying Duration of gonotrophic cycle Egg half-life Time from hatching to emergence Survivorship from larva first instar to pupa Sex ratio (male:female)	Field collection, [21, p. 121] [29] [20, 35, 29] 15 -30 days Estimation (to be determined) 95 -99% 55 -75 4 -7 days Lab data, [21, p. 104] 8 -11 days Lab data, [21, p. 106] 67 -69% 49% Production data (ILM)
	τ M τ F γ i	Adult male half-life Adult female half-life Mating competitiveness of sterilizing males	5 -9 days 15 -21 days 1	Lab data, [21, p. 50] Lab data, [21, p. 50] Lab [21, pp. 51-53], field [26]

  For later use, we introduce the usual relations ≪, < and ≤ on R d (where d ≥ 1) as the coordinatewise partial orders on R d induced by the cone R d + . More precisely, for x, y ∈ R d , • x ≤ y if and only if for all 1

	Symbol Name	Formula	Value interval
	b µ L ν L ν E ν E µ E µ M µ F	Effective fecundity Larva death rate Larva to adult transition rate Larval coefficient for effective hatching rate Egg death rate Adult male death rate Adult female death rate	r viable N eggs τ gono -log(r L ) τ L 1 τ L ν L ν L + µ L log(2) τ E log(2) τ M log(2) τ F	7.46 -14.85 0.034 -0.05 0.09 -0.125 0.64 -0.79 0.023 -0.046 0.077 -0.139 0.033 -0.046
	Table 2: Conversion of the biological parameter from Table 1 into mathematical parameters for
	systems (1) and (2)
	2 Theoretical study of the simplified model

Table 3

 3 Thus according to the values taken by ν E in Table 3, page 15, we have the following bounds for N : 29 ≤ N ≤ 105.

	, page 15, we have
	N = 117.5	ν E ν E + 0.03	.
	Remark 3.1

Table 4 :

 4 CPU times for the numerical simulations

Table 5 :

 5 Effort ratio ρ * = M crit i /M * (M i = 0) for various values of ν E . For this range of parameters, ρ * is practically independent on β ∈ [10 -4 , 1].

		005 0.010 0.020 0.030 0.050 0.100 0.150 0.200 0.250
	ρ *	16	30	48	60	76	93	101	106	108

  E \β 10 -4 10 -3 10 -2 10 -1 10 0

						10 -4 10 -3 10 -2 10 -1 10 0
	0.005	63	151	204	253 303	323	445	571	697	824
	0.010	93	180	232	281 331	361	475	592	708	825
	0.020 118	203	256	304 354	381	485	590	695	800
	0.030 130	215	267	315 365	391	488	587	685	783
	0.050 141	226	278	327 377	440	530	621	713	804
	0.100 152	236	289	337 387
	0.150 156	240	293	341 391
	0.200 158	242	295	343 393
	0.250 160	244	296	344 395

ν

  4 10 -3 10 -2 10 -1 10 0

	258	351	448	545	642
	286	374	464	553	643
	301	381	462	544	625
	307	383	461	538	615
	332	404	477	550	623

Table 6 :

 6 Left: under-estimation of the entrance time into the basin of 0 from the analytic formula[START_REF]Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control[END_REF]. Middle and right: over-estimation of the entrance time into the basin of 0 from formula (9) with ε = ∈ {1.2, 1.4, 1.6, 1.8, 2, 4, 8}. Results are shown in Table7for φ = 1.2, φ = 2 and φ = 8. E \β 10 -4 10 -3 10 -2 10 -1 10 0

	M * +
	M * + +φM crit
					10 -4 10 -3 10 -2 10 -1 10 0	10 -4 10 -3 10 -2 10 -1 10 0
	0.005 168	286	363	435 504	148	262	338	409 478	128	237	311	380 449
	0.010 200	305	376	441 505	180	283	352	417 480	160	258	326	391 454
	0.020 219	313	377	437 495	199	292	355	415 473	180	270	333	392 450
	0.030 225	314	375	434 492	207	295	355	413 471	188	274	334	392 450
	0.050 228	314	373	431 488	212	297	355	413 470	194	278	336	394 452
	0.100 231	314	372	430 488	215	298	356	414 472	200	282	340	398 456
	0.150 232	315	373	431 489	217	300	358	416 474	202	285	343	401 459
	0.200 233	316	375	433 491	219	302	360	418 476	205	287	345	403 462
	0.250 234	318	376	434 493	220	303	362	420 478	206	289	347	406 464

i

, when applicable, for φ = 8 (middle) and φ = 4 (right).

Then we compute numerically the entrance time for a range of releasing efforts. In details, computations were performed for M i = φM crit i with φ ν

  . Last but not least, ν E \β 10 -4 10 -3 10 -2 10 -1 10 0

				10 -4 10 -3 10 -2 10 -1	10 0	10 -4	10 -3	10 -2	10 -1	10 0
	0.005 399	680	863 1034 1199	587 1036 1338 1619	1893	2024	3749	4920	6027	7115
	0.010 854 1302 1603 1880 2154	1283 2009 2499 2962	3416	4548	7343	9257 11112 12912
	0.020 1513 2166 2603 3022 3423	2296 3367 4092 4782	5452	8317 12451 15331 18040 20736
	0.030 1950 2726 3253 3761 4264	2989 4260 5132 5976	6815	10862 15871 19319 22691 26040
	0.050 2482 3421 4059 4686 5315	3837 5381 6434 7483	8529	14058 20188 24395 28588 32774
	0.100 3100 4218 5000 5777 6553	4817 6675 7975 9268 10563	17891 25266 30457 35640 40813
	0.150 3383 4581 5434 6278 7122	5274 7268 8688 10095 11502	19651 27618 33285 38913 44541
	0.200 3545 4806 5694 6578 7461	5549 7638 9117 10588 12060	20757 29073 34979 40865 46750
	0.250 3649 4956 5869 6779 7689	5717 7884 9405 10922 12438	21443 30036 36123 42188 48254

Table 8 :

 8 Total effort ratio to get into the basin of 0 for various values of (ν E , β), with M

  357) 1112 (1, 519) 1253 (1, 585) 1386 (1, 647) 1504 (2, 771) 0.030 1125 (3, 363) 1371 (1, 510) 1538 (1, 572) 1696 (2, 693) 1839 (2, 752) 0.050 1383 (2, 379) 1669 (1, 496) 1875 (1, 556) 2066 (2, 672) 2238 (2, 730) 0.100 1655 (2, 370) 1997 (1, 480) 2238 (1, 539) 2458 (2, 650) 2678 (2, 708) 0.150 1772 (1, 388) 2134 (1, 473) 2394 (2, 583) 2632 (2, 641) 2871 (2, 699) 0.200 1834 (1, 384) 2213 (1, 470) 2482 (2, 578) 2731 (2, 636) 2979 (1, 738) 0.250 1873 (1, 382) 2263 (1, 468) 2531 (2, 575) 2787 (2, 633) 3043 (2, 692) 286) 11921 (10, 344) 14025 (10, 402) 15778 (10, 460) 8794 (10, 206) 12114 (10, 289) 14709 (10, 347) 17305 (10, 405) 19900 (10, 463) 9522 (10, 209) 13603 (10, 291) 15948 (10, 350) 18762 (10, 408) 21576 (10, 466) 10431 (10, 211) 14201 (10, 293) 17138 (10, 352) 19586 (10, 410) 22524 (10, 468) 10709 (10, 212) 14584 (10, 295) 17601 (10, 353) 20618 (10, 412) 23133 (10, 470)

	νE\β 0.005 282 (2, 287) 10 -4 0.010 547 (1, 344) 0.020 900 (1, 10 -4 10 -3 10 -2 10 -1 10 0 384 (2, 491) 448 (1, 608) 502 (1, 682) 554 (1, 752) 698 (2, 497) 796 (1, 602) 884 (1, 669) 969 (2, 805) 1095 (10, 135) 2317 (10, 168) 4139 (10, 188)	10 -3 1838 (10, 248) 3575 (10, 268) 6015 (10, 280)	10 -2 2450 (10, 323) 4536 (10, 337) 7573 (10, 343)	10 -1 2986 (10, 393) 5499 (10, 402) 8909 (10, 402) 10246 (10, 460) 10 0 3522 (10, 462) 6323 (10, 466)
	5448 (10, 196)	7829 (10, 283)	9506 (10, 343) 11183 (10, 402) 12581 (10, 460)
	7155 (10, 201)	9818 (10,

  8794) 289 (10, 12114) 347 (10, 14709) 405 (10, 17305) 463 (10, 19900) 0.150 209 (10, 9522) 291 (10, 13603) 350 (10, 15948) 408 (10, 18762) 466 (10, 21576) 0.200 211 (10, 10431) 293 (10, 14201) 352 (10, 17138) 410 (10, 19586) 468 (10, 22524) 0.250 212 (10, 10709) 295 (10, 14584) 353 (10, 17601) 412 (10, 20618) 470 (10, 23133)

	νE\β 0.005 135 (10, 1095) 10 -4 0.010 168 (10, 2317) 0.020 188 (10, 4139) 0.030 196 (10, 5448) 0.050 201 (10, 7155) 0.100 206 (10, 10 -4 10 -3 10 -2 10 -1 10 0 248 (10, 1838) 323 (10, 2450) 393 (10, 2986) 462 (10, 3522) 268 (10, 3575) 337 (10, 4536) 402 (10, 5499) 466 (10, 6323) 280 (10, 6015) 343 (10, 7573) 402 (10, 8909) 460 (10, 10246) 283 (10, 7829) 343 (10, 9506) 402 (10, 11183) 460 (10, 12581) 286 (10, 9818) 344 (10, 11921) 402 (10, 14025) 460 (10, 15778) 456 (1, 317) 528 (1, 629) 534 (1, 1012) 642 (1, 1179) 708 (1, 1300) 771 (1, 1414) 830 (1, 1522) 10 -3 10 -2 10 -1 10 0 667 (1, 420) 752 (1, 474) 826 (1, 521) 896 (1, 565) 661 (1, 749) 735 (1, 833) 803 (1, 909) 868 (1, 982) 527 (1, 1246) 627 (1, 1445) 690 (1, 1588) 749 (1, 1724) 807 (1, 1860) 514 (1, 1513) 605 (1, 1749) 666 (1, 1925) 724 (1, 2090) 782 (1, 2257) 494 (1, 1787) 581 (1, 2072) 640 (1, 2279) 698 (1, 2485) 755 (1, 2692)
	483 (1, 1896) 569 (1, 2200) 628 (1, 2428) 686 (1, 2652) 744 (1, 2877)
	477 (1, 1953) 563 (1, 2272) 622 (1, 2510) 680 (1, 2745) 738 (1, 2979)
	473 (1, 1988) 559 (1, 2317) 618 (1, 2562) 676 (1, 2802) 734 (1, 3043)

  2 xx f (w + , y) and ∂ 2 xx f (0, y). We write ∂ 2 xx f (x, y) > 0 ⇐⇒ e -(x+y) 2 + 2ψ -(4ψ + 1)x + ψx 2 > 2ψ Hence ∂ 2 xx f (0, y) > 0 if and only if y < log(1 + 1 ψ ). Similarly, ∂ 2 xx f (w + , y) < 0 if and only if

	y > log 1 +	1 ψ	-(2 +	1 2ψ	)w + +	1 2	w 2 + -w + .

Acknowledgments: This study was part of the "Phase 2 SIT feasibility project against Aedes albopictus in Reunion Island", jointly funded by the Regional Council of La Reunion and the European Regional Development Fund (ERDF) under the 2014-2020 Operational Program. YD was (partially) supported by the DST/NRF SARChI Chair M3B2 grant 82770. YD and MS also acknowledge partial supports from the STIC AmSud "MOSTICAW Project", Process No. 99999.007551/2015-00.

whose eigenvalues are -(ν E + µ E ), -µ M and -µ F .

At a non-zero steady state we recall that bF

so that

.

The characteristic polynomial of -N -1 0 M 0 is

which is equal to

and we rewrite it as

We find P (0) > 0 (since X 2 > 0) and

so that J is stable if and only if P (1) < 0. (P is increasing and then decreasing on (0, +∞)). This condition reads

Let us treat first the case when M i = 0. The stability condition rewrites

that is, for a nonzero steady state,

If M crit i > 0, we know that there are exactly two steady states between 0 and 1/λ for M i = 0, which we denote by

Lemma 9 We have the following lower bound:

However it must be expected that min(t E ♭ , t M ♭ ) > t F ♭ , and we can give explicit approximations of t E ♭ and t M ♭ .

A first upper bound. We compare the solution of (2) with the solution of the linear system

where

The following property follows from the fact that ( 2) is cooperative: Lemma 10 Solutions of (2) and [START_REF] Nadin | Hindrances to bistable front propagation: application to wolbachia invasion[END_REF] with initial data such that (E 0 , M 0 , F 0 ) ≤ (E 0 e , M 0 e , F 0 e ) satisfy:

∀t ≥ 0, (E(t), M (t), F (t)) ≤ (E e (t), M e (t), F e (t)).

We use the under-estimation of E -given by [START_REF] Li | Modelling releases of sterile mosquitoes with different strategies[END_REF], to define, for X = (X i ) i = (E, M, F ) and i ∈ {1, 2, 3},

Lemma 11 For any solution X e to (25) satisfying the assumption of Lemma 10, we have the upper bound on the entrance time: τ (M i ) ≤ max t E min , t M min , t F min . Analytic computations are made in Section B.3.

An second upper bound in two steps. Let ρ * := M i / M + be the under-estimated effort ratio. When using the above one-step approach, we conclude with a finite upper bound for τ (M i ) if and

Expanding upon the same idea as for the lower bound, we let ǫ = M + /( M + + M i ) so

If r 0 -< 0 the we can use the simple upper bound E e (t) ≤ r 0 + e κ + t . This condition reads

In this case, we know that E e (t) ≤ E -

Then, under the same condition we have s 0 ± > 0. By using the fact that

In addition, we have t E min > t F min if and only if

Remark B.2 For small ǫ, the previous estimations roughly show that

Finally, we need to compute the condition M e (t)

where

.

The condition we need to compute is therefore

We assume that the male half-life is shorter than that of the females and of the eggs, so that σ F , σ E > 1. Under the stronger assumptions that r -< 0 < r + and

we obtain that α > 0. We simply treat two subcases: first if µ M + κ + < 0 (small µ M ) then we obtain α + < 0 < α -and thus

Second, if µ M + κ -> 0 (large µ M ) then we obtain α -< 0 < α + and thus t M min := 1 -κ + log (N -1) α + α + ψ .

In the last case (when µ M is large), we can check that t M min > t E min is equivalent to α + α + > r 0 + , which holds since α > 0 and α + > r 0 + . In this case we obtain

.


Alexandr Klimchik

Damien Chablat

Anatol Pashkevich

Advancement of MSA-technique for stiffness modeling of serial and parallel robotic manipulators

Keywords: Stiffness modeling, matrix structural analysis, serial robots, parallel robots

Introduction

In many modern robotic applications, manipulators are subject to essential external loadings that affect the positioning accuracy and provoke non-negligible compliance errors [START_REF] Klimchik | Efficiency evaluation of robots in machining applications using industrial performance measure[END_REF]. For this reason, manipulator stiffness analysis becomes one of the most important issues in the design of robot mechanics and control algorithms. It allows the designer to achieve required balance between the dynamics and accuracy. However, to make the stiffness analysis efficient, it should rely on a simple and computationally reasonable method that is able to deal with flexible links and rigid connections, passive and elastic joints that are common for serial and parallel robots.

At present, there exist three main techniques in this area [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF]: Finite Element Analysis (FEA), Matrix Structural Analysis (MSA) and Virtual Joint Modeling (VJM). The most accurate of them is the FEA [START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF], which allows modeling links and joints with their true dimension and shape. However, it is usually applied at the final design stage because of the high computational expenses. The MSA is considered as a compromise technique, which incorporates the main ideas of the FEA, but operates with larger elements. This obviously leads to the reduction of the computational expenses, but it requires some non-trivial actions for including of passive and elastic joints in the related mathematical model. Some reviews of existing works on manipulator stiffness analysis can be found in [START_REF] Gonçalves | A comparison of stiffness analysis methods for robotic systems[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF][START_REF] Liu | Stiffness Modeling of Parallel Mechanisms at Limb and Joint/Link Levels[END_REF][START_REF] Yeo | Design and analysis of cable-driven manipulators with variable stiffness[END_REF][START_REF] Klimchik | Serial vs. quasi-serial manipulators: Comparison analysis of elasto-static behaviors[END_REF], where the authors addressed different aspects of the FEA, MSA and VJM techniques. Among recent contributions devoted to the MSA, it is worth mentioning the work of Cammarata [START_REF] Cammarata | Unified formulation for the stiffness analysis of spatial mechanisms[END_REF], who merged together flexibilities of the link and adjacent joints but did not avoid manual work of lines/columns merging in the global stiffness matrix. There are also several works that deal with the MSA application to particular parallel and serial manipulators [START_REF] Shi | Static stiffness modelling of EAST articulated maintenance arm using matrix structural analysis method[END_REF][START_REF] Azulay | Comparative analysis of a new 3×PPRS parallel kinematic mechanism[END_REF].

The most essential contribution to the robot-oriented modification of the MSA was done by Deblase et. all [START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF]. The authors proposed a general technique to take into account passive joints and rigid connections via the additional constraints. Nevertheless, some manual procedures of merging matrix components were not avoided, as well as preloadings and elastic connections were not treated. For this reason, this paper focuses on the advancement of the MSA technique for robotic applications allowing to reduce related manual operations.

MSA models of manipulator links

Modeling of the flexible link. If the link flexibility is non-negligible, the 2-node "free-free" stiffness model should be used, which can be obtained either from the approximation of the link by a beam or using the CAD-based technique [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF]. In both cases, the link is described by the linear matrix equation

T ij ij ij i j                    I D I 0 D 0I t d t (2)
where () [] ij  d denotes the 3 3  skew-symmetric matrix derived from the vector

( ) ( ) ( ) ( ) [ , , ] ij ij ij ij x y z T ddd  d
describing the link geometry from the i th to j th node.

The second group of equations describing the force equilibrium can be derived using rigid body static equations. Using the above-introduced definition for () ij D , this equation can be written in a form

() 6 6 6 1 ij T i j        W 0 W ID (3)
that is also convenient for aggregation of the model components.

MSA models of manipulator's joints

Modeling of the rigid joint. If two adjacent links are connected by means of a rigid joint, the stiffness model must include relations describing two principal rules of structural mechanics: (a) displacement compatibility, (b) force equilibrium. If the nodes corresponding to the adjacent link are denoted as , ij , the first of these rules can be expressed as

ij    tt yielding the following constraint   6 6 6 6 6 1 6 12 i j             t 0 t II (4)
The second rule can be written as ij  W W 0 leading to the second constraint  

6 6 6 6 6 1 6 12 i j         W I I 0 W (5)
Similar equations can be written for the rigid connection of three and more links.

Modeling of the passive joint.

If two adjacent links (ai), (jb) are connected by means of a passive joint, the displacement compatibility condition must be replaced by

  r ij i j      Λ t t 0, where r ij
Λ is a rank-deficient matrix that defines directions for the passive joint that do not admit free relative motions. This allows presenting the first constraint in the form

1 12 ** rr i r r j ij ij             Λ 0 t Λ t (6)
that contains a rectangular matrix * r ij Λ of the rank r . In the simple cases, the matrix * r ij Λ can be easily derived from r ij Λ by simple elimination of zero rows. The force equilibrium condition is presented in this case as

  r ij i j    Λ W W 0 that gives the second constraint in the following form 1 12 ** rr ij ij i r r j         W 0 W ΛΛ (7)
In addition, passive joints do not transmit the force/torque in the direction corresponding to zero columns in the matrix r ij Λ , which yields complementary equations

1 * 1 * ; pp ij i i j p j p      Λ W 0 Λ W 0 (8)
where * p ij Λ denotes the full-rank matrix of size 6 p  defining passive directions, and 6 pr . Modeling of the elastic joint. If two adjacent links are connected by an elastic joint, that can be treated as a passive joint with springs. In this case, the deflection compatibility condition remains the same. However, the force equilibrium must be slightly revised and presented as ** ;

( )

e e e i ij ij ij ji ij          W W 0 W Λ K Λ t t (9)
where

                       t 0 0 I I t K ΛK 0 Λ W Λ 0 W (10)
It should be stressed that the above expressions are valid for the so-called "no preloading case", when equal deflections do not generate the elastic forces, i.e. [START_REF] Azulay | Comparative analysis of a new 3×PPRS parallel kinematic mechanism[END_REF] and the right-hand side of (10) should be slightly modified.

ji e e e i ij ij ij ij i j        W W 0 Λ K Λ t t WW
Modeling of the actuated joint. The actuated joint ensures transmission of the force/torque between the manipulator links. In the frame of MSA technique, it can be presented either as a rigid, passive or elastic connection. In the first two cases, the actuating effort is transmitting in the direction corresponding to the vectors 1 ,..., r uu while in the last case transmission can be performed in the direction defined by the vectors 1 6 ,..., r  uu . Hence, the actuated joint may be included in the global stiffness model as a set of two linear matrix constraints.

Including of boundary conditions and loadings

Modeling of rigid support. The rigid connection of the link (jb) to the robot base can be presented as a special case of the rigid joint with an eliminated link (ai) and zero deflection j  t0 . This simplifies the deflection compatibility constraint as

  6 6 6 1 j       I0 t (12)
that contains 6 linear equations to be included in the global stiffness model. Corresponding reaction wrench at the rigid support may be computed as

() 12
jb jb   W K t ( 13
)
where the deflection b t is obtained from the solution of the global stiffness equations.

Modeling of passive support. The passive connection of the link (jb) to the robot base can be presented as a special case of the passive joint with an eliminated link (ai) and zero deflection in the non-passive directions * r ij j    Λ t 0 . This simplifies the deflection compatibility constraint down to

1 * 6 jr r ij r          t 0 Λ (14)
that contributes r linear equations to the global stiffness model. In addition, it is necessary to take into account that the wrench components for the passive direction are equal to zero

* p ij j  Λ W 0 (15)
that contributes another p equations to the global stiffness model. Totally, this gives 6 rp  independent linear constraints.

Modeling of elastic support. Similarly, the elastic connection of the link (jb) to the robot base can be presented as a special case of the elastic joint with an eliminated link (ai) and zero deflection in the non-elastic directions * r ij j    Λ t 0 . This allows us to use the same deflection compatibility constraint as above, which contributes r linear equations to the global stiffness model. In addition, it is necessary to consider that the wrench components corresponding to the non-rigid directions are produced by the elastic forces satisfying the Hook's law

              t K Λ 0 W Λ (16)
that contributes another e scalar equations to the global stiffness model. Totally, this gives 6 re  independent linear constraints.

Including the external loading. In robotic manipulators, both serial and parallel, there is at least one node that is not connected directly to the robot base. It corresponds to the end-effector that interacts with robot environment by applying the force/torque to the external objects. For the global stiffness model, the end-effector produces the boundary conditions that should include the vector of the external wrench e W , which is also necessary for computation of the Cartesian stiffness matrix.

To take into account the external loading, the global stiffness model must be completed by the linear constraint derived from the force equilibrium at the node e, i.e.

ije  W W W , which can be rewritten in the form

    6 6 6 6 6 12 6 1 12 1 i j e        W I I W W (17)
Similarly, it is possible to derive the constraints produced by the external wrenches applied to other nodes.

Aggregation of MSA model components

In contrast to large mechanical structures consisting of a huge number of flexible components, the robotic manipulator is rather simple for MSA analysis. Since the number of flexible links in manipulator is relatively small, the assembling stage can be simplified and the columns/rows merging operations can be avoided and replaced by adding relevant constraints.

As it was shown above, the MSA equations for robotic manipulator are derived from three main sources: (i) link models, (ii) joint models and (iii) boundary conditions. The first of them describes the force-displacement relations for all links (both flexible and rigid) yielding 12 scalar equations per link relating 24 variables { , , , }

, } { , } {} {} { , } {} { , } { } {} {} {} {} { } { } , r r j ij r r ij ij p ij r r i ij ij i e e e ij ij ij r ij p e e ij ij e ij ij r ij e e ij ij                                                                     Λ ΛΛ Λ t Λ Λ Λ Λ Λ Λ Λ Λ Λ Λ W II K K I K I * 0 0 { } { } { } {} {} { } ext i ext i ext i ij ij ext i                                                   0 0 0 0 W 0 W 0 0 W W 0 0 0 0 W W (18)
After rearranging the matrix rows and introducing relevant definitions for the blocks, the global stiffness model can be presented as and to obtain the desired Cartesian stiffness matrix.

Conclusion

The paper deals with the advancement of the Matrix structural analysis for robotics applications. In contrast to previous results, it is suitable for robots with flexible links, rigid connections, passive and elastic joints with external loadings and preloadings. The paper describes in detail required matrix transformations that allow the user to obtain desired force-deflection relation and torque/deflections in internal mechanical elements. Its allows user straightforward aggregation of stiffness model equations avoiding traditional column/row merging procedures in the extended stiffness matrix.

  W 0 . In the case when the springs are initially preloaded by the wrench 0 ij W , the force equilibrium equation must be presented in the form

   t t W W . The second group of equations ensures the displacement compatibility and force/torque equilibrium for each internal connection (both rigid, passive and elastic), it includes two types of relations written for { , ,...} ij  tt and { , ,...}ijWWseparately. Independent of the connection type, each joint provides 12 scalar equations if it connects two links, 18 scalar equations for the connection of three links, etc. The third group of equations is issued from the manipulator connections to the environment, they are presented as the force/displacement constraints for certain nodes that give up to 6 scalar equations for { , ,...} ij  tt or at least 6 scalar equations for { , ,...} ij WW depending on the connection type. To derive the global stiffness model of the considered manipulator, let us combine all displacements and all wrenches in single vectors and denote them as {} i t and {} i W . This allows us to present the desired model in the form of the following block-matrix equation of sparse structure

  the node variables {} i t in two groups corresponding internal nodes m t and to the end effector e t one can get the following linear system

Modeling of the rigid link. If the link flexibility is negligible, the above stiffness model should be replaced by a "rigidity constraint". Applying to this link the

				( ) 11 ( ) 21 ij ii ( ) 12 ( ) 22 1 12 2 ij ij ij jj                 Wt KK Wt KK	  	(1)
	where	,  tt are the deflections at the link ends,	, WW are the link end
		ij								ij
	wrenches, i and j are the node indices corresponding to the link ends,
	and ( ) 11 ij KKKK are 66  stiffness matrices. The linear relation (1) should ( ) ( ) ( ) 12 21 22 ,,, ij ij ij
	be represented in a global coordinate system.
											rigid
	body kinematic equation, one can get the following relations between the nodal
	displacements	[   p t	;	 φ	]	,	[   p t	;	 φ	]
			i	i		i			j	j	j
		( )	6 6					6 1	;	( )	33	[	()	]
											3 3	3 3	6	6

  the matrix *

		e Λ of size	6 e ,	e	6  corresponds to the non-rigid direc-r
		ij
	ij p tions of the joint (similar to * Λ ) and e ij K is	e e  stiffness matrix describing elastic
	properties of the joint. In the compact form, the letter equations can be presented as
							i
	6 6		6 6		6 6	6 6	j
	e	e	e	e		e	(6 ) 1 e
	ij	** ij ij i	* ji j	e	e 6 (6 ) 24	i
							i	2 41

Acknowledgments. The work presented in this paper was supported by the grant of Russian Science Foundation №17-19-01740.