Résoudre des équations différentielles à retard pour modéliser… des marmottes?

La version 6.3 du logiciel COMSOL Multiphysics^® permet désormais d’utiliser les solutions antérieures d’une étude temporelle au pas de temps actuel. En d’autres termes, la solution d’un temps antérieur peut être utilisée pour affecter la solution du temps courant. Cela a de nombreuses applications, notamment pour la modélisation des équations différentielles à retard. Examinons cette nouvelle fonctionnalité dans le contexte d’un exemple quelque peu fantaisiste.

Rêveries d’un randonneur solitaire

Il y a quelque temps, j’ai eu la chance de randonner un peu en montagne et j’ai aperçu une marmotte. Après avoir pris quelques photos, j’ai commencé à réfléchir au cycle de vie de cet adorable herbivore. J’ai rapidement songé à trouver une équation pour décrire l’évolution de la population de marmottes qui se régalent dans les prairies de montagne.

La charmante marmotte ayant inspiré cet article de blog.

J’ai repensé à certains de mes manuels universitaires sur les méthodes numériques et je me suis souvenu de toutes sortes d’équations qui pourraient être utiles, comme les équations de Lotka–Volterra et l’équation différentielle logistique. Inspirons nous de ces dernières et partons du principe que nous devons écrire deux équations différentielles ordinaires (EDO) : une pour le taux de variation du nombre de marmottes dans la prairie, M(t), et une autre pour le taux de variation de la quantité de biomasse comestible, B(t).

La première équation, correspondant au taux de variation de la population de marmottes, comporte des termes liés à la reproduction, R(t), et à la mortalité, D(t), que nous examinerons ensuite plus en détail. La deuxième équation fait intervenir trois constantes : G_0 est le taux de croissance de la végétation, B_{max} est la masse maximale de végétation pouvant pousser dans la prairie et E_0 est le taux de consommation des marmottes par jour.

En ce qui concerne l’équation de la population de marmottes, nous considérons pour commencer un terme de reproduction qui dépend linéairement de la population actuelle:

Le taux de mortalité des marmottes, quant à lui, doit être lié aux ressources disponibles. Nous définissons donc dans une étape intermédiaire la surface fourragère par jour:

où A_0 est la taille de la prairie. Nous supposons que plus la surface fourragère est grande, plus les chances de rencontrer un prédateur sont élevées. Nous allons également considérer qu’il y a un nombre fixe de prédateurs par zone, P(t) = P_0. Le taux de mortalité des marmottes est donc proportionnel au nombre actuel de marmottes, M(t), à la surface fourragère journalière de chaque marmotte et à la densité des prédateurs. Cela conduit à une équation pour le taux de mortalité dépendant des deux inconnues:

Nous supposons également que toute la colonie de marmottes sort d’hibernation en même temps, et résolvons ces équations pour prédire l’évolution de la population au cours de l’été, après quoi la colonie retourne dans son terrier pour l’hiver. Nous allons également inclure une condition d’arrêt pour arrêter le solveur si le nombre de marmottes ou la quantité de végétation tombe en dessous d’une valeur critique, impliquant l’effondrement de la colonie ou de l’environnement.

Capture d’écran montrant comment modéliser un ensemble d’EDO. Différentes EDO peuvent avoir des unités différentes. Une Condition d’arrêt est également ajoutée au Solveur temporel .

Ces types d’équations peuvent être résolus dans COMSOL Multiphysics^® à l’aide de l’interface EDO et EAD globales, comme le montre la capture d’écran ci-dessus. Cette interface définit les deux variables à résoudre, M(t) et B(t), ainsi que l’équation de chacune d’entre elles en termes de dérivées temporelles et de valeurs initiales, M_\text{init}=M(t=0) et B_\text{init}=B(t=0). Nous pouvons résoudre ces équations pour prévoir l’évolution de la population durant l’été.

Solutions des équations de la dynamique de la biomasse végétale (vert) et de la population de marmottes (brun), avec une annotation de la population à la fin de l’été. La population de marmottes augmente d’abord fortement, puis diminue en raison d’une prédation accrue.

Les résultats montrent que la population de marmottes augmente fortement dans un premier temps et consomme rapidement la nourriture disponible. Cela conduit à une réduction de la population due à l’augmentation de la prédation. Au bout d’un certain temps, la population de marmottes atteint un équilibre. (Une remarque pour le lecteur intéressé: déterminer s’il s’agit d’un équilibre stable peut s’avérer assez difficile.)

Tenir compte de facteurs supplémentaires avec les équations différentielles à retard

À ce stade, nous disposons de résultats qui semblent plausibles, mais nous pouvons également avoir quelques réserves. Il y a au moins deux facteurs importants qui devraient être pris en compte :

1. Les marmottes commencent à s’accoupler à la sortie de l’hibernation mais ne se reproduisent pas instantanément ; il y a une période de gestation.
2. Si la population de marmottes augmente pendant une période suffisamment longue, davantage de prédateurs seront attirés dans la prairie.

Voyons comment décrire mathématiquement ces deux phénomènes et comment les mettre en œuvre dans notre modèle. Pour traiter le premier point, nous devons faire du taux de reproduction une fonction de la population au moment de l’accouplement et considérer la période de gestation, T_g, au cours de laquelle le taux de reproduction est égal à zéro. Ainsi, notre équation pour le taux de reproduction devient:

où H(t-T_g) est la fonction de Heaviside. Le taux est également réduit par la fraction \frac{M\left( T_g \right)}{M_\text{init}}, c’est-à-dire le déclin total de la population au cours de la période de gestation. Nous supposons qu’une fois que les jeunes marmottes commencent à chercher de la nourriture, la prédation des marmottes plus âgées cesse. Par conséquent, cette fraction reste fixe après t=T_g.

Nous avons maintenant une équation différentielle à retard qui peut être saisie dans le logiciel en combinant une instruction if() avec l’opérateur at(time,variable). Nous pouvons écrire le terme de reproduction sous la forme suivante:

R0*if(t > Tg,at(t-Tg,M)*at(Tg,M)/Minit,0)

Le premier argument de l’opérateur at() correspond à l’instant antérieur auquel nous voulons évaluer le second argument. Pour activer cette fonctionnalité permettant de faire référence à des instants antérieurs en cours de résolution, nous devons modifier les réglages du solveur, comme le montre la capture d’écran ci-dessous. L’utilisation d’un pas de temps Strict ne dépassant pas le temps de retard permet de bien décrire le début de la saison des naissances. Pour le reste, le modèle est résolu comme précédemment.

Capture d’écran montrant comment activer l’utilisation des opérateurs temporels en cours de résolution.

Voyons comment évolue la population lorsqu’on tient compte de ce taux de reproduction non constant. Nous observons d’abord le déclin de la population pendant la période où aucune marmotte ne naît, puis une augmentation de la population qui commence à se stabiliser.

Population de marmottes (brun) et quantité de végétation (vert) lorsqu’un terme de retard est introduit dans le taux de natalité.

Pour finir, modifions le taux de mortalité pour tenir compte de l’augmentation des prédateurs qui se produit si le nombre moyen de marmottes augmente pendant un laps de temps suffisamment long. En d’autres termes, nous devons suivre la moyenne temporelle de la population de marmottes au cours des quelques jours précédents.

Plutôt que d’évaluer cette intégrale sur tous les pas de temps précédents de la période T_p, nous pouvons appliquer le théorème fondamental du calcul et écrire une équation pour le taux de variation de la population moyenne au cours des jours précédents:

Cette équation pourrait être résolue simplement en ajoutant une autre équation globale à notre modèle, mais elle serait évaluée à partir de t=0. Nous utiliserons donc une autre instruction if() pour définir M(t<0) = M_\text{init}. Nous utilisons la moyenne mobile pour mettre à l’échelle la densité de prédateurs P(t) = \left( P_0 /M_\text{init}\right) M_\text{ave}(t), ce qui affecte le taux de mortalité. En résolvant à nouveau, nous voyons l’effet sur la population à la fin de l’été.

Bonne nouvelle ! La population a augmenté et les marmottes peuvent maintenant retourner dans leurs terriers pour leur longue hibernation.

Résultats montrant la population de marmottes (en brun) et la quantité de biomasse (en vert), le taux de reproduction et le taux de mortalité dépendant des solutions précédentes.

Les observations expérimentales confirment une augmentation de la population de marmottes.

Quelques mots sur la modélisation de bilans de population

Il convient de souligner à ce stade que l’exemple que nous avons présenté ici vise à démontrer une nouvelle fonctionnalité intéressante dans le contexte d’un modèle facile à comprendre. En fait, les modèles réels de dynamique des populations peuvent être beaucoup plus sophistiqués que cela, utilisant souvent une approche compartimentée ou basée sur les intervalles, où différentes variables sont utilisées pour suivre différentes fractions de la population totale. Le modèle SEIR (Susceptible, Exposed, Infectious, Recovered and immune) utilisé en épidémiologie est un exemple représentatif de cette approche, comme évoqué dans notre article de blog “Modeling the Spread of COVID-19 with COMSOL Multiphysics^®”.

Cette approche basée sur les intervalles est particulièrement intéressante en génie chimique, où elle est souvent utilisée pour suivre la taille des gouttelettes et des précipités. Les modèles suivants mettent en évidence ce type de modélisation :

• Oil–Water Flow Through an Orifice — A Droplet Population Model
• Precipitation of Barium Sulfate
• Crystallization of Benzoic Acid in a Mixed Suspension, Mixed Product Removal Crystallizer

Mises en garde et remarques finales

Nous avons montré les principes fondamentaux de l’implémentation et de la résolution des équations différentielles à retard dans COMSOL Multiphysics^® version 6.3. Bien que ce modèle simple puisse produire des dynamiques très intéressantes, gardez à l’esprit que cet exemple est à visée pédagogique. Une fois que vous aurez compris cet exemple simple, vous serez prêt à vous attaquer à des modèles plus sophistiquées. Vous pourriez, par exemple, vouloir combiner ces équations différentielles à retard en temps avec des équations différentielles spatiales. Chaque interface du logiciel prend en charge ces nouveaux opérateurs. Quelle que soit la combinaison d’interfaces que vous souhaitez utiliser dans votre modélisation multiphysique, la méthode à suivre pour faire référence à des solutions précédentes sera exactement la même que celle présentée ici.

Vous voulez essayer le modèle présenté ci-dessus ? Téléchargez le fichier MPH correspondant dans la Bibliothèques d’applications: