Riječnik

Odaberite jednu od ključnih riječi s lijeve strane ...

Epidemic modelingLessons for the real world

Vrijeme čitanja: ~20 min

In this section we'll use our models as a lens through which to examine some important questions for real-world epidemic planning. It's worth repeating the caveat that such model-based conclusions should always be held skeptically and weighed rigorously against broader real-world considerations. However, this exercise can nevertheless be useful, because it can help build intuition about the basic mechanisms underlying important phenomena which are present in the models and known by experts to be present in real-world epidemics as well.

What difference does it make if we start with a handful of infectious individuals rather than just one? The difference is profound and has significant implications for real-world infectious diseases.

Exercise
Modify the code for the function SIR_simulation below to set the number of initially infectious individuals to 10 rather than 1. What happens to the resulting histogram?

Use your observation to draw a conclusion about the strategy of temporary reduction of R_0 in the context of this model.

using Plots
import Base.Iterators: countfrom
 
function SIR_simulation(population_size, infection_probability)
    statuses = fill("susceptible", population_size, 1)
    statuses[1, 1] = "infectious"
    transmissions = []
    for t in countfrom(2)
        n_infectious = sum(statuses[:, t-1] .== "infectious")
        if n_infectious == 0
            break
        end
        statuses = [statuses fill("susceptible", population_size)]
        for k in 1:population_size
            if statuses[k, t-1] == "recovered" || statuses[k, t-1] == "infectious"
                statuses[k, t] = "recovered"
            end
        end
        for j in 1:population_size
            if statuses[j, t-1] == "infectious"
                for k in 1:population_size
                    if statuses[k, t-1] == "susceptible"
                        if rand() < infection_probability
                            push!(transmissions, [(j, t-1),(k, t)])
                            statuses[k, t] = "infectious"
                        end
                    end
                end
            end
        end
    end
    statuses, transmissions
end
 
population_size = 20
infection_probability = 2 / population_size
statuses, transmissions = SIR_simulation(population_size, infection_probability)
statuses
 
n = population_size = 1000
p = infection_probability = 2/population_size
n_recovered(n, p) = sum(SIR_simulation(n, p)[1][:, end] .== "recovered")
histogram([n_recovered(n, p) for _ in 1:5_000], xlims = (0, 1000), 
          nbins = 100, label = "", xlabel = "final number recovered",
          ylabel = "number of runs")

Solution. What we see is that even with 10 initially infectious individuals, the disease spreads to a substantial percentage of the population with very high probability:

using Plots
import Base.Iterators: countfrom
 
function SIR_simulation(population_size, 
                        infection_probability,
                        n_initially_infected)
    statuses = fill("susceptible", population_size, 1)
    statuses[1:n_initially_infected, 1] .= "infectious"
    transmissions = []
    for t in countfrom(2)
        n_infectious = sum(statuses[:, t-1] .== "infectious")
        if n_infectious == 0
            break
        end
        statuses = [statuses fill("susceptible", population_size)]
        for k in 1:population_size
            if statuses[k, t-1] == "recovered" || statuses[k, t-1] == "infectious"
                statuses[k, t] = "recovered"
            end
        end
        for j in 1:population_size
            if statuses[j, t-1] == "infectious"
                for k in 1:population_size
                    if statuses[k, t-1] == "susceptible"
                        if rand() < infection_probability
                            push!(transmissions, [(j, t-1),(k, t)])
                            statuses[k, t] = "infectious"
                        end
                    end
                end
            end
        end
    end
    statuses, transmissions
end
 
n = population_size = 1000
p = infection_probability = 2/population_size
n_recovered(n, p) = sum(SIR_simulation(n, p, 10)[1][:, end] .== "recovered")
histogram([n_recovered(n, p) for _ in 1:5_000], xlims = (0, 1000), 
          nbins = 100, label = "", xlabel = "final number recovered",
          ylabel = "number of runs")

If R_0 is substantially larger than 1 and the initial number of infectious individuals is more than a few, then the probability of widespread infection is very high. In the 5000 runs used to produce the histogram above, all 5000 eventually reached approximately 80% of the population.

The final proportion infected depends on the R_0 value but not on the number of individuals who infectious initially.

While we derived this observation in the context of an SIR model, it does apply to real-life outbreaks as well: unless the disease is , R_0 matters far more than the number of individuals who are infectious to start.

When you heart the argument that an infectious disease should be ignored because the number of cases is still a small proportion of the population, it's analogous to arguing that an oven fire should be ignored until it at least engulfs the kitchen: flammability matters more than the size of the initial fire, unless the initial fire is nonexistent.

This observation also helps us see why a temporary reduction in R_0 does not prevent a spike in cases later once R_0 goes back up. Unless the conditions for a lower R_0 value are sustained long enough to eliminate the virus altogether, the post-reduction phase is essentially a fresh run of the model with large number of susceptible individuals, a high R_0 value, and enough initial infections to ensure widespread infection. To avoid a population-level outbreak, we need either a long-term reduction in R_0 or widespread immunity achieved through a vaccine.

Containment vs Mitigation

The SIR model may be made more realistic in a wide variety of ways. One of the most glaring omissions is the lack of in the model: in the real world, the virus spreads by the kind of physical contact which is much more likely between people who live close to one another. Let's see what happens when we incorporate spatial relationships into the model.

Let's arrange our individuals into a square grid, and we'll let each person transmit to its four neighbors:

using Plots

function infection_plot(statuses)
    p = plot(ratio = 1, legend = false, axis = false, grid = false)
    for (status, color) in (("susceptible", "gray"), 
                            ("infectious", "tomato"), 
                            ("recovered", "darkgreen"))
        pts = [Tuple(i) for i in CartesianIndices(statuses) if statuses[i] == status]
        if length(pts) > 0
            scatter!(p, pts, color = color, markersize = 6)
        end
    end
    p
end
 
"Return the four nodes adjacent to (i, j)"
function neighbors(i, j, sidelength, barrier = Set())
    filter(k -> all(1 .≤ k .≤ sidelength) && k ∉ barrier, [(i+1, j), (i-1, j), (i, j+1), (i, j-1)])
end
 
function spatial_simulation(sidelength, n_timesteps, R₀, barrier = Set())
    statuses = fill("susceptible", sidelength, sidelength, n_timesteps)
    statuses[1:3, 1:3, 1] .= "infectious"
    for t in 2:n_timesteps
        for i in 1:sidelength
            for j in 1:sidelength
                if statuses[i, j, t - 1] == "infectious"
                    nbs = neighbors(i, j, sidelength, barrier)
                    for (k,l) in nbs
                        if statuses[k, l, t-1] == "susceptible"
                            if rand() < R₀/length(nbs)
                                statuses[k, l, t] = "infectious"
                            end
                        end
                    end
                    statuses[i, j, t] = "recovered"
                elseif statuses[i, j, t - 1] == "recovered"
                    statuses[i, j, t] = "recovered"
                end
            end
        end
    end
    statuses
end
 
sidelength = 30
n_timesteps = 150
barrier_height = 25
barrier = Set([(10, k) for k in 1:barrier_height])
R₀ = 2.2
statuses = spatial_simulation(30, 150, R₀, barrier)
t = 150
infection_plot(statuses[:, :, t])
scatter!(collect(barrier), color = :black, markershape = :square)

Varying t, we obtain the following animation:

Exercise
Experiment with other variations on this model to explore the following idea: barriers don't seem to help much unless they manage to completely contain the virus, because community spread quickly becomes the dominant mode of transmission.

This exercise is more open-ended than the others in this course. You might want to do this exploration in a Jupyter notebook.

 
Bruno
Bruno Bruno