Modéliser la dynamique des populations animales : la prédation

Dans le document Systèmes dynamiques et équations différentielles, nous avons commencé à étudier un écosystème très simple, composé de deux populations, de lièvres et de lynx, jusqu’ici considérées comme isolées l’une de l’autre. Dans ces conditions, la population des lièvres croît exponentiellement et celle des lynx décroît exponentiellement. Mais les lynx sont des prédateurs des lièvres… C’est en capturant des lièvres et en s’en nourrissant qu’ils peuvent se développer. À l’inverse, la population des lièvres est directement affectée par ces captures. L’évolution de l’effectif des lynx et celle des lièvres sont ainsi liées. Comment rendre compte mathématiquement de cette interdépendance ?

Système d’équations différentielles

Plus il y a de proies, plus il est facile pour un prédateur d’en capturer une ; symétriquement, plus il y a de prédateurs, plus les proies sont susceptibles de les rencontrer avec une issue tragique pour elles. Le nombre de proies capturées à tout instant est donc proportionnel au produit N P, produit du nombre de proies par le nombre de prédateurs. Il est donc égal à ω N P, où la valeur positive du paramètre ω est fonction de l’« efficacité » des prédateurs. Dans l’équation servant à calculer la croissance des proies, à la suite du terme de croissance intrinsèque r N, il convient donc d’ajouter le terme – ω N P qui rend compte de la capture de proies. On peut ainsi compléter l’équation différentielle :

N’ = r N – ω N P

Cette équation peut aussi être comprise en la réécrivant :

N’ = (r – ω P) N

(r – ω P) est le taux de croissance de la population N ; ce taux décroît quand le nombre de prédateurs P augmente et il tend vers r si ce nombre décroît et tend vers 0.

Symétriquement, la capture des lièvres contribue au développement des lynx. On introduit le facteur e qui rend compte de l’efficacité de la prédation (0 < e < 1) et on obtient l’équation différentielle :

P’ = – m P + e ω N P.

On obtient ainsi deux équations différentielles couplées, car P intervient dans la première et N dans la seconde. Elles sont de plus non-linéaires par la présence du produit N P.

Une autre manière d’interpréter un tel système d’équations différentielles consiste à le réécrire en calculant l’intégrale de 0 à t des deux membres qui les composent.

On obtient ici :

\[N(t) = N_0 + \int_{0}^{t} (rN(u) \,- \omega N(u) P(u)) du\]
\[P(t) = P_0 + \int_{0}^{t} (-m P(u) \,+ e \omega N(u) P(u)) du\]

La première équation se lit alors : la valeur de N à l’instant t est égale à la somme algébrique de toutes les contributions, positives et négatives, pendant l’intervalle de temps [0, t], à l’effectif de la population initiale N(0).

On peut réécrire ces équations intégrales afin de bien distinguer les contributions positives et négatives :

\[N(t) = N_0 + \int_{0}^{t} rN(u) du \,- \int_{0}^{t} \omega N(u) P(u) du\]
\[P(t) = P_0 \,- \int_{0}^{t} mP(u) du \,+ \int_{0}^{t} e \omega N(u) P(u) du\]

Les valeurs N(t) et P(t) intègrent – résument en quelque sorte – l’histoire du système, de l’instant 0 à l’instant t : [N(t), P(t)] est l’état du système à l’instant t ; N et P sont appelées variables d’état.

Les modèles constitués d’équations différentielles dites « ordinaires », comme celles employées ici, sont déterministes : l’évolution en est entièrement déterminée, pour un jeu de valeurs des paramètres, par la donnée d’un état initial [N₀, P₀]. Symétriquement, la donnée de l’état [N(t), P(t)] à l’instant t résume le passé du système. Il détermine également son futur.

Résolution des équations

Le caractère non-linéaire de ce système d’équations différentielles rend difficile la recherche de solutions explicites, c’est-à-dire des expressions mathématiques des fonctions N(t) et P(t).

Si la résolution sur calculateur analogique est restée longtemps la seule disponible, elle est supplantée depuis plusieurs décennies par la résolution numérique sur ordinateur. Étudions-en le principe à travers l’exposé de la méthode la plus simple, la méthode d’Euler.

On sait, par définition même de la dérivée, que N’(t) est égal à la limite, quand dt tend vers 0, de (N(t+dt) – N(t)) /dt. Pour une valeur de dt suffisamment petite, soit Δt, on peut écrire l’approximation :

N(t+Δt) ≈ N(t) + Δt N’(t)

Et donc :

N(t+Δt) ≈ N(t) + Δt [r N(t) – ω N(t) P(t)]

De même :

P(t+Δt) ≈ P(t) + Δt [– m P(t) + e ω N(t) P(t)]

Ces expressions livrent un schéma itératif pour calculer les couples de valeurs [N(Δt), P(Δt)], [N(2Δt), P(2Δt)], [N(3Δt), P(3Δt)]… autrement dit les états successifs du système aux temps Δt, 2Δt, 3Δt, etc. à partir des valeurs initiales N(0) et P(0) et à condition, bien sûr, de fixer des valeurs aux paramètres r, m, ω et e.

