Vaccination

Can vaccinating the contacts of cases, or the whole population, contain an outbreak that isolation and contact tracing alone do not? This page covers ring vaccination of traced contacts, post-exposure prophylaxis, mass and group vaccination campaigns, multi-dose schedules and limited vaccine supply.

It uses the outbreak model from Isolation and contact tracing: a Poisson offspring distribution with mean 3 (R = 3), an exponential generation time with mean 5 days, and an incubation period of LogNormal(1.5, 0.5) (mean and standard deviation of the log, median about 4.5 days). Infectiousness starts at infection. All times are in days, counted from the infection of the first index case. Each scenario simulates a number of outbreaks, stopped at max_cases cases, and containment_probability is the proportion that died out before reaching that cap (see the glossary). If Julia syntax such as x -> ..., do blocks or c ? a : b is new, see Julia for R users.

using EpiBranch
using Distributions
using StableRNGs

clinical = clinical_presentation(incubation_period = LogNormal(1.5, 0.5))

scenario(interventions = AbstractIntervention[], attributes = clinical) =
    ModelSpec(BranchingProcess(Poisson(3.0), Exponential(5.0)); interventions, attributes)

iso = Isolation(onset_to_isolation_delay = Exponential(2.0), duration = Inf)
ct = ContactTracing(
    probability = 0.7, isolation_to_trace_delay = Exponential(1.0),
    action = Quarantine(duration = Inf)
)

iso isolates symptomatic cases a mean of 2 days after onset, for good. ct traces 70% of their contacts a mean of 1 day after the infector's isolation and quarantines them, also for good.

How a dose protects

A dose protects from delay_to_immunity days after it is given (0 by default). It has up to four effects, each a probability; two of them, onward_efficacy and post_exposure_efficacy, are only on RingVaccination:

Parameter	What it prevents	Acts if protection starts
efficacy	the vaccinated person being infected	before the exposure
severity_efficacy	a severe outcome, such as death, in a vaccinated person who is infected anyway	when the clinical step that reads it says so (before symptom onset in the example below)
onward_efficacy	the vaccinated person infecting others	before each of their transmissions
post_exposure_efficacy	an infection the person already has	before symptom onset (if before the exposure, it blocks the infection itself)

efficacy is leaky by default: each exposure after protection starts is blocked with that probability. mode = AllOrNothingMode() instead makes that proportion of vaccinated people fully protected and leaves the rest unprotected. In a branching process every exposure is a different person, and the two give the same results on average. They differ in models where a person can be exposed more than once, such as households and networks.

severity_efficacy is the one effect the model does not apply on its own. It reaches the outbreak only through a clinical step in the model's progression (see Clinical transitions), such as Death or Hospitalisation, whose probability reads severity_efficacy and immunity_time for the person:

Death(delay = LogNormal(2.5, 0.4),
      probability = (rng, ind) ->
          immunity_time(ind) <= onset_time(ind) ?
              0.7 * (1 - severity_efficacy(ind)) : 0.7)

Compare many replicates

Every scenario uses its own random draws. Two scenarios run once each can differ by chance, and adding a dose can even seem to lower containment. Compare scenarios over many simulated outbreaks, and check that a difference is larger than the Monte Carlo error: for a containment probability p estimated from n outbreaks the standard error is about sqrt(p * (1 - p) / n).

Ring vaccination

RingVaccination vaccinates traced contacts and needs ContactTracing in the same scenario:

rv = RingVaccination(efficacy = 0.8)

rng = StableRNG(42)
results = simulate(scenario([iso, ct]), 200; max_cases = 500, rng = rng)
println("Isolation + tracing: $(round(containment_probability(results), digits=3))")

rng = StableRNG(42)
results = simulate(scenario([iso, ct, rv]), 200; max_cases = 500, rng = rng)
println("Isolation + tracing + ring vaccination: $(round(containment_probability(results), digits=3))")

Isolation + tracing: 0.33
Isolation + tracing + ring vaccination: 0.33

Under this tracing, vaccination makes no difference

