Respecter la positivité d’un modèle

Modèle positif

En dynamique de populations, les variables des modèles représentent des densités. Un modèle sera bien posé s’il assure la non-négativité de toutes les composantes des trajectoires issues de conditions initiales non-négatives (système positif).

Il est assez facile de vérifier mathématiquement si un modèle est bien positif, en évaluant les dérivées de l’état sur les frontières de l’orthant positif \mathbb{R}^n_+. Cependant, même si le modèle est bien posé, des trajectoires peuvent se rapprocher très fortement des frontières de \mathbb{R}^n_+ et parfois entraîner une instabilité numérique et l’apparition de trajectoires simulées qui ne respectent pas la positivité. Certains schémas numériques permettent de respecter la contrainte de positivité par construction pour certaines classes de modèles (e.g. Blanes, Iserles, and Macnamara (2022)).

Avec DifferentialEquations.jl, il est possible d’ajouter une contrainte d’appartenance à un domaine (comme e.g. \mathbb{R}^n_+) que le résultat de l’intégration par le schéma numérique doit vérifier à chaque itération; si la contrainte n’est pas respectée, alors le pas d’intégration est diminué. Cette méthode est valable pour tous les schémas numériques disponibles.

Modèle de Bauch (2005)

Le Modèle de Bauch (2005) analyse les interactions entre dynamique des épidémies et comportement individuel de vaccination, sur la base d’un modèle de type SIR. Le modèle s’écrit: \left\{\begin{array}{l} \dot S = \mu (1-x) - \beta S I - \mu S,\\ \dot I = \beta S I - \gamma I - \mu I,\\ \dot x = k x (1-x)(\omega I - 1), \end{array}\right. avec S et I les densités d’invidus susceptibles et infectés; les paramètres \mu,~\beta,~\gamma sont associés respectivement aux processus de naissance/mortalité, transmission et guérison.

La variable x représente la proportion d’individus vaccinés à la naissance. Cette proportion suit une dynamique d’imitation, qui tend à augmenter la couverture vaccinale des nouveaux nés lorsque la prévalence de la maladie est plus forte, avec k la vitesse d’ajustement et \omega un paramètre de sensibilité à la dangerosité de la maladie.

Il est facile de vérifier la positivité de ce modèle (\dot z_i\geq 0,~\forall i), néanmoins les simulations présentent une instabilité numérique, susceptible de générer des trajectoires négatives.

Simulations

Simulation sans contrainte particulière sur l’état

Si l’on ne prend pas de précautions, les simulations peuvent générer des trajectoires négatives1.

1 cf. à ce propos la légende de la figure 3 dans Bauch (2005), (2.10) et (2.12) devant y être respectivement comprises comme (2.11) et (2.10).

Nous reprenons ici pour partie les paramètres de Bauch (2005), dans un objectif illustratif.

#
# 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(".")
# paramètres
μ = 1/(365*50)
γ = 0.1
k = 0.001
ω = 6000
β = 10*(μ + γ)
param = [μ, γ, k, ω, β]

# horizon temporel
tmax = 365*200
tspan = (0.0, tmax)
tstep = 0.1

Définissons les équations du modèle:

# modèle
function imit_vacc(etat, param, t)
    μ, γ, k, ω, β = param

    S = etat[1]
    I = etat[2]
    x = etat[3]

    dS = μ*(1-x) - β*S*I - μ*S
    dI = β*S*I - γ*I - μ*I
    dx = k*x*(1-x)*(ω*I-1)

    return [dS, dI, dx]
end

Et intégrons depuis une condition initiale, en utilisant la méthode usuelle de DifferentialEquations.jl (avec le shema numérique Rodas4P()).

using DifferentialEquations

# conditions initiales
S0 = 0.05
I0 = 0.0001
x0 = 0.95
etat0 = [S0, I0, x0]

# probleme aux conditions initiales
tosim = ODEProblem(
    imit_vacc,
    etat0,
    tspan,
    param,
    saveat = 1,
)

# intégration
tol = 1e-10
sol = solve(tosim, Rodas4P(), reltol = tol, abstol = tol)

Une représentation graphique montre que la positivité des variables n’est pas respectée par l’intégration numérique.

Code 
# représentation graphique
using CairoMakie

# figure setup
fig1 = Figure()

ax1 = Axis(
    fig1[1, 1],
    xlabel = "temps",
    ylabel = "densités de populations",
)

ax2 = Axis(
    fig1[1, 2],
    xlabel = "temps",
    ylabel = "taux de vaccination",
)

# plots
lines!(
    ax1,
    sol.t,
    sol[1,:];
    linewidth = 2,
    label = L"$S$",
)

lines!(
    ax1,
    sol.t,
    sol[2,:];
    linewidth = 2,
    label = L"$I$",
)

axislegend(ax1, position = :rc)

lines!(
    ax2,
    sol.t,
    sol[3,:];
    linewidth = 2,
    label = L"$x$",
)

axislegend(ax2, position = :rc)

# titre général en ligne 0
Label(fig1[0, :], "Simulation sans contrainte de positivité", fontsize = 20)

fig1

En fait, on constate que les simulations ne respectent pas le domaine de définition des variables; ici le nombre d’individus infectés I diverge vers -\infty tandis que le nombre d’individus susceptibles S diverge vers +\infty, alors que ces deux variables sont mathématiquement contraintes au domaine [0,1].

Simulation avec contrainte de positivité

Il est possible d’imposer aux solvers de DifferentialEquations.jl de respecter l’appartenance des variables à un domaine spécifique à chaque pas de temps, via un kwarg isoutofdomain passé à solve(). Cet argument est assigné à une fonction qui renvoie true ou false selon que l’état est ou pas dans le domaine souhaité. Si la contrainte n’est pas respectée, le pas de temps d’intégration du solver est réduit.

On commence par définir une fonction qui renvoie false quand l’état a toutes ses composantes positives, et true sinon. Cette fonction est réalisée sur la base de la fonction any() et d’une fonction anonyme qui teste la positivité:

function isnotpositive(u, par, t)
    return any(z -> z < 0, u)
end

Et on passe cet argument pour l’intégration numérique:

sol_constr = solve(
    tosim,
    Rodas4P(),
    reltol = tol,
    abstol = tol,
    isoutofdomain = isnotpositive, # contrainte sur l'état
)

Le résultat de l’intégration (même schéma numérique, mêmes tolérances, même horizon temporel) permet d’obtenir le cycle limite attendu dans ce cas (Bauch (2005)):

Code 
# représentation graphique
fig2 = Figure()

ax1 = Axis(
    fig2[1, 1],
    xlabel = "temps",
    ylabel = "densités de populations",
)

ax2 = Axis(
    fig2[1, 2],
    xlabel = "temps",
    ylabel = "taux de vaccination",
)

lines!(
    ax1,
    sol_constr.t,
    sol_constr[1,:];
    linewidth = 2,
    label = L"$S$",
)

lines!(
    ax1,
    sol_constr.t,
    sol_constr[2,:];
    linewidth = 2,
    label = L"$I$",
)

axislegend(ax1, position = :rc)

lines!(
    ax2,
    sol_constr.t,
    sol_constr[3,:];
    linewidth = 2,
    label = L"$x$",
)

axislegend(ax2, position = :rc)

# titre général en ligne 0
Label(fig2[0, :], "Simulation avec contrainte de positivité", fontsize = 20)

fig2

References

Bauch, C. T. 2005. “Imitation Dynamics Predict Vaccinating Behaviour.” Proceedings of the Royal Society B: Biological Sciences 272 (1573): 1669–75.
Blanes, S., A. Iserles, and S. Macnamara. 2022. “Positivity-Preserving Methods for Ordinary Differential Equations.” ESAIM: Mathematical Modelling and Numerical Analysis 56 (6): 1843–70.

Reuse