Peak size, peak time and final epidemic size

A model trajectory becomes useful when it is converted into clearly defined biological outcomes. This lesson calculates the epidemic peak and final size while explaining exactly what each quantity does—and does not—mean.

Biological scenario

A deterministic SIR simulation begins with 990 susceptible, 10 infectious and no recovered people. We solve it for 200 days using \(\beta=0.3\) and \(\gamma=0.1\) per day.

We want to answer three questions: How high does the infectious population rise? When does that maximum occur? How many people experience infection during the outbreak?

Define the outcomes first

OUTCOME 1Peak size

The largest value reached by the infectious population.

\[I_{\max}=\max_t I(t)\]

OUTCOME 2Peak time

The time at which the infectious maximum occurs.

\[t_{\mathrm{peak}}=\operatorname*{arg\,max}_t I(t)\]

OUTCOME 3Final epidemic size

The total number infected during the outbreak under a stated definition and sufficiently late end time.

Peak size and peak time in stored data

1 SEARCHnp.argmax(I)

Find the index of the largest stored infectious value.

2 PEAK SIZEI[peak_index]

Use that index in the infectious array.

3 PEAK TIMEtime[peak_index]

Use the same index in the time array.

4 PERCENTAGE100 * peak_I / N

Express the peak relative to population size.

The calculated peak time is the time of the largest stored value. A finer reporting grid normally estimates the continuous peak time more precisely.

“Final epidemic size” requires a definition

Suppose the simulation ends at time \(T\). Several related quantities can be reported:

New infections after time 0\[S(0)-S(T)\]

Initially infectious plus new infections\[I(0)+S(0)-S(T)\]

Recovered at the end\[R(T)\]

When are the last two equal?

If \(R(0)=0\) and the epidemic has effectively ended so that \(I(T)\approx0\), then

\[I(0)+S(0)-S(T)\approx R(T).\]

If initially recovered people are present, or infectious people remain at \(T\), simply calling \(R(T)\) the final epidemic size can be misleading.

Why the ending time matters

The final array entry is only the outcome at the chosen simulation end:

final_S = S[-1]
final_I = I[-1]
final_R = R[-1]

Index -1 means the last stored entry. It does not automatically mean that the epidemic is biologically finished.

We therefore check whether final_I is below a small biological threshold. This threshold is an interpretation choice, not a proof that the mathematical solution has reached exactly zero.

Build an outcome function

def calculate_outcomes(time, S, I, R):
    peak_index = np.argmax(I)

    peak_I = I[peak_index]
    peak_time = time[peak_index]

    final_S = S[-1]
    final_I = I[-1]
    final_R = R[-1]

    new_infections = S[0] - final_S
    total_infected = I[0] + new_infections

    return {
        "peak_I": peak_I,
        "peak_time": peak_time,
        "final_S": final_S,
        "final_I": final_I,
        "final_R": final_R,
        "new_infections": new_infections,
        "total_infected": total_infected
    }

The function uses the initial and final compartment values explicitly, making the chosen definitions visible.

Run the model and calculate the outcomes

Interactive PythonSIR epidemic outcomes

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

def sir_rhs(t, y, beta, gamma):
    """Return deterministic SIR derivatives."""
    
    S, I, R = y
    N = S + I + R
    
    infection_flow = beta * S * I / N
    recovery_flow = gamma * I
    
    return [
        -infection_flow,
        infection_flow - recovery_flow,
        recovery_flow
    ]

def calculate_outcomes(time, S, I, R):
    """Calculate peak and final outcomes from an SIR solution."""
    
    peak_index = np.argmax(I)
    
    peak_I = I[peak_index]
    peak_time = time[peak_index]
    
    final_S = S[-1]
    final_I = I[-1]
    final_R = R[-1]
    
    new_infections = S[0] - final_S
    total_infected = I[0] + new_infections
    
    return {
        "peak_index": peak_index,
        "peak_I": peak_I,
        "peak_time": peak_time,
        "final_S": final_S,
        "final_I": final_I,
        "final_R": final_R,
        "new_infections": new_infections,
        "total_infected": total_infected
    }

# Initial state and parameters
S0 = 990
I0 = 10
R0 = 0
N = S0 + I0 + R0

beta = 0.3
gamma = 0.1

