Stochastic modeling of the Monkeypox 2022 epidemic with cross-infection hypothesis in a highly disturbed environment

1.   Introduction

Due to the major world-wide impact of vector-borne diseases on human health, many strategies, integrating more or less the three main actors of transmission (pathogen, vectors, man) were developed to reduce their spread. In this work, we are interested in investigating strategies targeting only the vector (mosquito belonging to the Aedes genus) of viral diseases such as dengue, chikungunya and zika. To this end, mathematical modelling has an important role to play in studying and conceiving different scenarios.

In this work, we focus on two strategies that have already been implemented in the field (e.g. ^[1,2]): the sterile insect technique and population replacement. We mention that there are other strategies based on genetic manipulation like, for example, the release of insects carrying a dominant lethal (RIDL) ^[3,4,5], the driving of anti-pathogen genes into natural populations ^[6,7,8], or also methods that combine both reduce and replace strategies ^[9].

The sterile insect technique consists in releasing sterilized males massively into the wild to mate with females in order to reduce the size of the insect population. It has been first studied by R. Bushland and E. Knipling and experimented successfully in the early 1950's by nearly eradicating screw-worm fly in North America. Since then, this technique has been studied on different pests and disease vectors ^[10,11]. In particular, it is interesting to seek to control mosquito populations, which has been mathematically modelled and studied in several articles, see e.g. ^{[12,13,14,15,16,17,18,19,20]}.

Recently, there has been increasing interest in the biology of Wolbachia and its application as a control agent against mosquito vector populations, taking advantage of phenomena called cytoplasmic incompatibility (CI) and pathogen interference (PI) ^[21,22]. 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, in the first mitotic divisions ^[23]. This is the so-called cytoplasmic incompatibility (CI). The pathogen interference (PI) is characterized by the inability of Aedes mosquitoes infected with some Wolbachia strains to transmit viruses like dengue, chikungunya and zika ^[24].

Once released, the mosquitoes breed with wild ones. Over time and if the releases are large and long enough, it can be expected that the majority of mosquitoes will be carriers of Wolbachia, due to cytoplasmic incompatibility. As a result of PI, the mosquito population then has reduced vectorial competence, reducing the risk of dengue, chikungunya and zika outbreaks. The technique of releasing mosquitoes carrying Wolbachia to replace the wild population is a type of population replacement strategy. It has been modelled and studied in several works, see e.g. ^{[25,26,27,28,29,30]}.

In this article, we are interested in studying the optimization of the release protocol. More precisely, given the duration of the experiment and given a certain amount of mosquitoes (which will serve to control the vector population), what should be the temporal distribution of releases in order to be as close as possible to the objective to be achieved at the final moment of the experiment? To answer this question, we first define a cost functional that will mathematically represent the objective we are trying to achieve. For the sterile insect technique, since the objective is to reduce the size of the population, the quantity to be minimized will be defined as the number of females at the final time. For the replacement of the population by Wolbachia infected mosquitoes, the quantity to be minimized will be the distance (in the least square sense) at the final moment to the infected equilibrium, corresponding to the state where all mosquitoes carry the bacteria Wolbachia (the entire population is infected). Obviously, we include as a constraint in the optimal control problem that the number of mosquitoes released during the experiment is limited.

Similar optimization problems for sterile insect or population replacement techniques, with different cost functionals, have been proposed in e.g. ^[31,32,33]. The main difference here with previous works is that we only consider the state at the final time, which seems natural but induces several additional technical difficulties.

The main lines of the document are as follows. In the following section, we present mathematical models for each of the two strategies. Starting from a model that integrates the entire mosquito life cycle, we use several hypotheses to simplify this system and arrive at two simple systems that model the two techniques studied in this article. In section 3, we present the cost functions and describe the optimization problems to be solved. We state existence results and some optimal control properties. These results are illustrated in the section 4 where a few numerical simulations are provided. Next, we present a conclusion and a discussion of our results. Finally, an appendix is devoted to some technical results, in particular those concerning to the existence of optimal control.

2.   Mathematical modelling

2.1.   Mosquito life cycle

The life cycle of a mosquito (male or female) occurs successively in two distinct environments: it includes an aquatic phase (egg, larva, pupa) and an non-aquatic phase (adult). A few days after mating, a female mosquito can lay a few dozen eggs, possibly spread over several breeding sites. Once laid, eggs of some species can withstand hostile environments (including adverse weather conditions) up to several months, before hatching. This characteristic contributes to the adaptability of mosquitoes and has allowed them to colonize temperate regions. After stimulation (e.g. rainfall), the eggs hatch to give birth to larvae that develop in the water and reach the pupa state. This larval phase can last from a few days to a few weeks. Then, the insect undergoes its metamorphosis. The pupa (also called nymph) remains in the aquatic state for 1 to 3 days and then becomes an adult mosquito (or imago): it is the emergence and the beginning of the air phase. The lifespan of an adult mosquito is estimated to be of a few weeks.

In many species, egg laying is only possible after a blood meal, i.e. the female must bite a vertebrate before each egg laying. This behaviour, called hematophagy, can be exploited by infectious agents (bacteria, viruses or parasites) to spread, alternately from a vertebrate host (humans, for what we are interested in here) to an arthropod host (here, the mosquito).

In order to model this life cycle dynamics, we introduce the following quantities:

● $E(t)$ density of (viable) eggs at time $t$;

● $L(t)$ larvae density at time $t$;

● $P(t)$ pupa density at time $t$;

● $F(t)$ and $M(t)$ density of adult females and males, respectively, at time $t$.

We consider the parameters:

● $\beta_E > 0$ is oviposition rate for females;

● $\delta_E, \delta_L, \delta_P, \delta_F, \delta_M > 0$ are death rates for eggs, larvae, pupa, adult females, and males, respectively;

● $\tau_E$ hatching rate for eggs;

● $\nu$ the probability that a pupa gives rise to a female, therefore $(1-\nu)$ is the probability to give rise to a male ($0 < \nu < 1$);

● $\tau_L$ and $\tau_P > 0$ transition rates from larval phase to pupa and from pupa to adult;

● intra-specific competition is only expected to occur in the aquatic phase. This models the occupation of breeding sites that can only support a limited number of eggs and the limited access to resources for larvae. In the larval compartment, this competition is described by the introduction of a positive constant denoted $c$ and is assumed to depend on the concentration of larvae: the larger the number of larvae, the greater the competition to find the nutrients essential for larval development. The environmental capacity for eggs is denoted $K$. This amount can be interpreted as the maximum density of eggs that females can lay in breeding sites.

From the above considerations, we can determine the dynamics of the mosquito population and obtain the following dynamical system

It is important to note that this system is only valid if a fairly large number of mosquitoes are considered, since in this model, the probability of a female mating with a male is equal to 1.

Such an assumption seems relevant for high number of mosquitoes and is done in several models ^[34]. In more generality, one may consider that the rate $\beta_E$ depends on $M$ as a function $\beta_E(M)$, which complexifies the study performed here. However, it is important to take this dependence into account when the size of the population is significantly reduced. For example, we refer to ^[19,41] and the references therein where the functions $\beta_E(M)$ are chosen in exponential form: it has a linear growth compared to $M$ for a small value of $M$, and on the contrary, it is almost constant for a large value of $M$.

In order to further simplify this system of ODE, we assume that the time dynamics of the pupa compartment is fast. Then, denoting $\tilde{t} = \varepsilon t$ a new time variable, and $\tilde{P}(\tilde{t}) = P(t)$, we have

As $\varepsilon \to 0$, we deduce that we may replace the third equation in system (2.1) by

which implies the relation $P = \frac{\tau_L}{\delta_P+\tau_P} L$.

To reduce further this system of equations, we will use some assumptions on the larval compartment. We first consider that the competition at the larvae stage is negligible (i.e. $c\ll 1$). Moreover, in favourable conditions, the larval stage may be really fast. Then, by the same token as above, this compartment may be considered at equilibrium leading to the relation

Injecting this relation, system (2.1) reduces to

where we use the notation $\beta_F = \frac{\tau_P\tau_L\tau_E}{(\delta_P+\tau_P)(\tau_L+\delta_L)}$.

2.2.   Sterile insect technique

As explained above, the sterile insect technique consists in releasing sterile males to mate with females with the aim of reducing the size of the population. We denote by $M_s$ the density of sterile males. Only females mating fertile males will be able to lay viable eggs. Assuming an uniform repartition of the population of mosquitoes, the probability that a female mates with a fertile male is given by $\frac{M}{M+\gamma M_s}$. The parameter $\gamma$ accounts for the fact that females may have a preference for fertile males. Introducing the sterile male population into system (2.2) leads to

It is clear that the extinction state, where $E = F = M = M_s = 0$ is a steady state. However, an important observation for this system is that this steady state cannot be reached. Indeed, under suitable assumptions on the parameters, it is unstable as stated in the following proposition (whose proof is given in the Appendix).

Proposition 2.1. Assume that

Then the steady state $(0, 0, 0, 0)$ for system (2.3) is unstable.

The first hypothesis of (2.4) reflects the fact that the mortality rate of released sterile mosquitoes is higher than that of wild mosquitoes ^[19]. The second hypothesis implies that the egg laying rate is sufficiently high. It should be noted that with the values taken in the field, these assumptions are satisfied (see Section 4).

Due to the high number of equations in system (2.3), we will reduce this system by making the following assumption : The death rate for males and females is the same ($\delta_F = \delta_M$) and the probability that a pupa emerges to a female or a male is the same ($\nu = \frac 12$). Thanks to this assumption, male and female densities satisfy the same equation. Hence assuming that initially these quantities are equal, we will have that $F = M$.

In this case (which will be the setting of this paper), system (2.3) reduces to

2.3.   Introduction of the bacteria Wolbachia