The two numbers are identical, by construction. Tracing starts after the infector is isolated, and isolation already stops any further exposure of their contacts. A dose that protects against future exposure (efficacy) has nothing left to prevent, with or without quarantine.

efficacy matters only where a contact can still be exposed after being traced: under leaky isolation (post_isolation_transmission > 0), with a finite isolation duration that releases the infector while the contact is still susceptible, when tracing starts before the infector is isolated (for example eligibility = OnSymptomOnset()), or in a ring wider than direct contacts (depth > 1) that passes through people who keep transmitting after they are traced. To protect contacts who have already been exposed, use post_exposure_efficacy or onward_efficacy (see Protecting a contact who has already been exposed).

Tracing that starts at the infector's symptom onset, before isolation, and does not quarantine, leaves contacts exposed after the trace, and there the dose acts:

ct_onset = ContactTracing(probability = 0.7,
    isolation_to_trace_delay = Exponential(1.0),
    eligibility = OnSymptomOnset(), action = FlagOnly())

for (label, interventions) in (
        ("no vaccination", [iso, ct_onset]),
        ("ring vaccination", [iso, ct_onset, rv]))
    runs = simulate(scenario(interventions), 1000; max_cases = 500, rng = StableRNG(42))
    println(rpad(label, 20), round(containment_probability(runs), digits = 3))
end

no vaccination      0.191
ring vaccination    0.23

Delay to immunity

A vaccine that takes time to protect can be too late. Here protection starts 7 days after the dose, compared with no vaccination and same-day protection over more outbreaks, with the standard error of each:

rv_delayed = RingVaccination(efficacy = 0.8, delay_to_immunity = 7.0)

n_delay = 4000
for (label, interventions) in (
        ("no vaccination", [iso, ct_onset]),
        ("same-day protection", [iso, ct_onset, rv]),
        ("7-day delay", [iso, ct_onset, rv_delayed]))
    runs = simulate(scenario(interventions), n_delay; max_cases = 500, rng = StableRNG(42))
    p = containment_probability(runs)
    println(rpad(label, 22), round(p, digits = 3),
        " (standard error ", round(sqrt(p * (1 - p) / n_delay), digits = 3), ")")
end

no vaccination        0.216 (standard error 0.007)
same-day protection   0.227 (standard error 0.007)
7-day delay           0.214 (standard error 0.006)

Contacts infected after the trace are infected in the days between the trace and their infector's isolation. Here the trace comes a mean of 1 day after onset and the isolation a mean of 2 days after it, so that gap is longer than a week for only about 2% of contacts. A dose that protects 7 days later comes too late for almost all of them. Runs of 20,000 outbreaks put the 7-day result level with no vaccination, so any gap between the two above is simulation error.

Counting doses

A run stopped at max_cases leaves the last cases' contacts undecided and not yet vaccinated, so the run here is bounded by generations instead: max_generations = 5, with max_cases = nothing, draws no contacts for cases in the fifth generation and leaves every earlier case with all its contacts decided. condition = 50:200 keeps simulating until an outbreak stops with between 50 and 200 cases and returns it. count(is_vaccinated, state.individuals) then counts everyone vaccinated, like sum(is_vaccinated(x)) in R:

rng = StableRNG(42)
state = simulate(scenario([iso, ct, rv]); condition = 50:200,
    max_cases = nothing, max_generations = 5, rng = rng)
n_vaccinated = count(is_vaccinated, state.individuals)
n_infected = count(is_infected, state.individuals)
# Cases in the fifth generation drew no contacts to vaccinate
n_complete = count(ind -> is_infected(ind) && ind.generation < 5, state.individuals)
println("Vaccinated: $n_vaccinated, Infected: $n_infected")
println("Doses per case: $(round(n_vaccinated / n_complete, digits = 1))")

Vaccinated: 152, Infected: 122
Doses per case: 2.2

The doses per case divide by the cases before the fifth generation, the ones whose contacts could be traced and vaccinated.

Rings beyond direct contacts

By default tracing reaches a case's direct contacts. depth = 2 also traces the contacts of those contacts, the second ring that Ebola ring vaccination protocols vaccinate around a confirmed case. An infected contact who meets the tracing trigger starts a ring of their own. Contacts who do not start a ring of their own (those not infected, or infected without meeting the trigger, such as asymptomatic contacts) are followed one step further so the ring can extend past them. The same RingVaccination vaccinates everyone in the ring:

