Populations en interaction (1)

La tordeuse du bourgeon de l’épinette (suite)

Pour débuter cette partie sur les populations en interactions, nous reprenons le modèle de tordeuse de Ludwig, Jones, and Holling (1978) en supposant que la population d’oiseaux réagit (lentement) à la démographie des tordeuses, via la prédation.

Le changement principal ici est la dimension du modèle (dimension 2) : les tailles de populations de tordeuses x et d’oiseaux y varient toutes deux au cours du temps en s’influençant l’une l’autre, avec une population d’oiseaux qui varie lentement (d’où le paramètre \varepsilon supposé petit) par rapport à celle des tordeuses.

Le modèle prend la forme : \left\{ \begin{array}{l} \displaystyle \dot x = rx\left(1-\frac{x}{K}\right) - \frac{\alpha x^2}{h^2+x^2}\ y \\[.3cm] \displaystyle \dot y = \varepsilon \left(\frac{n \alpha x^2}{h^2+x^2}\ y -m y\right) \end{array} \right. \tag{1}

Il faut ajuster un peu la manière de coder pour prendre en compte ces deux dimensions. Commençons par les paramètres divers.

using DifferentialEquations
using Plots
using DataFrames

# paramètres
r = 5.0      # natalité
K = 10.0     # mortalité
α = 1.0      # taux max de prédation
h = 0.5      # constante de demi-saturation

ϵ = 0.01     # timescale, \epsilon + TAB
n = 5.0      # gain à la prédation
m = 3.0      # mortalité

par_tord_ois = [r, K, α, h, ϵ, n, m]

# temps d'intégration
tspan = (0.0, 400.0)
tstep = 0.02

# conditions initiales
x0 = 1.0    # tordeuses
y0 = 2.5    # oiseaux
etat0 = [x0, y0]

Puis le modèle :

function tord_ois(u, p, t)
    r, K, α, h, ϵ, n, m = p
    x = u[1]
    y = u[2]

    dx = r*x*(1 - x/K) - α*x^2/(h^2 + x^2)*y
    dy = ϵ*(n*α*x^2/(h^2 + x^2)*y - m*y)

    return [dx, dy]
end

Problème d’intégration et simulation :

prob_tord_ois = ODEProblem(
      tord_ois,
      etat0,
      tspan,
      par_tord_ois;
      saveat = tstep,
)

sol_tord_ois = solve(prob_tord_ois, reltol = 1e-6)

sol_tord_ois = DataFrame(sol_tord_ois)
rename!(sol_tord_ois, :timestamp => :time, :value1 => :x, :value2 => :y)

Représentation graphique contre le temps

# color definitions
mygreen = palette(:tab10)[3]
myorange = palette(:tab10)[2]
myblue = palette(:tab10)[1]
myred = palette(:tab10)[4]

# initialisation graphique et solution simulée
fig1 = plot(
    sol_tord_ois.time,
    sol_tord_ois.x;
    color = myblue,
    linewidth = 2,
    label = "tordeuses \$x\$",
    xlabel = "temps",
    ylabel = "densités de populations",
    title = "Dynamiques des tordeuses avec\n population d'oiseaux variable",
)

plot!(
    fig1,
    sol_tord_ois.time,
    sol_tord_ois.y ./ 2,  # scale by 2 for aesthetics
    color = myorange,
    linewidth = 2,
    label = "oiseaux \$y/2\$",
)

On observe ici des bifurcations dynamiques avec le passage de la population de tordeuse d’une branche d’équilibre stable à l’autre, qui créé un comportement de type cycle d’hysteresis.

La situation se comprend bien sur le diagramme de bifurcations (y, x). On retrace le diagramme de bifurcations (cf. la page sur les populations exploitées).

Code

using Symbolics
using Polynomials

@variables X
D = Differential(X)

# lieu des équilibres positifs
Y = r/(α*X)*(1-X/K)*(h^2+X^2)

# dénominateur de la dérivée
dnumerator = Symbolics.arguments(Symbolics.value(simplify(expand_derivatives(D(Y)))))[1]

# on récupère les coefficients X^k du polynôme
coefs_dict = Symbolics.value(dnumerator).dict
dd = Dict(Symbolics.degree(first(kv)) => kv[2] for kv ∈ coefs_dict)
# on rajoute dans le dictionnaire le coefficient constant
dd[0] = substitute(dnumerator, Dict(X=>0))

# on définit le polynôme à partir du dictionnaire
dnumpoly = SparsePolynomial(dd, :X)

# on calcule les racines en filtrant les racines entre 0 et K via une fonction anonyme s-> K > s > 0
droots = filter(s -> K > s > 0, roots(dnumpoly))

# vecteur pour le tracé du diagramme de bifurcation
xplot1 = 0.08:.01:droots[1]
xplot2 = droots[1]:.01:droots[2]
xplot3 = droots[2]:.01:K
yeq1 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot1]
yeq2 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot2]
yeq3 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot3]