To model the strategy consisting of releasing Wolbachia infected mosquitoes to replace the wild population, we introduce the infected population into (2.2). Let us denote $E_i$, $F_i$, $M_i$ the infected by Wolbachia eggs, female and male compartments.$E_u$, $F_u$, $M_u$ correspond to the uninfected compartments.

Assuming an uniform repartition of the population of mosquitoes, then the probability for a female to mate with an infected male is equal to the proportion of infected males into the population, i.e. $\frac{M_i}{M_u+M_i}$. Similarly, the probability to mate with a uninfected male is $\frac{M_u}{M_u+M_i}$. To model the cytoplasmic incompatibility, we introduce a parameter denoted $s_h$, corresponding to the fraction of uninfected female's eggs fertilized by infected males which will not hatch (and will thus not be counted in the egg compartments which represent only viable eggs). We have $0 < s_h\leq 1$, the case $s_h = 1$ correspond to the perfect cytoplasmic incompatibility (see e.g. ^[35] for a discussion on cytoplasmic incompatibility). From system (2.2), we construct the following system taking into account infected and uninfected mosquitoes:

In this system, we have introduced the following two parameters: $\eta < 1$ modelling the fecundity reduction of infected females with respect to uninfected females, $\delta > 1$ modelling the increase of mortality for infected mosquitoes.

As above, for the sterile insect technique, we make use of the same set of assumptions ($\nu = \frac 12$ and $\delta_F = \delta_M$) to reduce the system by considering that the quantity of males and females is the same : $M_u = F_u$ and $M_i = F_i$.

Under these assumptions, system (2.6) for the Wolbachia strategy reduces to

3.   Towards optimization problems

We introduce in this section the optimization problems considered in this work. Since for both strategies, the idea consists in releasing mosquitoes (sterile males or infected by Wolbachia), the control is on the release function which will be denoted $u$. We assume that the release occurs in a time interval $[0, T]$ for $T > 0$ given. Obviously some constraints should be satisfied by the release function. We assume that there exists $C\geq 0$ and $\overline{U}\geq 0$ such that $0\leq u\leq \overline{U}$ a.e. and $\int_0^T u(t)dt \leq C$. The first bound means that $u$, the instantaneous rate of mosquito release (number of mosquitoes per unit of time) is bounded by a constant $\overline{U}$ all along the period $[0, T]$; the second means that the total number of released mosquitoes is bounded by a positive constant $C$. Both assumptions are natural considering that one cannot produce an infinite number of mosquitoes to release nor release them at an infinite rate.

Before the experiments begin, it is assumed that the systems are in equilibrium. That is why, for each system, we first determine the equilibria. We present successively the optimization problem considered for the sterile insect technique and for the population replacement by shape{Wolbachia} infected mosquitoes.

3.1.   Sterile insect technique

Let us first consider the system (2.5) for the sterile insect technique. The following lemma gives the equilibria.

Lemma 3.1. Under the assumption (2.4), there are two equilibria for system (2.5): the extinction equilibria given by $(E_1^*, F_1^*, M_s^*) = (0, 0, 0)$, and the non-extinction equilibria $(E_2^*, F_2^*, M_s^*) = \left(\bar E_2, \frac{\nu \beta_F}{\delta_F}\bar E_2, 0\right)$ with $\bar E_2 = K\left(1-\frac{(\tau_E+\delta_E)\delta_F}{\nu \beta_E \beta_F}\right)$. Moreover, the non-extinction equilibrium is linearly asymptotically stable.

Let us denote $u$ the release function of sterile male mosquitoes. Then system (2.5) reads

In our minimisation problem, we want to find the release function $u$ under the above mentioned physical constraints for which the solution of the final time is the closest possible to the extinction equilibrium. More precisely, let us introduce the cost functional

We want to solve the problem

We insist that in this problem of optimal control, we want to minimize the number of eggs as well as the number of females. This is due to the complex life cycle of mosquitoes during which eggs can stay long before they hatch. In other words, even if the number of females is greatly reduced, there may be a resurgence of mosquitoes after a certain time if the egg stock is not negligible. Similarly, if we only minimize the number of eggs, females can lay a large number of eggs. It is therefore important to minimize both the number of eggs and the number of females.

The following result, whose proof is postponed to the Appendix, gives the existence of a solution to this problem.

Proposition 3.2. Under the assumption (2.4), problem (3.2) has a solution $u^*$. Moreover, assuming that

the optimal control strategy uses the maximal amount of mosquitoes, in other words

and there exists $T_0\in (0, T)$ such that $u^* = 0$ on $(T_0, T)$.

3.2.   Population replacement

Let us consider the reduced model (2.7) for the introduction of the bacteria Wolbachia. The following Lemma gives the equilibria for this system:

Lemma 3.3. Let us consider that $1 < \delta$, $\eta < 1$, $0 < s_h\leq 1$. We denote $b = \frac{\nu\beta_F\beta_E}{(\tau_E+\delta_E)}$. Assume moreover that

Then there are four distinct non-negative equilibria:

