Epidemic modelingLessons for the real world

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.

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

Continue

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 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 value but not on the number of individuals who infectious initially.

Continue

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 fully eradicatedsignificantly reduced, 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.

Continue

This observation also helps us see why a temporary reduction in does not prevent a spike in cases later once goes back up. Unless the conditions for a lower 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 value, and enough initial infections to ensure widespread infection. To avoid a population-level outbreak, we need either a long-term reduction in or widespread immunity achieved through a vaccine.

Continue

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 geographytransportationRNA 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 , we obtain the following animation: