OPTIMAL RELEASE OF MOSQUITOES TO CONTROL DENGUE TRANSMISSION

. In order to prevent the propagation of human diseases transmitted by mosquitoes (as dengue or zika), a solution is to release mosquitoes infected by Wolbachia . In this study, we model the release and the propagation over time and space of such infected mosquitoes in a population of uninfected ones. The aim of this study is to investigate the best location in space of the release to ensure invasion by the infected mosquitoes. R ´ESUM´E. Aﬁn de pr´evenir la propagation de maladies transmises `a l’homme par les moustiques (comme la dengue ou le zika), une solution consiste `a relˆacher des moustiques infect´es par la bact´erie Wolbachia . Dans cette ´etude, nous mod´elisons le relˆacher et la propagation dans le temps et l’espace de ces moustiques infect´es dans une population de moustiques hˆotes non infect´es. Le but de cette ´etude est d’´etudier le meilleur emplacement dans l’espace des relˆachers aﬁn d’assurer l’invasion par les moustiques infect´es.


Introduction
Aedes aegypti is the main vector transmitting dengue viruses.This mosquito can also transmit chikungunya, yellow fever and Zika infection.According to the World Health Organization, 390 million people are infected by dengue every year and 3.9 billion people, in 128 countries, are at risk of infection by dengue viruses.As there is no treatment for dengue fever, the current method of preventing dengue 1 Sorbonne Université, CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France.luis@ann.jussieu.fr,martin.strugarek@gmail.com.
The authors acknowledge supports by the Project Analysis and simulation of optimal shapesapplication to life-sciences of the Paris City Hall.N.V. and J.P.Z.acknowledge support from the France-Brasil network for mathematics.
virus transmission and epidemics is to target the vector, i.e. the mosquito.Beyond preventing mosquitoes from accessing egg-laying habitats by environmental management and modification, one of the most promising control techniques is to transform mosquito population with a virus-suppressing Wolbachia bacteria.The idea of using Wolbachia for disease control was first proposed in the 1960s [8] but applying it to Aedes aegypti population is very recent.Wolbachia bacterium strains were isolated from Drosophila melanogaster in laboratory just before 2000 [7,10] but were introduced into Aedes aegypti embryos only on 2009 [11].The capability of this bacteria to suppress dengue virus and other pathogens transmission by Aedes aegypti was shown in laboratory around 2010 [12,2,21].It was also shown that this bacteria shortens life span [22] and most of the infected adults do not reach the infectious stage.But the most important modification induced by the bacteria is cytoplasmic incompatibility (CI) [11].Cytoplasmic incompatibility is used by the bacteria to spread rapidly into natural population [18] by producing non-viable eggs when uninfected females mate with infected males.Reproduction between infected males and females lead to infected eggs.As this bacteria is vertically transmitted (from mother to off-springs), uninfected males mating with infected females give rise only to infected eggs.
Formally, a proportion 1 − s h of uninfected female's eggs fertilized by infected males actually hatch.Cytoplasmic incompatibility is complete when s h = 1.We denote by b 1 , respectively b 2 , the net fecundity rate of uninfected females, respectively infected females.Death rate for uninfected mosquitoes is denoted d 1 .As Wolbachia decreases lifespan, death rate of infected mosquitoes d 2 verifies d 2 > d 1 .Is is also observed that Wolbachia infected mosquitoes tend to have reduced fertility, then b 2 ≤ b 1 .Finally, we denote κ the carrying capacity.Cytoplasmic incompatibility and vertical transmission drive the spatial spread of the infected population producing a bistable dynamic of Wolbachia [19].If the infected population is installed above a sufficient threshold frequency Θ compared to the uninfected population, it will spread and tend to increase to 1, otherwise it will tend to decline to zero.
We are interested on optimizing the release of Wolbachia-infected mosquitoes into a wild host population of mosquitoes.We denote u the release function.
For fixed maximal time T > 0 and domain Ω, the system of equation that we consider is the following: (1) in Ω.
The system of equations (1) models the propagation across time and space of the infected-mosquito population n 2 into an already existing uninfected mosquito population n 1 .The two coupled equation driving the dynamics of n 1 and n 2 are classical bi-stable reaction-diffusion equations.Note that in the reaction term of the first equation the term − n2 n1+n2 stands for the vertical transition of the disease whereas the coefficient s h models that this vertical transmission may or not be perfect because of the cytoplasmic incompatibility.The diffusion coefficient is denoted D; it is assumed to be the same for both population since both populations belongs to the same genus of mosquitoes.The last term of the second equation +u stands here to model the releases of infected mosquitoes developed in laboratory: it is on this control that we will act upon.More precisely, a question we want to address in this work is to know what should be the shape of the release function u to be as close as possible to the total invasion of the infected population into the domain.
The outline of this paper is the following.In the next section, we introduce the optimal control problem and prove the existence of an optimum for this problem.In Section 2, we consider a toy problem, which is a very simplified version of the full problem, for which we can solve explicitly the optimal problem and find the optimum.In Section 3, we investigate numerically the optimization of the spatial releases of mosquitoes.Finally, we end this paper with a conclusion and perspective for future works.An appendix is devoted to recalling the reduction of system (1).

Optimal Control Problem
We are going to simplify the problem.Instead of studying the coupled equations (1), we are going to follow the proportion of mosquitoes p(t, x) = n2(t,x) n1(t,x)+n2(t,x) as in [15].This reduction is clearly justified in the limit of large population in [15] (see also [1,Section 2.3]).The formal computation is explained in the Appendix.In order to simplify the reading, we perform the scaling x = x √ D not to keep the diffusion coefficient along the computations.Obviously, for the numerical simulations performed in Section 3, we have to keep in mind this scaling.
Denoting by p the proportion of infected mosquitoes, and u the release function, the dynamics is governed by the reaction-diffusion equation ( 2) The general optimal control problem we want to investigate involves the leastsquares functional J defined by ( 4) which models that one aims at steering the system as close as possible to the target state.In some sense, it stands for the research of a control strategy ensuring the persistence of infected mosquitoes at the time horizon T .Of course, it is relevant from the biological point of view to impose several constraints on the control function u.Indeed, the production of Wolbachia-infected mosquitoes is limited, which imposes that the total number of mosquitoes released is bounded.Hence, the control function u is assumed to belong to the set (5) modeling an upper limit on the instantaneous number of Wolbachia-infected individuals released at time t, as well as on the total number of released mosquitoes.
We then deal with the following optimal control problem: Since this problem involves the minimization over function depending on time and space variables, it is difficult to study.Then, we will reduce it to a simpler one by assuming that the time distribution of the control function is given.
1.1.Modeling of the optimal control problem.In order to weaken the difficulty of Problem (P full ), we introduce a simpler, although still relevant, problem by assuming that: • releases are done periodically in time (for instance every week) and are impulses in time1 ; • at each release, the largest allowed amount of mosquitoes is released, corresponding to the maximal production capacity per week (which is relevant, according to the comparison principle).As a consequence, we will be interested in determining the optimal way of releasing spatially the infected mosquitoes.Let us denote by t 0 = 0 < t 1 < . . .< t N = T , t i = i∆T , the release times.Rewriting the L 1 constraint on the control as u, 1 D ,D((0,T )×Ω) ≤ C, the control function reads where the pointwise constraint is modified into 0 The new optimal design problem reads It is possible to recast System (2) without source measure terms, coming from the specific form of the control functions.For the sake of simplicity, we provide here a naive formal analysis, but claim that this can be proven rigorously by using a standard variational analysis.Let us approximate the Dirac measure at t = t i by the function 1 ε 1 [ti,ti+ε] .Making the change of variable t = t i + τ , and introducing p given by p(τ, x) = p(t, x), one gets from system (2) that p solves Letting formally ε go to 0 and denoting, with a slight abuse of notation, still by p the formal limit of the system above yields Let us denote G the anti-derivative of 1 g vanishing at 0, namely Then, by a direct integration of ( 6) on [0, 1], we obtain Coming back on the function p yields Hence we arrive at the system (7) and the optimization problem reads where p is the solution of (7).In the next Section, we investigate the existence of solutions for this problem.
Theorem 1. Problem (P reduced ) has a solution.
Proof.For the sake of readability, we only provide the proof in the case N = 2. Indeed, there is no additional difficulty to deal with the general case whose proof follows exactly the same lines.The proof is divided into several steps.Let ) N be a minimizing sequence for Problem (P reduced ).
Notice that, since u belongs to V T,C,M and G −1 takes its value in [0, 1[, we infer from the maximum principle that 0 ≤ p(t, •) < 1 for a.e.t ∈ [0, T ] so that one has for all u ∈ V T,C,M It follows that inf u∈V T ,C,M J(u) belongs to (0, |Ω| 2 ) and, in particular, is finite.
Step 1: Convergence of the minimizing sequence.Let p n be the solution to (7) associated to the control function u n and let us introduce . By induction, one easily shows that v n is uniformly bounded in L ∞ .Since the class V T,C,M is closed for the L ∞ weak-star topology, there exists v ∞ ∈ V T,C,M such that, up to a subsequence, v n converges weakly-star to v ∞ in L ∞ .Here and in the sequel, we will denote similarly with a slight abuse of notation a given sequence and any subsequence.
Multiplying the main equation of ( 7) by p n and integrating by parts, we infer from the above estimates the existence of a positive constant C such that for every n ∈ N, which also reads for every n ∈ N.
By using the pointwise bounds on p n , it follows that p n is uniformly bounded in L 2 ([0, T ], H 1 (Ω)).Furthermore, by using (7), one gets that the sequence ∂ t p n is uniformly bounded in L 2 ([0, T ], W −1,1 (Ω)).According to the Aubin-Lions theorem (see [14]) we infer that p n converges (up to a subsequence) to and weakly in L 2 ([0, T ], H 1 (Ω)).Passing to the limit in (7) yields that p ∞ is a weak solution to (8) It is standard that any solution to this bistable reaction-diffusion equation is continuous in time. Introducing .) by passing to the limit as n → +∞ in the variational formulation on p n , using adapted test-functions belonging to

This is a consequence of the weak convergence of p
, which is a consequence of the continuity and convexity 2 of G for p ∈ [0, 1).
Step 2: Conclusion.Let us first show that u ∞ belongs to V T,C,M .Since the derivative of G is 1/g which is positive, G is increasing and therefore, one has 0 ≤ u ∞ ≤ m a.e. in Ω.For the integral condition (namely, Ω u ≤ C), let us distinguish between two cases: Case 1: if m|Ω| ≤ C, the conclusion follows immediately. 2Indeed, one has G (p) = κb 2  (1 − s h p 2 ) ((p − 1)(s h p − 1)) 2 .which is positive whenever p belongs to [0, 1].
Case 2: if m|Ω| > C, let us use that G is, as aforementioned, lower semicontinuous for the weak-star topology of L ∞ .Thus, we deduce that It follows that u ∞ belongs to V T,C,M and one concludes by using the Fatou Lemma: We finally infer that u ∞ solves Problem (P reduced ).
Remark 1.The uniqueness issue remains open, even for simple domain.It is likely that symmetries of the release domain play an important role.
It is interesting to notice that, in a very particular case, we have an explicit expression of the minimizer for this problem.Let p M be the solution of (2) associated to u M identically.Since G −1 is an increasing function by the comparison principle we have for all time t ∈ [0, T ] and a.e.
Evaluating this expression at time t = T , the expected conclusion follows by noting that the constant function equal to M on (0, T ) × Ω belongs to U T,C,M .
1.3.Computation of derivatives.As a preliminary remark, we claim that for any element u of the set V T,C,M and any admissible perturbation h, the mapping , where p denotes the unique weak solution of (7), is differentiable in the sense of Gâteaux at u in the direction h.Indeed, proving such a property is standard in calculus of variations and rests upon an application of the implicit function theorem.In the sequel, and with no confusion possible, we will denote by ṗ the Gâteaux-differential of p at u in direction h and by dJ(u), h the Gâteaux-differential of J at u in direction h, namely Let us make the cone of admissible perturbations precise.We call "admissible perturbation" any element of the tangent cone T u,V T ,C,M to the set V T,C,M at u. Definition 1.The cone T u,V T ,C,M is the set of N -tuples h = (h 0 , . . ., h N −1 ) ∈ (L ∞ (Ω)) N such that, for any i ∈ {0, . . .N − 1} and for any sequence of positive real numbers ε n decreasing to 0, there exists a sequence of functions h n i ∈ L ∞ (0, T ) converging to h i as n → +∞, and u i + ε n h n i ∈ V T,C,M for every n ∈ N (see e.g.[5]).
where q is the unique solution of the backward problem Proof.By using the preliminary discussion, one has ( 9) where ṗ denotes the unique solution of the system (10) Let us multiply the main equation of this system by q and then integrate by parts with respect to the variables t and x.By using in particular the Green formula, we get successively that and therefore, yielding the desired conclusion.
Remark 2. For practical purposes, it may be useful to notice that where q denotes the solution of the initial boundary value problem

A toy Problem
This section is devoted to investigating a simpler version of (P full ) corresponding to the case N = 1 with f = 0.More precisely, let p be the solution of ( 11) Then, the optimization toy problem reads is the unique solution of Equation (11).Note that the equation ( 11) has to be understood in a weak sense, since u 0 ∈ L ∞ (Ω) ⊂ L 2 (Ω) (see for example [17,Section 10.7]).
For this simple problem, we are able to solve explicitely the optimization problem : Theorem 2. Problem (11) has a unique solution u 0 , which is constant and equal to min 1, M, C |Ω| .Proof.First, note that Problem (P toy ) has a solution.Indeed, it is standard that the mapping L 2 (Ω) u 0 → p ∈ C 0 ([0, T ], L 2 (Ω)) is continuous.Therefore, so is Ĵ by composition of continuous mappings.The conclusion follows by observing that V T,C,M is a compact subset of L 2 (Ω).
The proof relies on a well-adapted rewriting of the criterion Ĵ.For that purpose, let us introduce the Neumann operator −∆ N on Ω defined on According to the spectral theorem, there exists an orthonormal family (φ j ) j≥1 consisting of (real-valued) eigenfunctions of −∆ N , associated with the non-decreasing sequence positive eigenvalues (λ j ) j≥1 .Moreover, by setting λ 0 = 0 and φ 0 = 1 √ |Ω| , the sequence (φ j ) j≥0 is a Hilbert basis of L 2 (Ω) and any solution p of (11) can be expanded in a unique way in L 2 (Ω) as ( 12) u 0j e −λj t φ j (x), with u 0j = u 0 , φ j L 2 (Ω) .By expanding the square in the definition of Ĵ, we then infer that

Gaussian Releases
From a practical point of view, not all controls u ∈ V T,C,M correspond to a release that could actually be conducted, as for example the constant solution of the toy problem of the previous section.To guarantee a solution that could be implemented, we restrict here the admissible controls to more accurately model the way mosquitoes are released in practice.
We thus consider that there are K ∈ N simultaneous releases and that each one results in a Gaussian distribution of mosquitoes centered around the position of the release x k ∈ Ω for k = 1, ..., K.Then, the feasible controls are of the form (14) u where the constants m and σ are chosen such that u K (•, x 1 , ..., x K ) ∈ V T,C,M .In particular, we choose to saturate the constraint on the total number of mosquitoes released, i.e. we take Ω u K (x)dx = C.The goal is then to find the best position of the releases and the optimization problem becomes where p ∈ C 0 ([0, T ], L 2 (Ω)) is the unique solution of ( 7) with control u K (•, x 1 , ..., x K ).

Remark 3.
Since Ω is a bounded domain in R 2 , the question of the existence of a minimizer is trivial.But, the uniqueness is still a challenging problem.
Proposition 3. Let (x 1 , ..., x K ) ∈ Ω K .For k ∈ {1, ..., K}, one has where q is the unique solution of the backward problem Proof.It is an easy application of the chain rule.First, we notice that Next, using Proposition 2, we find thanks to the chain rule that for all k ∈ {1, ..., K} We deduce the result from the last equality.
3.1.Numerical Resolution.We now present the computation of the numerical solution of (P K ).For this we use a direct method which consists in carrying out a discretization of Equation ( 7) and of the control in order to obtain a finite dimensional optimization problem with constraints.We can then compute an approximation of a local minimizer of (P K ) with a numerical optimization solver.
Our results were obtained with the finite element toolbox FreeFem++ [6] which contains an implementation of the optimization routine Ipopt [20].We therefore consider a finite element basis of functions (ϕ i ) i that allows us to discretize the control as u h (x, x 1 , ..., x K ) = i u i ϕ i (x) and the proportion of infected mosquitoes as p h (t, x) = i p i (t)ϕ i (x), the finite element approximation of the solution of the PDE (7) with initial condition G −1 (u h (x, x 1 , ..., x K )).The cost function can be computed with numerical integration as J h (x 1 , ..., x K ) = Ω (1 − p h (T, x)) 2 dx.In addition, Ipopt requires the gradient of the cost function and thanks to Proposition 3 we have where q h (0, x) is the finite element approximation of the solution of the backwards PDE.
Remark 4. Because of Proposition 1, we were interested in the case M > C |Ω| .In addition, we have fixed C such that the constant solution u = C |Ω| leads to extinction (as T tends to +∞) but there exists R ∈]0, C πM [ such that the function u(x) = M × 1 B(0,R) (x) belongs to V T,C,M and leads to invasion (as T tends to +∞).
We now present numerical simulations for the parameters given in Table 1.The birth and death rates are given per day, whereas the unit of the carrying capacity is per m 2 and the diffusion coefficient is given per m 2 per day.The numerical values are taken from [1] and references therein.We consider a square domain of 1 hectare, a final time of 200 days and we set the total amount of mosquitoes released such that C < G(θ)|Ω|.In Figure 1 we show the control u K (•, x 1 , ..., x K ) for K = 3, 4, 5, 6 releases and for each case the same total amount of mosquitoes is released.For the case of 6 releases we display in Figure 2 the time dynamics of the proportion of infected mosquitoes p(t, •).As expected, it leads to the total invasion of the domain.

Conclusion
We investigate in this work the optimization of the release of Wolbachia-infected mosquitoes into a host population in the aim to replace the wild population by a Wolbachia-infected population unable to transmit several diseases to human.To conduct this study, we first reduce the optimal problem under investigation by assuming that the time distribution is given.Then we obtain existence of a minimum for this latter problem.Finally, reducing again the control problem by considering that the releases are modelled by Gaussian distributions, some numerical computations are performed.
Optimisation strategies for release protocols of mosquitoes have been investigated by several authors [16,4,3].However, in these papers, only the time optimization  of the releases is investigated.Up to our knowledge, this work is the first attempt in optimizing spatially the releases, which is of great interest for experiments in the field.The preliminary results obtained in this paper should be continued.In particular, the optimality conditions for the system (P K ) should be studied in a future work in the aim to find properties of the optimal solution.The numerical simulations should also be continued to better representation of what is observed in the field.

Proposition 1 .
Let N ∈ N * and M ≤ C |Ω| .Then u = M is the unique solution of Problem (P reduced ).Proof.It is a direct application of the comparison principle.Let u * be a solution of Problem (P reduced ).By contradiction, let us assume that u * is not identically equal to M a.e. in Ω.Then, let t i be a release time for which the associated control function u * i is not identically equal to M in Ω. Recall that u * i ≤ M .Let us denote by p * the solution of the problem (2) associated to the control function u * .Let u M be the control function defined by u M i = M and u M j = u * j for all j ∈ {0, . . ., N − 1}\{i}.