# Fine reporting grid for outcome estimation
end_time = 200
reporting_times = np.linspace(0, end_time, 4001)

# Solve the SIR model
solution = solve_ivp(
    sir_rhs,
    t_span=(0, end_time),
    y0=[S0, I0, R0],
    t_eval=reporting_times,
    args=(beta, gamma),
    rtol=1e-9,
    atol=1e-11
)

if not solution.success:
    raise RuntimeError(solution.message)

S, I, R = solution.y
outcomes = calculate_outcomes(solution.t, S, I, R)

# Determine whether the epidemic is effectively finished
infectious_threshold = 0.01
finished = outcomes["final_I"] < infectious_threshold

# Display outcomes
print("SIR epidemic outcomes")
print("------------------------------------------------")
print(f"Population:                    {N}")
print(f"Peak infectious:               {outcomes['peak_I']:.2f}")
print(f"Peak percentage:               {100 * outcomes['peak_I'] / N:.2f}%")
print(f"Peak time:                     {outcomes['peak_time']:.2f} days")
print(f"Final susceptible S(T):        {outcomes['final_S']:.2f}")
print(f"Final infectious I(T):         {outcomes['final_I']:.6f}")
print(f"Final recovered R(T):          {outcomes['final_R']:.2f}")
print(f"New infections after t=0:      {outcomes['new_infections']:.2f}")
print(f"Total infected incl. initial:  {outcomes['total_infected']:.2f}")
print(f"Below infectious threshold:    {finished}")

# Check the relationship when R(0)=0
difference = abs(
    outcomes["total_infected"] - outcomes["final_R"]
)
print(f"Difference from final R:       {difference:.6f}")

# Plot and mark the peak
plt.figure(figsize=(8, 4.8))
plt.plot(solution.t, S, label="Susceptible")
plt.plot(solution.t, I, label="Infectious")
plt.plot(solution.t, R, label="Recovered")

plt.scatter(
    outcomes["peak_time"],
    outcomes["peak_I"],
    color="black",
    zorder=5,
    label="Infectious peak"
)
plt.axvline(
    outcomes["peak_time"],
    color="grey",
    linestyle=":"
)

plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("SIR trajectory and infectious peak")
plt.grid(alpha=0.25)
plt.legend()
plt.tight_layout()

Output

Run the code to see the result.

Interpret each reported quantity

Output	Question answered	Important qualification
Peak infectious	What is the largest simultaneous infectious population?	It is not the total number infected over the outbreak.
Peak time	When is the stored infectious maximum reached?	Precision depends partly on reporting-time resolution.
New infections	How many initially susceptible people become infected?	Calculated as \(S(0)-S(T)\).
Total infected including initial	How many people participate in the outbreak as infectious cases, including those infectious at time 0?	Calculated as \(I(0)+S(0)-S(T)\) for this SIR setup.
Final recovered	How many people are in \(R\) at time \(T\)?	Equals total outbreak infections only under the stated initial and end conditions.

Peak prevalence is not incidence

Peak infectious population

The maximum number infectious at the same time. This is a prevalence-type quantity.

New infections through time

The flow \(\beta SI/N\) measures the instantaneous rate at which new infections occur. This is related to incidence.

A high daily infection flow and a high simultaneous infectious population are related but are not the same outcome.

Biological interpretation

The peak size indicates the largest simultaneous infectious burden, which may be relevant to healthcare demand and workforce absence. Peak time indicates when that burden is expected under the model.

Final epidemic size describes cumulative outbreak involvement rather than simultaneous burden. Its interpretation must state whether initially infectious or initially recovered people are included.

Important reporting rules

Define every outcome before calculating it.
Report whether outcomes are counts or percentages.
State the simulation end time.
Check that infectious prevalence is sufficiently small before calling an outcome “final.”
State whether initial infectious people are included in final epidemic size.
Do not confuse peak prevalence with cumulative infections or incidence.
Remember that deterministic outputs need not be whole numbers.

Try ending the simulation too early

Change end_time = 200 to end_time = 40 and change the reporting grid accordingly:

reporting_times = np.linspace(0, end_time, 801)

Check whether the infectious threshold condition is satisfied.
Compare final_R with total_infected.
Explain why the day-40 result should not be called the final epidemic size.
Restore the 200-day simulation and compare again.