● Wolbachia invasion $(E^*_{uW}, F^*_{uW}, E^*_{iW}, F^*_{iW}) : = \left(0, 0, K\Big(1-\frac{\delta\delta_F}{b\eta}\Big), K\Big(\frac{\nu\beta_F}{\delta\delta_F}-\frac{\nu\beta_F}{b\eta}\Big)\right)$ is stable;

● Wolbachia extinction $(E^*_{uE}, F^*_{uE}, E^*_{iE}, F^*_{iE}) : = \left(K\Big(1-\frac{\delta_F}{b}\Big), K\Big(\frac{\nu\beta_F}{\delta_F}-\frac{\nu\beta_F}{b}\Big), 0, 0\right)$ is stable;

● co-existence steady state $(E^*_{uC}, F^*_{uC}, E^*_{iC}, F^*_{iC})$ is unstable, with $F^*_{uC} = \frac{\nu\beta_F}{\delta_F}E^*_{uC}$, $F^*_{iC} = \frac{\nu\beta_F}{\delta\delta_F}E^*_{iC}$ and

● extinction $(0, 0, 0, 0)$ is unstable.

Notice that the first assumption in (3.4) boils down to consider that the birth rate is larger than the death rate and is generically satisfied for mosquito populations. Since $s_h$ is expected to be close to $1$ (the case $s_h = 1$ being the perfect cytoplasmic incompatibility case), the second inequality may be seen as a condition on $K$ to be large enough.

As above, we denote $u$ the release function of Wolbachia-infected mosquitoes. Assume that the system is initially at the Wolbachia free equilibrium, we want to determine an optimal release function $u$ which brings the system as close as possible to the Wolbachia invasion equilibrium. More precisely, let us consider $(E_u, F_u, E_i, F_i)$ solution to the following Cauchy problem:

We introduce the following cost function

with the standard notation for the positive part $X_+ = \max \{X, 0\}$ for $X\in \mathbb{R}$. The cost functional considered in this study models that one wants to minimize the distance (in a least square sense) of the solution at final time $T$, to the steady state $(E^*_{uW}, F^*_{uW}, E^*_{iW}, F^*_{iW})$ corresponding to Wolbachia invasion (with the notation in Lemma 3.3). We investigate the following optimization problem

The following result gives the existence of a solution. Its proof is given in the Appendix.

Proposition 3.4. Under the same assumptions as in Lemma 3.3, problem (3.6) has a solution.

System (2.7) may even be simplified further by assuming a fast dynamics of the aquatic phase and a large fertility, in the spirit of ^[36]. This leads to a simple differential equation on the proportion of infected female mosquitoes, for which the optimization problem has been studied in detail in ^[37]. In particular, the optimality conditions have been derived by using the so-called Pontryagin Maximum Principle (PMP, see Lemmas 4, 5, 6). It has been proved for this simplified system that the optimal strategy uses the maximal amount of mosquitoes, in other words that $\int_0^T u^*(t)\, dt = C$ whenever $\overline{U}T > C$. It is likely that the same property holds true when considering the more realistic system (3.5) even if it should be more tedious to show it. This can be observed in the simulations shown in the next section.

4.   Numerical simulations

We will now give some solutions of the optimal control problems (3.2) and (3.6). For this purpose, we will use the opensource optimization routine GEKKO (see ^[38]). Differential algebraic equations are implicitly solved by using orthogonal collocation over finite elements. The optimization problem is solved by using the IPOPT library which implements a primal-dual interior point method ^[39]. Following experimental results reported in the field (see e.g. ^[1,2]), where releases are made from several weeks to three months, we choose $T = 80$ days to perform our numerical simulations.

4.1.   Sterile insect technique

In this section, we will give some illustrations of the optimal strategy given by the optimal control problem (3.2). We will use the parameter values of Table 1 coming from [19,Table 1-3]. We recall that $\nu$ is assumed to be equal to $0.5$.

As in ^[19], in order to get results relevant for an island of $74\, ha$ (hectares) with an estimated male population of about $69\, ha^{-1}$, the total number of males is equal to $M^* = F^* = 69\times 74 = 5106.$ If we assume that $M_s = 0$, we get

Thus, the value of $K$ is given by the expression

Figures 1 and 2 show some numerical simulations for the sterile insect technique. For each figure, we display the dynamics of eggs $E$ (left plot), females $F$ (middle plot), and optimal release strategy $u$ (right plot). In Figure 1, we provide some optimal strategies for the problem (3.2) with $T = 80$, $\gamma = 1$, a total quantity of sterile mosquitoes $C = 150000$ and different values for the maximum instantaneous release allowed $\overline{U}$ : $\overline{U} = 5000$ (1st line), $\overline{U} = 10000$ (2nd line), and $\overline{U} = 20000$ (3rd line). In Figure 2, we focus on the influence of the mating competitiveness parameter $\gamma$. The choice $\gamma = 1$ models that females have equal probability to mate with fertile or sterile males. Reducing $\gamma$ allows us to model the fact that the probability that females mate with sterile males is smaller than that of females mating with fertile males. In the numerical simulations, the parameter values are $\gamma = \frac 13$ (1st line), $\gamma = \frac 23$ (2nd line), and $\gamma = 1$ (3rd line).

