Contrôle d’épidémies par lâchers de moustiques

Piste rouge Le 26 octobre 2018  - Ecrit par  Luis Almeida, Yannick Privat, Martin Strugarek, Nicolas Vauchelet Voir les commentaires

Le moustique est l’un des animaux les plus dangereux pour l’être humain. Il existe une technique surprenante permettant de limiter sa prolifération dans certaines régions du globe, consistant à organiser des lâchers de moustiques infectés par la bactérie Wolbachia. Nous expliquons comment élaborer une stratégie efficace en utilisant les outils issus de la théorie du contrôle et de l’optimisation.

Remarque : la piste devient carrément noire au niveau des équations (NDLR).

Introduction

L’animal le plus dangereux pour l’être humain est, de loin, le moustique. On estime, en effet, que les maladies transmises par les moustiques sont responsables de plus de 700 000 décès par an parmi la population humaine (voir par exemple les statistiques compilées sur le blog de Bill Gates [Gates]).
Il existe près de 3000 espèces de moustiques dont la plupart n’ont qu’un rôle pollinisateur et ne s’intéressent pas à l’être humain.
Cependant plusieurs centaines d’espèces ont besoin du sang de mammifères (dont l’être humain) dans leur développement. En effet, les femelles moustiques utilisent les protéines et le fer contenus dans le sang pour déclencher le fonctionnement de leurs ovaires.
En se nourrissant, elles peuvent transmettre à l’être humain des agents pathogènes : on dit que le moustique est un vecteur de nombreuses maladies. Par exemple, des moustiques anophèles sont des vecteurs de la malaria, des moustiques du genre Aedes sont des vecteurs de la dengue, du chikungunya, du zika, de la fièvre jaune, d’autres moustiques du genre Culex transmettent la fièvre du Nil, ...
De plus, avec le réchauffement climatique et l’activité humaine, le territoire de ces moustiques s’étend rapidement (on pourra par exemple consulter les cartes de propagation d’Aedes albopictus en France métropolitaine sur le site du ministère [SanteGouv]).

Lutter contre ces maladies est devenu un sujet majeur de santé publique.
Malheureusement, pour la plupart de ces maladies, il n’existe à l’heure actuelle aucun vaccin.
Par conséquent, une stratégie de lutte efficace consiste à chercher à contrôler la population de moustiques.
Pour cela plusieurs techniques sont développées. Certaines techniques consistent à limiter les sites de reproduction des moustiques (par exemple les eaux stagnantes). Si les insecticides sont toujours largement utilisés, il est désormais bien connu que non seulement cela a un impact négatif sur l’environnement, mais également que les moustiques développent rapidement des résistances qui réduisent l’efficacité des produits. Citons aussi la technique de l’insecte stérile (voir par exemple [le site de l’IAEA]) et l’utilisation de la bactérie Wolbachia. C’est cette dernière stratégie que nous proposons d’étudier dans cet article.

Modéliser la dynamique d’une population est une branche importante des mathématiques. Les outils principalement utilisés sont les systèmes d’équations différentielles.
Dans cet article, nous proposons de montrer comment ces outils mathématiques permettent d’optimiser certaines stratégies de contrôle de population.
Dans un premier temps, à partir de la description du cycle de vie des moustiques, nous proposons une façon de modéliser mathématiquement la dynamique d’une telle population. Ensuite, sous certaines hypothèses, nous simplifions le système dynamique obtenu pour avoir un modèle plus facile à traiter mathématiquement.
Pour ce modèle simplifié, nous proposons alors une étude mathématique et écrivons le problème de contrôle pour lequel nous énonçons les résultats obtenus.

Modélisation

Le cycle de vie des moustiques

Le cycle de vie d’un moustique (mâle ou femelle) se déroule successivement dans deux milieux bien distincts : il comporte une phase aquatique (œuf, larve, pupe), et une phase adulte aérienne.
On dit qu’il s’agit d’un insecte holométabole, car il subit une métamorphose complète pour passer du stade de larve à celui d’adulte.

Quelques jours après avoir été fécondée, une femelle moustique peut déposer quelques dizaines d’œufs éventuellement répartis entre plusieurs endroits de ponte.
Une fois déposés, les œufs de certaines espèces (dont le fameux « moustique tigre », Aedes albopictus, désormais bien connu en France métropolitaine, voir [SanteGouv]) peuvent résister jusqu’à plusieurs mois et aussi à des conditions climatiques défavorables avant d’éclore. Cette caractéristique contribue à la capacité d’adaptation des moustiques, et leur a permis de coloniser des régions tempérées.
Après stimulation (par exemple une pluie), les œufs éclosent pour donner naissance à des larves qui vont se développer dans l’eau et atteindre l’état de pupes. Cette phase larvaire peut durer de quelques jours à quelques semaines. Ensuite l’insecte fait sa métamorphose. La pupe (aussi appelée « nymphe ») reste à l’état aquatique de 1 à 3 jours puis devient un moustique adulte (ou « imago ») : c’est l’émergence, puis le début de la phase aérienne. Grossièrement, la durée de vie d’un moustique adulte est estimée à quelques semaines.

Chez de nombreuses espèces, la ponte n’est possible qu’après un repas de sang, c’est-à-dire que la femelle doit piquer un vertébré avant chaque ponte. Ce comportement, qu’on appelle hématophagie, peut être mis à profit par des agents infectieux (comme des bactéries, virus ou parasites) pour se propager, passant alternativement d’un hôte vertébré (l’homme, pour ce qui nous intéresse ici) à un hôte arthropode (ici, le moustique). Lorsque cet agent infectieux cause une maladie chez le vertébré, on la qualifie de « maladie vectorielle », dont l’arthropode hématophage est le « vecteur ».

