# 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

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

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 DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, Plots

# Parameters, variables, and derivatives
@parameters t x
@parameters dS dI brn ϵ
@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 => 0.5, dI => 0.1, brn => 3, ϵ => 0.1])

# Method of lines discretization
# Need a small dx here for accuracy
dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization);
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 198-element Vector{Float64}:
0.9062790519529313
0.9125333233564304
0.9187381314585725
0.9248689887164855
0.9309016994374948
0.9368124552684678
0.9425779291565073
0.9481753674101716
0.9535826794978997
0.9587785252292473
⋮
0.18443279255020154
0.18763066800438638
0.190482705246602
0.19297764858882513
0.19510565162951538
0.1968583161128631
0.19822872507286887
0.1992114701314478
0.19980267284282716

### Solving time dependent SIS epidemic model

# Solving SIS reaction diffusion model
sol = solve(prob, Tsit5(), saveat=0.2);

# Retriving 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

steadystateprob = SteadyStateProblem(prob)
steadystate = solve(steadystateprob, DynamicSS(Tsit5()))
u: 198-element Vector{Float64}:
0.3317836646213583
0.3317732755347774
0.33178911139346345
0.33175474345790923
0.33179088720238725
0.3317337314324058
0.33179111779166187
0.3317124422636641
0.33179191510139755
0.3316930492687938
⋮
0.6560941213709633
0.6561639796910776
0.6562347277330004
0.6563041441217199
0.6563698199311431
0.6564293278195016
0.6564800894280143
0.6565194894492582
0.6565447998587352

### 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}).$

function episize!(dS, dI)
newprob = remake(prob, p=[dS, dI, 3, 0.1])
episize!(exp(1.0),exp(0.5))
0.6656214868938777