It can first be observed that in both situations, as provided for in Proposition 3.2, the optimal control strategy uses the maximal amount of mosquitoes ($\int_0^T u(t)\, dt = C$) and does not act at the end of the time interval. In addition, it seems preferable to concentrate most releases towards the end of the time interval.

We notice in Figure 1 that the region where the constraint $u = \bar{U}$ is saturated seems to disappear as $\bar{U}$ becomes larger. Moreover, it seems that the remaining number of eggs and females is not strongly impacted by the value of $\bar{U}$. We observe in Figure 2 that the parameter $\gamma$ plays a role in the optimal strategy. In particular, it seems that the optimal strategy is to act on a wider time interval for larger $\gamma$. Moreover, the final number of eggs and females is smaller for larger $\gamma$. Indeed, it seems natural that by increasing the competitiveness of infertile men, we increase the effectiveness of the strategy.

4.2.   Population replacement

The section is devoted to illustrate the optimal strategy produced by the optimal control problem (3.6) with system (3.5) as constraint. The parameters values for system (3.5) are given in Tables 1 and 2. The expression of $b$ is given in Lemma 3.3. The value of the cytoplasmic incompatibility parameter $s_h$ (corresponding to the fraction of eggs from uninfected females fertilize by infected males which will not hatch) comes from ^[40]. The fecundity reduction $\eta$ of infected females with respect to uninfected females and the increase of mortality $\delta$ for infected mosquitoes have been fixed following ^[28].

We assume that Wolbachia has not been introduced before, and thus, that system (3.5) was initially at the equilibrium $(E_{uE}^*, F_{uE}^*, 0, 0)$. As for the sterile insect strategy, the initial density of mosquitoes will be equal to $F_{uE}^* = 5106$. We can deduce the value of $K$ thanks to the expressions of the equilibrium $(E^*_{uE}, F^*_{uE}, E^*_{iE}, F^*_{iE})$ in Lemma 3.3. The numerical results are displayed in Figure 3, when we take a total amount of mosquitoes $C = 150000$, and in Figure 4 for a total amount of mosquitoes $C = 1000$. For each figure, we display in dashed lines the position of the coexistence equilibria. We know in particular that once infected quantities and uninfected quantities have reached values which are respectively above and below the values of the coexistence steady state, then they are in the basin of attraction of the invasion steady state. We first notice that, as in the case of the sterile insect technique, there exists a time after which the control function $u$ vanishes.

We can also draw some conclusions by comparing the Figures 3 and 4. Indeed, we observe that when the amount of mosquitoes is large enough, it is better to act at the beginning of the process. However, if the amount of mosquitoes is low, it seems preferable to release them later. To interpret this observation, we recall that on the one hand the dynamic system (3.5) with $u = 0$ has two stable steady states ("Wolbachia invasion" and "Wolbachia extinction"), on the other hand the cost function that we want to minimize is the distance to the stable steady state "Wolbachia invasion" at the final time $T$. Then there are two possibilities:

● Either the quantity of mosquitoes to be released is large enough to reach the basin of attraction of the Wolbachia invasion in steady state. In this case, it is better to go as quickly as possible to this basin of attraction. Indeed, once the dynamical system solution enters this attraction basin, it converges to the desired steady state, even if $u = 0$. This is the situation depicted in the Figure 3.

● Or the amount of mosquitoes to be released is too small to reach the basin of attraction of the Wolbachia invasion in steady state. In this situation, once the releases stopped, the solution of the dynamical system (3.5) with $u = 0$ moves away from the desired stable state. It therefore seems preferable to release them as late as possible. This is observed in the Figure 4.

This observation should be related to the threshold phenomenon which has been observed in ^[37], where the authors have approximated the model by a scalar equation whose main unknown is the proportion of infected adult mosquitoes. For this very simple system, the authors proved that the control is bang-bang and that there is a threshold on the total amount of mosquitoes below which it is preferable to act at the end of the time interval and above which the action occurs at the beginning of the interval. For the model (3.5) with six equations, we remark, in the figures 3 and 4, that the optimal strategy is more complex. However, the threshold phenomenon seems still to be true.

The probability of cytoplasmic incompatibility $s_h$ plays a key role in the Wolbachia strategy. In Figure 5, we compare several values for $s_h$. For a given amount of mosquitoes infected with Wolbachia disease to be released, if the probability of cytoplasmic incompatibility is too low, we do not enter in the attraction basin. For a large values for $s_h$ (resp. for a small value), we recover the structure of the optimal release observed in Figure 3 (resp. Figure 4). In addition, we notice that if we take a perfect cytoplasmic incompatibility ($s_h = 1$) instead of the value given in ^[40], we have a similar behaviour of the solution and optimal release.

5.   Discussion and conclusion