ct2 = ContactTracing(
    probability = 0.7, isolation_to_trace_delay = Exponential(1.0),
    action = Quarantine(duration = Inf), depth = 2
)

rng = StableRNG(42)
state = simulate(scenario([iso, ct2, rv]); condition = 50:200,
    max_cases = nothing, max_generations = 5, rng = rng)
doses_depth2 = count(is_vaccinated, state.individuals)
n_complete2 = count(ind -> is_infected(ind) && ind.generation < 5, state.individuals)
println("Doses with a second ring: $doses_depth2, ",
    "per case: $(round(doses_depth2 / n_complete2, digits = 1))")

Doses with a second ring: 136, per case: 3.7

The second ring around a case in the fourth generation would reach the sixth, past the bound, so the doses per case here are an underestimate. A wider ring reaches more people and uses more doses, but the extra doses protect only people still exposed after they are traced. Here every infected ring member with symptoms is isolated and starts a ring of their own. Their contacts are traced after that isolation, when the second ring's doses have nothing left to protect against. The second ring prevents infections when it passes through people who keep transmitting after they are traced, such as asymptomatic contacts under tracing without quarantine. Compare dose counts and containment across depths for your own setting.

Protecting a contact who has already been exposed

By the time a contact is traced under the default tracing, their infector is already isolated. The contact faces no new exposure from that infector, and whether they were infected before is already decided. efficacy then has nothing to block, but the dose can still act on an infection the contact already has:

• post_exposure_efficacy stops the infection, with that probability, if protection starts before the contact's symptom onset. The contact transmits as usual until then and not at all afterwards. They never develop symptoms, and nothing else in their clinical course (hospitalisation, death, recovery) happens after that point. They still count as a case, and appear in the line list with a date_infection_aborted and no onset date.

• onward_efficacy leaves the infection and its disease alone, and blocks each of the contact's transmissions after protection starts with that probability, whenever their onset falls.

A dose can set both: it then stops some infections and makes the rest less infectious. Both act only on transmission after protection starts. Neither does much for a contact who has already infected most of their own contacts. Quarantine from the trace already stops that later transmission, and the examples here trace without quarantine. The clearest measure of the effect is the mean number of secondary cases per traced case, an effective reproduction number among traced cases. A run stopped at max_cases leaves the last cases' contacts undecided, so the function below counts cases from runs bounded at five generations instead, below that bound, where every contact was decided:

ct_noquarantine = ContactTracing(probability = 0.7,
    isolation_to_trace_delay = Exponential(1.0), action = FlagOnly())

function secondary_per_traced(runs)
    mean(count(id -> is_infected(s.individuals[id]), ind.secondary_case_ids)
    for s in runs for ind in s.individuals
    if is_infected(ind) && is_traced(ind) && ind.generation < 5)
end
traced_runs(interventions) = simulate(scenario(interventions), 400;
    max_cases = nothing, max_generations = 5, rng = StableRNG(42))

for (label, vaccination) in [
    ("no vaccine", nothing),
    ("efficacy = 0.9", RingVaccination(efficacy = 0.9)),
    ("onward_efficacy = 0.9",
        RingVaccination(efficacy = 0.0, onward_efficacy = 0.9)),
    ("post_exposure_efficacy = 0.9",
        RingVaccination(efficacy = 0.0, post_exposure_efficacy = 0.9)),
]
    interventions = vaccination === nothing ? [iso, ct_noquarantine] :
        [iso, ct_noquarantine, vaccination]
    runs = simulate(scenario(interventions), 400; max_cases = 500, rng = StableRNG(42))
    println(rpad(label, 30), "secondary cases per traced case ",
        round(secondary_per_traced(traced_runs(interventions)), digits = 2),
        ", containment ", round(containment_probability(runs), digits = 3))
end