La valeur Δt ne possède ici aucune interprétation vis-à-vis du système de modélisation. Elle est choisie suffisamment petite pour que l’approximation soit satisfaisante, mais pas trop pour ne pas alourdir les calculs.

Il existe de nombreuses autres méthodes de résolution numérique (on parle aussi d’intégration numérique) de systèmes d’équations différentielles. Elles mettent en œuvre des formules plus complexes d’approximation des dérivées. Certaines ajustent automatiquement le pas d’intégration Δt en fonction de l’erreur commise.

Plutôt que d’afficher les valeurs de N et de P au cours du temps sous la forme de colonnes de nombres fastidieuses et peu parlantes, il est préférable de faire tracer les graphes des fonctions N(t) et de P(t). La valeur de Δt pouvant être fixée arbitrairement petite, il est possible d’obtenir par calcul de très nombreux points et les courbes affichées, appelées trajectoires, apparaîtront continues.

Dans ce cas particulier de deux variables, il est aussi possible de tracer, dans le plan (N, P) la courbe reliant, dans l’ordre
chronologique, tous les points de coordonnées [N(t), P(t)] produits par une résolution numérique.

C’est en modifiant les conditions initiales et les valeurs des paramètres et en relançant la résolution numérique du système d’équations différentielles qu’il est possible d’accroître notre compréhension du comportement du système dynamique étudié. Le modèle a été élaboré terme par terme, équation par équation ; la simulation permet de déterminer le comportement qui résulte de l’interaction de ces termes et équations ; elle en effectue en quelque sorte la synthèse.

Pour accéder à l’applet Java, autorisez les applets du domaine https://interstices.info. Si votre navigateur n’accepte plus le plug-in Java, mais que Java est installé sur votre ordinateur, vous pouvez télécharger le fichier JNLP, enregistrez-le avant de l’ouvrir avec Java Web Start.

Applet réalisée par Hamdi Ben Abdallah.

Système proies – prédateurs. L’évolution des effectifs des populations des proies et des prédateurs apparaît sous deux formes. À gauche, les courbes des effectifs en fonction du temps. Les variables ont été normées de façon à varier entre 0 et 1, et les valeurs des paramètres des équations modifiées en conséquence. On observe les oscillations périodiques des effectifs. Un examen plus approfondi des courbes montre un décalage entre le pic de l’effectif d’une population et le pic de l’autre. À droite, le comportement du système à deux variables proies – prédateurs est représenté dans le plan : en abscisse, l’effectif des proies, en ordonnée, celui des prédateurs. Dans ce plan, le comportement oscillatoire se traduit par une courbe fermée.

Analyse qualitative

Bien entendu, le changement de valeur d’un ou plusieurs paramètres ou d’une ou plusieurs conditions initiales demande que la résolution numérique soit relancée. Il s’agit là de l’inconvénient majeur de la résolution numérique par rapport à la résolution mathématique explicite.

Cependant, sans même résoudre analytiquement ce système d’équations différentielles, il est possible d’en faire une étude générale sans affecter de valeurs numériques aux paramètres. On parle alors d’étude qualitative. Il est ainsi possible de discuter de l’existence d’un ou de plusieurs états d’équilibre, c’est-à-dire de valeurs N* et P* de N et de P pour lesquelles les dérivées N’ et P’ sont nulles : le système ne peut quitter de lui-même cet état et les valeurs de N et de P restent constantes dans le temps.

Il est de plus possible de caractériser un point d’équilibre, qui peut être stable si une petite perturbation de N ou de P est suivie d’un retour au point [N*, P*], ou instable dans le cas contraire. Il est enfin possible de tracer le « portrait » du système dynamique en délimitant, dans le plan (N, P), les régions où les signes de N’ et P’ sont constants, en traçant les lieux des points où la dérivée P’ ou N’ est nulle, etc.

Cet exemple très simple nous a permis de montrer comment un modèle se construit à partir d’hypothèses sur le système qui simultanément se voit précisé dans ses frontières (quelles entités retenir ?) et dans les interactions qu’entretiennent les entités entre elles.

Extensions du modèle

Nous pouvons poursuivre cette démarche de modélisation. Reconsidérons tout d’abord l’hypothèse, bien peu réaliste, que les proies disposent de ressources non limitées et reprenons l’équation différentielle de la croissance des proies en l’absence des prédateurs : N’ = r N. Nous souhaitons exprimer le fait que cette croissance est limitée par les ressources disponibles et qu’au fur et à mesure que la population croît, elle consomme ces ressources (nourriture, espace) dont la réduction affecte en retour la croissance. Nous introduisons un terme de freinage dans l’équation de croissance de la population des proies non soumises à prédation :

N’ = r (1 – N / K) N, où K est un paramètre à valeur positive, appelé « capacité limite ».