In this article, we study the problem of optimizing a dissemination protocol for a population replacement strategy and for the sterile insect technique applied to the control of the genus Aedes mosquito population. In our approach, we are looking for a control function $u$ minimizing the distance to the desired equilibrium (replacement or extinction of the wild population) at the final treatment stage. In particular, we show the existence of such optimal control and, after establishing some properties, we have illustrated them with numerical simulations.

As mentioned in the introduction, our choice of cost functional differs from other cost functionals that can be found in the literature, since we have chosen to minimize the distance from the desired equilibrium at the final time while other works consider the integral over the whole time of release ^[31,32]. As a consequence, the calculated release functions differ from those obtained in these other works. But some similarities can be observed. Indeed, for the Wolbachia strategy with both cost functionals, it seems preferable to act mainly at the beginning of the time interval when the number of available mosquitoes to release allows us to achieve the desired equilibrium.

We discuss now some limitations of our models which will be addressed in a future work ^[41]. First, in both situations, we use a cost functional that measures the distance, in the least square sense, from the desired equilibrium (extinction $F = E = 0$ for the sterile insect technique, and the Wolbachia invasion equilibrium). Thanks to this choice, we are able to provide the time distribution of the release function in order to be as close as possible to the desired steady state. This approach is entirely justified when one considers a given number of mosquitoes to be released during a given period of time and we want to optimize the release protocol. However, we may be interested in a different approach; for example, we may want to minimize the number of mosquitoes to use, since the production of such mosquitoes can be costly, financially speaking. Or we may be interested in reaching the desired state of equilibrium at the end of the release period. Indeed, with the cost functionals considered in our work, we cannot guarantee the success of the strategy. To answer this interesting question, we must first determine the basin of attraction for the desired steady state and then use a different cost functional for which the study in this article should be adapted.

Second, for the mathematical modelling, we have made several assumptions in the aim to simplify the system in order to derive models that could be tractable for a mathematical study. It is likely that some of these assumptions could be weakened. In particular, as already mentioned, the pertinence of system (2.1) for a small population is not so clear. Since the population of mosquitoes is usually high, the use of such models is often justified. However, when we aim at eradicating this population by using the sterile insect technique, the behaviour of the system close to the extinction steady state plays an important role. Therefore, a model that more accurately describes the dynamics near the extinction steady state is needed. To this end, one strategy is to use a birth rate of $\beta_E$ depending on the male density.

For instance, in ^[19], a function $\beta_E$ depending exponentially of the male density has been considered, taking into account an Allee effect which guarantees the stability of the extinction steady state.

In addition, it has been observed that the mortality rate of males may be higher than that of females. Thus, assuming that the number of male mosquitoes is the same as the number of female mosquitoes is a very strong hypothesis that should be weakened.

Nevertheless, the current work, and the rigorous mathematical results that we have been able to demonstrate in this simplified framework, should be a useful step towards a future understanding of more general and realistic models.

Acknowledgements

The authors were partially supported by the Project "Analysis and simulation of optimal shapes - application to life sciences" of the Paris City Hall.

Conflict of interest

The authors declare that they have no conflict of interest.

A.   Proof of Proposition 2.1

Thanks to Assumption (2.4), there exists $\varepsilon > 0$ small enough such that

Using that $M_s(t) = e^{-\delta_st}M_{s0}$, $M(t)\geqslant e^{-\delta_Mt}M_{0}$ for $t\geqslant0$ and (2.4), we deduce that there exists $t^* > 0$ such that

Let us assume by contradiction that the extinction steady state is stable. Then, we place ourselves in a neighbourhood of this equilibrium in which system (2.3) is monotonous. We deduce from standard comparison principle for monotonous system that, for all $t\geqslant t^*$

where, $(E_1, F_1, M_1)$ solves the following system, for $t\geqslant t^*$,

complemented with initial data $(E_1, F_1, M_1)(t^*) = (E, F, M)(t^*)$. We may study the stability of the extinction steady state for this later system. The Jacobian of this system at $0$ is given by

The characteristic polynomial for this matrix is given by

We have that $P_0(x)\to -\infty$ as $x\to +\infty$, and under assumption (A.1), $P_0(0) > 0$, then $P_0$ admits a positive root. Hence, Jac$(0)$ admits a positive eigenvalue, and the extinction state for system (A.3) is unstable. From the comparison in (A.2), we conclude the proof.

B.   Proof of Proposition 3.2

For the analysis of the optimal control problem (3.2), it will be useful to notice that the solutions of System (3.1) remain bounded.

Lemma B.1. Let $u\in {\mathcal U}_{C, \overline{U}}$ and $(E, F, M_s)$ be the solution of System (3.1) associated to the control function choice $u$. For every $t\in [0, T]$, one has

Proof. Notice as a preliminary remark that a standard barrier argument ensures the positiveness of solutions to System (3.1). Let us show the right inequality on $E(\cdot)$. One has $E(0) < K$. Assume by contradiction the existence of $t_0\in (0, T]$ such that $E(t_0) = K$. Without loss of generality, we assume that $t_0$ is the first solution of the equation $E(t) = K$ on $(0, T]$. Since $\frac{dE}{dt}(t_0) = -(\tau_E+\delta_E)K < 0$, we infer that $E(t) > K$ for $t < t_0$, close enough to $t_0$ whence a contradiction.

The left-hand side inequality on the function $E(\cdot)$ follows directly from the observation that the right-hand side of the first equation of System (3.1) is bounded by below by $-(\tau_E+\delta_E)E$ and a Gronwall argument.

Regarding now the inequalities on $F$, we claim that the left inequality follows from the positiveness of $E$. Moreover, by using that $E(\cdot) < K$ and the expression of $b$, we get

Existence of an optimal control.  Let us consider a minimizing sequence $(u_n)_{n\in \mathbb{N}}$ and denote by $(E_n, F_n, {M_s}_n)_{n\in \mathbb{N}}$ the corresponding solution to System (3.1). Noting that the class ${\mathcal U}_{C, \overline{U}}$ of admissible controls is compact for the $L^\infty$ weak-star topology, we infer the existence of $u^*\in {\mathcal U}_{C, \overline{U}}$ such that $(u_n)_{n\in \mathbb{N}}$ converges up to a subsequence to $u^*$ for the $L^\infty$ weak-star topology. Since

the sequence $({M_s}_n)_{n\in \mathbb{N}}$ converges in $H^1(0, T)$ up to a subsequence to $M_s^*$ given by

According to Lemma B.1, the triple $(E_n, F_n, {M_s}_n)$ is uniformly bounded on $[0, T]$. By using this boundedness property, one easily gets that $(\frac{dE_n}{dt})_{n\in \mathbb{N}}$ and $(\frac{dF_n}{dt})_{n\in \mathbb{N}}$ are bounded in $C^0([0, T])$ and therefore, $(E_n, F_n)_{n\in \mathbb{N}}$ is bounded in $W^{1, \infty}(0, T)$. According to the Ascoli theorem, the sequence $(E_n, F_n)_{n\in \mathbb{N}}$ converges to some $(E^*, F^*)\in W^{1, \infty}(0, T)$ in $C^0([0, T])$. As a consequence, according to (3.1), $(\frac{dE_n}{dt}, \frac{dF_n}{dt}, \frac{d{M_s}_n}{dt})_{n\in \mathbb{N}}$ is bounded in $L^2(0, T)$ and therefore, $(E_n, F_n, {M_s}_n)_{n\in \mathbb{N}}$ also converges (up to a subsequence) to $(E^*, F^*, M_s^*)$ in $H^1(0, T)$. We then infer from all the considerations above that $(E^*, F^*, M_s^*)$ satisfies System (3.1) and that $(J(u_n))_{n\in \mathbb{N}}$ converges, up to a subsequence, to $J(u^*)$. The existence follows.

First order optimality conditions. Let $u^*$ be an optimal control for Problem (3.2) and $(E^*, F^*, M_s^*)$ be the corresponding trajectories, solutions of (3.1) for $u = u^*$. To write the first order optimality conditions, we will use the Pontryagin Maximum Principle (PMP). To take into account the integral constraint on $u$, it is convenient to introduce a new state variable $y$ solving the o.d.e.

in such a way that the constraint $\int_0^T u(t)\, dt\leq C$ rewrites as the terminal condition $y(T)\leq C$. Under this form, it is more standard to write the Hamiltonian $\mathcal{H}$. Hence, it suffices to add the term $\lambda u$ in the definition of $\mathcal{H}$, where $\lambda$ is the Lagrange multiplier associated to the constraint $\int_0^T u(t)\, dt \leq C$.

Let us introduce the function $f_E$ defined by

as well as the Hamiltonian of Problem (3.2), given by

According to the Maximum Principle (see, e.g. ^[42]), there exists an absolutely continuous mapping $p:[0, T]\to \mathbb{R}^3$ called adjoint vector such that the so-called extremal $((E^*, F^*, M_s^*, y^*), (p_1^*, p_2^*, p_3^*, \lambda^*), u^*)$ satisfies a.e. in [0, T]:

● Adjoint equations:

and in addition, ${\lambda^*}' = 0$ which implies that $\lambda^*$ is a constant (still denoted $\lambda^*$ with a slight abuse of notation).

● Minimality condition:

and therefore, one has

● Transversality conditions: we impose the terminal conditions

on the adjoint state, where $\xi\in \mathbb{R}_+$ satisfies moreover the complementary condition $\xi (y(T)-C) = 0$.

We infer from (B.3) that $\lambda^*\geq 0$.

The total number of mosquitoes is used. Let us start with a preliminary lemma, whose proof is postponed at the end of this section.

Lemma B.2. Let us assume that $E(0) < K$. Then, the solution $(E, F, M_s)$ of System (3.1) satisfies

Let us argue by contradiction, considering $u^*$ a solution of Problem (3.2) and $(E^*, F^*, M_s^*)$ the associated trajectory. If $\int_0^T u^*(t)\, dt < C$, then one has necessarily $\xi = 0$ or equivalently $\lambda^* = 0$. We will reach a contradiction by showing that one has $p_3^* < 0$ on $(0, T)$. Indeed, if $p_3^* < 0$ on $(0, T)$, then one has necessarily $u^* = \overline{U}$ on $(0, T)$ according to (B.2) since $\lambda^* = 0$. But $\overline{U}$ is not feasible according to condition (3.3), yielding a contradiction.

