GitHub

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.2
 0.19980267284282716
 0.1992114701314478
 0.19822872507286887
 0.1968583161128631
 0.19510565162951538
 0.19297764858882516
 0.19048270524660196
 0.18763066800438638
 0.1844327925502015
 ⋮
 0.8518246325898284
 0.8574220708434928
 0.8631875447315321
 0.8690983005625053
 0.8751310112835144
 0.8812618685414275
 0.8874666766435696
 0.8937209480470687
 0.9
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.0

Solving 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.33573429289878326
 0.33573935504356833
 0.33574723403815415
 0.3357573868690472
 0.33576928643422044
 0.33578242253361645
 0.33579630283220613
 0.3358104537894754
 0.335824421549485
 0.33583777278594046
 ⋮
 0.6764808507305107
 0.6764571776237701
 0.6764313104627986
 0.676404576570814
 0.6763783176148898
 0.6763538845209214
 0.6763326324173349
 0.6763159156296987
 0.6763050827479459

The 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.6589652754981223

References:

  • 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.