Steady state of SIS (suspected-infected-suspected) reaction-diffusion model
Considering the following SIS reaction diffusion model:
\[\left\{\begin{array}{l} S_{t} = d_{S} S_{x x}-\beta(x) \frac{S I}{S+I}+\gamma(x) I=0, \quad 0<x<1 \\ I_{t} = d_{I} I_{x x}+\beta(x) \frac{S I}{S+I}-\gamma(x) I=0, \quad 0<x<1 \\ S_{x}=I_{x}=0, \quad x=0,1, \end{array}\right.\]
where $\int_{0}^{1} S(x,t)+I(x,t)dx = 1$. $S(x,t)$ and $I(x,t)$ denote the density of susceptible and infected populations at location $x$ and time $t$, $d_{S}$ and $d_{I}$ represent the diffusion coefficients for susceptible and infected individuals, and $\beta(x)$, $\gamma(x)$ are transmission and recovery rates at $x$, respectively.
We want to solve the steady state problem (same notations for convenience):
\[\left\{\begin{array}{l} d_{S} S_{x x}-\beta(x) \frac{S I}{S+I}+\gamma(x) I=0, \quad 0<x<1 \\ d_{I} I_{x x}+\beta(x) \frac{S I}{S+I}-\gamma(x) I=0, \quad 0<x<1 \\ S_{x}=I_{x}=0, \quad x=0,1, \end{array}\right.\]
where $\int_{0}^{1} S(x)+I(x)dx = 1$.
Note here elliptic problem has condition $\int_{0}^{1} S(x)+I(x)dx = 1$.
using OrdinaryDiffEq, SteadyStateDiffEq, ModelingToolkit, MethodOfLines,
DomainSets, Plots
# Parameters, variables, and derivatives
@parameters t x
@parameters dS=0.5 dI=0.1 brn=3 ϵ=0.1
@variables S(..) I(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
# Define functions
function γ(x)
y = x + 1.0
return y
end
function ratio(x, brn, ϵ)
y = brn + ϵ * sin(2 * pi * x)
return y
end
# 1D PDE and boundary conditions
eq = [
Dt(S(t, x)) ~
dS * Dxx(S(t, x)) -
ratio(x, brn, ϵ) * γ(x) * S(t, x) * I(t, x) / (S(t, x) + I(t, x)) +
γ(x) * I(t, x),
Dt(I(t, x)) ~
dI * Dxx(I(t, x)) +
ratio(x, brn, ϵ) * γ(x) * S(t, x) * I(t, x) / (S(t, x) + I(t, x)) -
γ(x) * I(t, x)]
bcs = [S(0, x) ~ 0.9 + 0.1 * sin(2 * pi * x),
I(0, x) ~ 0.1 + 0.1 * cos(2 * pi * x),
Dx(S(t, 0)) ~ 0.0,
Dx(S(t, 1)) ~ 0.0,
Dx(I(t, 0)) ~ 0.0,
Dx(I(t, 1)) ~ 0.0]
# Space and time domains
domains = [t ∈ Interval(0.0, 10.0),
x ∈ Interval(0.0, 1.0)]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [S(t, x), I(t, x)], [dS, dI, brn, ϵ])
# Method of lines discretization
# Need a small dx here for accuracy
dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx], t)
# Convert the PDE system into a DAE problem
prob = discretize(pdesys, discretization);DAEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 202-element Vector{Float64}:
0.9
0.9062790519529313
0.9125333233564304
0.9187381314585725
0.9248689887164855
0.9309016994374948
0.9368124552684678
0.9425779291565073
0.9481753674101716
0.9535826794978997
⋮
0.18763066800438638
0.190482705246602
0.19297764858882513
0.19510565162951538
0.1968583161128631
0.19822872507286887
0.1992114701314478
0.19980267284282716
0.2
du0: 202-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
⋮
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0Solving time-dependent SIS epidemic model
# Solving SIS reaction diffusion model
sol = solve(prob; saveat = 0.2);
# Retrieving the results
discrete_x = sol[x]
discrete_t = sol[t]
S_solution = sol[S(t, x)]
I_solution = sol[I(t, x)]
p = surface(discrete_x, discrete_t, S_solution)
display(p)Solving steady state problem
Change the elliptic problem to steady state problem of reaction diffusion equation.
See more solvers in Steady State Solvers · DifferentialEquations.jl
sys, tspan = symbolic_discretize(pdesys, discretization)
odeprob = ODEProblem(mtkcompile(sys), nothing, tspan)
steadystateprob = SteadyStateProblem(odeprob)
steadystate = solve(steadystateprob, DynamicSS(FBDF()))retcode: Success
u: 198-element Vector{Float64}:
0.33573429289930035
0.3357393550440854
0.3357472340386713
0.3357573868695643
0.33576928643473747
0.3357824225341335
0.335796302832723
0.33581045378999214
0.33582442155000164
0.3358377727864569
⋮
0.6764808507313231
0.6764571776245812
0.676431310463608
0.6764045765716223
0.6763783176156971
0.6763538845217277
0.6763326324181405
0.6763159156305036
0.6763050827487503The effect of human mobility on endemic size
Set the endemic size $f(d_{S},d_{I}) = \int_{0}^{1}I(x;d_{S},d_{I}).$
# Get the discretized I variables from the system
I_vars = filter(s -> contains(string(s), "I("), unknowns(odeprob.f.sys))
function episize!(dS_val, dI_val)
newprob = remake(odeprob, p = [dS => dS_val, dI => dI_val, brn => 3, ϵ => 0.1])
steadystateprob = SteadyStateProblem(newprob)
steadystate = solve(steadystateprob, DynamicSS(FBDF()))
y = sum(steadystate[v] for v in I_vars) * dx
return y
end
episize!(exp(1.0), exp(0.5))0.6589652754621023References:
- Allen L J S, Bolker B M, Lou Y, et al. Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model[J]. Discrete & Continuous Dynamical Systems, 2008, 21(1): 1.