## Słownik

Wybierz jedno ze słów kluczowych po lewej stronie…

# Epidemic modelingDeterministic SIR models

Czas czytania: ~20 min

One difficulty with the stochastic SIR model in the previous section is that the results are random. This can make certain kinds of analysis more difficult, and it means that we have to run many simulations to get a sense of the results. It can be helpful to get a single deterministic result which represents the typical behavior of the random system. We can do this by replacing each random incremental change with a deterministic step equal to the average behavior for that step.

For example, consider the first step. We have one infectious individual and susceptible individuals, so the expected number of individuals which get infected . Rather than choosing this number randomly, we'll just directly change the value by . This number is not necessarily an integer, , because it represents an average, not an exact count of a set of individuals.

For the next step, we'll have infectious individuals, susceptible individuals, and recovered individual. On average, there are transmission opportunities, and each one results in infection with probability . Therefore, the expected number of new infections at this time step is approximately . Then we perform the same update using the value as the number of new infections.

Let's continue this calculation in code:

using Plots
n = 1000
R₀ = 1.5
s = [1.0 * (n - 1)]
i = [1.0]
r = [0.0]
for _ in 1:40
n_infected = R₀ * s[end] * i[end] / n
push!(s, s[end] - n_infected)
push!(i, i[end] - i[end] + n_infected )
push!(r, r[end] + i[end])
end

p = plot([s i r], label = ["s(t)" "i(t)" "r(t)"],
xlabel = "time", ylabel = "count")

To see how this is exercise is valuable, let's plot a simulation of the stochastic model from the previous section alongside the results from this non-random approach:

include("data-gymnasia/sir-simulation.jl")
population_size = 1000
infection_probability = 1.5 / population_size
statuses, transmissions = SIR_simulation(population_size, infection_probability)
p = plot([s i r], label = ["s(t)" "i(t)" "r(t)"],
xlabel = "time", ylabel = "count")
for status in ("susceptible", "infectious", "recovered")
inf_totals = sum(statuses .== status, dims = 1)[:]
scatter!(p, inf_totals, label = "")
end
p

We can see that the deterministic approach plays the role of the mean in the law of large numbers: it's the steady value around which the curve for the random system fluctuates.

## Connection to differential equations

In mathematical notation, the system we simulated above is called a system of difference equations. Writing the change in as , the equations might be written in a math textbook as

If we want to obtain even smoother behavior by making time continuous as well, we should use a time increment of rather than 1. Since the amount by which or or changes in a given short interval of length is approximately times as much as it would change over one unit of time, we should modify the equation above by multiplying the right-hand side by :

If we divide both sides by and let , we obtain the following system of differential equations:

While this system does not have an analytical solution, we can solve it numerically and plot the results using Julia's DifferentialEquations package:

using DifferentialEquations, Plots
import ParameterizedFunctions: @ode_def

sir_ode = @ode_def SIRModel begin
ds = -p*s*i
di = p*s*i - i
dr = i
end p

n = 1000
initialvalues = 1.0 * [n - 1, 1, 0]
timerange = (0.0, 40.0)
p = 1.5 / n

sir_problem = ODEProblem(sir_ode, initialvalues, timerange, [p])

sir_solution = solve(sir_problem)
plot(sir_solution)
plot!(ylabel = "count")

One advantage of the differential equations approach is that it allows us to make a rigorous statement about how the behavior of the random, discrete system is very close to the behavior of the deterministic, continuous system (in the sense that the infection curves are close with high probability). If you'd like to see a theorem to this effect, check out this 2015 paper by Ekkehard Beck and Benjamin Armbruster of Northwestern University.

Exercise
You've probably seen graphs labeled Flattening the Curve on social media. The idea is to reduce the transmission rate so that the the peak (of the graph of number of infections versus time) is lower and comes later. This is valuable for the health care system, because it reduces the required capacity and gives time to expand capacity.

Even for the simple SIR model, we can see numerically how flattening the curve is directly related to . Make a plot of for two values to see that the curve is flatter for the smaller one.

Note: we supply the keyword argument vars = (0, 2) to the plot method for the differential equation solution object to say that we want to plot time on the horizontal axis (variable 0) and on the vertical axis (variable 2).

using DifferentialEquations, Plots

sir_ode = @iode_def SIRModel begin
ds = -p*s*i
di = p*s*i - i
dr = i
end p

n = 1000
function plot_infections!(plt, n, R₀)
initialvalues = 1.0 * [n - 1, 1, 0]
timerange = (0.0, 40.0)
p = R₀ / n
sir_problem = ODEProblem(sir_ode, initialvalues, timerange, [p])
sir_solution = solve(sir_problem)
plot!(plt, sir_solution, vars = (0, 2), label = string(R₀))
end

plt = plot(ylabel = "number infected")
plot_infections!(plt, n, #= insert R₀ value here! =#)
plot_infections!(plt, n, #= insert R₀ value here! =#)
plt

Solution.

plt = plot(ylabel = "number infected")
plot_infections!(plt, n, 1.5)
plot_infections!(plt, n, 2.0)
plt
Bruno