no vaccine                    secondary cases per traced case 1.94, containment 0.15
efficacy = 0.9                secondary cases per traced case 1.94, containment 0.15
onward_efficacy = 0.9         secondary cases per traced case 1.66, containment 0.232
post_exposure_efficacy = 0.9  secondary cases per traced case 1.64, containment 0.248

onward_efficacy and post_exposure_efficacy both reduce the secondary cases per traced case, from 1.94 to 1.66 and 1.64, while efficacy changes nothing. Containment rises with them, from 0.15 to 0.232 and 0.218. At 400 outbreaks per scenario each containment estimate has a standard error of about 0.02, so these gains are three to four times the Monte Carlo error. The gap between onward_efficacy and post_exposure_efficacy is within it, and another set of runs can place those two in either order.

Post-exposure prophylaxis

There is no separate post-exposure prophylaxis (PEP) intervention: antivirals or antibiotics given to traced contacts are a RingVaccination dose with delay_to_immunity = 0 (the default) and a post_exposure_efficacy. How much it achieves depends on speed, because protection has to start before symptom onset, and the incubation period here averages about five days:

for days in [0.0, 2.0, 5.0, 21.0]
    let pep = RingVaccination(efficacy = 0.0, post_exposure_efficacy = 0.9,
            delay_to_immunity = days)
        results = traced_runs([iso, ct_noquarantine, pep])
        cases = [ind for s in results for ind in s.individuals
                 if is_infected(ind) && is_traced(ind)]
        stopped = count(ind -> isfinite(EpiBranch.infection_aborted_time(ind)), cases)
        println("Protection after $(lpad(Int(days), 2)) days: ",
            round(Int, 100 * stopped / length(cases)), "% of traced cases stopped, ",
            "secondary cases per traced case ",
            round(secondary_per_traced(results), digits = 2))
    end
end

Protection after  0 days: 44% of traced cases stopped, secondary cases per traced case 1.64
Protection after  2 days: 26% of traced cases stopped, secondary cases per traced case 1.82
Protection after  5 days: 8% of traced cases stopped, secondary cases per traced case 1.92
Protection after 21 days: 0% of traced cases stopped, secondary cases per traced case 1.95

The share of infections stopped, and the reduction in secondary cases, fall quickly as the delay grows; with protection three weeks after the trace nothing is left.

Do not set efficacy and post_exposure_efficacy together

Where protection is already in place at the exposure, post_exposure_efficacy blocks the infection itself. It covers every contact that efficacy would protect, and setting both counts that protection twice. For an asymptomatic contact, who has no onset, blocking at exposure is all it can do. post_exposure_efficacy needs an incubation period, set by clinical_presentation.

Contacts of stopped infections are not traced

Under quarantine neither post_exposure_efficacy nor onward_efficacy has transmission left to block, and post_exposure_efficacy can even lower containment. A contact whose infection is stopped never develops symptoms, so they are never isolated and the people they infected before their dose are not traced from them. A programme that relies on symptom onset to trigger tracing loses that trigger for treated contacts. Tracing with depth = 2 reaches them instead, as contacts of contacts of the original case.

Clustered vaccine refusal

If each traced contact decides independently whether to accept a dose, every community ends up with about the same coverage. In practice refusal clusters by household or community: the people who avoid tracing also tend to decline vaccination.

vaccine_acceptance gives each community its own acceptance probability, drawn from a distribution, and everyone in the community shares it. Communities come from groups: list groups(...) before vaccine_acceptance(...) in the population characteristics. coverage on the dose is the probability a traced contact accepts it; the function (rng, ind) -> ind.state[:vaccine_acceptance] reads each person's community value:

community = groups(20)  # 20 communities
acceptance = vaccine_acceptance(propensity = Beta(2, 2))  # mean 0.5, varies by community

rv_clustered = RingVaccination(efficacy = 0.8,
    coverage = (rng, ind) -> ind.state[:vaccine_acceptance])
rv_independent = RingVaccination(efficacy = 0.8, coverage = 0.5)

The function below takes each community with at least 5 traced members, in each simulated outbreak, and records the share of them vaccinated:

function coverage_by_community(states)
    out = Float64[]
    for state in states, g in 1:20
        members = filter(state.individuals) do ind
            is_traced(ind) && get(ind.state, :group, 0) == g
        end
        length(members) >= 5 || continue
        push!(out, count(is_vaccinated, members) / length(members))
    end
    return out
end

rng = StableRNG(42)
clustered = coverage_by_community(simulate(
    scenario([iso, ct, rv_clustered], [clinical, community, acceptance]), 50;
    max_cases = 300, rng = rng))

rng = StableRNG(42)
independent = coverage_by_community(simulate(
    scenario([iso, ct, rv_independent], [clinical, community]), 50;
    max_cases = 300, rng = rng))

share_extreme(x) = count(c -> c < 0.2 || c > 0.8, x) / length(x)
println("Mean community coverage: clustered $(round(mean(clustered), digits=3)), independent $(round(mean(independent), digits=3))")
println("Communities under 20% or over 80% covered: clustered $(round(share_extreme(clustered), digits=3)), independent $(round(share_extreme(independent), digits=3))")

Mean community coverage: clustered 0.484, independent 0.499
Communities under 20% or over 80% covered: clustered 0.256, independent 0.0

Average coverage is about the same either way. Clustering spreads it out, leaving some communities almost fully covered and others almost untouched.

If acceptance is exactly 0 or 1 in each community, as with (rng, ind) -> Float64(rand(rng, Bernoulli(0.5))), each community accepts or refuses as a whole. With a Beta distribution, as here, communities differ but people within one still decide individually.

A community's acceptance lasts the whole outbreak, so a community that refuses keeps refusing in later generations. Give GroupVaccination the same coverage function to cluster refusal within the groups it vaccinates. MassVaccination has no coverage; its eligibility_time returns a day instead, and Inf means never vaccinated, so (rng, ind) -> rand(rng) < ind.state[:vaccine_acceptance] ? 30.0 : Inf vaccinates each person on day 30 with their community's acceptance as the probability.

Clustered refusal and containment

Under ct, ring vaccination cannot change whether an outbreak is contained (see the warning in Ring vaccination). Measuring the effect of clustering on containment needs a scenario where vaccination acts. Here tracing does not quarantine, and onward_efficacy reduces a vaccinated contact's own onward transmission. Every scenario uses the same population characteristics (clinical, community, acceptance), so vaccination is the only difference in how they are set up:

rv_clustered_onward = RingVaccination(efficacy = 0.8, onward_efficacy = 0.8,
    coverage = (rng, ind) -> ind.state[:vaccine_acceptance])
rv_independent_onward = RingVaccination(efficacy = 0.8, onward_efficacy = 0.8,
    coverage = 0.5)

n_outbreaks = 10000
containment(interventions) = containment_probability(simulate(
    scenario(interventions, [clinical, community, acceptance]), n_outbreaks;
    max_cases = 200, rng = StableRNG(42)))

for (label, interventions) in (
        ("no vaccination", [iso, ct_noquarantine]),
        ("clustered refusal", [iso, ct_noquarantine, rv_clustered_onward]),
        ("independent refusal", [iso, ct_noquarantine, rv_independent_onward]))
    p = containment(interventions)
    println(rpad(label, 20), round(p, digits = 3),
        " (standard error ", round(sqrt(p * (1 - p) / n_outbreaks), digits = 3), ")")
end

no vaccination      0.212 (standard error 0.004)
clustered refusal   0.227 (standard error 0.004)
independent refusal 0.219 (standard error 0.004)

Vaccinating about half the traced contacts raises containment a little, from 0.212 to 0.226 with clustered refusal and 0.224 with independent refusal. The standard error of a difference between two scenarios is about 1.4 times that of each, about 0.006 here, so the gain from vaccination is around two and a half standard errors. The 0.002 between clustered and independent refusal is well within the noise of these runs.

Clustering has little effect here because groups(20) puts each person in a community independently of who infected them. A case's contacts are spread across communities. A community's acceptance rarely lines up with the people that case goes on to infect. Clustering matters more when communities follow transmission, as households or a contact network do: a community that refuses then keeps transmitting within itself, which can lower containment at the same average coverage. Measure the effect in the model you are using.

