Experimentarium Digitale

Version stochastique du modèle de Lotka-Volterra

Chazottes Jean-René , Monticelli Marc — 2018-07-13T14:58:10Z

$$ \begin{cases} \dot{x}=x(1-x-y)\\ \dot{y}=\beta(x-\alpha)y \end{cases} $$

Version stochastique du modèle logistique

Chazottes Jean-René , Monticelli Marc — 2018-07-13T14:56:42Z

On considère le modèle logistique déterministe

$$ \dot{x}= x\left(\gamma-cx\right) $$

où $\dot{x}$ signifie la dérivée de $x$ par rapport au temps $t$. Si on se donne $x(0)=x_0\geq 0$, la solution $x(t)$ au temps $t$ peut représenter la densité d'une population.
On a deux paramètres. Le paramètre $\gamma=\lambda-\mu$, avec $\lambda,\mu>0$ qu'on va définir plus bas, et on suppose que $\lambda>\mu$. Le paramètre $c>0$ mesure l'intensité de la compétition intraspécifique : la croissance de $x(t)$ qui serait exponentielle si $c$ valait $0$ est contrecarrée par le terme $-cx^2$, d'autant plus fortement que $c$ est grand. Il y a deux points fixes : $x=0$, qui est instable (répulsif), et $x_*=\gamma$, qui est stable (attractif). Si $x(0)>0$, $x(t)$ va tendre vers $x_*$ quand $t\to+\infty$.

Parallèlement, on introduit une version stochastique du modèle logisitique déterministe sous la forme d'un processus (de Markov) de naissance et mort $N(t)$ défini par les taux de saut suivants où $i$ est un entier $>0$ :

$i\leadsto i+1$ avec un taux $\lambda i$
$i\leadsto i-1$ avec un taux $\mu i+ci^2$.

On se donne un nombre d'individus $N_0$ au temps $t=0$.
Si au temps $t$ on a $i$ individus ($i>0$) alors soit l'horloge qui suit une loi exponentielle de paramètre $\lambda i$ sonne la première et on passe à $i+1$ individus, soit c'est l'horloge suivant une loi exponentielle de paramètre $\mu i$ qui sonne la première et on passe à $i-1$ individus. Et ainsi de suite.
On constate que si $i=0$, on reste dans cet état, ce qui modélise l'extinction dont on peut démontrer qu'elle a lieu avec probabilité un. Elle est donc certaine mais le temps d'extinction est aléatoire.

Comment relier ces deux modèles ? L'idée est « renormaliser » les sauts avec un paramètre d'échelle $K>0$ pour qu'ils soient de plus en plus petits quand $K$ augmente, autrement dit on renormalise $N(t)$ par $K$, ce qui fait que $N(t)/K$ est une densité. On va donc avoir des trajectoires du processus qui vont tendre vers une trajectoire continue, sans sauts, quand $K\gg 1$. A cause du terme non-linéaire $ci^2$ modélisant la compétition, on se rend compte qu'il faut aussi renormaliser $c$ par $K$. Autrement dit, on a maintenant :

$\frac{i}{K}\leadsto \frac{i+1}{K}$ avec un taux $\lambda i=\lambda K \frac{i}{K}$
$\frac{i}{K}\leadsto \frac{i-1}{K}$ avec un taux $\mu i+\frac{c i^2}{K}= K \frac{i}{K}\big(\mu+\frac{c}{K}i\big)$.

Enfin, n'oublions pas que résoudre l'équation différentielle du modèle logisitique déterministe nécessite qu'on se donne une condition initiale $x_0$ et qu'il faut la relier à $N_0$. Ici on prend tout simplement $x_0=\frac{N_0}{K}$. Dans l'expérience numérique interactive suivante, on peut observer ce qui se passe. On a pris $c=1$ puisqu'il ne joue par de rôle. On a fixé $x_0=0,5$.

On observe ce à quoi on s'attend intuitivement : quand $K$ est petit, disons égal à $10$, le processus fluctue beaucoup autour de la solution du modèle déterministe et on observe facilement son extinction qui a lieu souvent dans la fenêtre de temps qu'on a fixée. Quand $K$ devient grand, c'est l'inverse.

Pour terminer, précisions que ces observations peuvent être transformées en énoncés mathématiques. On renvoie à cet article pour plus de détails.

Version stochastique du modèle de compétition entre 2 espèces

Chazottes Jean-René , Monticelli Marc — 2018-07-13T14:56:13Z

Version stochastique du modèle de Allee

Chazottes Jean-René , Monticelli Marc — 2018-07-13T14:51:01Z

$$ \dot{x}=x\left( \frac{x}{\gamma_0}-1\right)\left( 1-\frac{x}{\gamma}\right) $$

Version stochastique du modèle proie-prédateur de Rosenzweig-McArthur

Chazottes Jean-René , Monticelli Marc — 2018-04-05T09:24:48Z

Nous avons déjà présenté le modèle proie-prédateur de Rosenzweig-McArthur

$$ \begin{cases} \dot{x} = x \left(1-\frac{x}{\gamma}\right)- \frac{xy}{1+x} \\ \dot{y} =\beta y\left(\frac{x}{1+x}-\alpha\right) \end{cases} $$

où $x(t)$ est la densité des proies, $y(t)$ celle des prédateurs, et $\alpha,\beta, \gamma$ sont des paramètres positifs. Nous montrons ici comment il est lié à un processus markovien de saut $(N_x(t), N_y(t))_{t\geq 0}$, convenablement normalisé par un paramètre $K$ qui mesure l'échelle des deux populations. On note $N_x(t)$ le nombre de proies au temps $t$ et $N_y(t)$ le nombre de prédateurs au temps $t$, qui sont des variables aléatoires. Supposons qu'au temps $t$ il y ait $n_x\geq 1$ proies et $n_y\geq 1$ prédateurs (c-à-d, $N_x(t)=n_x$ et $N_y(t)=n_y$). Entre les temps $t$ et $t+h$ ($h$>0) la probabilité qu'il y ait
– la naissance d'une proie vaut

$$ n_x h + o(h) \, ; $$

– la mort d'une proie dûe à la compétition entre proies vaut

$$ \frac{n_x^2}{\gamma K} h+ o(h) \, ; $$

– la mort d'une proie dûe à la prédation vaut

$$ \frac{n_x n_y}{K+n_x} h+ o(h) \, ; $$

– la naissance d'un prédateur vaut

$$ \frac{\beta n_x n_y}{K+n_x}+ o(h) \, ; $$

– la mort d'un prédateur vaut

$$ \alpha \beta n_y h+ o(h)\,. $$

La probabilité qu'aucun de ces événements n'ait lieu entre $t$ et $t+h$ vaut

$$ 1-n_x \left(1+\frac{n_x}{\gamma K} +\frac{n_y}{K+n_x} \right) h - \beta n_y \left(\frac{n_y}{K+n_x} +\alpha \right) h + o(h). $$

Enfin, la probabilité qu'il y ait tout autre événement (comme la naissance de deux proies par exemple) vaut $o(h)$.
Voici une expérience numérique où vous pouvez comparer, en fonction de $K$, les trajectoires déterministes (en noir) du modèle de Rosenzweig-McArhtur avec celles du processus markovien de saut (en rouge).

Contentons-nous ici d'énoncer le théorème de Kurtz qui est le résultat de base pour comprendre ce qui se passe. Soit $T>0$ un horizon temporel. Supposons $N_x(0)$ et $N_y(0)$ donnés. Alors, pour tout $\varepsilon>0$

$$ \lim_{K\to+\infty} \mathbb{P}\left( \sup_{0\leq t\leq T}\mathrm{dist}\left(\left(\frac{N_x(t)}{K},\frac{N_y(t)}{K}\right),(x(t),y(t)) \right)\ge \varepsilon\right)=0 $$

où $(x(t),y(t))$ est la solution du modèle de Rosenzweig-McArthur démarrant de $x(0)=N_x(0)/K$, $y(0)=N_y(0)/K$, et $\mathrm{dist}$ est la distance euclidienne. Ajoutons que si on fixe $K$ et qu'on fait tendre $T$ vers l'infini, alors l'extinction a lieu avec probabilité égale à un.

A model of spread of infectious diseases

Chazottes Jean-René , Monticelli Marc — 2017-01-09T09:46:00Z

Publié dans : « Differential Equations. An invitation through embedded visual interactive digital experiments » p 84. 2017

We study one of the simplest models used in epidemics. We divide a given population into three disjoint groups. The population of susceptible individuals is denoted by $S$, the infected population by $I$, and the recovered population by $R$. Of course, each of these is a function of time. We make the following assumptions:

