Cambiar idioma

English

Iniciar sesión en Mathigon

o

Nueva cuenta     Restablecer la contraseña     

Epidemic modeling
by Samuel S. Watson
Introduction
Stochastic SIR models
Deterministic SIR models
Lessons for the real world
Links

Compartir

Restablecer progreso

Esto eliminará su progreso y los datos de chat de todos los capítulos de este curso y no se puede deshacer.

Glosario

Seleccione una de las palabras clave de la izquierda ...

Epidemic modelingDeterministic SIR models

Tiempo de leer: ~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 n-1 susceptible individuals, so the expected number of individuals which get infected p(n-1). Rather than choosing this number randomly, we'll just directly change the value by p(n-1). 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 i_1 = p(n-1) infectious individuals, s_1 = n - 1 - p(n-1) susceptible individuals, and r_1 = 1 recovered individual. On average, there are transmission opportunities, and each one results in infection with probability R_0/n. Therefore, the expected number of new infections at this time step is approximately s_1i_1R_0/n. Then we perform the same update using the value s_1i_1R_0/n 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")

The evolution of the numbers of susceptible, infectious, and recovered individuals over time, where each incremental change is made according to the average behavior of the corresponding change in the random system.

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

The deterministic equations describe the average behavior of the random system, and the random system sometimes follows the deterministic trajectory quite closely.

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 s as \Delta s, the equations might be written in a math textbook as

\begin{align*}\Delta s_t &= -\frac{s_ti_tR_0}{n} \\ \Delta i_t &= \frac{s_ti_tR_0}{n} - i_t \\ \Delta r_t &= i_t. \end{align*}

If we want to obtain even smoother behavior by making time continuous as well, we should use a time increment of \Delta t rather than 1. Since the amount by which s or i or r changes in a given short interval of length \Delta t 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 \Delta t:

\begin{align*}\Delta s_t &= -\frac{s_ti_tR_0}{n}\Delta t \\ \Delta i_t &= \left(\frac{s_ti_tR_0}{n} - i_t\right) \Delta t \\ \Delta r_t &= i_t \Delta t. \end{align*}

If we divide both sides by \Delta t and let \Delta t \to 0, we obtain the following system of differential equations:

\begin{align*}\frac{ds}{dt} &= -s_ti_tR_0/n \\ \frac{di}{dt} &= s_ti_tR_0/n - i_t \\ \frac{dr}{dt} &= i_t. \end{align*}

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")

The differential equation solution is very similar to the difference equation solution.

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 R_0. Make a plot of i_t for two R_0 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 i 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

Decreasing XEQUATIONX638XEQUATIONX decreases the maximum number of simultaneous cases while also putting that peak further into the future. There are also fewer individuals infected overall when XEQUATIONX639XEQUATIONX is smaller.

Bruno
Bruno Bruno