TangenteX.com

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 !

Le modèle de Lotka-Volterra est un système dynamique autonome non-linéaire. Son analyse mathématique aboutit très rapidement à des détours 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.

Pour ma part, je vous propose d'aborder ce modèle sous sa forme continue et d'en aborder les grandes lignes à travers une simulation.

Poser le problème et le mettre en équations

Le probléme 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 a également pour effet de les adimensionner! On écrira donc x(t) = X(t)/X0  et y(t) = Y(t)/Y0, avec X0 et Y0 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 \( \Delta \)x(t)/x(t) = \( \Delta \)X(t)/X(t). Pour établir la dérivée, il faut passer à la limite de ce rapport lorsque \( \Delta \)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 équations 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 dûe 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 x0 et y0 les populations initiales de proies et de prédateurs (c'est un problème de Cauchy) : c'est le système de Lotka-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 présente 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 la calcule à partir des valeurs propres de la matrice jacobienne du système). 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 TangenteX. Nous allons écrire un programme qui trace l'évolution des deux populations dans le temps ainsi que le diagramme de phases, c'est à dire l'évolution de y(t) en fonction de x(t).

Pour programmer la chose, je vous propose deux implémentations : l'une en SciLab (ou MatLab au choix...) et l'autre en Python. Ce programme est aisément portable en Maple et même en C/C++ (en utilisant la routine RK4 vue ici).

Implémentation avec SciLab

Affichage de la courbe d'évolution et d'un diagramme de phases

Ce premier programme, LotkaVolterra-1.sce, téléchargeable dans la bibliothèque de codes de TangenteX, 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');

Lignes de champ et point d'équilibre

Le programme LotkaVolterra-2.sce, téléchargeable dans la bibliothèque de codes de TangenteX, 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'équilibre 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
 end
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.

Les résultats

Affichage des courbes d'évolution des populations

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

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 diagramme de phases, 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 :

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 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 est convergent.

Implémentation avec Python

Courbe d'évolution des populations

Le script Python

L'implémentation en Python du modèle de LotkaVolterra est très similaire à son implémentation en SciLab. Nous y retrouvons la définition du système différentiel :

def LotkaVolterra(X, t=0):
    return array([ a*X[0] - b*X[0]*X[1] ,
                   c*X[0]*X[1] - d*X[1] ])

qui ne présenta aucune originalité ! Le vecteur X[0] contient la population des proies et le vecteur X[1] celle des prédateurs.

L'intégration du système est réalisée en utilisant la fonction odeint() :

X = integrate.odeint(LotkaVolterra, X0, t)
proies, predateurs = X.T

Par commodité, j'extrais la population des proies et celle des prédateurs du vecteur solution X.

En enfin j'affiche les deux courbes en utilisant les fonctions classiques de MatPlotLib.

Le script LotkaVolterra1.py est téléchargeable dans la bibliothèque de codes de TangenteX

Le résultat

En utilisant les mêmes paramètres que dans l'implémentation en SciLab, j'obtiens les courbes suivantes :

Dynamique des populations

Sans surprise, elles sont identiques à celles obtenues avec Scilab.

Diagramme de phases

Le script Python

Le script Python suit la même logique que le programme Scilab pour tracer les différentes orbites dans le plan de phase. Ce qui se traduit par le code Python suivant :

values  = linspace(0.3, 0.9, 5)
for v in values:
    X0 = v * X_f1                             
    X = integrate.odeint(LotkaVolterra,X0,t)
    plot( X[:,0], X[:,1],'r-')

Nous retrouvons la boucle qui permet de construire plusieurs orbites à partir d'un point d'origine qui varie. Ici, je ne le fait pas varier selon un point choisi par l'utilisateur, mais parmi les valeurs du vecteur values. Le nombre d'orbites est fixé par le vecteur values, que vous pouvez adapter si besoin.

Pour tracer les lignes de champ du système dynamique, matérialisées par des flèches orientées comme avec Scilab, j'utilise la fonction quiver() de MatPlotLib. Le code est un peu plus long que sous Scilab :

x = linspace(0, xmax, nb_points)
y = linspace(0, ymax, nb_points)
X1 , Y1  = meshgrid(x, y)                      
DX1, DY1 = LotkaVolterra([X1, Y1])
M = (hypot(DX1, DY1))    # calcul de la norme d'une flèche        
M[ M == 0] = 1.
DX1 /= M                 # normalisation de la norme d'une flèche
DY1 /= M
Q = quiver(X1, Y1, DX1, DY1, M, pivot = 'middle')

Le script LotkaVolterra2.py est téléchargeable dans la bibliothèque de codes de TangenteX.

Le résultat

Le diagramme de phases obtenu est très similaire à celui obtenu avec Scilab, seules les conditions initiales sont différentes :

Diagramme de phase

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 croît 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'interactions 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)

Contenu et design par Dominique Lefebvre - www.tangenteX.com octobre 2019    Licence Creative Commons    Contact :