the total population is constant (no births or deaths), so that $\dot{P}=0$, where $P(t)=S(t)+I(t)+R(t)$. We denote by $P_0$ this constant;
the rate of transmission of the disease is proportional to the number of encounters between susceptible and infected individuals. The easiest way to characterize this assumption mathematically is to put $\dot{S}=-\beta SI$ for some constant $\beta>0$ which represents the average number of contacts by an infected individual;
The constants $\nu>0$ and $\mu>0$ represent the recovery rate and rate of loss of immunity, respectively.

We can schematically represent this model by the following diagram.

We get the following equations: $$ \begin{cases} \dot{S}= -\beta SI + \mu R \\ \dot{I}=\beta SI-\nu I\\ \dot{R}=\nu I -\mu R \end{cases} $$

At first glance, this is a three-dimensional system. In fact, it reduces to a two-dimensional one. Indeed, we can for instance eliminate $R$ since $$ R=P_0-S-I. $$ Thus, if we determine $S(t)$ and $I(t)$, we then derive $R(t)$ for free. It results that the initial model boils down to the two-dimensional system $$ \begin{cases} \dot{S}= -\beta SI + \mu (P_0-S-I) \\ \dot{I}=\beta SI-\nu I\\ \end{cases} $$ with four (positive) parameters $\beta,\mu,\nu$ and $P_0$.
We can check easily that there are at most two fixed points, one at $(P_0,0)$ and the other at $$ (\bar{S},\bar{I})=\left( \frac{\nu}{\beta},\frac{\mu\big(P_0-\frac{\nu}{\beta}\big)}{\nu+\mu}\right). $$ The first fixed point clearly corresponds to no disease whatsoever in the population. The second fixed point only exists when $$ P_0\geq \frac{\nu}{\beta}. $$ When $P_0=\nu/\beta$, we have a bifurcation as the two fixed points coalesce at $(P_0,0)$. The number $\nu/\beta$ is called the threshold level for the disease.
We can also find the $I$-nullclines which are given by $I=0$ and $S=\nu/\beta$, and the $S$-nullcline which is given by the graph of the function $$ I(S)=\frac{\mu(P_0-S)}{\beta S+\mu}. $$

Notice that we are only interested by the triangular region of the plane $(S,I)$ given by $I,S\geq 0$ and $S+I\leq P_0$. Indeed, $S$ and $I$ represent population densities, and $P_0=S+I+R\geq S+I$. Therefore, the model is consistent if this region is invariant, in the sense that if we start in it, we stay there forever. Note that the $I$-axis is not invariant, while on the $S$-axis, solutions increase up to the fixed point $(P_0,0)$. One can check that the vector field along the boundary of this triangular region is either tangent or points inward.

Let's give more details. The Jacobian matrix at a point $(S,I)$ reads $$ \begin{pmatrix} -\beta I-\mu & -\beta P_0-\mu \\ \beta I & \beta S-\nu \end{pmatrix}. $$ When $P_0>\nu/\beta$, the determinant of the Jacobian matrix at $(P_0,0)$ is equal to $\mu(\nu-\beta P_0)<0$. At the other fixed point, one checks that the determinant is equal to $\beta \bar{I}(\beta \bar{S}+\mu)>0$, and the trace is equal to $-\beta\bar{I}-\mu<0$. Hence, $(P_0,0)$ is a saddle, whereas $(\bar{S},\bar{I})$ is a sink.
When $P_0>\nu/\beta$, the only fixed point is $(P_0,0)$. The determinant of the Jacobian matrix at this point is $>0$, and the trace is $<0$. We thus have a sink.

Modèle proie-prédateur en temps discret

Chazottes Jean-René , Monticelli Marc — 2015-07-10T09:41:38Z

$$ \begin{cases} x_{n+1} =r x_n(1-x_n-\alpha y_n) \\ y_{n+1} =y_n(\beta x_n -d) \end{cases} $$

Deux proies & un prédateur : illustration de la sensibilité aux conditions initiales

Chazottes Jean-René , Monticelli Marc — 2014-10-24T11:42:00Z

On considère un modèle où deux populations, l'une de densité $x$, l'autre de densité $y$, sont en compétition (voir cette expérience numérique interactive). Il y a en plus une troisième population, de densité $z$, qui exerce sa prédation sur ces deux populations.