Mass vaccination

MassVaccination offers vaccine to every contact the outbreak generates, whether traced or not; index cases are not vaccinated. eligibility_time is the day a person can be vaccinated, and protection starts delay_to_immunity days later; a person is protected against an exposure only if both have passed by then. Here everyone becomes eligible on day 30 and protection starts 14 days later:

mv = MassVaccination(efficacy = 0.85, eligibility_time = 30.0,
    delay_to_immunity = 14.0)

rng = StableRNG(42)
results = simulate(scenario(), 200; max_cases = 500, rng = rng)
println("No interventions: $(round(containment_probability(results), digits=3))")

rng = StableRNG(42)
results = simulate(scenario([mv]), 200; max_cases = 500, rng = rng)
println("Mass vaccination from day 30: $(round(containment_probability(results), digits=3))")

No interventions: 0.04
Mass vaccination from day 30: 0.04

With R = 3, an outbreak that takes off reaches the case cap well before protection starts on day 44. Containment stays level with the baseline without interventions, at 0.04 for both. At 200 outbreaks the standard error on each estimate is about 0.015.

For a gradual rollout in which each person becomes eligible at a random time, pass a distribution. Exponential(60.0) makes the mean wait 60 days:

mv_random = MassVaccination(efficacy = 0.85,
    eligibility_time = Exponential(60.0),
    delay_to_immunity = 14.0)

For a rollout by age, pass a function of the individual. This needs ages: add demographics to the population characteristics. The rule ind.state[:age] >= 65 ? 30.0 : 90.0 reads "day 30 if aged 65 or over, otherwise day 90":

mv_age = MassVaccination(
    efficacy = 0.85,
    eligibility_time = (rng, ind) -> ind.state[:age] >= 65 ? 30.0 : 90.0,
    delay_to_immunity = 14.0,
)

attrs = [clinical, demographics(age_distribution = Uniform(0, 90))]
rng = StableRNG(42)
results = simulate(scenario([mv_age], attrs), 200; max_cases = 500, rng = rng)
println("Age-stratified rollout: $(round(containment_probability(results), digits=3))")

Age-stratified rollout: 0.055

Containment is at the no-intervention baseline again: the older group is protected from day 44 and everyone else from day 104, and an outbreak that takes off has reached the case cap well before either. The printed 0.055 and the baseline 0.04 are two noisy estimates of the same quantity, as the standard error above suggests; runs of 3,000 give 0.055 and 0.059, a gap of under one standard error of the difference.

Efficacy that varies between people

efficacy can be one number for everyone, a distribution (people differ at random, with one draw per person, as for a varied immune response), or a function of the individual such as one that depends on age. Here older people respond less well:

mv_heterogeneous = MassVaccination(
    efficacy = (rng, ind) -> ind.state[:age] >= 65 ? 0.7 : 0.9,
    eligibility_time = (rng, ind) -> ind.state[:age] >= 65 ? 30.0 : 90.0,
    delay_to_immunity = 14.0,
)

Multi-dose schedules

For prime and boost (or longer) schedules, list one MassVaccination or RingVaccination per dose, each with its own dose_label. A dose protects against an exposure only if its protection has started by then. :prime and :boost are labels; the leading colon makes them Julia symbols.

prime = MassVaccination(efficacy = 0.6, eligibility_time = 30.0,
    delay_to_immunity = 14.0, dose_label = :prime)
boost = MassVaccination(efficacy = 0.5, eligibility_time = 60.0,
    delay_to_immunity = 14.0, dose_label = :boost)

rng = StableRNG(42)
results = simulate(scenario([prime, boost]), 200; max_cases = 500, rng = rng)
println("Two-dose rollout: $(round(containment_probability(results), digits=3))")

Two-dose rollout: 0.04

This is the baseline again, for the same reason: the prime protects from day 44 and the boost from day 74, by which time the outbreak has reached the cap.

Enter the efficacy each dose adds

