Vaccination and time-dependent interventions

An intervention changes either the population state or the model parameters. We represent pre-outbreak vaccination by reducing the initial susceptible population, and a timed contact intervention by changing the transmission parameter during the epidemic.

Biological scenario

A population of 1,000 begins with 10 infectious people. We compare:

no intervention;
40% vaccination coverage with 80% efficacy before the outbreak;
the same vaccination plus a day-30 intervention reducing \(\beta\) from 0.30 to 0.12 per day.

The aim is to compare infectious peaks, peak times and infections occurring after the simulation begins.

Two interventions, two mathematical mechanisms

Vaccination before the outbreak

Changes the initial conditions by moving effectively immunised people out of \(S(0)\).

\[V_{\mathrm{immune}}=c\varepsilon S_{\mathrm{eligible}},\]

where \(c\) is coverage and \(\varepsilon\) is vaccine efficacy.

Timed contact intervention

Changes transmission during the simulation through a time-dependent parameter \(\beta(t)\).

\[\beta(t)=\begin{cases}\beta_1,&t<T_L,\\\beta_2,&t\geq T_L.\end{cases}\]

Vaccination changes the initial state

There are \(N-I(0)\) initially non-infectious people eligible to be susceptible in this simplified scenario:

\[S_{\mathrm{eligible}}=N-I(0).\]

With coverage \(c=0.4\) and efficacy \(\varepsilon=0.8\), the effectively immunised fraction is \(c\varepsilon=0.32\), not 0.40.

\[V_{\mathrm{immune}}=0.4\times0.8\times990=316.8.\]

\[S(0)=990-316.8=673.2.\]

The standard SIR symbol \(R\) is better interpreted here as removed or immune, because it contains vaccine-induced immunity as well as people who later recover from infection.

Program a time-dependent transmission parameter

Before day 30
\(\beta(t)=0.30\)

→

From day 30
\(\beta(t)=0.12\)

def beta_at_time(
    t,
    beta_before,
    beta_after,
    intervention_day
):
    if t < intervention_day:
        return beta_before
    else:
        return beta_after

t is the time supplied by the ODE solver.
if selects the pre-intervention value before the switching day.
else selects the post-intervention value from that day onwards.
The function returns one value of \(\beta(t)\) for each solver evaluation.

Use \(\beta(t)\) inside the SIR equations

beta = beta_at_time(
    t,
    beta_before,
    beta_after,
    intervention_day
)

infection_flow = beta * S * I / N

The equations retain the same structure, but the infection flow now uses the transmission value appropriate to the current time.

Scenarios compared

No intervention

No vaccine immunity and \(\beta=0.30\) throughout.

Vaccination

32% of initially eligible susceptible people effectively immunised; \(\beta=0.30\) throughout.

Vaccination + timed intervention

Same vaccine immunity; \(\beta\) changes to 0.12 on day 30.

Run the intervention comparison

Interactive PythonVaccination and timed intervention

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def beta_at_time(
    t,
    beta_before,
    beta_after,
    intervention_day
):
    """Return the transmission parameter at time t."""
    
    if t < intervention_day:
        return beta_before
    else:
        return beta_after

def sir_intervention_rhs(
    t,
    y,
    beta_before,
    beta_after,
    intervention_day,
    gamma
):
    """Return SIR derivatives with time-dependent transmission."""
    
    S, I, R = y
    N = S + I + R
    
    beta = beta_at_time(
        t,
        beta_before,
        beta_after,
        intervention_day
    )
    
    infection_flow = beta * S * I / N
    recovery_flow = gamma * I
    
    return [
        -infection_flow,
        infection_flow - recovery_flow,
        recovery_flow
    ]

def simulate_scenario(
    name,
    population,
    I0,
    coverage,
    efficacy,
    beta_before,
    beta_after,
    intervention_day,
    gamma,
    end_time=200
):
    """Solve one vaccination and transmission scenario."""
    
    if not 0 <= coverage <= 1:
        raise ValueError("coverage must lie between 0 and 1.")
    if not 0 <= efficacy <= 1:
        raise ValueError("efficacy must lie between 0 and 1.")
    if intervention_day < 0:
        raise ValueError("intervention_day cannot be negative.")
    
    eligible_susceptible = population - I0
    immune_from_vaccine = (
        coverage * efficacy * eligible_susceptible
    )
    
    S0 = eligible_susceptible - immune_from_vaccine
    R0 = immune_from_vaccine
    
    reporting_times = np.linspace(0, end_time, 2001)
    
    solution = solve_ivp(
        sir_intervention_rhs,
        t_span=(0, end_time),
        y0=[S0, I0, R0],
        t_eval=reporting_times,
        args=(
            beta_before,
            beta_after,
            intervention_day,
            gamma
        ),
        rtol=1e-8,
        atol=1e-10,
        max_step=0.5
    )
    
    if not solution.success:
        raise RuntimeError(solution.message)
    
    S, I, R = solution.y
    peak_index = np.argmax(I)
    
    return {
        "name": name,
        "time": solution.t,
        "S": S,
        "I": I,
        "R": R,
        "S0": S0,
        "I0": I0,
        "vaccine_immune": immune_from_vaccine,
        "peak_I": I[peak_index],
        "peak_time": solution.t[peak_index],
        "outbreak_infections": I0 + S0 - S[-1],
        "beta_before": beta_before,
        "beta_after": beta_after,
        "intervention_day": intervention_day
    }