# diagramme de bifurcations
figbif = plot(
    yeq1,
    xplot1;
    linewidth = 2,
    color = mygreen,
    label = "équilibres stables",
    legend = :left,
    xlabel = "population d'oiseaux \$y\$",
    ylabel = "population de tordeuses \$x\$",
    title = "Diagramme de bifurcations pour le modèle de tordeuses",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

plot!(
    figbif,
    yeq2,
    xplot2;
    linewidth = 2,
    color = myred,
    label = "équilibres instables",
)

plot!(
    figbif,
    yeq3,
    xplot3;
    linewidth = 2,
    color = mygreen,
    label = "",
)

plot!(
    figbif,
    [0, maximum(yeq1)],
    [0, 0];
    color = myred,
    lw = 2,
    label = "",
)

On trace sur ce diagramme la trajectoire simulée plus haut :

plot!(
    figbif,
    sol_tord_ois.y,
    sol_tord_ois.x;
    color = myblue,
    linewidth = 2,
    linealpha = 0.5,
    label = "trajectoire",
)

Le modèle proie-prédateur de Lotka et Volterra

Nous considérons le modèle de dynamique de populations de Lotka (1925) et Volterra (1926) :

\left\{\begin{array}{l} \dot x = rx - c xy,\\ \dot y = bxy - m y. \end{array}\right. \tag{2}

Avec x la population de proies et y la population de prédateurs.

Dynamiques

Il n’y a pas de difficulté particulière à la simulation par rapport au modèle de la tordeuse du bourgeon de l’épinette avec population d’oiseaux variables.

Code

# conditions initiales
x0 = 1.0
y0 = 1.95
etat0 = [x0, y0]

# paramètres
r = 1.0
c = 1.0
b = 1.0
m = 1.0
par_lovo = [r, c, b, m]

# temps d'integration
tspan = (0.0, 30.0)
tstep = .01

# définition du modèle
function lovo(u, par, t)
    r, c, b, m = par
    x = u[1]
    y = u[2]
    dx = r*x - c*x*y
    dy = b*x*y - m*y
    return [dx, dy]
end

# problème
prob_lovo = ODEProblem(lovo, etat0, tspan, par_lovo, saveat = tstep)

# intégration
sol_lovo = solve(prob_lovo, reltol = 1e-6)

# dataframe
sol_lovo = DataFrame(sol_lovo)
rename!(sol_lovo, :timestamp => :time, :value1 => :x, :value2 => :y)

figlv = plot(
    sol_lovo.time,
    sol_lovo.x;
    linewidth = 2,
    color = myblue,
    label = "proies",
    xlabel = "temps",
    ylabel = "densité de populations",
    title = "Modèle de Lotka Volterra",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

plot!(
    figlv,
    sol_lovo.time,
    sol_lovo.y;
    linewidth = 2,
    color = myorange,
    label = "prédateurs",
)

Espace d’état

Commençons par tracer les isoclines nulles ainsi que l’orientation du champs de vecteur dans l’espace d’état. Ce dernier utilise la fonction quiver().

Caution

Ce qui suit fonctionne avec Plots.jl, néanmoins la façon de tracer le champs de vecteurs semble vraiment peu naturelle. La librairie graphique Makie que nous verrons un peu plus loin permet de faire ces représentations bien plus facilement, comme nous le voyons ici.

# fonctions qui renvoient les composantes de la dérivée mises à l'échelle
scale = 10
der_x(x, y) = lovo([x y], par_lovo, 0)[1]/scale
der_y(x, y) = lovo([x y], par_lovo, 0)[2]/scale

# valeur de x et y formant une grille sur laquelle évaluer le champs de vecteurs
xrange = range(0., 2, length=11)'  # note the quote : broadcast pour créer une grille via der_x.() et der_y.()
yrange = range(0. ,2, length=11)

# champs de vecteurs
figplan = quiver(
    xrange,
    yrange,
    quiver = (der_x.(xrange, yrange), der_y.(xrange, yrange));
    ylim = (-0.05, 2),
    xlim = (-0.05, 2),
    color = :lightgray,
    grid = false,
    title = "Modèle de Lotka Volterra",
    xlabel = "proies",
    ylabel = "prédateurs",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

On rajoute les isoclines nulles et les équilibres :

Code

# isoclines nulles
xplot = 0:2
yplot = 0:2

# dot x = 0
plot!(
    figplan,
    xplot,
    ones(length(xplot)).*r./c;
    linewidth = 2,
    color = mygreen,
    label = "\$\\dot x = 0\$",
)
plot!(
    figplan,
    zeros(length(yplot)),
    yplot;
    linewidth = 2,
    color = mygreen,
    label = "",
)

# dot y = 0
plot!(
    figplan,
    ones(length(yplot)).*m./b,
    yplot;
    linewidth = 2,
    color = myred,
    label = "\$\\dot y = 0\$",
)
plot!(
    figplan,
    xplot,
    zeros(length(xplot));
    linewidth = 2,
    color = myred,
    label = "",
)

# équilibres
plot!(
    figplan,
    (0, 0);
    markershape = :circle,
    color= myred,
    label = "",
)
plot!(
    figplan,
    (m/b, r/c);
    markershape = :circle,
    color= myorange,
    label = "",
)

Et enfin la trajectoire :

Code

# trajectoire
plot!(
    figplan,
    sol_lovo.x,
    sol_lovo.y;
    color = myblue,
    linewidth = 2,
    label = "trajectoire",
)

display(figplan)

Figure 1: Plan de phase du modèle de Lotka Volterra

Intégrale première

Nous illustrons en 3D l’intégrale première du modèle et superposons la trajectoire simulée plus haut sur ce graphique. L’intégrale première s’écrit : H(x,y)= -r\log(y)+c y-m\log(x) + bx

# l'intégrale première
function int_prem(x, y, par = par_lovo)
    r, c, b, m = par
    return -r*log(y) + c*y - m*log(x) + b*x
end

Pour tracer le graphique en 3D, nous utilisons le backend plotly.

plotly()

# x, y
xsurf = 0.2:0.1:3.0
ysurf = 0.2:0.1:3.0

# calcul de la surface via une compréhension de liste
zsurf = [int_prem(x, y, par_lovo) for x in xsurf, y in ysurf]
# et H(x0, y0)
zplane = [int_prem(x0, y0, par_lovo) for x in xsurf, y in ysurf]

# on trace H(x0, y0)
figsurf = plot(
    xsurf,
    ysurf,
    zplane; # array de taille length(xsurf) x length(ysurf)
    st = :surface,
    color = myorange,
    alpha =.5,
    label = "",
)

# l'intégrale première H(x,y)
plot!(
    figsurf,
    xsurf,
    ysurf,
    zsurf; # array de taille length(xsurf) x length(ysurf)
    st=:surface,
    alpha = .6,
    camera = (30, 10),
    color = :viridis, # colormap
    xlabel = "proies",
    ylabel = "prédateurs",
    title = "Intégrale première du modèle de Lotka Volterra",
)

# la trajectoire
# un simple plot avec 3 coordonnees: x et y de la solution, et z = H(x,y)
plot!(
    figsurf,
    sol_lovo.x,
    sol_lovo.y,
    ones(length(sol_lovo.x)).*int_prem(x0, y0, par_lovo);
    color = myred,
    lw=4,
    label = "",
)

Allons plus loin en étudiant le modèle de Rosenzweig MacArthur, par ici.

References

Lotka, A. J. 1925. Elements of Physical Biology. Williams; Wilkins.

Ludwig, D., D. D. Jones, and C. S. Holling. 1978. “Qualitative Analysis of Insect Outbreak Systems: The Spruce Budworm and Forest.” Journal of Animal Ecology 47: 315–32.

Volterra, V. 1926. “Variazioni e Fluttuazioni Del Numero d’individui in Specie Animali Conviventi.” Memoria Della Reale Accademia Nazionale Dei Lincei 2: 31–113.

Reuse