Populations isolées

Le modèle de Malthus

Nous considérons le modèle proposé par Malthus (1798) : \dot x = (n-m)x, \tag{1} avec n le taux de natalité, et m le taux de mortalité.

Préliminaires

Pour simuler ce modèle, c’est à dire intégrer numériquement les solutions au problème de Cauchy correspondant à l’Equation 1 avec x(0)=x_0\geq0, nous allons utiliser les routines de DifferentialEquations.jl avec la méthode par défaut. Nous utiliserons aussi le package Plots.jl pour les représentations graphiques¹.

¹ si vous n’avez jamais installé un package, il faudra commencer par le faire avant de pouvoir l’utiliser…

using DifferentialEquations
using Plots

Conditions initiales, paramètres et temps

Nous définissions la condition initiale de la simulation :

x0 = 0.1

Nous définissons les paramètres, et les encapsulons dans un vecteur de paramètres :

n = 3.0     # natalité
m = 2.0     # mortalité
par_malthus = [n, m] # packing in an array

Enfin, les propriétés du temps d’intégration :

tspan = (0.0, 10.0)
tstep = 0.1

Système dynamique

On définit le modèle comme une fonction renvoyant la dérivée de x en fonction de l’état (u est la notation pour DifferentialEquations.jl), de paramètres (p), et du temps (t). L’ordre des arguments est important, c’est sous cette forme que les routines d’intégration attendent le modèle.

function malthus(u, p, t)
    n, m = p                # unpacking
    x = u                # use x notation
    return dx = (n-m)x      # return derivative
end

Note

L’ordre des arguments u, p, t est important pour le solver de DifferentialEquations.jl.

Intégration

On commence par définir le problème de Cauchy à intégrer, comme un ODEProblem avec arguments: le modèle, la condition initiale, les bornes d’intégration, les paramètres ainsi que l’option saveat = tstep permettant de récupérer la solution tout les tstep pas de temps (il y a plein d’options de sortie de l’intégration, dont même une fonction du temps !)

prob_malthus = ODEProblem(
    malthus,          # modèle
    x0,               # condition initiale
    tspan,            # tspan
    par_malthus;      # paramètres
    saveat = tstep,   # keyword argument, option de sortie
)

On intègre le modèle via solve, défini par DifferentialEquations.jl pour des struct de type ODEProblem :

sol_malthus = solve(prob_malthus)

Le type de solution renvoyé par le solveur est assez complexe et comprend de nombreux champs informatifs sur le calcul. On peut accéder au temps de simulation via sol_malthus.t :

first(sol_malthus.t, 3)

3-element Vector{Float64}:
 0.0
 0.1
 0.2

ainsi qu’aux valeurs de la variable x calculées le long du temps via sol_matlhus.u :

first(sol_malthus.u, 3)

3-element Vector{Float64}:
 0.1
 0.1105170918098962
 0.12214028021690636

Même si ce n’est pas indispensable, il est possible de transformer la solution renvoyée facilement en DataFrame, qui peut permettre des manipulations plus faciles.

using DataFrames

sol_malthus = DataFrame(sol_malthus)
rename!(sol_malthus, :timestamp => :time, :value => :x)

Note

remarquez la fonction rename! qui modifie en place le dataframe. Par convention les fonctions dont le nom finit par ! modifient leur argument en place.

si bien que :

first(sol_malthus, 3)

3×2 DataFrame

Row	time	x
	Float64	Float64
1	0.0	0.1
2	0.1	0.110517
3	0.2	0.12214

Représentation graphique

On peut représenter graphiquement la simulation de la croissance de la population au cours du temps (ici via le dataframe).