Pour certaines valeurs des paramètres, la dynamique de ces trois populations conduit à un attracteur complexe, du même genre que celui du modèle de Lorenz. L'attracteur n'est ni un point, ni un cycle limite mais un ensemble fractal.
Voici les équations différentielles de ce modèle :

$$ \begin{cases} \dot{x} & = x(1-x-y-10z)\\ \dot{y} & = y(1-1.5x-y-z)\\ \dot{z} & = z(-1+5x+0.5y-0.01z) \end{cases} $$

Vous pouvez manipuler l'attracteur en 3D sur cette page.
L'expérience numérique ci-dessous illustre la sensibilité aux conditions initiales. On prend deux conditions initiales proches sur l'attracteur. Pendant un certain temps, les deux trajectoires correspondantes sont proches mais elles se séparent nettement, avant de se rapprocher à nouveau à cause du fait que les trajectoires sont confinées dans l'attracteur (sans jamais se croiser à cause du déterminisme). Et ainsi de suite.
Le phénomène remarquable est qu'on ne peut pas prédire quand les deux trajectoires vont se rapprocher à nouveau : il n'y a pas de périodicité, les mouvements ont l'air erratique alors que le système est déterministe : c'est ce qu'on appelle le chaos déterministe.
Pour plus de clarté, on projette la dynamique dans le plan $(Oxy)$ et on visualise également $x(t)$.

Dans l'expérience suivante, on lance 3000 conditions initiales très proches les unes des autres.

Modèle deux proies - un prédateur

Chazottes Jean-René , Monticelli Marc — 2014-09-25T15:46:51Z

MAJ 3/06/2015 : nécessite WebGL

Pour certaines valeurs des paramètres, la dynamique de ces trois populations conduit à un attracteur complexe, du même genre que celui du modèle de Lorenz. Voici les équations :

$$ \begin{cases} \dot{x} & = x(1-x-y-10z)\\ \dot{y} & = y(1-1.5x-y-z)\\ \dot{z} & = z(-1+5x+0.5y-0.01z) \end{cases} $$

Vous pouvez manipuler en 3D l'attracteur :

Le modèle proie-prédateur de Volterra avec compétition entre proies

Chazottes Jean-René , Monticelli Marc — 2013-10-11T15:21:31Z

Le modèle que vous pouvez expérimenter numériquement ci-après est une variante du modèle proie-prédateur de Loka-Volterra. Il s'écrit

$$ \begin{cases} \dot{x}=x(1-x-y)\\ \dot{y}=\beta(x-\alpha)y \end{cases} $$

où $\alpha$ et $\beta$ sont des paramètres positifs. La différence avec le modèle proie-prédateur de Loka-Volterra est qu'en l'absence de prédateurs, la densité des proies suit l'équation $\dot{x}=x(1-x)$ au lieu de l'équation $\dot{x}=x$. La conséquence est que $x(t)$ ne peut plus tendre vers l'infini mais tend vers $1$, comme le montre l'expérience numérique ci-dessous.

ENI équation logistique à mettre.

Le lecteur peut expérimenter numériquement le modèle en manipulant l'expérience numérique interactive ci-dessous et constater que les solutions se comportent très différemment de celles du modèle proie-prédateur de Loka-Volterra.

On constate que l'équilibre $(x^*,y^*)=(\alpha,1-\alpha)$, qui correspond à la coexistence entre les proies et les prédateurs, attire toutes les solutions dont les conditions initiales sont des densités de proies et de prédateurs positives. Cela s'interprète comme un retour à l'équilibre si on a perturbé le système. On parle de coexistence stable.

On voit aussi que cet équilibre n'existe pas si $\alpha$ devient supérieur à $1$ car ce n'est plus un point du quadrant positif et qu'il n'a donc plus de sens biologique. Quand $\alpha$ vaut $1$, cet équilibre se confond avec l'équilibre $(1,0)$ qui est celui vers lequel tend la densité de proies s'il n'y a pas de prédateurs.

La conclusion est qu'en prenant en compte la compétition entre proies, qui empêche qu'elles prolifèrent sans borne en l'absence de prédateurs, on a complètement détruit le caractère oscillatoire périodique des solutions du modèle proie-prédateur de Loka-Volterra. Pour retrouver de telles oscillations, il faut en fait modifier le terme d'interaction : c'est l'objet du modèle de Rosenzweig-McArthur.

Code HTML pour intégrer cette simulation dans vos pages :