Optimization of spatial control strategies for population replacement, application to Wolbachia

This article is dedicated to our friend Enrique Zuazua on the occasion of his 60 th birthday. Abstract In this article, we are interested in the analysis and simulation of solutions to an optimal control problem motivated by population dynamics issues. In order to control the spread of mosquito-borne arboviruses, the population replacement technique consists in releasing into the environment mosquitoes infected with the Wolbachia bacterium, which greatly reduces the transmission of the virus to the humans. Spatial releases are then sought in such a way that the infected mosquito population invades the uninfected mosquito population. Assuming very high mosquito fecundity rates, we ﬁrst introduce an asymptotic model on the proportion of infected mosquitoes and then an optimal control problem to determine the best spatial strategy to achieve these releases. We then analyze this problem, including the optimality of natural candidates and carry out ﬁrst numerical simulations in one dimension of space to illustrate the relevance of our approach.

However, the suppression of one population of insects might have consequences on the environment. Then, other approaches aim at replacing the wild population of mosquitoes by another population inoffensive to human. One strategy under investigation consists in using the bacteria Wolbachia taking advantage of phenomena called cytoplasmic incompatibility (CI) and pathogen interference (PI) [7,30]. In key vector species such as Aedes aegypti, if a male mosquito infected with Wolbachia mates with a non-infected female, the embryos die early in development [38]. This is the so-called cytoplasmic incompatibility (CI). Moreover, it has been observed that Aedes mosquitoes infected with some Wolbachia strains are not able to transmit viruses like dengue, chikungunya and zika [36], this is the pathogen interference (PI). Then, one may release mosquitoes artificially infected by Wolbachia to mate with wild ones. Over time and if the releases are large and long enough, it can be expected that the majority of mosquitoes will carry Wolbachia, due to cytoplasmic incompatibility. As a result of PI, the mosquito population then has reduced vectorial competence.
In this paper, we focus on the Wolbachia strategy and investigate the question of optimizing the spatial distribution of the releases. Several mathematical models have been proposed for the Wolbachia technique, see e.g. [13,14,19,28]. In these papers, the authors model the time dynamics of the mosquitoes population. Then, the question of optimizing the time of releases has been investigated e.g. in [1,2,6,8]. However the spatial distribution of mosquitoes may have an impact on the success of the strategy. It is therefore relevant to add spatial dependence in mathematical models, which makes the study much more complicated.
In order to have a model simple enough to be tractable from a mathematical point of view, the authors in [4] introduce a model focusing only on the proportion of Wolbachia-infected mosquitoes, denoted p in the sequel: p := n in n in + n un where n in is the density of Wolbachia-infected mosquitoes and n un the density of uninfected mosquitoes. This quantity solves a scalar reaction-diffusion equation where D is a diffusion coefficient and f is a bistable function. 1 For this model, the conditions to initiate the spatial spread are well-known [33]. It has been proved later in [32] that this model may be rigorously derived from a more general system governing the dynamics of Wolbachia-infected and Wolbachia-uninfected mosquitoes by performing a large fecundity asymptotics.
In this study, we are investigating the question of the best spatial strategy for mosquito release, i.e., giving a certain amount of mosquitoes, we are trying to determine optimal locations to release them in order to ensure the invasion of the environment by Wolbachia-infected mosquitoes. If we denote u the release function, then the above model is modified into x) = f (p(t, x)) + u(t, x)g(p(t, x)), t ∈ (0, T ), x ∈ Ω, ∂ ν p(t, x) = 0, t ∈ (0, T ), x ∈ ∂Ω, p(0, x) = 0, x ∈ Ω, (1.1) where Ω is an open bounded connected subset of R d with a regular boundary ∂Ω. The function g is positive and vanishes when p = 1. The derivation of (1.1) will be detailed in Section 2.
Let us summarize the main assumptions on f and g we will use in the sequel.
(H f,g ) A first study was carried out in [3], giving rise to the first very simple numerical experiments. In the present article, we seek to complete the results of this study, by analyzing qualitatively the solutions and by proposing adapted numerical strategies. Let us mention that a problem of the same nature has been investigated in [24], mainly from a numerical point of view. Authors characterize optimal vaccination strategies to minimizes the costs associated with infections by the Zika virus and vaccines in the state of Rio Grande do Norte in Brazil.
In [25], an optimal control problem close to the one investigated hereafter is tackled. The authors consider a population whose evolution is driven by a reaction-diffusion equation and look at determining initial data submitted to L 1 and L ∞ constraints, maximizing the total size of the population. In our article, we choose to deal with a least square criterion instead of the average criterion considered in [25] (and more recently in [23]). For reasons that will appear later, constant solutions are natural candidate to solve the considered optimal control problem. We show, as in [25], that for certain families of parameters, the constant functions are local minimizers of the optimal control problem. On the other hand, we complete this first analysis and also manage to show for our model, that these same functions can be or not be global minimizers depending on the considered range of parameters. These results are also illustrated numerically. Finally, it is worth mentioning that exact controllability issues for similar reaction-diffusion systems have been investigated in [20,22].
The outline of the paper is the following. In Section 2, we present the derivation of system (1.1) and present the optimal control problem we are looking at. Section 3 contains the main mathematical results of this paper. Their proofs are given in Section 4. More precisely, the rigorous derivation of system (1.1) is explained in Section 4.1 and Section 4.2 is devoted to the mathematical study of the optimal control problem. Finally, numerical illustrations with the description of the numerical algorithm are provided in Section 5.

