Epidemic modelingDeterministic SIR models
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,
For the next step, we'll have infectious individuals, susceptible individuals, and recovered individual. On average, there are
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
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
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
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.
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
plt = plot(ylabel = "number infected") plot_infections!(plt, n, 1.5) plot_infections!(plt, n, 2.0) plt