Let us show that $p_3^* < 0$ on $(0, T)$. To this aim, we introduce

Then, the first two equations of the adjoint system (B.1) read

Let $v:t\mapsto p_1^*(t)e^{\int^t A}$. The last system becomes

Therefore, $v$ solves the ODE

Notice that

Introduce the change of variable $s = \int_0^t e^{\int^z(\delta_F+A)}\, dz$, $\tilde T = \int_0^T e^{-\int^z(\delta_F+A)}\, dz$ and the function $w$ defined on $[0, \tilde T]$ by $w(s) = v(t)$. Then, the function $w$ satisfies the Cauchy system

where $t$ has to be understood as a function of $s$ in this system. We then infer that $w"(T) > 0$ and therefore $w$ is convex in a neighbourhood of $\tilde T$. Since it is also positive and decreasing according to the terminal conditions, it follows that $w$ cannot vanish on $[0, \tilde T]$. We successively infer that $v$ is positive on $[0, T]$ and so is $p_1^*$.

Recall that $p_3^*$ satisfies the equation

and therefore

according to the reasoning above and Lemma B.2. It follows that $t\mapsto p_3^*(t)e^{-\delta_s t}$ increases on $[0, T]$ and vanishes only at $T$. Thus, $p_3^* < 0$ on $[0, T)$ and we are done.

Structure of the control. We have shown that $\lambda^*\neq0$ and therefore, $\lambda^* < 0$. According to the first order optimality conditions (and (B.2) in particular), since $p_3^*(T) = 0$ and $p_3^*$ is continuous, we infer that $u^* = 0$ in a neighbourhood of $T$.

Proof of Lemma B.2. Using Lemma B.1, $E(t) < K$ for all $t\in[0, T]$. After some computations, we thus obtain

C.   Proof of Proposition 3.4

This proof is very similar to the one of Proposition 3.2. For the sake of completeness but to avoid redundancies, we only provide a sketch of proof. Let us consider a minimizing sequence $(u_n)_{n\in \mathbb{N}}$ and denote by $(E_u^n, F_u^n, {E_i}^n, F_i^n)_{n\in \mathbb{N}}$ the corresponding solution to System (3.5).

● Since $(u_n)_{n\in \mathbb{N}}$ is uniformly bounded in $L^\infty(0, T)$, it converges to some element $u^*\in {\mathcal U}_{C, \overline{U}}$.

● The 4-tuple $(E_u^n, F_u^n, E_i^n, F_i^n)_{n\in \mathbb{N}}$ is bounded in $H^1(0, T)$. First, observe that a standard barrier argument shows that each element of this 4-tuple is positive. Moreover, given $n\in \mathbb{N}$, one has $E_i^n(t)+E_u^n(t) < K$. Indeed, one has $E_i^n(0)+E_u^n(0)\leq K$. Assuming by contradiction that $\max_{t\in [0, T]}E_i^n(t)+E_u^n(t) > K$, let $t_0$ be the first time in $(0, T)$ such that $E_i^n(t_0)+E_u^n(t_0) = K$. Then, one computes

yielding a contradiction. It follows that $(E_u^n)_{n\in \mathbb{N}}$ and $(E_i^n)_{n\in \mathbb{N}}$ are bounded in $C^0([0, T])$. Since

and

it follows that $(F_u^n)_{n\in \mathbb{N}}$ and $(F_i^n)_{n\in \mathbb{N}}$ are also bounded in $C^0([0, T])$. Since $\frac{d}{dt} F_u^n \geq - \delta_{F} F_u^n$, a Gronwall inequality yields $F_u^n(t)\geq K\left(\frac{\nu\beta_F}{\delta_F}-\frac{\nu\beta_F}{b}\right)e^{-\delta_F t}$. Therefore, $F_u^n+F_i^n$ cannot vanish on $[0, T]$ and we finally get the expected conclusion.

● By boundedness of $(E_u^n, F_u^n, {E_i}^n, F_i^n)_{n\in \mathbb{N}}$ in $(H^1(0, T))^4$, there exists $(E_u^*, F_u^*, E_i^*, F_i^*)\in (H^1(0, T))^4$ such that $(E_u^n, F_u^n, {E_i}^n, F_i^n)_{n\in \mathbb{N}}$ converges up to a subsequence to $(E_u^*, F_u^*, E_i^*, F_i^*)\in (H^1(0, T))^4$, weakly in $H^1(0, T)$ and strongly in $L^2(0, T)$. Standard variational arguments show not only that $(E_u^*, F_u^*, E_i^*, F_i^*)$ satisfies (3.5) associated to the control function $u^*$, but also that $(J(u_n))_{n\in \mathbb{N}}$ converges to $J(u^*)$. The existence follows.