# Shared settings
N = 1000
I0 = 10
gamma = 0.1

# Run three scenarios
scenarios = [
    simulate_scenario(
        name="No intervention",
        population=N,
        I0=I0,
        coverage=0,
        efficacy=0.8,
        beta_before=0.30,
        beta_after=0.30,
        intervention_day=30,
        gamma=gamma
    ),
    simulate_scenario(
        name="Vaccination",
        population=N,
        I0=I0,
        coverage=0.40,
        efficacy=0.80,
        beta_before=0.30,
        beta_after=0.30,
        intervention_day=30,
        gamma=gamma
    ),
    simulate_scenario(
        name="Vaccination + timed intervention",
        population=N,
        I0=I0,
        coverage=0.40,
        efficacy=0.80,
        beta_before=0.30,
        beta_after=0.12,
        intervention_day=30,
        gamma=gamma
    )
]

# Display outcome table
print("Scenario                         Vaccine immune   Peak I   Peak day   Outbreak infections")
print("------------------------------------------------------------------------------------------")

for result in scenarios:
    print(
        f"{result['name']:32s} "
        f"{result['vaccine_immune']:14.1f} "
        f"{result['peak_I']:8.2f} "
        f"{result['peak_time']:10.1f} "
        f"{result['outbreak_infections']:19.2f}"
    )

# Prepare transmission schedules for plotting
plot_time = scenarios[0]["time"]
beta_no_intervention = np.full_like(plot_time, 0.30)
beta_combined = np.where(plot_time < 30, 0.30, 0.12)

# Plot infectious curves and beta schedule
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

for result in scenarios:
    axes[0].plot(
        result["time"],
        result["I"],
        label=result["name"]
    )

axes[0].axvline(30, color="grey", linestyle=":")
axes[0].set_xlabel("Time (days)")
axes[0].set_ylabel("Infectious population")
axes[0].set_title("Infectious trajectories")
axes[0].grid(alpha=0.25)
axes[0].legend(fontsize=8)

axes[1].plot(
    plot_time,
    beta_no_intervention,
    label="No timed change"
)
axes[1].plot(
    plot_time,
    beta_combined,
    label="Timed intervention"
)
axes[1].axvline(30, color="grey", linestyle=":")
axes[1].set_xlabel("Time (days)")
axes[1].set_ylabel("Transmission parameter")
axes[1].set_title("Transmission schedules")
axes[1].set_ylim(0, 0.35)
axes[1].grid(alpha=0.25)
axes[1].legend()

fig.tight_layout()

Output

Run the code to see the result.

Understand the new code

Code	Technical meaning	Modelling use
`0 <= coverage <= 1`	A chained comparison requiring a proportion between 0 and 1.	Prevents impossible vaccination coverage.
`np.where(condition, a, b)`	Creates an array choosing a where a condition is true and b otherwise.	Creates the displayed piecewise transmission schedule.
`np.full_like(time, 0.30)`	Creates an array matching the shape of time and fills it with 0.30.	Displays constant transmission.
`max_step=0.5`	Prevents the solver from taking an internal step longer than half a day.	Encourages resolution around the abrupt day-30 parameter change.
`outbreak_infections`	Uses \(I(0)+S(0)-S(T)\).	Excludes people immune from vaccination who never enter the outbreak.

Biological interpretation

Vaccination lowers the number initially susceptible to infection. The timed intervention reduces the infection flow from day 30 onward. Both mechanisms can lower infectious burden, but they act through different parts of the model.

The effect of a timed intervention depends strongly on when it begins. If it begins after the natural infectious peak, it cannot reduce the earlier peak retrospectively.

Assumptions and limitations

Assumption	Consequence
Vaccine protection is immediate.	No delay between vaccination and immunity is represented.
Effective protection is all-or-nothing.	Successfully protected people are placed directly in the immune compartment.
Protection does not wane.	Vaccinated immune people remain immune throughout the simulation.
\(\beta\) changes instantly on day 30.	Behavioural change is represented as an abrupt step rather than a gradual transition.
Parameters are otherwise constant.	Changing adherence and other time-varying factors are omitted.

This is a scenario analysis under model assumptions. It does not establish that an observed intervention caused a particular real-world outcome.

Try changing intervention timing

In the combined scenario, change intervention_day=30 to 15 and then to 60. Also change the plotted schedule condition and vertical line to match.

Compare peak sizes and peak times.
Compare outbreak infections.
Determine whether day 60 occurs before or after the uncontrolled peak.
Explain why timing changes the outcome even when \(\beta_2\) is unchanged.