plot(
    sol_malthus.time,                   # abscisses
    sol_malthus.x;                      # ordonnées
    palette = :tab10,                   # palette de couleurs
    linewidth = 2,
    title = "\n Modèle de Malthus \$n=$n, m=$m\$",
    label = "population \$x\$",
    ylabel = "densité de population \$x(t)\$",
    xlabel = "temps \$t\$",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

Note

Les titres et autres chaînes de caractères peuvent utiliser des formules LaTeX via \$, e.g. "\$x\$".

La notation "$m" accède à la valeur du paramètre m et la renvoit dans la chaine de caractères.

Le modèle logistique

Nous considérons ici le modèle “logistique” proposé par Verhulst (1838) :

\dot x = r x \left(1-\frac{x}{K}\right), \tag{2} avec r le taux de croissance intrinsèque de la population et K la capacité de charge de l’environnement.

Il n’y a pas de difficulté particulière par rapport aux simulations précedentes. On va réutiliser la condition initiale et les paramètres de temps définis précédemment.

Il faut cependant définir les paramètres du modèle Equation 2 :

Code

r = 1.0      # natalité
K = 10.0     # mortalité
par_logistic = [r, K] # packing

ainsi que le nouveau système dynamique :

Code

function logistic(u, p, t)
    r, K = p            # unpacking
    x = u            # use x notation
    return dx = r*x*(1-x/K)    # return derivative
end

et le problème de Cauchy correspondant :

Code

prob_logistic = ODEProblem(
    logistic,         # modèle
    x0,               # condition initiale
    tspan,            # tspan
    par_logistic;     # paramètres
    saveat = tstep,   # option de sortie
)

On simule et on transforme la solution en dataframe :

Code

sol_logistic = solve(prob_logistic)

sol_logistic = DataFrame(sol_logistic)
rename!(sol_logistic, :timestamp => :time, :value => :x)

On trace la solution, en rajoutant les équilibres stable (x=K) et instable (x=0) :

# solution
plot(
    sol_logistic.time,
    sol_logistic.x;
    palette = :tab10,
    linewidth = 2,
    title = "\n Modèle logistique \$r=$r, K=$K\$",
    label = "population \$x\$",
    legend = :right,
    ylabel = "densité de population \$x(t)\$",
    xlabel = "temps \$t\$",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

# équilibre 0
plot!(
    sol_logistic.time,
    zeros(length(sol_logistic.time));
    color = "red",
    linewidth = 2,
    linestyle = :dash,
    linealpha = .5,
    label = "équilibre instable",
)

# équilibre K
plot!(
    sol_logistic.time,
    ones(length(sol_logistic.time)).*K;
    color = "green",
    linewidth = 2,
    linestyle = :dash,
    linealpha = .5,
    label = "équilibre stable",
)

Note

Remarquez l’utilisation de plot! qui modifie en place la première figure en rajoutant des éléments (ici les équilibres et leurs labels).

Effets Allee

On s’intéresse à un modèle de dynamique de population avec “effets Allee forts”, souvent attribué à Gruntfest, Arditi, and Dombrovsky (1997)² :

² mais de nombreuses variations de cette forme polynomiale existent dans la littérature depuis Bazykin (1985)

\dot x = r x \left(\frac{x}{\epsilon}-1\right)\left(1-\frac{x}{K}\right), \tag{3} avec r le taux de croissance intrinsèque de la population (par analogie avec la logistique), K la capacité de charge de l’environnement et \epsilon le seuil en dessous duquel la population n’est pas viable (‘seuil de Allee’).

Nous souhaitons illustrer la bistabilité caractéristique du comportement de ce modèle à partir de la simulation depuis plusieurs conditions initiales.

Commençons par définir le modèle et les paramètres (en conservant les r et K précédents):

Code

epsilon = 2                 # seuil de Allee
par_allee = [r, K, epsilon] # packing
tspan_allee = (0.0, 3.0)

function allee(u, p, t)
    r, K, epsilon = p        # unpacking
    x = u                 # use x notation
    return dx = r*x*(x/epsilon - 1)*(1 - x/K)     # derivative
end

Nous définissons une fonction qui réalise une simulation en un seul appel depuis une condition initiale x_0 :

function int_allee(x0, tspan = tspan_allee, param = par_allee)
    prob_allee = ODEProblem(
        allee,            # modèle
        x0,               # condition initiale
        tspan,            # tspan
        param;            # paramètres
        saveat = tstep,   # option de sortie
        )

    sol_allee = solve(prob_allee)
    sol_allee = DataFrame(sol_allee)
    rename!(sol_allee, :timestamp => :time, :value => :x)

    return sol_allee
end

Note

La fonction int_allee comporte un nombre d’arguments variable: int_allee(x0) renverra la solution avec les arguments tspan et param à leurs valeurs par défaut (ce sont des vararg); c’est le même appel que int_allee(x0, tspan_allee, par_allee).

Nous définissons un vecteur de conditions initiales différentes :

x0step = 1.35
x0vec = x0step:x0step:K     # range de valeur

Finalement on réalise la figure. La stratégie diffère un peu de ce qui a été vu ci-dessus. Nous commençons par initier un graphique fig, et faisons une boucle for pour tracer chacune des simulations correspondant aux différentes conditions initiales³ :

³ la palette de couleur est accessoire, juste pour l’esthétique

# custom color palette
init_cgrad = palette([:steelblue, :lightblue], length(x0vec))

# initialisation du graphique
fig = plot(;          # on initie uniquement le graphique
    palette = init_cgrad,
    title = "\n Modèle à effets Allee \$r=$r, K=$K\$, \$ϵ=$epsilon\$",
    legend = :right,
    ylabel = "densité de population \$x(t)\$",
    xlabel = "temps \$t\$",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
    )

# boucle de plot avec intégration pour differentes conditions initiales
for x0 in x0vec         # x0 parcourt l'array x0vec
    plot!(
        fig,              # on modifie fig
        int_allee(x0).time, # abscisses
        int_allee(x0).x;    # ordonnées
        linewidth = 2,      # keyword arguments après ;
        label = "",
    )
end

# équilibre 0
plot!(
    fig,
    int_allee(0).time,
    zeros(length(int_allee(0).time));
    color = "green",
    linewidth = 2,
    linestyle = :dash,
    linealpha = .5,
    label = "équilibre stable",
)

# équilibre epsilon
plot!(
    fig,
    int_allee(0).time,
    ones(length(int_allee(0).time)).*epsilon;
    color = "red",
    linewidth = 2,
    linestyle = :dash,
    linealpha = .5,
    label = "équilibre instable",
)

# équilibre K
plot!(
    fig,
    int_allee(0).time,
    ones(length(int_allee(0).time)).*K;
    color = "green",
    linewidth = 2,
    linestyle = :dash,
    linealpha = .5,
    label = "",
)

display(fig)      # actually shows the plot P

On peut sauvegarder la figure dans différents formats (e.g. png, pdf) :

savefig(fig, "dyn_allee.pdf")
savefig(fig, "dyn_allee.png")

Passons aux populations exploitées.

References

Bazykin, A. D. 1985. Mathematical Biophysics of Interacting Populations. Nauka, Moscow (in Russian).

Gruntfest, Y., R. Arditi, and Y. Dombrovsky. 1997. “A Fragmented Population in a Varying Environment.” Journal of Theoretical Biology 185: 539–47.

Malthus, T. R. 1798. An Essay on the Principle of Population. J. Johnson.

Verhulst, P. F. 1838. “Notice Sur La Loi Que La Population Poursuit Dans Son Accroissement.” Correspondance Mathématique Et Physique 10: 113–21.

Reuse

CC BY-NC 4.0

--- title: "Populations isolées" engine: julia ---  ```{julia} #| include: false #| eval: true # # this cell does not appear in the rendering, but is executed # # for reproducibility purposes, we load the local julia environment/project, containing # [13f3f980] CairoMakie v0.11.6 # [a93c6f00] DataFrames v1.6.1 # [864edb3b] DataStructures v0.18.16 # [0c46a032] DifferentialEquations v7.12.0 # [e9467ef8] GLMakie v0.9.6 # [91a5bcdd] Plots v1.40.0 # [f27b6e38] Polynomials v4.0.6 # [90137ffa] StaticArrays v1.9.1 # [0c5d862f] Symbolics v5.16.1 # # to share the environment, copy Project.toml and Manifest.toml files in some directory # `activate` the local environment # if necessary `instantiate` to get the correct package versions # # to check if some PackageName.jl is used from the local environment ## Pkg.status("PackageName") using Pkg Pkg.activate(".") ``` ## Le modèle de Malthus {#sec-malthus} Nous considérons le modèle proposé par @malthus : $$ \dot x = (n-m)x, $$ {#eq-malthus} avec $n$ le taux de natalité, et $m$ le taux de mortalité. ### Préliminaires Pour simuler ce modèle, c'est à dire intégrer numériquement les solutions au problème de Cauchy correspondant à l'@eq-malthus avec $x(0)=x_0\geq0$, nous allons utiliser les routines de `DifferentialEquations.jl` avec la méthode par défaut. Nous utiliserons aussi le package `Plots.jl` pour les représentations graphiques^[si vous n'avez jamais installé un package, il faudra commencer par le faire avant de pouvoir l'utiliser...]. ```{julia} using DifferentialEquations using Plots ``` ### Conditions initiales, paramètres et temps Nous définissions la condition initiale de la simulation : ```{julia} x0 = 0.1 ``` Nous définissons les paramètres, et les encapsulons dans un vecteur de paramètres : ```{julia} n = 3.0 # natalité m = 2.0 # mortalité par_malthus = [n, m] # packing in an array ``` Enfin, les propriétés du temps d'intégration : ```{julia} tspan = (0.0, 10.0) tstep = 0.1 ``` ### Système dynamique On définit le modèle comme une fonction renvoyant la dérivée de $x$ en fonction de l'état (`u` est la notation pour `DifferentialEquations.jl`), de paramètres (`p`), et du temps (`t`). L'ordre des arguments est important, c'est sous cette forme que les routines d'intégration attendent le modèle. ```{julia} function malthus(u, p, t) n, m = p # unpacking x = u # use x notation return dx = (n-m)x # return derivative end ``` ::: {.callout-note} L'ordre des arguments `u, p, t` est important pour le solver de `DifferentialEquations.jl`. ::: ### Intégration On commence par définir le problème de Cauchy à intégrer, comme un `ODEProblem` avec arguments: le modèle, la condition initiale, les bornes d'intégration, les paramètres ainsi que l'option `saveat = tstep` permettant de récupérer la solution tout les `tstep` pas de temps (il y a plein d'options de sortie de l'intégration, dont même une fonction du temps !) ```{julia} prob_malthus = ODEProblem( malthus, # modèle x0, # condition initiale tspan, # tspan par_malthus; # paramètres saveat = tstep, # keyword argument, option de sortie ) ``` On intègre le modèle via `solve`, défini par `DifferentialEquations.jl` pour des `struct` de type `ODEProblem` : ```{julia} sol_malthus = solve(prob_malthus) ``` Le type de solution renvoyé par le solveur est assez complexe et comprend de nombreux champs informatifs sur le calcul. On peut accéder au temps de simulation via `sol_malthus.t` : ```{julia} #| output: true first(sol_malthus.t, 3) ``` ainsi qu'aux valeurs de la variable $x$ calculées le long du temps via `sol_matlhus.u` : ```{julia} #| output: true first(sol_malthus.u, 3) ``` Même si ce n'est pas indispensable, il est possible de transformer la solution renvoyée facilement en `DataFrame`, qui peut permettre des manipulations plus faciles. ```{julia} using DataFrames sol_malthus = DataFrame(sol_malthus) rename!(sol_malthus, :timestamp => :time, :value => :x) ``` ::: {.callout-note} remarquez la fonction `rename!` qui modifie en place le dataframe. Par convention les fonctions dont le nom finit par `!` modifient leur argument en place. ::: si bien que : ```{julia} #| output: true first(sol_malthus, 3) ``` ### Représentation graphique On peut représenter graphiquement la simulation de la croissance de la population au cours du temps (ici via le dataframe). ```{julia} #| output: true plot( sol_malthus.time, # abscisses sol_malthus.x; # ordonnées palette = :tab10, # palette de couleurs linewidth = 2, title = "\n Modèle de Malthus \$n=$n, m=$m\$", label = "population \$x\$", ylabel = "densité de population \$x(t)\$", xlabel = "temps \$t\$", margin = .5Plots.cm, topmargin = 1Plots.cm, ) ``` ::: {.callout-note} Les titres et autres chaînes de caractères peuvent utiliser des formules LaTeX via `\$`, e.g. `"\$x\$"`. La notation `"$m"` accède à la valeur du paramètre `m` et la renvoit dans la chaine de caractères. ::: ## Le modèle logistique {#sec-logistic} Nous considérons ici le modèle "logistique" proposé par @verhulst1838 : $$ \dot x = r x \left(1-\frac{x}{K}\right), $$ {#eq-logistic} avec $r$ le taux de croissance intrinsèque de la population et $K$ la capacité de charge de l'environnement. Il n'y a pas de difficulté particulière par rapport aux simulations précedentes. On va réutiliser la condition initiale et les paramètres de temps définis précédemment. Il faut cependant définir les paramètres du modèle @eq-logistic : ```{julia} #| code-fold: true r = 1.0 # natalité K = 10.0 # mortalité par_logistic = [r, K] # packing ``` ainsi que le nouveau système dynamique : ```{julia} #| code-fold: true function logistic(u, p, t) r, K = p # unpacking x = u # use x notation return dx = r*x*(1-x/K) # return derivative end ``` et le problème de Cauchy correspondant : ```{julia} #| code-fold: true prob_logistic = ODEProblem( logistic, # modèle x0, # condition initiale tspan, # tspan par_logistic; # paramètres saveat = tstep, # option de sortie ) ``` On simule et on transforme la solution en dataframe : ```{julia} #| code-fold: true sol_logistic = solve(prob_logistic) sol_logistic = DataFrame(sol_logistic) rename!(sol_logistic, :timestamp => :time, :value => :x) ``` On trace la solution, en rajoutant les équilibres stable ($x=K$) et instable ($x=0$) : ```{julia} #| code-fold: false #| output: true # solution plot( sol_logistic.time, sol_logistic.x; palette = :tab10, linewidth = 2, title = "\n Modèle logistique \$r=$r, K=$K\$", label = "population \$x\$", legend = :right, ylabel = "densité de population \$x(t)\$", xlabel = "temps \$t\$", margin = .5Plots.cm, topmargin = 1Plots.cm, ) # équilibre 0 plot!( sol_logistic.time, zeros(length(sol_logistic.time)); color = "red", linewidth = 2, linestyle = :dash, linealpha = .5, label = "équilibre instable", ) # équilibre K plot!( sol_logistic.time, ones(length(sol_logistic.time)).*K; color = "green", linewidth = 2, linestyle = :dash, linealpha = .5, label = "équilibre stable", ) ``` ::: {.callout-note} Remarquez l'utilisation de `plot!` qui modifie en place la première figure en rajoutant des éléments (ici les équilibres et leurs labels). ::: ## Effets Allee {#sec-allee} On s'intéresse à un modèle de dynamique de population avec "effets Allee forts", souvent attribué à @gruntfest1997^[mais de nombreuses variations de cette forme polynomiale existent dans la littérature depuis @bazykin1985] : $$ \dot x = r x \left(\frac{x}{\epsilon}-1\right)\left(1-\frac{x}{K}\right), $$ {#eq-allee} avec $r$ le taux de croissance intrinsèque de la population (par analogie avec la logistique), $K$ la capacité de charge de l'environnement et $\epsilon$ le seuil en dessous duquel la population n'est pas viable ('seuil de Allee'). Nous souhaitons illustrer la bistabilité caractéristique du comportement de ce modèle à partir de la simulation depuis plusieurs conditions initiales. Commençons par définir le modèle et les paramètres (en conservant les $r$ et $K$ précédents): ```{julia} #| code-fold: true epsilon = 2 # seuil de Allee par_allee = [r, K, epsilon] # packing tspan_allee = (0.0, 3.0) function allee(u, p, t) r, K, epsilon = p # unpacking x = u # use x notation return dx = r*x*(x/epsilon - 1)*(1 - x/K) # derivative end ``` Nous définissons une fonction qui réalise une simulation en un seul appel depuis une condition initiale $x_0$ : ```{julia} function int_allee(x0, tspan = tspan_allee, param = par_allee) prob_allee = ODEProblem( allee, # modèle x0, # condition initiale tspan, # tspan param; # paramètres saveat = tstep, # option de sortie ) sol_allee = solve(prob_allee) sol_allee = DataFrame(sol_allee) rename!(sol_allee, :timestamp => :time, :value => :x) return sol_allee end ``` ::: {.callout-note} La fonction `int_allee` comporte un nombre d'arguments variable: `int_allee(x0)` renverra la solution avec les arguments `tspan` et `param` à leurs valeurs par défaut (ce sont des *vararg*); c'est le même appel que `int_allee(x0, tspan_allee, par_allee)`. ::: Nous définissons un vecteur de conditions initiales différentes : ```{julia} x0step = 1.35 x0vec = x0step:x0step:K # range de valeur ``` Finalement on réalise la figure. La stratégie diffère un peu de ce qui a été vu ci-dessus. Nous commençons par initier un graphique `fig`, et faisons une boucle `for` pour tracer chacune des simulations correspondant aux différentes conditions initiales^[la palette de couleur est accessoire, juste pour l'esthétique] : ```{julia} #| output: true # custom color palette init_cgrad = palette([:steelblue, :lightblue], length(x0vec)) # initialisation du graphique fig = plot(; # on initie uniquement le graphique palette = init_cgrad, title = "\n Modèle à effets Allee \$r=$r, K=$K\$, \$ϵ=$epsilon\$", legend = :right, ylabel = "densité de population \$x(t)\$", xlabel = "temps \$t\$", margin = .5Plots.cm, topmargin = 1Plots.cm, ) # boucle de plot avec intégration pour differentes conditions initiales for x0 in x0vec # x0 parcourt l'array x0vec plot!( fig, # on modifie fig int_allee(x0).time, # abscisses int_allee(x0).x; # ordonnées linewidth = 2, # keyword arguments après ; label = "", ) end # équilibre 0 plot!( fig, int_allee(0).time, zeros(length(int_allee(0).time)); color = "green", linewidth = 2, linestyle = :dash, linealpha = .5, label = "équilibre stable", ) # équilibre epsilon plot!( fig, int_allee(0).time, ones(length(int_allee(0).time)).*epsilon; color = "red", linewidth = 2, linestyle = :dash, linealpha = .5, label = "équilibre instable", ) # équilibre K plot!( fig, int_allee(0).time, ones(length(int_allee(0).time)).*K; color = "green", linewidth = 2, linestyle = :dash, linealpha = .5, label = "", ) display(fig) # actually shows the plot P ``` \ On peut sauvegarder la figure dans différents formats (e.g. png, pdf) : ```{julia} savefig(fig, "dyn_allee.pdf") savefig(fig, "dyn_allee.png") ``` Passons aux [populations exploitées](pop_exploitees.qmd).