Lorsque N croît à partir d’une valeur N(0) < K, N/K croît également, puisque N’ est positive, et tend vers 1 ; (1 – N / K) tend donc vers 0 et N’ également, tout en restant positive ; N continue donc à croître, mais toujours plus lentement en approchant de la valeur K. On observe une croissance « en S », appelée croissance logistique.

En introduisant cette formulation dans notre système d’équations précédent, nous obtenons :

N’ = r (1 – N / K) N – ω P N

P’ = – m P + e ω N P

Le comportement de ce système peut être lui aussi étudié en lançant des résolutions numériques en modifiant les valeurs des paramètres r, K, ω, m et e, ainsi que les conditions initiales N(0) et P(0). La simulation fait apparaître des comportements oscillatoires amortis des deux populations.

Pour accéder à l’applet Java, autorisez les applets du domaine https://interstices.info. Si votre navigateur n’accepte plus le plug-in Java, mais que Java est installé sur votre ordinateur, vous pouvez télécharger le fichier JNLP, enregistrez-le avant de l’ouvrir avec Java Web Start.

Applet réalisée par Hamdi Ben Abdallah et Régis Monte.

Système proies – prédateurs avec terme logistique. L’introduction d’un terme logistique dans l’équation qui régit la dynamique de l’effectif des proies produit un comportement oscillatoire amorti des effectifs des deux populations, d’autant plus amorti que la valeur de K est petite. Les valeurs de N et de P tendent vers un équilibre (N*, P*). En ce point, les dérivées dN/dt et dP/dt sont nulles et les variables N et P gardent les valeurs N* et P*. Sur la droite, la courbe n’est plus fermée ; elle converge vers le point d’équilibre (N*, P*).

Des populations non homogènes

Dans les modèles décrits jusqu’à présent, une seule variable est associée à une population. L’équation différentielle correspondante rend ainsi compte du comportement en moyenne des individus de la population. En particulier, l’âge des individus n’intervient pas. Ainsi, tout individu de la population des prédateurs est-il susceptible de chasser ; à l’inverse, aucune proie n’est plus facile à capturer qu’une autre et toute capture procure le même bénéfice à la population des prédateurs. Si cette hypothèse se révèle inadéquate, il peut se révéler nécessaire de scinder la population en plusieurs tranches d’âge et associer une variable, et donc une équation différentielle, à l’effectif de chaque tranche.

Soit ainsi N_i(t) l’effectif d’une tranche d’âge i de la population étudiée. La variation dans le temps de cet effectif résulte des individus vieillissants en provenance de la tranche inférieure i-1, de la mortalité dans la tranche d’âge i, et du vieillissement des individus qui quittent de ce fait cette tranche d’âge et passent dans la tranche supérieure.

Nous avons vu précédemment que dans l’équation différentielle dx/dt= – x / λ, λ s’interprétait comme la durée de vie moyenne d’un individu de la population d’effectif x (t). Autrement dit, λ est le temps de résidence moyen d’un individu dans la population x(t).

Il s’en suit, en l’absence de toute autre interaction (absence de prédation), que la population d’une tranche d’âge i évolue selon l’équation différentielle :
\[\frac{dN_i(t)}{dt} = \frac{N_{i-1}(t)}{\lambda_{i-1}} \,- m_i N_i(t) \,- \frac{N_i(t)}{\lambda_i}\]
où λ_i est la « largeur » de la tranche d’âge i et m_i son taux de mortalité.

L’équation relative à la première tranche d’âge est quelque peu différente, puisque c’est la seule dont l’accroissement est dû aux naissances, dont le nombre peut être calculé en appliquant un taux de natalité spécifique à chacune des tranches d’âge. En l’absence ici aussi de toute autre interaction, l’évolution de l’effectif de la première tranche d’âge N₁ est ainsi régie par l’équation :
\[\frac{dN_1}{dt} = \sum_{i=1}^{n} r_i N_i \,- \frac{N_1}{m_1}\]
où r_i est le taux de natalité de la tranche i et n est le nombre de tranches d’âge.

Une grande diversité de modèles

L’examen de la littérature spécialisée dévoile quantités de variantes et d’extensions de ces premiers modèles présentés ici. Il s’agit en particulier d’ étudier des écosystèmes dans lesquels un grand nombre d’espèces interagissent au sein de réseaux trophiques. Une espèce y est simultanément la proie d’une ou plusieurs autres et le prédateur, qui peut consommer une ou plusieurs espèces, est lui-même la proie d’une ou plusieurs autres. À raison d’une variable d’état, et donc d’une équation différentielle, par espèce, la taille des systèmes différentiels croît, et avec elle, le nombre de paramètres dont il faut déterminer la valeur, ainsi que le nombre de simulations à effectuer pour mieux appréhender le comportement de ces systèmes complexes.

Au sein des populations humaines, les systèmes d’équations différentielles sont mis en œuvre pour étudier des phénomènes de propagation, en particulier d’épidémie.