Doses protect independently: an exposure is blocked unless every dose whose protection has started fails to block it. With efficacies e1 and e2, overall protection is 1 - (1 - e1) * (1 - e2). The boost's efficacy is therefore the protection it adds among people the prime left unprotected. For 60% protection after one dose and 80% after two, set the prime to 0.6 and the boost to (0.8 - 0.6) / (1 - 0.6) = 0.5, as above. Setting the boost to 0.8 would give 92% overall.

is_vaccinated(ind; dose_label = :boost) tells whether a person had a given dose, and linelist has one date column per dose, such as date_vaccination_prime and date_immunity_boost.

Two doses in a ring

Ring doses are given at the trace. A second dose sets dose_delay, the days from the trace to that dose, and names the dose it follows with requires_dose. Only contacts who have had the earlier dose by the time the later one is due receive it. The boost's coverage is then the share retained between doses. List each dose after the dose it requires; listing a boost before its prime is an error.

prime_ring = RingVaccination(efficacy = 0.6, delay_to_immunity = 21.0,
    coverage = 0.8, dose_label = :prime)
boost_ring = RingVaccination(efficacy = 0.5, dose_delay = 28.0,
    delay_to_immunity = 14.0, coverage = 0.9,
    requires_dose = :prime, dose_label = :boost)

rng = StableRNG(42)
results = simulate(scenario([iso, ct, prime_ring, boost_ring]), 200;
    max_cases = 500, rng = rng)
println("Two-dose ring: $(round(containment_probability(results), digits=3))")

Two-dose ring: 0.325

A protocol that gives the boost four to six weeks after the trace uses a distribution, dose_delay = Uniform(28.0, 42.0), drawn once per contact. delay_to_immunity, post_exposure_efficacy and onward_efficacy also take a distribution or a function of the individual, each drawn once per contact when the dose is given. A contact keeps the same protection against every exposure.

The mean number of each dose per simulated outbreak shows what the second dose costs:

doses(label) = mean(count(i -> is_vaccinated(i; dose_label = label), s.individuals)
    for s in results)
println("Primed: $(round(doses(:prime), digits = 1)), ",
    "boosted: $(round(doses(:boost), digits = 1))")

Primed: 567.2, boosted: 510.1

A later dose is recorded when it falls due, even if the contact was infected in the meantime. The boost count is a count of doses scheduled.

Under this tracing neither dose changes the outbreak: no ring member is infected after being traced (see the warning in Ring vaccination). Even under ct_onset, where infections do follow the trace, they follow within days. Both doses here start protecting weeks after the trace, too late for those infections (see Delay to immunity). The 0.325 above is not comparable with the 0.33 printed for isolation and tracing alone: each dose's coverage is drawn, which shifts the random stream, and setting both efficacies to zero gives 0.325 as well. Removing the boost and re-running is not a controlled comparison either. The two runs use different random draws, and their containment can differ by chance in either direction (see the note on comparing many replicates above).

Repeat campaign visits

GroupVaccination vaccinates every member of a group, such as a household or community, once a case in it is detected (lab-confirmed by default), dose_delay days later. A campaign that visits each group several times, and reaches a person missed on one visit on a later one, can be set up with its coverage and dose_delay. People who refuse outright are a separate population characteristic.

With three visits, each reaching a willing person with probability c, the probabilities that the first, second and third visit is the first to reach them are c, (1-c)*c and (1-c)^2*c. Their sum is the probability of ever being reached, and the coverage among willing people. Rescaled to sum to one, they give the distribution of the delay from the group's trigger to the visit that reaches them:

visit_times = [0.0, 7.0, 14.0]
reach_per_visit = 0.6
first_reached = [reach_per_visit * (1 - reach_per_visit)^(i - 1)
                 for i in eachindex(visit_times)]
campaign_reach = sum(first_reached)
visit_delay = DiscreteNonParametric(visit_times, first_reached ./ campaign_reach)

# Each person's willingness is drawn once, when they are created, and stored
# as :willing; it is separate from whether a visit reaches them.
prob_willing = 0.9
willingness = (rng, ind) -> (ind.state[:willing] = rand(rng) < prob_willing)
repeated_campaign = GroupVaccination(efficacy = 0.8,
    coverage = (rng, ind) -> ind.state[:willing] ? campaign_reach : 0.0,
    dose_delay = visit_delay)