Modelling
In the whole article, we will consider a given bounded connected open domain Ω of R d assumed to have a Lipschitz boundary. Let T > 0 denote a fixed horizon of time.

Model with two compartments
In order to justify the introduced model on the proportion of Wolbachia-infected mosquitoes, we first explain how to derive it. Let us denote n in the density of infected mosquitoes and n un the density of uninfected mosquitoes. The dynamics of these quantities is governed by the reaction-diffusion system complemented by initial conditions n in (t = 0, x) = n init in (x) 0, n un (t = 0, x) = n init un (x)> 0, where the following notations are used: -u: instantaneous releases of Wolbachia infected mosquitoes. It is on this control that we will act upon. At this step, we do not make the admissible space of controls precise, this will be done in what follows; -d un , d in = δd un with δ > 1: death rates, respectively for uninfected and infected mosquitoes. We assume that d in > d un since Wolbachia decreases lifespan; -F un , F in = (1 − s f )F un : net fecundity rates, respectively for uninfected and infected mosquitoes. We assume that F in < F un since Wolbachia reduces fecundity; -ε : parameter without dimension quantifying the fecundity, we assume ε 1 meaning that the fecundity is considered to be large; -s h ∈ (0, 1): cytoplasmic incompatibility parameter (fraction of uninfected females' eggs fertilized by infected males which will not hatch). Formally, a proportion 1 − s h of uninfected female's eggs fertilized by infected males actually hatch. Cytoplasmic incompatibility is perfect when s h = 1; -K: carrying capacity; -D: dispersal coefficient.
All the constants above are assumed to be positive. Existence and uniqueness of solutions for such reactiondiffusion system is by now well-known see e.g. [12,27]. The equations driving the dynamics of n in and n un are bistable and monostable reaction-diffusion equations, respectively. Note that in the reaction term of the second equation, the term − nin nin+nun 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.
In accordance with [2], we will assume moreover that the relation holds true. It is notable that such a parameters choice is relevant since for Wolbachia-infected Aedes mosquitoes, and more precisely in the case of wMel strain, CI is almost perfect in these species-strain combination (see [10]) meaning that s h is close to 1. Furthermore, such mosquitoes typically have a slightly reduced fecundity. In that particular case, one has s f 0.1, δ 1.1 and s h 0.9 so that (2.2) holds true. To model optimal strategies with an adapted optimal control problem, it is convenient to introduce the Wolbachia-infected equilibrium (n * in,W , 0) for the uncontrolled system, defined by that is (n * in,W , 0) is a stationary solution of (2.1a)-(2.1b). A possible approach hence consists in looking for controls steering the system as close as possible to the target state (n * in,W , 0). In some sense, it stands for the research of a control strategy ensuring the persistence of infected mosquitoes at the time horizon T . This leads to define the least squares functional J T given by where (n in , n un ) denotes the unique solution to the reaction-diffusion system (2.1a). Here, we use the notation x + = max{x, 0}. Observe that the presence of this maximum in the definition of J T does not induce nondifferentiability since the mapping x → x 2 + from R to R is C 1 .

Reduction for large fecundity
When the fecundity is large compared to other parameters, it is relevant to consider the asymptotics ε → 0, which allows us to reduce system (2.1a)-(2.1b). This reduction is inspired by [2] where the authors consider a differential system. We first explain formally how to reduce this system and state the main result, the rigorous approach is postponed to Section 4.1. Since n in and n un will depend on ε, we use the notation n ε in and n ε un .
We introduce the notation F for the antiderivative of f , In what follows, we will assume: This assumption is necessary to guarantee that invasion of the infected population may occur in space by local release. We will check that this assumption is satisfied for the particular choice of parameters we will consider for the numerical experiments in Remark 5.1.
We consider system (2.1a)-(2.1b) with Neumann boundary conditions to model that the boundary acts as a barrier, and initial conditions satisfying We assume also that the initial conditions are well-prepared, i.e. (2.14) A typical example of initial conditions is when the system is as the Wolbachia-free equilibrium for which n init,ε in = 0 and n init,ε un = K(1 − εdun Fun ). In this case, assumption (2.14) is obviously satisfied.
Convergence result. Following the ideas in [32], where a similar asymptotic limit is performed, we derive an asymptotic model on the proportion of infected mosquitoes, as the fecundity rates tend to +∞.
Let us now define the least squares functional J ε T given by where (n ε in , n ε un ) have been introduced at the very beginning of Section 2.2 and n * in,W in (2.3). As a corollary of the convergence result above, let us make the asymptotic behavior of the functional J ε T precise. Corollary 2.3. Under the assumptions (2.13)-(2.14) on the initial data, let (u ε ) be a sequence converging towards u 0 , weakly star in L ∞ ((0, T ) × Ω) as ε 0. Then, p ε (T, ·) converges towards p 0 (T, ·) strongly in L 2 (Ω), and moreover, J ε T (u ε ) converges towards J 0 T (u 0 ) defined by where p denotes the solution of (2.8), as ε 0.
The proof of this result is postponed to Section 4.1. 3.
In what follows, we will rather deal with the proportion p 0 to model optimal releases strategies. The following section is dedicated to modeling issues about the optimal control problem we will deal with.

Toward an optimal control problem
In this section, we will introduce an optimal control problem modeling optimal mosquito releases. For this purpose, we assume all fecundity rates large, which legitimates the use of the asymptotic model (2.8)-(2.9) introduced in Section 2.2. We will focus on time-pulsed releases, which will lead us to further simplify the problem.
In order not to cumulate all the difficulties related to the search for release distributions in time and space, we will suppose that one release, which is an impulse in time, 2 is done at the beginning of the experiment, i.e. u(t, x) will be assimilated to a particular approximation of a Dirac impulse in time, namely u 0 (x)δ {t=0} . More precisely, we will consider as choice of release term, the function where u 0 ∈ L ∞ (Ω) will be given. Making the change of variable t = τ η, and introducingp given byp(τ, x) = p 0 (t, x), one gets from system (2.8) thatp solves We now provide a purely formal argument to justify the optimal control problem we will deal with. Letting formally η go to 0 and denoting, with a slight abuse of notation, still byp 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 (2.17) on [0, 1], we obtain Hence we arrive at the system where f and g are given by (2.9).
To take into account biological constraints on the release procedure, we will moreover assume that the release function is such that: -the local release of mosquitoes is bounded : 0 ≤ u 0 ≤ M a.e. in Ω with M > 0; -the total number of used mosquitoes is bounded (production limitation), reading with C ∈ (0, M T ). Note that it is relevant to choose the parameter C strictly lower than M T . In the converse case, it would mean that the choice u 0 (·) = M is admissible, so that the local maximal number of mosquitoes can be released (almost) everywhere in Ω. Since producing infected mosquitoes has an important cost, it is reasonable from a biological point a view to assume that such a release is not possible.
This leads to introduce the admissible set V C,M given by The goal is to be as near as possible to the equilibrium p = 1 at time T . Let us denote (with a slight abuse of notation) by J T , the least squares functional defined by Observe that coincides, up to a positive multiplicative constant, with the asymptotic functional J 0 T given by (2.16). The optimization problem thus reads where p is the solution of (2.18). From now on and without loss of generality, we will assume in what follows that the diffusion coefficient D is equal to 1.

Main results
Constant solutions are natural candidates to solve Problem (P reduced ). Indeed, it has been observed in ( [3], Thm. 2.1) that in the very simple case where f (·) = 0 and G : x → x, Problem (P reduced ) has a unique solution u 0 , which is constant and equal to min 1, M, C |Ω| . Furthermore, as stated in the following result, constant solutions equal to M are optimal for a given range of the parameters. We show moreover that, outside of this range, constant functions remain critical points and show that they are still local minimizers whenever C is small enough. We also comment on the sharpness of this result by highlighting that for certain parameters, constant functions may not be global minimizers for Problem (P reduced ).
According to Corollary 4.8, it is enough to concentrate on the constant function equal to C/|Ω|.
Let us assume that f and g satisfy (H f,g ). Problem (P reduced ) has a solution.
(i) For every M ∈ (0, C/|Ω|], the constant function u M equal to M is the unique solution to Problem (P reduced ). (ii) Let us assume that M |Ω| > C. The constant function u(·) = C/|Ω| is a critical point for Problem (P reduced ) (meaning that it satisfies the first order optimality conditions stated in Proposition 4.7). Furthermore, if C ≤ |Ω|G(θ), there exists K T > 0 such that for every h ∈ L 2 (Ω), the second order differential of J T at u satisfies and it follows that the function u is a local minimizer for Problem (P reduced ).
Let us comment on the sharpness of Theorem 3.1. As will be emphasized hereafter and in Section 5, we do not expect that u solves Problem (P reduced ) for all values of C ∈ (|Ω|G(θ), M |Ω|). In some case, this will be confirmed numerically, by using u as starting point of optimization algorithms and obtain at convergence a nonconstant minimizerũ such that J T (ũ) < J T (u).
Actually, even in the case C < min{G(θ), M }|Ω|, where we know from Theorem 3.1 that the constant solution u is a local minimizer, under some conditions on |Ω| and C, we may construct non constant initial date u 0 such that J T (u 0 ) < J T (u), as stated below in Proposition 4.13.
Recalling that θ c ∈ (θ, 1) is defined by θc 0 f (p) dp = 0. We assume The following result shows that under some conditions on Ω, M and C, the constant function u is not a global minimum of the optimization problem (P reduced ).
Proposition 3.2. Let us assume (3.1) and that: -C is large enough; Graphs of the function f , its second order derivative f and the function g by using the data from [15] (see Table 1).
-the inradius 3 of Ω is large enough; -T is large enough.
Then the constant solution u is not a global minimum for Problem (P reduced ).
A proof of this result is provided in Section 4.2.4.
Remark 3.3. The conditions stated in Proposition 3.2 are not sharp; the obtention of necessary and sufficient condition for constant solution to be a global minimizers seems to be intricate and we let it open. A related problem concerns the issue of finding sufficient and necessary conditions guaranteeing invasion in a bistable reaction-diffusion system that is, up to our knowledge still open, and we refer to [25] for partial answers in this direction.
Remark 3.4. It is notable that, for the sets of parameters from [15] below, the functions f and g satisfy (H f,g ). Indeed, the function f vanishes once on (0, 1) and its root z satisfies θ < z 0.466, while θ c 0.582 (see Fig. 1).

Proofs
This section will be devoted to prove Theorems 2.2 and 3.1.

Model reduction
In this section, we will give a proof of Theorem 2.2, allowing us to reduce the system (2.1a)-(2.1b) to a scalar reaction-diffusion equation for the proportion as the parameter ε goes to 0. It is inspired by [2] in which the authors use a model composed by two differential equations.

Uniform a priori estimates
We first establish some uniform bounds with respect to ε > 0.
Lemma 4.1. Assume the assumptions of Theorem 2.2 hold. Let u ε be given in V C,M , ε ∈ (0, 1) and (p ε , n ε ) be the unique solution of (2.5)-(2.6). Then, 3 In other words, the radius of the largest ball inscribed in Ω.
Proof. By nonnegativity of n init in and n init un , it is standard to deduce the nonnegativity of n ε in and n ε un (Indeed 0 is a subsolution for (2.1a) and for (2.1b), see e.g. [9]). Moreover, since R + is invariant for the equation of n ε un and n init un is non identically equal to zero, we deduce that n ε un > 0 on Ω × (0, T ] (see e.g. [37, th. 2]). Therefore, Consider the function h defined in (2.7). We remark that the denominator is positive.
Lemma 4.2. Under above assumptions, for ε > 0 small enough, we have the uniform estimate and for some nonnegative constants C and C 0 .
Proof. On the one hand, multiplying equation (2.5) by εn ε and integrating on Ω, we get Since from Lemma 4.1, we know that n ε and p ε are uniformly bounded in L ∞ ([0, T ] × Ω), we deduce (4.2) after an integration in time.
On the other hand, we fix ε 0 > 0 small enough such that, for all ε ε 0 , we have |n ε | ≤ C 1 < 1 ε on [0, T ] × Ω for some constant C 1 > 0 (which is always possible thanks to Lem. 4.1). Then, we multiply by p ε the equation satisfied by p ε (2.6) and integrate over Ω, we deduce for some nonnegative constant C 3 . Then, using a Cauchy-Schwarz inequality, we get From (4.2) and the well-known inequality 2ab ≤ a 2 + b 2 , we deduce after an integration in time where we recall that u ε ∈ V C,M and p ∈ [0, 1]. Taking ε small enough, we get the desired estimate.

Compactness result and proof of Theorem 2.2
We first recall the following compactness result (see [29]).
Step 2. Passing to the limit. We now pass to the limit in the weak formulation of equations (2.5) and (2.6).
From the weak formulation of (2.5), we deduce that for any test function φ ∈ C ∞ c ((0, T ) × Ω), we have From the L ∞ -bound of Lemma 4.1, we deduce that the term of the left hand side and the last term of the right hand side converge to 0 as ε → 0. For the first term of the right hand side, we may pass into the limit thanks to the weak convergence of n ε , the strong convergence of p ε , and the weak convergence of u ε . We obtain, for As a consequence (2.7) is verified almost everywhere.
We are left to pass into the limit in the weak formulation of (2.6).
From the above convergence it is straightforward to pass into the limit into the first two terms of the left hand side. For the third term, we use estimate (4.2), and a Cauchy-Schwarz inequality to get as ε → 0, thanks to Lemma 4.1 and 4.2. We may pass into the limit for the first term of the right hand side of (4.3) since (p ε ) converges strongly and a.e., and (n ε ) converges weakly. Then, for the last term of the right hand side of (4.3), we verify .
From the strong L 2 convergence of (p ε ) and the L ∞ bounds in Lemma 4.1, we deduce that the last two terms go to 0 as ε → 0. For the first term, we write It is then straightforward to conclude the convergence towards 0 of the first term of the right hand side of (4.4).
Finally, passing into the limit ε → 0 into (4.3), we obtain We conclude by using the fact that (n 0 , p 0 , u 0 ) verifies the relation (2.7).

Proof of Corollary (2.3)
First, observe that p ε (T, ·) converges weakly-star to p 0 (T, ·) in L 2 (Ω). Indeed, the proof is standard. Consider the variational formulation (4.3) on p ε where test functions φ are now chosen in The variational formulation (4.3) is then modified by the addition of in the left-hand side term. Since p ε (T, ·) is bounded in L 2 (Ω) it converges weakly to some limit in L 2 (Ω) up to a subsequence. Passing to the limit in the variational formulation and using Theorem 2.2 allows us to identify the closure point of p ε (T, ·) as p 0 (T, ·). Finally, by uniqueness of the closure point, we infer that the whole sequence p ε (T, ·) converges to p 0 (T, ·). Since the L 2 (Ω)-norm is lower semicontinuous for the weak-topology, we get that Let us multiply the equation (2.6) by p ε and then, integrate it on (0, T ) × Ω, we obtain By assumption on the initial data, we have the convergence of the first term of the right hand side. For the second term of the right hand side, the weak convergence in L 2 ((0, T ), H 1 (Ω)) guarantee that this term is upper semi-continuous, then Due to the strong convergence in L 2 ((0, T ) × Ω) of the sequence (p ε ) ε and using also the uniform bound of the sequence (n ε ) ε (see Lem. 4.1), we deduce the convergence of the third term of the right hand side of (4.5). Using also the strong convergence of (p ε ) ε and the weak convergence of (n ε ) ε , we get the convergence of the fourth term in the right hand side of (4.5). Finally, from a Cauchy-Schwarz inequality we have Thanks to the estimates in Lemma 4.2, we get that this latter term goes to 0 as ε → 0. Finally, we have proved that It follows that p ε (T, ·) L 2 (Ω) converges to p 0 (T, ·) L 2 (Ω) as ε 0, and since p ε (T, ·) converges to p 0 (T, ·) weakly in L 2 (Ω), it follows that this convergence is in fact strong, whence the claim.

4.2.
Analysis of the optimal control problem (P reduced )

Existence of an optimal control
As a preliminary remark, note that existence of an optimal control has been shown in ( [3], Thm. 1.1) in a more general setting. To make this article self-contained, we recall the argument hereafter. The analysis to follow is valid under the assumption (H f,g ) on f and g. It is not restricted to the particular choice of functions f and g given by (2.9). Proof. In what follows, we will denote by p u0 the solution to Problem (2.18) associated to the control choice u 0 . Let (u n 0 ) n∈N ∈ (V C,M ) N be a minimizing sequence for Problem (P reduced ). Notice that, since u n 0 belongs to V C,M and the range of G −1 is included in [0, 1[, we infer from the maximum principle that 0 ≤ p u n 0 (t, ·) < 1 for a.e. t ∈ [0, T ] so that (J T (u n 0 )) n∈N is bounded and inf u0∈V C,M J T (u 0 ) is finite. Since the class V C,M is closed for the L ∞ weak-star topology, there exists u ∞ 0 ∈ V C,M such that, up to a subsequence, u n 0 converges weakly-star to u ∞ 0 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 (2.18) by p u n 0 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 By using the pointwise bounds on p u n 0 , one gets that p u n 0 is uniformly bounded in L 2 (0, T ; H 1 (Ω)). Furthermore, according to (2.18), the sequence ∂ t p u n 0 is uniformly bounded in L 2 (0, T ; W −1,1 (Ω)). The Aubin-Lions theorem (see [29] and Lem. 4.3) yields that p u n 0 converges (up to a subsequence) to p ∞ ∈ L 2 (0, T ; H 1 (Ω)), strongly in L 2 (0, T ; L 2 (Ω)) and weakly in L 2 (0, T ; H 1 (Ω)). Furthermore, using that the sequence ∂ t p u n 0 is uniformly bounded in L 2 (0, T ; W −1,1 (Ω)) also yields that ∂ t p ∞ belongs to L 2 (0, T ; W −1,1 (Ω)). Furthermore, reproducing the standard variational argument used in the proof of Corollary 4.1.3 to show the weak convergence of p ε (T, ·) to p 0 (T, ·) in L 2 (Ω) as ε 0, one shows that for all t ∈ [0, T ], p u n 0 (t, ·) also converges weakly, up to a subsequence, to p ∞ (t, ·) in L 2 (Ω).
Passing to the limit in (2.18) yields that p ∞ is a weak solution to

It is standard that any solution to this bistable reaction-diffusion equation is continuous in time.
It remains to show that u ∞ 0 = G(p ∞ (0 + , .)). Note first that G is convex since g is decreasing on (0, 1) under assumption (H f,g ). According to the convergence results above, since p u n 0 (0, ·) converges weakly (up to a subsequence) to p ∞ (0, ·) in L 2 (Ω) and since u n 0 := G(p u n 0 (0 + , .)), we get that G(p ∞ (0, ·)) = u ∞ 0 and hence, p ∞ = p u ∞ 0 by passing to the limit as n → +∞ in the variational formulation on p u n 0 . Finally, let us show that u ∞ 0 belongs to V C,M . Since the derivative of G is 1/g which is positive, G is increasing and therefore, one has 0 ≤ u ∞ 0 ≤ M a.e. in Ω. Moreover, since Ω u n 0 ≤ C rewrites u n 0 , 1 L ∞ ,L 1 ≤ C, we immediately get that the integral condition is satisfied by u ∞ 0 . Therefore, u ∞ 0 solves Problem (P reduced ).

First and second order optimality conditions
We now state first and second order optimality conditions. The objective is twofold: first, we will analyze the optimality of constant solutions, and second, we will use them to derive adapted numerical algorithms.
where q denotes the adjoint state, solving the backward p.d.e.
and p denotes the solution to (2.18) associated to the control choice u 0 . Furthermore, the second order derivative of J T at u 0 reads for every admissible perturbation h, whereṗ denotes the solution to the linear system x ∈ Ω. (4.7) Proof. As a preliminary remark, we claim that for any element u of the set V C,M and any admissible perturbation h, the mapping u ∈ V C,M → p ∈ L 2 (0, T ; H 1 (Ω)), where p u denotes the unique weak solution of (2.18), is differentiable in the sense of Gâteaux at u in direction h. Indeed, proving such a property is standard in calculus of variations and rests upon the implicit function theorem. Let u ∈ V C,M . Let h denote an admissible perturbation. Observing that p u+εh solves the system Letṗ denote the derivative of ε → p(u + εh) at ε = 0. Standard computations yield thatṗ solves the linearized reaction-diffusion system x ∈ Ω. (4.8) Furthermore, according to the chain rule, one has Let us multiply the main equation of (4.8) by q u , and integrate then two times by parts on (0, T ) × Ω. One thus gets Similarly, let us multiply the main equation of (4.6) byṗ, and integrate then by parts on (0, T ) × Ω. We obtain (4.10) By comparing (4.9) and (4.10), we infer that leading to the following duality identity: By using (4.8) and (4.6), we rewrite the expression above as Thus the desired expression of the derivative follows.
Let us now compute d 2 J T (u 0 ). Since J T is two times differentiable, one has whereq is given byq A standard reasoning enables us to prove thatq solves the linear p.d.e.
x ∈ Ω, (4.11) withṗ, the solution of the linear p.d.e. (4.7). One has By using the main equation in Systems (4.11) and (4.7), one gets The Green formula finally yields whence the expected expression for the second order derivative.
Let us now derive first and second order optimality conditions for this problem.
where q solves the adjoint system (4.6) associated to the control choice u 0 . Let u * 0 be a solution to Problem (P reduced ). Then, there exists λ ∈ [0, +∞) such that (called necessary first order optimality condition) or equivalently, the function Λ defined by vanishes identically in Ω. Moreover, one has λ Ω u * 0 (x) dx − C = 0 (slackness condition). Moreover, the second order optimality conditions for this problem read: Proof. Let us introduce the Lagrangian functional associated to Problem (P reduced ), given by According to Proposition 4.6, and denoting by d un the differential operator with respect to the variable u, the Euler inequation associated to Problem (P reduced ) reads: d un L(u, λ) · h ≥ 0 for all admissible perturbation h of u * 0 in {u 0 ∈ L ∞ (0, T ), 0 ≤ u 0 ≤ M a.e. in Ω}. This can be rewritten for all functions h as above. The analysis of such optimality condition is standard in optimal control theory (see for example [21]) and yields: Moreover, one has λ Ω u * 0 (x) dx − C = 0 (slackness condition). It remains to show that such conditions also rewrite Λ(·) = 0 in Ω. It is straightforward that if the optimality conditions above are satisfied, then Λ(·) = 0 in Ω. Let us examine the converse sense, assuming that Λ(·) = 0 in Ω. Then, for a.e. x ∈ {u * 0 = 0}, one has We infer from this result that either the pointwise or the integral constraint is saturated by every minimizer u * 0 . Corollary 4.8. Let u * 0 be a solution to Problem (P reduced ). Then, one has necessarily Proof. Let us first assume that M ≥ C/|Ω|. Let us argue by contradiction, assuming that Ω u * 0 < C. Let p (resp. q) denote the solution to the direct problem (2.18) (resp. the adjoint problem (4.6)) associated to the control choice u * 0 . According to Theorem 4.7 and its proof, the slackness condition implies that λ = 0. Recall that one has p(t, x) ∈ (0, 1) for a.e. (t, x) ∈ (0, T ) × Ω, as highlighted in Section 2.3, and therefore q(T, ·) ∈ (−1, 0) a.e. in Ω. A simple comparison argument yields that q is negative in (0, T ) × Ω (see e.g. [9]). Since G is bijective and increasing, so is G −1 and we infer that ψ is negative in Ω. By using Theorem 4.7, we get that necessarily, u * 0 (·) = M , which is in contradiction with the assumption above on M and C. The case where M < C/|Ω| is solved hereafter, in the proof of Theorem 3.1.

Optimality of constant solutions
This section is devoted to the proof the our main results, that is Therem 3.1. Let us first show (i). The proof rests upon a simple comparison argument: one shows more precisely that u M solves Problem (P reduced ) as soon as it belongs to V C,M which is equivalent to the condition above on the parameters.
Let u ∈ V C,M . Let p and p M denote the solutions to System (2.18) corresponding respectively to the control choices u and u M .
Since u belongs to V C,M and G −1 is increasing, one has G −1 (u(x)) ≤ G −1 (M ) for a.e. x ∈ Ω, meaning that p(0 + , ·) ≤ p M (0 + , ·) on Ω. According to the parabolic comparison principle, we infer that p(t, ·) ≤ p M (t, ·) on Ω, for all t ∈ [0, T ), so that one gets in particular that p * (T, ·) ≤ p M (T, ·) in Ω, and therefore, J T (u M ) ≤ J T (u). Uniqueness follows from the monotonicity of G and the comparison principle, since 0 u M a.e. in Ω.
Let us now prove (ii). Set c = C/|Ω|. According to the optimality conditions (4.12), since c < M , the function u identically equal to the constant c satisfies the first order optimality conditions if, and only if, there exists λ ∈ R + such that ψ(·) = −λ in Ω. Since (G −1 ) (u(·)) is constant in Ω, this is equivalent to say that q(0 + , ·) is constant in Ω.
First, observe that, by uniqueness of the solutions to the reaction-diffusion system (2.18), the associated solution p is constant in space. Moreover, writing p(t, ·) = p(t) with a slight abuse of notation, one easily sees that p solves the ODE (4.13) Standard uniqueness arguments coming from the Cauchy-Lipschitz theorem show that if p(0 + ) / ∈ {0, θ, 1} (the set of roots of f ), then f (p(·)) does not vanish on [0, T ] and has hence a constant sign.
Proceeding similarly for the solution q to System (4.6) associated to p = p drives us to look for constant solutions with respect to the space variable. Let q denote such a solution (whenever it exists). Hence, it solves q (t) = −f (p(t))q(t), t ∈ [0, T ], q(T ) = p(T ) − 1 and therefore, By uniqueness of the solution to (4.6), it follows that q solves (4.6). Now, if p(0 + ) = θ, meaning that u = G(θ), then p(·) = θ and one has q(t) = (θ − 1)e (T −t)f (θ) for all t ∈ [0, T ]. All in all, we get that q(0 + , ·) is constant on Ω and the switching function ψ, which is constant, reads by using that p(t) ∈ (0, 1) for all t ∈ [0, T ] and that G is bijective and increasing. We infer that the first order optimality conditions are satisfied by u.
To investigate the second order optimality conditions, it is convenient to introduce the Hilbert basis {w n } n∈N * of L 2 (Ω) made of the Neumann-Laplacian eigenfunctions defined by: α n w n with α n = h, w n L 2 (Ω) for all n ∈ N * .
Using that the solution p to (2.18) does not depend on the space variable, it is standard to expandṗ aṡ where v n solves the o.d.e. v n (t) = (−λ n + f (p(t))) v n (t) and v n (0) = ( According to Proposition 4.6, one thus computes Using that {w n } n∈N * is orthonormal in L 2 (Ω), we finally get the following diagonalized expression of the second order derivative The signature of d 2 J T (u)(h, h) seen as an infinite quadratic form with respect to h is then directly given by the sign of the coefficients δ n . Notice that for all n ∈ N * , one has Let us first assume that C ≤ |Ω|G(θ), meaning that (G −1 )(c) ≤ θ. In that case, since p solves (4.13), and that the three roots of f are 0, θ and 1, one infers that p is a decreasing function and that f (p) remains negative all along (0, T ). Furthermore, on (0, θ), the function f is positive. Finally, one computes (G −1 ) (c) = (G −1 ) (c)g (G −1 (c)) which is negative since so is g on (0, 1). Combining all these facts, we infer that and since p(T ) < 1, it follows that for every n ∈ N * . Therefore, by setting we get that for every admissible perturbation h, one has Expanding J T at the second order at u, it is then standard that this condition implies that u is a local minimizer for the functional J T .
In other words, this lemma states the existence of radially symmetric steady-states to the stationary equation associated to (2.18). We then deduce the existence of stationary subsolutions for System (2.18) that are positive and compactly supported, provided the domain contains a large enough ball, in other words that the inradius of Ω be large enough.
Using that w α is a subsolution, we deduce the following comparison result.
Let us introduce Notice that the family of subsolutions (w α ) α have already been used to provide a sufficient condition on the release function to initiate propagation of infected mosquitoes [33]. Remark 4.12. It is worth mentioning that in the one dimensional case, the expressions for R α and C α are completely explicit: We are now in position to prove Proposition 3.2 that we rewrite more precisely using the notations above.
Proposition 4.13. Let us assume (3.1). Assume moreover the existence of α ∈ (θ c , G −1 (M )] such that Ω contains strictly a ball of radius R α , and C α ≤ C. Then the constant solution u := C |Ω| is not a global minimizer of the optimization problem P reduced whenever T is large enough.
Proof. From assumption (3.1), we have G −1 (u) < θ, hence we have already seen in Section 4.2.3 that the solution, denoted p, of (2.18) with initial data G −1 (u) is constant in space and decreasing with respect to time. More precisely, it solves the ODE Hence, when t → +∞, p(t) decays to 0. For any α ∈ (θ c , G −1 (M )] satisfying the assumptions above, the subsolution w α defined in Corollary 4.10 is such that G(w α ) ∈ V C,M . From Corollary 4.11, if we take u 0 ∈ V C,M such that G −1 (u 0 ) ≥ w α , then for all t ≥ 0, the corresponding solution to (2.18) verifies p(t, ·) ≥ w α . Hence J T (u 0 ) ≤ 1 2 Ω (1 − w α (x)) 2 dx. Moreover, since p(T ) → 0 as T → +∞, we have that for T large enough Hence, u is not a global minimum of J T at time T since J T (u 0 ) < J T (u) = 1 2 Ω (1 − p(T )) 2 dx.
Remark 4.14. If G −1 (M ) > θ c , if the inradius of Ω is large enough and if C is large enough, it is always possible to find α satisfying the assumptions of Proposition 4.13. For instance, it suffices to choose α = G −1 (M ) and to take C ≥ C G −1 (M ) and the inradius of Ω large enough so that (3.1) holds and a ball of radius R G −1 (M ) be included in Ω.

Numerical experiments
In this section, we provide some numerical approximations of solutions for the optimal control problem (P reduced ).
The parameter values are given in Tables 1 and 2. We will assume that Ω is an interval (0, L), i.e. d = 1. From these tables, we deduce that s f = 0.1, δ = 4 3 , and thus θ = 36 . System (2.18) will be discretized with an explicit Euler scheme in time and a standard finite difference approximation of the Laplacian. In all simulations, the number of steps in space and time will be fixed to 20 and 200 respectively (in order to satisfy the CFL condition). The solution of the optimal control problem will be obtained by testing and combining two approaches: -a Uzawa type algorithm, based on the gradient computation of Proposition 4.6. It consists in alternating at each iteration a step of minimization of the Lagrangian associated with the problem with respect to the primal variable (u 0 ) and a step of maximization with respect to the Lagrange multiplier associated with the integral constraint. The minimization step is performed with a projected gradient type method, where L ∞ constraints on u 0 are taken into account by means of a projection operator. -the opensource optimization routine GEKKO (see [5]) solving the optimization problem using the IPOPT (Interior Point OPTimizer) library, a software package for large-scale nonlinear problems by an interiorpoint filter line-search algorithm (see [35]). This algorithm has been initialized with the previous control obtained by using the aforementioned Uzawa type algorithm.
Remark 5.1. According to Remark 2.1, the assumption (H f,g ) is satisfied for the particular choices of functions f and g given by (2.9) under (2.10) and (2.11) which hold true for the values of the parameters in Table 1. Furthermore, it is easy to check numerically that the assumption (2.12) is satisfied for the values of the parameters taken from the case at hand (see e.g. [33]). Indeed, for the parameter values, f < 0 on (0, θ), which implies that F < 0 on (0, θ), and moreover F (1) > 0 (see Fig. 1).
Let us distinguish between two cases: Case C/|Ω| > M .
In Figure 2, the local minimizers of Problem (P reduced ) for C = 1.2 and M = 0.02 (left) (resp. M = 0.03 (right)) obtained by using the aformentioned Uzawa and Gekko algorithms are reported. We observe the extinction (resp. the invasion) of the population. One recovers the theoretical result stated in item (i) of Theorem 3.1, in other words that the constant function equal to M solves Problem (P reduced ) whenever C ≥ M |Ω| (see Table 3). In this situation, the space dependency has no impact on the time dynamics, i.e. the dynamics is the same as if there is no diffusion. Then, since it is a bistable dynamics, when M < G(θ) there is extinction of the population, whereas there is invasion when M > G(θ).
This situation is illustrated in Figures 3 and 4 with Gekko algorithm and Figure 5 with Uzawa algorithm for C ∈ {0.5, 0.8} and M ∈ {0.04, 0.4}. We can see in Figure 3 that when the number total of mosquitoes released is too low (when C = 0.5), then the infected population decreased until the extinction of this population. On the contrary, if the number total of mosquitoes released is higher (when C = 0.8), then we obtain an invasion of the infected mosquitoes. The simulation with the Uzawa algorithm in Figure 5 recovers the fact that C/L is a local minimizer for Problem (P reduced ). Indeed this algorithm seems to converges always to this constant solution. Nevertheless, it is not a global minimum since Gekko provides a better control as it is illustrated thanks to the values of J T (u) reported in Table 3. Moreover, we see on Figure 3 that invasion of the infected population seems  to occurs whereas the infected population seems to go to extinction in Figure 5. This is also in concordance with the result stated in Proposition 4.13.

Perspectives
In a near future, we foresee to investigate a more involved model, closer to practical experiments, where one aims at determining release distributions in time and space, assuming that: -releases are done periodically in time (for instance every week) and are impulses in time 4 ; -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. Considering N releases, we denote by t 0 = 0 < t 1 < . . . < t N −1 < T , t i = i∆T , the release times.
As done in this article, System (1.1) can be recast without source measure terms, coming from the specific form of the control functions. In a second time, we will also look at dropping the assumption on the frequency of releases and determine optimal times of releases (in the spirit of [2], where a simpler ODE model were considered).
Another interesting question is also raised by the spatial heterogeneities. Indeed, in field experiments the environment is not homogeneous in space. Then an important issue, from an experimental point of view, is to determine how to adapt the releases with respect to the spacial heterogeneities to optimize the success of the replacement strategies.