En raison de l’impact majeur des maladies vectorielles (au premier rang desquelles le paludisme) sur la santé humaine dans le monde, de nombreuses stratégies, intégrant plus ou moins les trois principaux acteurs de la transmission (agent pathogène, vecteurs, homme) ont été élaborées pour réduire leur circulation. Nous nous intéressons ici à des stratégies ciblant uniquement le vecteur (moustique appartenant au genre Aedes) de maladies virales telles que la dengue et le chikungunya.

Afin de modéliser mathématiquement ces stratégies de lutte, il nous faut dans un premier temps décrire dans un modèle la dynamique en temps de la population de moustiques. Pour cela, nous considérons un modèle compartimental représentant chaque phase du cycle de vie des moustiques tel que décrit ci-dessus. Nous introduisons les quantités suivantes :

  • $E(t)$ la densité des œufs à l’instant $t$ ;
  • $L(t)$ la densité larvaire à l’instant $t$ ;
  • $P(t)$ la densité des pupes à l’instant $t$ ;
  • $F(t)$ la densité des femelles adultes fécondées à l’instant $t$ ; dans ce modèle simplifié, on suppose que toutes les femelles sont fécondées immédiatement après leur émergence.
  • $M(t)$ la densité des mâles adultes à l’instant $t$.

Nous introduisons les paramètres suivants :

  • $\beta_E>0$ le taux de ponte des femelles ;
  • $\delta_E, \delta_L, \delta_P, \delta_F, \delta_M >0$ les taux de mort des œufs, des larves, des pupes, et des adultes femelles et mâles, respectivement ;
  • $\tau_E$ le taux d’éclosion des œufs ;
  • $\nu$ la probabilité que la pupe donne une femelle, et donc $(1-\nu)$ celle d’avoir un mâle ($0<\nu<1$) ;
  • $\tau_L$ et $\tau_P >0$ les taux de transition de la phase larvaire vers la pupe et de la pupe vers la phase adulte ;
  • la compétition intraspécifique est supposée n’intervenir que dans la phase aquatique : pour les œufs et les larves. Celle-ci traduit d’une part l’occupation des sites de pontes qui ne peuvent accueillir qu’un nombre limité de femelles, d’autre part, l’accès limité aux ressources pour les larves.
    Dans le compartiment larvaire, cette compétition est décrite par l’introduction d’une constante positive notée $c$ et est supposée dépendre linéairement de la concentration des larves : plus le nombre de larves est grand plus la compétition pour trouver les nutriments essentiels à la maturation des larves est importante.
    On note $K$ la capacité environnementale d’occupation des sites de pontes pour les femelles. Cette quantité peut être interprétée comme la densité maximale de femelles pondeuses que l’environnement est capable de supporter.

La dynamique générale peut-être représentée par le graphique ci-dessous dans lequel chaque état est représenté par une boîte reliée aux autres par des flèches décrivant la transition entre les états ; les flèches pointant vers le bas représentent la mort :

À partir des considérations ci-dessus, nous pouvons déterminer la dynamique de chaque compartiment. Ainsi, pour le compartiment des oeufs, nous obtenons
\[ \frac{d}{dt} E = \underbrace{\beta_E F}_\textrm{naissance} \underbrace{\left(1-\frac{F}{K}\right)}_\textrm{compétition pour les sites} - \underbrace{\tau_E E}_\textrm{éclosion} - \underbrace{\delta_E E}_\textrm{mort}. \]
Pour le compartiment larvaire, la dynamique modélisée dans le graphique s’écrit
\[ \frac{d}{dt} L = \underbrace{\tau_E E}_\textrm{éclosion} - \underbrace{c L^2}_\textrm{compétition intraspécifique} - \underbrace{\tau_L L}_\textrm{transition vers l'état de pupe} - \underbrace{\delta_L L}_\textrm{mort}. \]
Pour les pupes, la dynamique s’écrit de la même façon :
\[ \frac{d}{dt} P = \underbrace{\tau_L L}_\textrm{transition des larves} - \underbrace{\tau_P P}_\textrm{émergence} - \underbrace{\delta_P P}_\textrm{mort}. \]
En ajoutant les compartiments décrivant les adultes mâles et femelles, nous obtenons le système différentiel suivant :
\[ (S_1)\qquad \left\{ \begin{aligned} \frac{d}{dt} E &= \beta_E F \left(1-\frac{F}{K}\right) - E \big( \tau_E + \delta_E \big), \\ \frac{d}{dt} L &= \tau_E E - L \big( c L + \tau_L + \delta_L \big), \\ \frac{d}{dt} P &= \tau_L L - (\delta_P + \tau_P) P, \\ \frac{d}{dt} F &= \nu \tau_P P - \delta_F F, \\ \frac{d}{dt} M &= (1-\nu)\tau_P P - \delta_M M. \end{aligned} \right. \]
Afin de simplifier ce système d’équations différentielles, nous considérons que la population de pupes est à l’équilibre.
En effet, la durée de cette phase étant beaucoup plus courte que celle des autres phases du cycle de vie du moustique, nous pouvons supposer que sa dynamique temporelle est très rapide.

Par conséquent, en écrivant la troisième équation à l’équilibre, on obtient
\[ 0 = \tau_L L - \delta_P P - \tau_P P, \]
ce qui donne la relation
$P=\frac{\tau_L}{\tau_P+\delta_P} L$.

En injectant cette relation dans le système $(S_1)$, celui-ci se réduit au système suivant :
\[ (S_2)\qquad \left\{ \begin{aligned} \frac{d}{dt} E &= \beta_E F\left(1-\frac{F}{K}\right) - (\tau_E+\delta_E) E, \\ \frac{d}{dt} L &= \tau_E E - (cL+\tau_L+\delta_L) L, \\ \frac{d}{dt} F &= \frac{\nu \tau_P\tau_L}{\tau_P+\delta_P} L - \delta_F F, \\ \frac{d}{dt} M &= \frac{(1-\nu) \tau_P\tau_L}{\tau_P+\delta_P} L - \delta_M M, \\ \end{aligned} \right. \]

Introduction de la bactérie Wolbachia

La bactérie Wolbachia est un endo-symbionte (elle se trouve nécessairement à l’intérieur d’une cellule) présent dans beaucoup d’espèces d’arthropodes.
Elle n’est pas présente naturellement dans les moustiques vecteurs du genre Aedes, mais ceux-ci peuvent être infectés artificiellement.
Wolbachia se transmet verticalement de la femelle à sa descendance et induit chez ces moustiques un phénotype particulier : l’incompatibilité cytoplasmique. Ceci signifie que la fécondation de femelles « saines » par des mâles infectés par la bactérie est rendue impossible.
Ce qui a motivé encore davantage l’intérêt d’étudier cette bactérie est l’observation que lorsqu’elle infecte des moustiques (du moins chez certaines espèces infectées par certaines souches de Wolbachia), ceux-ci deviennent nettement moins susceptibles de transmettre des virus de la dengue, du chikungunya et du zika.
Par conséquent, des lâchers de moustiques infectés par cette bactérie sont en cours d’étude et de réalisation [1] dans divers pays dans l’objectif de modifier la population existant à l’état sauvage en l’infectant par la bactérie Wolbachia, réduisant ainsi considérablement leur capacité vectorielle et donc la circulation de maladies vectorielles virales.

Image au microscope de Wolbachia (en vert)

Bien évidemment, l’utilisation efficace d’une telle stratégie nécessite l’emploi d’outils mathématiques.
Nous proposons une approche déterministe et mécaniste, où l’on traduit l’objectif recherché (le remplacement de population) par un problème de contrôle optimal.

Afin de modéliser mathématiquement ce processus, il convient d’introduire la population infectée par Wolbachia dans le modèle $(S_2)$ présenté ci-dessus. Appelons donc $F_i$ les compartiments de femelles infectées par Wolbachia et $F_s$ les compartiments de femelles saines.

Compte tenu de la complexité intrinsèque du modèle biologique, il est à la fois utile et pertinent de le simplifier en utilisant certaines propriétés biologiques caractéristiques des espèces de moustiques en jeu. Cela nous conduit à considérer le système
\[ \left\{ \begin{aligned} \frac{d}{dt} F_s &= b F_s \left(1-s_h \frac{F_i}{F_s+F_i}\right)\left(1-\frac{F_s+F_i}{K}\right) - \delta_{F} F_s, \\ \frac{d}{dt} F_i &= \eta b F_i\left(1-\frac{F_s+F_i}{K}\right) - \gamma \delta_{F} F_i, \end{aligned} \right. \]
avec la notation $b=\frac{\tau_P\tau_L\beta_E}{(\tau_P+\delta_P)(\tau_L+\delta_L)}$.

Dans ce système, nous avons introduit les deux paramètres suivants :
$\eta<1$ modélisant la réduction de la fécondité des femelles infectées par rapport aux femelles saines,
$\gamma>1$ traduisant l’augmentation du taux de mortalité des individus infectés.

Pour aller plus loin.

Voici comment il est possible de simplifier le système complet. Appelons $E_i$, $L_i$, $F_i$, $M_i$ les compartiments d’œufs, de larves, de femelles et de mâles infectés par Wolbachia et $E_s$, $L_s$, $F_s$, $M_s$ les compartiments sains.
En supposant une répartition uniforme de la densité de moustiques, la probabilité pour une femelle de rencontrer un mâle infecté est égale à la proportion de mâles infectés dans la population totale, soit $\frac{M_i}{M_s+M_i}$. De même la probabilité de rencontrer une mâle sain est $\frac{M_s}{M_s+M_i}$.
Pour caractériser l’incompatibilité cytoplasmique, nous introduisons
un paramètre noté $s_h$, correspondant à la fraction d’œufs des femelles saines fécondés par des mâles infectés qui n’éclot pas. Dans le cas général, on a $0< s_h\leq 1$, le cas de l’incompatibilité cytoplasmique parfaite correspondant à $s_h=1$.
En prenant en compte les croisements entre moustiques infectés et moustiques sains, le système $(S_2)$ est modifié en :
\[ (S_3)\qquad \left\{ \begin{aligned} \frac{d}{dt} E_s &= \beta_E F_s\left(\frac{M_s}{M_s+M_i}+(1-s_h)\frac{M_i}{M_s+M_i}\right)\left(1-\frac{F_s+F_i}{K}\right) - (\tau_E+\delta_E) E_s, \\ \frac{d}{dt} L_s &= \tau_E E_s - \big(c(L_s+L_i)+\tau_L+\delta_L\big) L_s, \\ \frac{d}{dt} F_s &= \frac{\nu \tau_P\tau_L}{\tau_P+\delta_P} L_s - \delta_{F} F_s, \\ \frac{d}{dt} M_s &= \frac{(1-\nu) \tau_P\tau_L}{\tau_P+\delta_P} L_s - \delta_{M} M_s, \\ \frac{d}{dt} E_i &= \eta \beta_E F_i\left(1-\frac{F_s+F_i}{K}\right) - (\tau_E+\delta_E) E_i, \\ \frac{d}{dt} L_i &= \tau_E E_i - \big(c(L_s+L_i)+\tau_L+\delta_L\big) L_i, \\ \frac{d}{dt} F_i &= \frac{\nu \tau_P\tau_L}{\tau_P+\delta_P} L_i - \gamma \delta_{F} F_i, \\ \frac{d}{dt} M_i &= \frac{(1-\nu) \tau_P\tau_L}{\tau_P+\delta_P} L_i - \gamma \delta_{M} M_i. \end{aligned} \right. \]

Compte tenu du nombre d’équations mis en jeu, ce système est assez complexe à étudier, nous allons donc le réduire en faisant quelques hypothèses.
Tout d’abord, on suppose que les taux de mortalité des mâles et des femelles sont les mêmes ($\delta_F=\delta_M $) et que la probabilité que les pupes donnent des mâles ou des femelles est la même ($\nu=\frac 12$).
Sous ces hypothèses, les mâles et les femelles sont solutions des mêmes équations différentielles, donc si les densités sont initialement égales, ce que nous allons supposer, elles restent égales à tout instant.
Dans la suite, nous ferons donc l’hypothèse $M_s=F_s$ et $M_i=F_i$.
Une conséquence de cette hypothèse est que les équations sur les densités des mâles ne sont pas utiles pour l’étude du modèle. On se réduit ainsi à un système de $6$ équations différentielles.

Il est raisonnable aussi de supposer que le taux d’oviposition $\beta_E$ et le taux d’éclosion des œufs $\tau_E$ sont beaucoup plus grands que leur taux de mortalité (en raison de la longue survie potentielle des œufs chez les espèces appartenant au genre Aedes).
Il en résulte que nous supposerons que la dynamique des œufs est à l’équilibre, ce qui conduit à
\[ E_s \sim \frac{\beta_E}{\tau_E} F_s \left(1-s_h \frac{F_i}{F_s+F_i}\right)\left(1-\frac{F_s+F_i}{K}\right), \qquad E_i \sim \frac{\eta \beta_E}{\tau_E} F_i\left(1-\frac{F_s+F_i}{K}\right). \]
Le système $(S_3)$ se réduit donc, dans un premier temps, à un système ne portant plus que sur la densité larvaire et la densité de femelles :
\[ (S_4)\qquad \left\{ \begin{aligned} \frac{d}{dt} L_s &= \beta_E F_s \left(1-s_h \frac{F_i}{F_s+F_i}\right)\left(1-\frac{F_s+F_i}{K}\right) - \big(c(L_s+L_i)+\tau_L+\delta_L\big) L_s, \\ \frac{d}{dt} F_s &= \frac{\nu \tau_P\tau_L}{\tau_P+\delta_P} L_s - \delta_{F} F_s, \\ \frac{d}{dt} L_i &= \eta \beta_E F_i\left(1-\frac{F_s+F_i}{K}\right) - \big(c(L_s+L_i)+\tau_L+\delta_L\big) L_i, \\ \frac{d}{dt} F_i &= \frac{\nu \tau_P\tau_L}{\tau_P+\delta_P} L_i - \gamma \delta_{F} F_i. \end{aligned} \right. \]

En négligeant la compétition à l’état larvaire ($c(L_s+L_i)\ll\tau_L+\delta_L$) et en supposant le compartiment larvaire à l’équilibre, on obtient
\[ L_s = \frac{\beta_E }{\tau_L+\delta_L} F_s \left(1-s_h \frac{F_i}{F_s+F_i}\right)\left(1-\frac{F_s+F_i}{K}\right), \qquad L_i = \frac{\eta \beta_E }{\tau_L+\delta_L} F_i\left(1-\frac{F_s+F_i}{K}\right) \]
En injectant ces deux expressions dans les équations sur la densité des femelles dans le système $(S_4)$, on aboutit au système de dimension $2$.

Ce dernier système, que nous allons munir d’un contrôle pour représenter l’action humaine, est l’objet de l’étude mathématique que nous proposons d’aborder dans cet article.

Stratégie de remplacement de population et optimisation

Traduction dans un système sans dimensions et valeurs des paramètres

Puisque les moustiques infectés par la bactérie Wolbachia sont de mauvais vecteurs de maladies telles que la dengue, le chikungunya, le zika,
et comme cette bactérie se transmet de mère à enfants, une stratégie de lutte contre ces maladies consiste à relâcher des moustiques artificiellement infectés par cette bactérie dans l’objectif de remplacer la population de vecteurs existante par une population de « moins bons vecteurs », incapable de propager ces maladies.

Afin de modéliser ces lâchers, nous introduisons dans le modèle mathématique ci-dessus le flux $\widetilde{u}(t)$ décrivant la densité de moustiques (femelles adultes) infectés que l’on relâche par unité de temps :
\[ \left\{ \begin{aligned} \frac{d}{dt} F_s &= b F_s \left(1-s_h \frac{F_i}{F_s+F_i}\right)\left(1-\frac{F_s+F_i}{K}\right) - \delta_{F} F_s, \\ \frac{d}{dt} F_i &= \eta b F_i\left(1-\frac{F_s+F_i}{K}\right) - \gamma \delta_{F} F_i + \widetilde{u}. \end{aligned} \right. \]

Nous pouvons illustrer sur ce dernier système le principe de l’adimensionnement : il est préférable, en utilisant des changements d’inconnues dans ce problème, de réécrire ce système sous la forme la plus simple possible, en éliminant tous les paramètres superflus. Cela simplifie d’autant les études théoriques ou numériques à venir et l’écriture mathématique sera la plus concise. Introduisons $S(t) = F_s(\delta_F t)/K$, $I(t) = F_i (\delta_F t)/K$, $u(t) = \widetilde{u} (\delta_F t)/K$ et $\alpha = b / \delta_F$, le système ci-dessus se réécrit alors
\[ (S_6)\qquad \left\{ \begin{aligned} \frac{d}{dt} S &= \alpha S \big(1-s_h \frac{I}{S + I}\big) \big(1-(S+I) \big) - S, \\ \frac{d}{dt} I &= \eta \alpha I \big(1-(S+I) \big) - \gamma I + u, \end{aligned} \right. \]
Ce dernier paramètre peut être compris comme le nombre de reproduction de base de la population saine (sans Wolbachia), qui caractérise la croissance de la population pour le système linéarisé en $0$ (voir par exemple [le lien]).
On a donc procédé à un changement de l’échelle de temps (la nouvelle variable $t$ est « sans dimension », c’est-à-dire qu’on compte le temps relativement à la demi-vie d’un moustique femelle adulte, soit $1/\delta_F$) en même temps qu’à une renormalisation des variables (par division par la capacité environnementale). On a obtenu un système adimensionné, et mis en évidence les 4 seuls paramètres (sans unité !) déterminants de ce modèle lorsque $u=0$, rassemblés dans la table ci-dessous.

Paramètre Interprétation Gamme de valeurs
$\alpha$ Nombre de reproduction de base $2 - 100$
$s_h$ Taux d’incompatibilité cytoplasmique $0.7 - 1$
$\eta$ Facteur d’impact sur la fécondité $0.5 - 1$
$\gamma$ Facteur d’impact sur la mortalité $1 - 1.5$

Les résultats obtenus sur le système adimensionné peuvent ensuite être retranscrits aux échelles du modèle biologique, en répétant les mêmes opérations : rétablissement de l’unité pour l’échelle de temps et remise à l’échelle des populations.

États stationnaires du système dynamique

Intéressons-nous, dans un premier temps, aux trajectoires du système $(S_6)$ en l’absence d’un terme de contrôle, c’est-à-dire en posant $u=0$.
Sur la figure ci-après, on a représenté quelques trajectoires du système $(S_6)$, autrement dit le lieu des points $(S(t),I(t))_{t\in [0,+\infty)}$ associés à divers choix de conditions initiales. L’observation de cette figure révèle l’existence de quatre zones distinctes (notées A, B, C, D sur la figure) correspondant à des jeux de conditions initiales dont la trajectoire associée converge quand $t$ tend vers $+\infty$ vers un état stationnaire du système, autrement dit une solution constante du système $(S_6)$. Chaque état stationnaire (représenté par un point noir sur la figure ) correspond à une des situations suivantes :

  1. état stationnaire n° 1 (en bas à gauche sur la figure : la population totale de moustiques (sains et infectés) est éteinte ;
  2. état stationnaire n° 2 (en haut à gauche sur la figure, noté $ES_{\rm cible}$) : la population de moustiques sains est éteinte et il ne reste plus que des moustiques infectés ;
  3. état stationnaire n° 3 (au milieu de la figure : il y a co-existence de moustiques sains et infectés ;
  4. état stationnaire n° 4 (en bas à droite sur la figure, noté $ES_{\rm init}$) : la population de moustiques infectés est éteinte et il ne reste plus que des moustiques sains.

La détermination de ces états stationnaires par le calcul se fait très simplement.

En effet, il suffit de chercher les solutions constantes du système $(S_6)$.
Par conséquent, les états stationnaires sont les couples $(\overline{S},\overline{I})$ solutions du système
\[ \begin{array}{ll} 0 &= \alpha \overline{S} \left(1-s_h \frac{\overline{I}}{\overline{S}+\overline{I}}\right)\left(1-(\overline{S}+\overline{I})\right) - \overline{S}, \\ 0 &= \eta \alpha \overline{I}\left(1-(\overline{S}+\overline{I})\right) - \gamma \overline{I}. \end{array} \]
La première équation implique
\[ \overline{S} = 0 \quad \mbox{ ou } \quad s_h \frac{I}{\overline{S}+\overline{I}} + \overline{S} + (1-s_h) \overline{I} = 1-\frac{1}{\alpha} \]
La seconde équation donne
\[ \overline{I} = 0 \quad \mbox{ ou } \quad \overline{S} + \overline{I} = 1-\frac{\gamma}{\eta \alpha}. \]
On en déduit donc les quatre états stationnaires ci-dessus.
L’état stationnaire n° 1 correspond à la solution $(\overline{S},\overline{I}) = (0,0)$.
L’état stationnaire n° 2 correspond à la solution $(\overline{S},\overline{I}) =\left(0,1-\frac{\gamma}{\eta \alpha}\right)$.
L’état stationnaire n° 3 correspond à la solution
\[ (\overline{S},\overline{I}) = \left(1-\frac{1}{s_h}\left(1-\frac{\eta}{\gamma}\right)\left(1-\frac{\gamma}{\eta \alpha}\right),\frac{1}{s_h}\left(1-\frac{\eta}{\gamma}\right)\left(1-\frac{\gamma}{\eta \alpha}\right)\right). \]
Enfin, l’état stationnaire n° 4 correspond à la solution $(\overline{S},\overline{I}) =\left(1-\frac{1}{\alpha},0\right)$.

Vers un problème de contrôle optimal

Une analyse plus poussée permet de montrer que les états stationnaires n° 2 et 4 sont asymptotiquement stables, autrement dit si l’on considère une répartition initiale de moustiques proche de l’état stationnaire, alors les densités de moustiques sains et infectés convergeront vers lui en temps long. À l’inverse, les états stationnaires n° 1 et 3 sont instables, ce qui signifie qu’il ne sont pas stables, ainsi que cela peut être observé sur la figure ci-dessous.
Un tel système ayant deux états stables est souvent appelé bistable.

Cette figure représente les trajectoires du système $(S_6)$ pour les choix de paramètres : $\alpha = 5$, $\eta=0.8$, $\gamma = 1.1$, $s_h=0.8$. Les quatre points noirs représentent les états stationnaires du système.

Dans l’espoir de contrôler la population de moustiques sains, nous avons cherché à mettre en œuvre une stratégie de contrôle, autrement dit nous nous sommes demandés comment effectuer des lâchers de moustiques infectés par la bactérie Wolbachia permettant de réduire la taille de la population de moustiques sains, voire de la remplacer par une population de moustiques infectés par la bactérie Wolbachia dont la capacité vectorielle pour certains virus sera fortement réduite.

L’étude mathématique a notamment permis de déterminer s’il était préférable d’effectuer des relâchers ponctuellement dans le temps ou régulièrement, de façon continue. Décrivons en quelques mots notre approche ainsi que les premiers résultats présentés dans [ASPV].

La modélisation par un problème d’optimisation consiste à choisir une fonction coût que l’on veut minimiser durant le processus de contrôle, ainsi que des contraintes mathématiques à respecter, qui assurent la pertinence biologique des solutions.

Choix de la fonction coût.

Fixons un horizon de temps $T>0$ (le temps de l’expérience). On recherche une fonction de contrôle $u$ (le flux de moustiques à relâcher) permettant d’amener le couple $(S,I)$ le plus près possible de l’état stationnaire $ES_{\rm cible}$ (correspondant à l’état où il ne reste que des moustiques infectés) à l’issue du temps $T$. S’il existe de nombreux critères de persistance des espèces en biomathématiques, ceux-ci sont en général utilisés sur des systèmes plus simples (car monostables) et reposent souvent sur la comparaison de solutions avec des sur- et sous-solutions, mettant en jeu les valeurs propres d’opérateurs linéarisés. Malheureusement, dans le cas du système $(S_6)$, en raison de la présence de plusieurs états stationnaires de natures variées (stables et instables), il est beaucoup plus ardu de déterminer un critère simple et naturel garantissant la survie de la population de moustiques infectés. Cette observation nous a conduits à choisir une fonction coût assez naïve, de type moindres carrés, mesurant la distance de l’état final à l’état stationnaire de remplacement de population $ES_{\rm cible}$. On introduit donc le critère $J$ défini par
\[ \boxed{J(u) = \frac 12 S(T)^2 + \frac 12 ((\overline{I}-I(T))_+)^2} \]

  • on a utilisé la notation $X_+=\max \{X,0\}$ pour $X\in \mathbb{R}$,
  • on désigne par $\overline{I}=1-\frac{\gamma }{\eta \alpha}$ la quantité de moustiques infectés dans l’état stationnaire $ES_{\rm cible}$, autrement dit $ES_{\rm cible}=(0,\overline{I})$,
  • $(S, I)$ est la solution du système $(S_6)$ associé au pire choix de donnée initiale $(S(0),I(0)) = ES_{\rm init} = \left(1-\frac{1}{\alpha},0\right)$.

Choix de l’ensemble des contraintes.

Trois contraintes ont été imposées sur le contrôle :

  • puisque $u(t)$ correspond au flux de moustiques infectés par Wolbachia que l’on souhaite lâcher, on doit avoir $u(t)\geq 0$ ponctuellement,
  • pour des raisons pratiques, il n’est pas possible de créer un nombre arbitrairement grand de moustiques infectés par Wolbachia, ni de les relâcher de façon instantanée. C’est pourquoi il est raisonnable d’imposer une contrainte ponctuelle de type $u(t)\leq M$ pour tout $t \in [0, T]$. Cette contrainte modélise le fait qu’un lâcher de moustiques est nécessairement distribué en temps. Lorsqu’on relâche cette contrainte ($M \to +\infty$), à l’inverse, on autorise des lâchers impulsionnels.
  • on suppose enfin que le nombre total de moustiques utilisés dans l’expérience est borné a priori, dû aux contraintes de production des moustiques infectés en laboratoire, autrement dit que
    \[ \int_0^T u(t) dt \leq C, \]
    où $C>0$ est une constante fixée. Sans cette contrainte sur le problème, on obtiendrait une solution biologiquement sans intérêt, consistant à relâcher en chaque temps le flux maximum autorisé : $u(t) \equiv M$.

En résumé, le problème de contrôle optimal s’écrit

$(\mathcal{P})\qquad$ minimiser le critère $J$ par rapport à $u$, sous contrainte que $u$ satisfasse les trois contraintes ci-dessus.

Résultats théoriques et numériques

Analyse théorique du modèle.

Le problème de contrôle optimal $(\mathcal{P})$ semble difficile à résoudre de façon analytique.
Observant que, de façon typique, le taux d’éclosion des œufs de moustiques sains/infectés est très largement supérieur au taux de mort des individus (voir tableau ci-dessous)
nous nous sommes intéressés à un problème asymptotique, obtenu en faisant tendre les taux d’éclosion des deux populations de moustiques vers l’infini.
Comme nous allons l’expliquer ci-dessous, cette limite asymptotique permet de réduire le problème de contrôle du système $(S_6)$ à un problème de contrôle d’une équation scalaire qui peut être résolu explicitement.
La solution de ce problème asymptotique sera alors proche (en un certain sens) de la solution du problème complet $(\mathcal{P})$, dans le cas où $\alpha \gg 1$.

Afin d’effectuer cette étude analytique, nous introduisons un petit paramètre $\varepsilon\ll 1$ de telle sorte que $\alpha=\frac{1}{\varepsilon}$.
En introduisant la proportion $p$ de femelles infectées et la quantité totale de femelles $N$, définies par
\[ p(t)=\frac{I(t)}{S(t)+I(t)}, \qquad N(t) = S(t)+I(t), \]
En faisant tendre $\varepsilon$ vers $0$, nous montrons que $p$ satisfait l’équation
\[ \begin{equation}\label{eq:p} \frac{d}{dt}p=f(p)+ug(p), \end{equation} \]
avec
\[ \begin{equation}\label{eq:fg} f(p) = p(1-p) \frac{\eta - \gamma(1-s_h p) }{(1-p)(1-s_h p)+\eta p}, \quad \mbox{ et } \quad g(p) = \frac{(1-p)(1-s_h p)}{(1-p)(1-s_h p) + \eta p}. \end{equation} \]
L’équation $\eqref{eq:p}$ est complétée par la donnée initiale $p(t=0)=0$.

Ce calcul formel peut être montré rigoureusement, ce qui prouve que le système $(S_6)$ peut être réduit au système $\eqref{eq:p}$ ne portant que sur la proportion de la population infectée dans la population totale.

L’interprétation de cette « limite singulière » est la suivante : à chaque instant, la population totale évolue sur une échelle de temps rapide afin de se maintenir à son niveau d’équilibre (égal à $1$ par renormalisation). Puisque la population totale peut être considérée comme constante dès lors que cette échelle de temps est rapide devant l’échelle à laquelle on agit sur le système, il est naturel que l’état de la population puisse alors être décrit de façon univoque par la proportion d’individus infectés.

Cette réduction va bien évidemment simplifier le problème de contrôle $(\mathcal{P})$.
En effet, nous constatons tout d’abord, que
\[ \lim_{\varepsilon\to 0} J(u) = \frac {1}{2} (1-p(T))^2 + \frac{1}{2} (1-p(T))^2. \]
Nous introduisons donc le critère
\[ J_0 = \frac 12 (1-p(T))^2. \]
Le problème de contrôle optimal réduit s’écrit alors

$(\mathcal{P}_{\mbox{réduit}})\qquad $ minimiser $J_0$ par rapport à $u$, sous contrainte que $u$ satisfasse les trois contraintes précédemment mentionnées.

Pour le problème réduit $(\mathcal{P}_{\mbox{réduit}})$, nous avons démontré en utilisant des techniques de calcul des variations que la stratégie optimale de contrôle consiste à effectuer des lâchers d’amplitude maximale de façon continue, soit au début (cas favorable où les contraintes $C$ et $M$ autorisent effectivement à établir l’infection par Wolbachia) soit à la fin de l’expérience (cas défavorable où l’infection ne peut être établie durablement), voir fig. ci-dessous (cliquer pour agrandir).

PNG - 50.2 ko
PNG - 55.1 ko
PNG - 53.8 ko

Note sur le calcul des variations.

Il s’agit d’une branche active des Mathématiques qui s’intéresse à l’optimisation de fonctions de fonctions (appelées parfois fonctionnelles), qui fournit un bon cadre ainsi que de multiples outils pour résoudre le problème $(\mathcal{P}_{\mbox{réduit}})$.

Certains problèmes se résolvant maintenant par des outils de calcul des variations sont connus depuis l’Antiquité. Une légende raconte qu’en 814 avant J.-C., la reine Didon débarqua sur les côtes de l’actuelle Tunisie. Le chef local lui fournit une peau de boeuf et lui proposa de lui offrir le territoire qu’elle pourrait délimiter à l’aide de cette peau. La reine Didon découpa alors la peau en de très fines lanières. Le long de la côte, elle délimita le futur État de Carthage en formant un demi-cercle avec les lanières.

Le mathématicien polonais Hermann Schwarz en 1875, puis le mathématicien allemand Adolf Hurwitz en 1901, sont parvenus à fournir une réponse rigoureuse à un problème assez similaire bien que légèrement plus simple : « Quelle est la figure de plus grande aire que l’on peut délimiter si son périmètre est fixé ? ». Cette question préoccupait depuis longtemps les mathématiciens sans que l’on soit parvenu, jusqu’alors, à trouver le bon cadre de résolution.

La réponse n’est pas triviale. Elle est fournie par une inégalité appelée inégalité isopérimétrique,
Le problème qu’a résolu la reine Didon était encore un peu plus compliqué car il obéissait à une contrainte supplémentaire : la mer devait délimiter une partie du futur territoire de la reine. Les techniques du calcul des variations en sont venues à bout dans le courant du XXe siècle.

Bien que le résultat de contrôle soit complet dans ce cas, on peut se poser la question de l’intérêt d’une telle démarche, puisque dans les cas réels, un taux d’éclosion ne peut pas être infini (et donc $\varepsilon > 0$). Il est en réalité possible de montrer que la stratégie optimale de contrôle sera, en un certain sens, extrêmement proche des stratégies évoquées précédemment, dès lors que les taux d’éclosion sont suffisamment grands. Cela signifie en particulier que, sans mener d’analyse supplémentaire, on peut prétendre que la stratégie de contrôle proposée ci-dessus est « quasi » optimale dans un certain régime de paramètres du modèle.

Par ailleurs, nous avons également réalisé des simulations numériques sur le problème $(\mathcal{P})$ qui nous ont permis d’obtenir une idée précise de la meilleure façon d’effectuer les relâchers de moustiques infectés dans tous les cas.

Simulations numériques.

Afin de déterminer la stratégie optimale de lâcher de moustiques, nous avons cherché à estimer la/les solution(s) du problème $(\mathcal{P})$ à l’aide d’un ordinateur. Étant donné que les ordinateurs ne peuvent traiter qu’un nombre fini de données et que le contrôle $u$ cherché est une fonction, dont le graphe est constitué d’une infinité de points, nous avons en réalité résolu une approximation du problème $(\mathcal{P})$ ne mettant en jeu qu’un nombre fini de variables. Plus précisément :

  • la première étape consiste à choisir une durée caractéristique $\Delta t$ aussi petite que possible et à découper l’intervalle de temps de l’expérience en un grand nombre d’intervalle de temps d’amplitude $\Delta t$, autrement dit à écrire
    \[ t_0=0<\Delta t<2\Delta t<3\Delta t < ...< N_t \Delta t=T, \]
    avec $N_t$, un entier très grand. On dit que l’on a discrétisé l’intervalle de temps $[0,T]$.
  • au lieu de chercher à déterminer $u(t)$ en tout temps $t$, nous recherchons une fonction $u^{N_t}$ telle que $u^{N_t}$ est constante sur chaque intervalle de temps $[0,\Delta t]$, $[\Delta t, 2\Delta t]$, ..., $[(N_t-1)\Delta t,N_t\Delta t]$. On dit que $u^{N_t}$ est une approximation de $u$ constante par morceaux.
  • On adapte alors le problème $(\mathcal{P})$, c’est-à-dire qu’on le modifie en un problème dont l’inconnue est $u^{N_t}$ et non plus $u$, de telle façon que la solution de ce problème soit une bonne approximation de $u$. On obtient alors un problème dit discrétisé. Cette étape n’est en général pas simple à réaliser et nécessite de mettre en œuvre des techniques dites d’analyse numérique.
  • On utilise alors un algorithme d’optimisation en dimension finie, c’est-à-dire une méthode systématique, qui, si elle est répétée un grand nombre de fois, est susceptible de fournir une approximation de la solution du problème discrétisé.

Quelques résultats obtenus sont retranscrits sur la figure ci-dessous (cliquer pour agrandir).

PNG - 39.9 ko
PNG - 46.8 ko
PNG - 50.2 ko

 

PNG - 47.5 ko
PNG - 33.6 ko
PNG - 35.6 ko

Ces figures traduisent l’évolution des moustiques et stratégie de contrôle. De gauche à droite, $C$ (la quantité renormalisée de moustiques utilisée) prend des valeurs de plus en plus grandes.
En haut : dynamique en temps du système. La population de moustiques sains est représentée en bleu et celle des moustiques infectés par Wolbachia en rouge. En bas : stratégie optimale de contrôle (l’axe des abscisses représente le temps et l’axe des ordonnées représente l’intensité des lâchers). Dans les deux derniers cas ($C$ assez grand), on note que pour la solution du problème d’optimisation, le flux maximal autorisé $u = M$ n’est utilisé qu’au tout début et à la toute fin de l’expérience.

Conclusion

Partant d’une description simple du cycle de vie des moustiques, nous avons transformé une question opérationnelle de lutte anti-vectorielle (contre la transmission de maladies telles que la dengue) en un problème d’optimisation, par la modélisation mathématique.

Le problème mathématique que nous avons ainsi obtenu peut être attaqué tant du point de vue théorique que numérique, nous conduisant à formuler des recommandations qualitatives au sujet de la lutte anti-vectorielle.

En affinant la modélisation (par exemple, en commençant par prendre en compte les variations spatio-temporelles de la population), nous souhaitons maintenant appliquer le même schéma et, à terme, utiliser effectivement les mathématiques pour préciser certaines stratégies de lutte et par là, aider à la réduction de la circulation des maladies vectorielles.

Bibliographie

[Gates] https://www.gatesnotes.com/Health/Mapping-the-End-of-Malaria.

[APSV] L. Almeida, Y. Privat, M. Strugarek, N. Vauchelet, Optimal releases for population replacement strategies, application to Wolbachia, Preprint Hal-01807624 (2018).

[ElimDengue] www.eliminatedengue.com, consulté le 3 août 2018.

[SanteGouv] http://solidarites-sante.gouv.fr/sante-et-environnement/risques-microbiologiques-physiques-et-chimiques/especes-nuisibles-et-parasites/article/cartes-de-presence-du-moustique-tigre-aedes-albopictus-en-france-metropolitaine, consulté le 3 août 2018.

Post-scriptum :

Les auteurs remercient chaleureusement les relecteurs Michele Triestino, Baptiste Mélès, Gérard Grancher et PseudoNeo pour leurs remarques et commentaires qui ont permis d’améliorer et de clarifier cet article de façon significative.

Article édité par Stéphane Labbé

Notes

[1Les seuls projets de ce type mis en œuvre à l’heure actuelle (2018) sont menés par le consortium « Eliminate Dengue - World Mosquito Program » [ElimDengue], issu de l’Université Monash (Australie) et présent dans plus de 10 pays.

Partager cet article

Pour citer cet article :

Luis Almeida, Yannick Privat, Martin Strugarek, Nicolas Vauchelet — «Contrôle d’épidémies par lâchers de moustiques» — Images des Mathématiques, CNRS, 2018

Commentaire sur l'article

Laisser un commentaire

Forum sur abonnement

Pour participer à ce forum, vous devez vous enregistrer au préalable. Merci d’indiquer ci-dessous l’identifiant personnel qui vous a été fourni. Si vous n’êtes pas enregistré, vous devez vous inscrire.

Connexions’inscriremot de passe oublié ?

Suivre IDM