campaign_model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    attributes = [clinical, groups(3), willingness],
    interventions = [Isolation(onset_to_isolation_delay = Exponential(1.0), duration = Inf),
        repeated_campaign])

println("Reached among willing people: $(round(campaign_reach, digits = 3)); ",
    "expected overall coverage: $(round(prob_willing * campaign_reach, digits = 3))")

Reached among willing people: 0.936; expected overall coverage: 0.842

Overall coverage is the share willing times the share of willing people reached. Those are the shares the parameters aim for, and what a run delivers is close to them, up to the sampling noise of a finite number of outbreaks:

campaign_runs = simulate(campaign_model, 200; max_cases = 300, rng = StableRNG(42))
people = [ind for s in campaign_runs for ind in s.individuals]
willing_people = filter(ind -> get(ind.state, :willing, false), people)
println("Measured: among willing people ",
    round(count(is_vaccinated, willing_people) / length(willing_people), digits = 3),
    ", overall ", round(count(is_vaccinated, people) / length(people), digits = 3))

Measured: among willing people 0.923, overall 0.833

A person with :willing == false refuses throughout the campaign, and a willing person left unvaccinated was never reached. For a reach probability of zero, use coverage = 0.0 directly. A single visit gives the usual campaign with one opportunity. Use a different dose_label for another dose; willingness can be shared across doses, as here, or stored separately for each.

Limitations of this campaign recipe

Each person's eventual dose and its date are decided once, and protection starts the usual delay_to_immunity after that date. Visits are not simulated one by one. The campaign cannot respond to limited team capacity or to how the outbreak develops, and failed visits are not recorded. Visits are assumed to reach people independently, with a fixed probability.

Limited vaccine supply

The scenarios above vaccinate everyone they reach. Real capacity is finite: a fixed number of teams, a daily dose limit, a stockpile. CapacityConstrained limits how many people a ring, group or mass vaccination reaches: at most budget_per_period people every period days. People beyond the limit are not vaccinated, and by default those traced first are served first:

rv_capped = CapacityConstrained(rv; budget_per_period = 5.0, period = 1.0)

state_uncapped = simulate(scenario([iso, ct, rv]); condition = 50:200, max_cases = 200, rng = StableRNG(42))
state_capped = simulate(scenario([iso, ct, rv_capped]); condition = 50:200, max_cases = 200, rng = StableRNG(42))

println("Doses without a limit: $(count(is_vaccinated, state_uncapped.individuals))")
println("Doses with at most 5 people admitted a day: $(count(is_vaccinated, state_capped.individuals))")

Doses without a limit: 421
Doses with at most 5 people admitted a day: 146

The simulation offers one case's traced contacts at a time, in the order it settles cases, but checks each proposed dose against the allowance of the day it is actually dated, whatever day the case making the trace was infected. A dose for one day can be admitted before a case settled later asks for an earlier day; the allowance for each day still holds across both.

capacity_usage reports how much of the allowance a simulated outbreak had used, and how much was available, by the time it stopped:

capacity_usage(rv_capped, state_capped)

(used = 84, available = 85.0)

priority decides who is served when more people are waiting than the allowance allows: it is a function of the individual and the simulation state, and people with lower values are served first. The default orders by trace time. A lottery among the people waiting:

rv_lottery = CapacityConstrained(rv; budget_per_period = 5.0, period = 1.0,
    priority = (ind, state) -> rand(state.rng))

period = Inf (the default) makes the budget a single stockpile for the whole outbreak:

rv_stockpile = CapacityConstrained(rv; budget_per_period = 200.0)

With a finite period, carry_over = true (the default) adds a day's unused allowance to the next day's; carry_over = false discards it, so only that day's own budget_per_period is ever available. With carry_over = true admissions can bunch: a quiet spell builds an allowance that a later day spends at once, well above budget_per_period in a day.

Which models support capacity limits

On the network and household models, capacity limits work for group vaccination, ring vaccination and mass vaccination, with any eligibility_window: a dose the race withdraws for falling outside the window frees the capacity it used. To limit an intervention of your own, see Capacity limits.