Le modèle de Lotka-Volterra

Un modèle de dynamique des populations

Un problème pratique de pêche

Vito Volterra, mathématicien et physicien italien du début du XXeme siècle fut confronté à un problème curieux de dynamique des populations. On lui exposa le phénomène suivant: entre 1905 et 1923, les effectifs des populations de sardines et de leurs prédateurs en mer Adriatique semblaient osciller, et qui plus est, en décalage de phase.

Volterra, spécialiste des systèmes dynamiques et des équations différentielles qui les décrivent, a eu vite fait de construire un modèle d'évolution de ces populations, qu'il publia en 1926. Un mathématicien américain, Alfred James Lotka, publia indépendamment en 1924 (ou 25 les dates varient selon les sources), un modèle équivalent. On le désigne aujourd'hui sous le nom de "modèle de Lotka-Volterra" pour ne pas faire de jaloux... Il a été étendu à plusieurs domaines de la dynamique des populations et est utilisé aujourd'hui très largement en biologie, écologie et même chimie. Pour ne citer que le cours d'écologie de Ricklefs et Miller, tout un chapitre (le 9) est consacré à l'étude des écosystèmes et au modèle de Lotka-Volterra. En passant, les auteurs, américains, citent Lotka et pas Volterra, ce qui n'est pas très élégant! Et je parle bien sur d'un cours d'écologie scientifique, pas d'écologie politique.

Le modèle de Lotka-Volterra est un système dynamique autonome non-lnéaire. Son analyse mathématique aboutit très rapidement à des détours profonds dans la théorie des équations différentielles. Aussi, ai-je été assez surpris de voir surgir le modèle de Lotka-Volterra dans le programme d'enseignement de spécialité en mathématique de terminale S. Il s'agit d'une approche discrète, dans le prolongement du cours sur les suites, et dans le cadre d'un TD sans doute..

A la réflexion, ce n'est pas idiot, à condition qu'on en reste au niveau d'un système discret, construit sur la manipulation de suites. Et encore, je ne suis pas persuadé que les élèves aient un bagage et une motivation pour les maths suffisants pour explorer ce modèle.

Pour ma part, je vous propose d'aborder ce modèle sous sa forme continue et d'en explorer les grandes lignes à travers une simulation, comme il est de coutume sur TangenteX.com.

Poser le problème et le mettre en équation

La problématique est relativement simple : étudier l'évolution dans le temps de deux populations, celle des proies que je note X(t) et celle des prédateurs que je note Y(t). Il n'aura échappé à personne que ces deux fonctions sont des fonctions de \( \mathbb{R} \)+ vers \( \mathbb{N} \), le nombre d'individus d'une population étant entier, jusqu'à preuve du contraire. Il n'est pas évident de travailler sur ce genre de fonctions, pour lesquelles nous ne disposons pas d'outils simples. Par contre, pour les fonctions de \( \mathbb{R} \)+ dans \( \mathbb{R} \), supposées continues et dérivables sur le domaine, nous disposons de tout l'arsenal des équations différentielles, sous une réserve toutefois...

Pour transformer nos X(t) et Y(t) en fonctions réelles, nous allons les normaliser, c'est à dire les transformer en rapport, ce qui également pour effet de les adimensionner! On écrira donc x(t) = X(t)/X₀ et y(t) = Y(t)/Y₀, avec X₀ et Y₀ respectivement la population à l'instant t0 des proies et des prédateurs.

Ce qui nous intéresse, c'est d'étudier les variations des populations pour un intervalle de temps très petit. Pour la population de proies (resp. prédateurs), cette variation s'écrit Δx(t)/x(t) = ΔX(t)/X(t). Pour établir la dérivée (matière fondamentale d'une équation différentielle!), il faut passer à la limite de ce rapport lorsque Δt tend vers 0. Oui, mais ce n'est réaliste que si les populations X(t) et Y(t) sont suffisamment grandes! C'est une limite du modèle: pas question de l'appliquer sur un groupe de 10 lapins et 5 renards (dans les simulations, les unités sont quelconques et les fonctions continues et dérivables par hypothèse)!

On supposera que les conditions mentionnées ci-dessus sont remplies et on va tenter de mettre en équation notre problème!

Pour ce faire, je vais mettre en équation les hypothèses tirées de l'écologie du système:

s'il n'existe aucun prédateur, le taux de croissance de la population de proies est constant, soit x'(t) = a*x(t), en nommant a ce taux de reproduction. Vous noterez que je ne fixe aucune limite à la croissance de la population de proies, ce qui est irréaliste sur le plan de l'écologie! Les ressources sont limitées (eau, nourriture) et il existe une compétition pour le territoire entre les proies. Un point à perfectionner...
s'il n'existe aucune proie, les prédateurs meurent de faim! On traduit cela par l'équation y'(t) = -c*y(t). c est le taux de mortalité des prédateurs. Notez le signe - qui traduit la diminution de la population.
on supposera que le taux de mortalité des proies du aux prédateurs est proportionnel au nombre de prédateurs. Ce qui nous donne une version améliorée de l'équation en x'(t), soit x'(t) = (a - b*y(t))x(t). Le coefficient b étant le taux de mortalité des proies en présence de prédateurs
Enfin, autre supposition, le taux de croissance de la population de prédateurs est proportionnel à la population des proies (il y a plus à manger!!), d'où y'(t) = (-c + d*x(t))y(t). Le coefficient d est le taux de reproduction des prédateurs en présence de proies.

J'obtiens donc le système différentiel autonome, c'est à dire dans lequel le temps n'apparait pas explicitement, ci-dessous:

x'(t) = x(t)*(a - b*y(t))
y'(t) = y(t)*(-c + d*x(t))

avec x₀ et y₀ les populations initiales de proies et de prédateurs (c'est un problème de Cauchy..). C'est le système de Lokta-Volterra.

Il n'existe pas de solution analytique à ce système, mais on peut au moins en dire quelques mots:

le théorème de Cauchy-Lipschitz nous dit qu'il existe des solutions locales et globales
on en déduit aussi que les solutions sont positives. C'est heureux...
Un autre théorème (Lyapunov) nous permet de démontrer que le système a deux états d'équlibre : (x,y) = (0,0) et (x,y) = (d/c, a/b).
les solutions du système de Lotka-Volterra sont périodiques
Sur une période, on sait calculer la population moyenne de proies <x> = d/c et <y> = a/b
Le calcul de la période n'est pas vraiment simple. On se contentera de la mesurer sur les données simulées.

Utiliser les équations pour construire une simulation

Muni du système différentiel, l'écriture d'un programme de simulation n'est pas très compliqué et nous l'avons déjà souvent fait sur notre site. Nous allons écrire un programme qui trace l'évolution des deux populations dans le temps ainsi que le portrait de phase, c'est à dire l'évolution de y(t) en fonction de x(t).

Pour programmer la chose, j'utiliserai Scilab (ou son équivalent commercial MatLab), langage déjà largement utilisé sur TangenteX.com.

Les programmes

Affichage de la courbe d'évolution et d'un portait de phase

Ce premier programme, LotkaVolterra-1, permet de tracer les courbes d'évolution des populations de proies et de prédateurs et un portait de phase. Il permettra d'illustrer le caractère périodique des solutions et leur déphasage.
Quelque commentaires sur le code, et d'abord la partie caractéristique de ce genre de programme, c'est à dire la fonction qui définit le système différentiel:

function [w] = LotkaVolterra(t,y)

w(1) = a*y(1) - b*y(1)*y(2);

w(2) = c*y(1)*y(2) - d*y(2);

endfunction

Vous noterez que le programme reprend exactement le système défini dans le chapitre précédent, dans lequel y1 représente la population des proies et y2 la population des prédateurs.
Rien à dire de particulier sur la programmation, qui est semblable à tous les autres exemples de programmes Scilab qui résolvent des EDO.
Je vous passe les initialisations diverses, comme celle des paramètres a,b,c et d de simulation. Bien sur, vous pourrez changer les valeurs choisies dans le source pour faire vos expériences, et même prévoir une saisie par l'utilisateur des valeurs de ces paramètres.

La résolution du système prend une ligne :

y = ode(y0,t0,t,LotkaVolterra);

Là aussi, rien de plus classique!
Et enfin, l'affichage. J'ai séparé la fenêtre graphique en deux avec l'instruction subplot pour afficher les deux courbes:

subplot(2,1,1);plot2d(t,y(1,:),style = 2);xtitle('Evolution des populations','Temps','Population');

subplot(2,1,1);plot2d(t,y(2,:), style = 3);

subplot(2,1,2);plot2d(y(1,:),y(2,:), style = 5);xtitle('Portrait de phase','Proies','Prédateurs');

Le source est téléchargeable ici.

Lignes de champ et point d'équilibre

Le programme LotkaVolterra-2 trace tout d'abord les lignes de champ figurées par des flèches tangentes aux orbites. Pour information, on appelle "orbite" la trajectoire du système y2 = f(y1) du système dans le plan de phase.
En cliquant sur un point quelconque du plan de phase, vous obtiendrez l'orbite correspondant aux conditions initiales x0,y0 correspondantes aux coordonnées du point cliqué. Cela permet de tracer plusieurs orbites et d'étudier ainsi la position du point d'équlibre non trivial (différent de (0,0))

LotkaVolterra-2 utilise bien sur le même système différentiel que le programme précédant. Il utilise aussi les mêmes données d'initialisation et les paramètres de simulation. Les points caractéristiques du code sont peu nombreux:

le code d'affichage des lignes de champ
la boucle de calcul et d'affichage des orbites en fonction des conditions initiales sélectionnées par l'utilisateur

Le code d'affichage des lignes de champ:

X = 0:dx:xmax;

Y = 0:dy:ymax;

fchamp(LotkaVolterra,0,X,Y);

Le code de la boucle de traitement des orbites :

xtitle("Modèle proies-prédateurs - Diagramme de phase");

while %t

[cbouton,x0,y0] = xclick();

if cbouton == 5 then

break

else

if x0 >= 0 & x0 <= xmax & y0 >= 0 & y0 <= ymax then

[y] = ode([x0;y0],t0,t,LotkaVolterra);

plot2d(y(1,:),y(2,:),style = 5);

end

Pour sortir de la boucle, cliquez sur le bouton droit de votre souris. Si vous voulez plus d'information sur les fonctions fchamp et xclick de Scilab, consultez l'aide intégrée de Scilab.

Le source du programme est téléchargeable ici.

Les résultats

Affichage des courbes d'évolution des populations

Voici les courbes produites par LotkaVolterra-1, avec les paramètres du code:

Evolution des populations proies - prédateurs selon le modèle de Lokta-Volterra

Sur la courbe du haut, en vert l'évolution de la population des prédateurs et en bleu celle des proies.

On constate plusieurs choses:

les amplitudes maxima sont constantes dans les deux populations. C'est écologiquement irréaliste, notre modèle est un peu simpliste.
la période est identique pour les deux populations, et vaut un peu moins de 3, lu sur la courbe. Le calcul donne environ 2.6...
le déphasage entre les maxima d'amplitude des populations est constant, ce qui est encore écologiquement irréaliste. Toujours la rusticité de notre modèle...

Le portait de phase, tracé ici pour une seule solution, n'a pour seule ambition que d'illustrer l'allure de la variation relative des deux populations. Il s'agit de la courbe y2 = f(y1), c'est à dire la population des prédateurs en fonction de celle des proies. La courbe est fermée, ce qui laisse deviner la stabilité du système, mais nous verrons ça ci-dessous avec les courbes du programme LotkaVolterra-2.

Affichage du champ de vecteurs et recherche du point d'équilibre

Voici le résultat restitué par LotkaVolterra-2:

Modèle de Lokta-Volterra - champ de vecteurs

On note d'abord les lignes de champ, matérialisées par les flèches sur le plan de phase. Vous pouvez noter que toutes les flèches sont tangentes à une orbite, ce qui est assez normal pour une ligne de champ.. Le tracé de ce champ permet d'avoir une bonne idée de la dynamique du système sans calculer sa solution. En particulier, on devine le point d'équilibre en (2,3), qui correspondant à la valeur théorique.

En rouge, apparaissent les différentes orbites tracées en choisissant des points initiaux de l'orbite au hasard dans le plan de phase. Plus le point initial se rapproche du point d'équilibre et plus l'orbite est petite et apparait elliptique. Les courbes sont fermées, ce qui montrent que notre schéma de calcul (celui de Scilab) est convergent.

Pour aller plus loin

Comme je l'ai écris plus haut, certaines des hypothèses retenues pour établir le système différentiel sont plus ou moins réalistes du point de vue de l'écologie des populations. Par exemple:

en absence de prédateurs, on considère que la population des proies croit exponentiellement, selon le taux de reproduction a. C'est ne pas tenir compte de la limitation en ressources, qui régule l'exponentielle! On pourrait introduire une contrainte, en s'inspirant de l'équation de la logistique, par exemple.
les fonctions d'interaction proies/prédateurs (figurées par les constantes b et d) peuvent être améliorées. J'utilise ici la "loi d'action de la masse", c'est à dire que chaque prédateur mange toutes les proies qu'il rencontre, le nombre de proies étant proportionnel à la population de proies. Ce n'est pas très réaliste! Le prédateur ne peut manger qu'un certain nombre de proies, quelque que soit la population disponible, car il sera vite rassasié! Pour en tenir compte, on remplace les termes constants par des fonctions comme la fonction de Holling type I ou II, ou d'autres fonctions (voir par exemple users.aims.ac.za/.../AIMS-Prey-predator-Models.ppt ou encore http://lama.univ-savoie.fr/~labart/generated/fichiers/MATH705/Pop_video.pdf)