Parameter comparisons and sensitivity experiments

A model outcome depends on its parameter values. Parameter comparison shows how selected scenarios differ, while sensitivity analysis measures how strongly an outcome responds to small parameter changes around a stated baseline.

Biological scenario

For an SIR epidemic in a population of 1,000, the baseline parameters are \(\beta=0.30\) and \(\gamma=0.10\) per day.

We want to investigate how transmission and recovery affect the infectious peak and final outbreak size.

Start with an exactly understandable sensitivity

Before using simulation, consider the basic reproduction number for the simple SIR model:

\[\mathcal{R}_0=\frac{\beta}{\gamma}.\]

At \(\beta=0.30\) and \(\gamma=0.10\), we have \(\mathcal{R}_0=3\).

Increase transmission by 1%

If \(\beta\) increases by 1% while \(\gamma\) remains fixed, \(\mathcal{R}_0\) increases by exactly 1%. Its normalized sensitivity to \(\beta\) is +1.

Increase recovery by 1%

If \(\gamma\) increases by 1%, \(\mathcal{R}_0\) decreases by approximately 1% for a small change. Its normalized local sensitivity to \(\gamma\) is −1.

The signs are intuitive: greater transmission increases epidemic potential, while faster recovery decreases it. This simple example gives meaning to a sensitivity coefficient before we calculate one numerically.

Four complementary questions

QUESTION 1Exact formula sensitivity

How does \(\mathcal{R}_0=\beta/\gamma\) respond algebraically?

QUESTION 2One parameter at a time

How do full epidemic curves change across selected \(\beta\) values?

QUESTION 3Two-parameter response

How does peak size vary across combinations of \(\beta\) and \(\gamma\)?

QUESTION 4Intervention timing

How does changing the day of a transmission reduction alter peak and final size?

QUESTION 5Local coefficient

What is the derivative-like response near the baseline?

QUESTION 6Robustness

Does the conclusion change with the perturbation size or parameter range?

Design a fair comparison

Population1,000

Initial state\(S_0=990,I_0=10,R_0=0\)

Time interval0–200 days

Numerical settingssame solver and tolerances

Principle	Reason
State the baseline.	A parameter effect is interpreted relative to a reference model.
Change only intended quantities.	Otherwise the source of an outcome difference becomes unclear.
Use biologically meaningful ranges.	Extreme values may answer an irrelevant question.
Use identical numerical settings.	Differences should not arise from inconsistent computation.
Choose outcomes before running.	This avoids selecting only results that appear interesting afterwards.

One-at-a-time comparison

Choose several transmission values:

\[\beta\in\{0.20,0.25,0.30,0.35,0.40\},\]

while fixing \(\gamma=0.10\). This isolates the model response to \(\beta\) along that particular line through parameter space.

One-at-a-time comparison is easy to interpret, but it does not reveal interactions that occur when several parameters change together.

Two-parameter experiment

A grid evaluates every selected combination of \(\beta\) and \(\gamma\). Results can be stored in a two-dimensional array:

peak_grid = np.zeros(
    (len(gamma_values), len(beta_values))
)

The row index identifies a \(\gamma\) value and the column index identifies a \(\beta\) value.

for row, gamma in enumerate(gamma_values):
    for column, beta in enumerate(beta_values):
        result = simulate_sir(beta, gamma)
        peak_grid[row, column] = result["peak_I"]

A loop inside another loop is a nested loop.
enumerate(values) supplies both an index and its corresponding value.
peak_grid[row, column] selects one cell in the two-dimensional array.

Local normalized sensitivity

For parameter \(p\) and outcome \(Y\), we estimate a local sensitivity using a small increase and decrease around the baseline:

\[S_p^Y\approx\frac{Y(p+\Delta p)-Y(p-\Delta p)}{2\Delta p}\frac{p}{Y(p)}.\]

This is dimensionless. A value near 1 means that, locally, a 1% parameter increase is associated with approximately a 1% outcome increase. A negative value means the outcome moves in the opposite direction.

This interpretation is local and approximate. It applies near the chosen baseline and for the selected perturbation size.

Absolute versus normalized sensitivity

Measure	Approximate form	Interpretation
Absolute sensitivity	\(\Delta Y/\Delta p\)	Outcome-unit change per parameter-unit change; its numerical size depends on units.
Normalized sensitivity	\((\Delta Y/Y)/(\Delta p/p)\)	Percentage outcome change per percentage parameter change; dimensionless and easier to compare.

Program the finite-difference calculation

change_fraction = 0.05
delta_beta = change_fraction * beta_baseline

outcome_plus = model(beta_baseline + delta_beta)
outcome_minus = model(beta_baseline - delta_beta)

derivative = (
    outcome_plus - outcome_minus
) / (2 * delta_beta)

sensitivity = (
    derivative * beta_baseline / outcome_baseline
)

A 5% change is used on each side of the baseline. The central difference uses both sides and is generally more balanced than a one-sided comparison.

Run the complete multi-scenario sensitivity experiment

Interactive PythonSIR parameter experiments

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 timed_sir_rhs(
    t,
    y,
    beta_before,
    beta_after,
    intervention_day,
    gamma
):
    """Return SIR derivatives with a timed beta change."""
    
    if t < intervention_day:
        beta = beta_before
    else:
        beta = beta_after
    
    return sir_rhs(t, y, beta, gamma)

def outcomes_from_solution(solution):
    """Extract shared outcomes from a solved SIR trajectory."""
    
    S, I, R = solution.y
    peak_index = np.argmax(I)
    
    return {
        "time": solution.t,
        "S": S,
        "I": I,
        "R": R,
        "peak_I": I[peak_index],
        "peak_time": solution.t[peak_index],
        "final_size": 10 + 990 - S[-1]
    }

def simulate_sir(beta, gamma):
    """Solve the shared constant-parameter SIR scenario."""
    
    reporting_times = np.linspace(0, 200, 2001)
    
    solution = solve_ivp(
        sir_rhs,
        t_span=(0, 200),
        y0=[990, 10, 0],
        t_eval=reporting_times,
        args=(beta, gamma),
        rtol=1e-8,
        atol=1e-10
    )
    
    if not solution.success:
        raise RuntimeError(solution.message)
    
    return outcomes_from_solution(solution)

def simulate_intervention(intervention_day):
    """Solve an SIR scenario with beta reduced on a chosen day."""
    
    reporting_times = np.linspace(0, 200, 2001)
    
    solution = solve_ivp(
        timed_sir_rhs,
        t_span=(0, 200),
        y0=[990, 10, 0],
        t_eval=reporting_times,
        args=(0.30, 0.12, intervention_day, 0.10),
        rtol=1e-8,
        atol=1e-10,
        max_step=0.5
    )
    
    if not solution.success:
        raise RuntimeError(solution.message)
    
    return outcomes_from_solution(solution)

def normalized_sensitivity(
    parameter_name,
    beta_baseline,
    gamma_baseline,
    outcome_name,
    change_fraction=0.05
):
    """Estimate one local normalized sensitivity."""
    
    baseline = simulate_sir(
        beta_baseline,
        gamma_baseline
    )[outcome_name]
    
    if parameter_name == "beta":
        parameter_baseline = beta_baseline
        delta = change_fraction * parameter_baseline
        
        outcome_plus = simulate_sir(
            beta_baseline + delta,
            gamma_baseline
        )[outcome_name]
        
        outcome_minus = simulate_sir(
            beta_baseline - delta,
            gamma_baseline
        )[outcome_name]
    
    elif parameter_name == "gamma":
        parameter_baseline = gamma_baseline
        delta = change_fraction * parameter_baseline
        
        outcome_plus = simulate_sir(
            beta_baseline,
            gamma_baseline + delta
        )[outcome_name]
        
        outcome_minus = simulate_sir(
            beta_baseline,
            gamma_baseline - delta
        )[outcome_name]
    
    else:
        raise ValueError("parameter_name must be beta or gamma.")
    
    absolute_sensitivity = (
        outcome_plus - outcome_minus
    ) / (2 * delta)
    
    normalized = (
        absolute_sensitivity
        * parameter_baseline
        / baseline
    )
    
    return normalized

# Baseline and ranges
beta_baseline = 0.30
gamma_baseline = 0.10

beta_values = np.array([0.20, 0.25, 0.30, 0.35, 0.40])
gamma_values = np.array([0.07, 0.10, 0.13])

# Exact R0 sensitivity example
R0_baseline = beta_baseline / gamma_baseline
R0_beta_sensitivity = 1.0
R0_gamma_sensitivity = -1.0

# One-at-a-time beta results
beta_results = [
    simulate_sir(beta, gamma_baseline)
    for beta in beta_values
]

# Two-dimensional peak grid
peak_grid = np.zeros(
    (len(gamma_values), len(beta_values))
)

for row, gamma in enumerate(gamma_values):
    for column, beta in enumerate(beta_values):
        peak_grid[row, column] = simulate_sir(
            beta, gamma
        )["peak_I"]

# Local normalized sensitivities
sensitivity_labels = [
    "beta → peak",
    "gamma → peak",
    "beta → final",
    "gamma → final"
]

sensitivity_values = [
    normalized_sensitivity(
        "beta", beta_baseline, gamma_baseline, "peak_I"
    ),
    normalized_sensitivity(
        "gamma", beta_baseline, gamma_baseline, "peak_I"
    ),
    normalized_sensitivity(
        "beta", beta_baseline, gamma_baseline, "final_size"
    ),
    normalized_sensitivity(
        "gamma", beta_baseline, gamma_baseline, "final_size"
    )
]

# Extra scenario: intervention timing
intervention_days = np.array([10, 20, 30, 40, 50, 60, 70])
intervention_results = [
    simulate_intervention(day)
    for day in intervention_days
]

intervention_peaks = np.array([
    result["peak_I"]
    for result in intervention_results
])

intervention_final_sizes = np.array([
    result["final_size"]
    for result in intervention_results
])

# Display exact example
print("Exact sensitivity example")
print("---------------------------------------------")
print(f"Baseline R0:                  {R0_baseline:.2f}")
print(f"Normalized sensitivity to beta:  {R0_beta_sensitivity:+.1f}")
print(f"Normalized sensitivity to gamma: {R0_gamma_sensitivity:+.1f}")

# Display one-at-a-time table
print("\nOne-at-a-time beta experiment")
print("beta   Peak I   Peak day   Final epidemic size")
print("------------------------------------------------")

for beta, result in zip(beta_values, beta_results):
    print(
        f"{beta:4.2f}  {result['peak_I']:7.2f}  "
        f"{result['peak_time']:8.1f}  "
        f"{result['final_size']:19.2f}"
    )

# Display intervention-timing table
print("\nIntervention timing experiment")
print("Day   Peak I   Peak day   Final epidemic size")
print("------------------------------------------------")

for day, result in zip(
    intervention_days,
    intervention_results
):
    print(
        f"{day:3d}  {result['peak_I']:7.2f}  "
        f"{result['peak_time']:8.1f}  "
        f"{result['final_size']:19.2f}"
    )

# Display local sensitivities
print("\nLocal normalized sensitivities")
print("---------------------------------------------")

for label, value in zip(
    sensitivity_labels,
    sensitivity_values
):
    print(f"{label:20s} {value:8.3f}")

# Create four complementary plots
fig, axes = plt.subplots(2, 2, figsize=(12, 9))

# Panel 1: beta trajectories
for beta, result in zip(beta_values, beta_results):
    axes[0, 0].plot(
        result["time"],
        result["I"],
        label=f"beta={beta:.2f}"
    )

axes[0, 0].set_title("A. One-at-a-time beta comparison")
axes[0, 0].set_xlabel("Time (days)")
axes[0, 0].set_ylabel("Infectious population")
axes[0, 0].grid(alpha=0.25)
axes[0, 0].legend(fontsize=8)

# Panel 2: two-parameter heatmap
image = axes[0, 1].imshow(
    peak_grid,
    origin="lower",
    aspect="auto",
    cmap="viridis"
)

axes[0, 1].set_xticks(range(len(beta_values)))
axes[0, 1].set_xticklabels(
    [f"{value:.2f}" for value in beta_values]
)
axes[0, 1].set_yticks(range(len(gamma_values)))
axes[0, 1].set_yticklabels(
    [f"{value:.2f}" for value in gamma_values]
)
axes[0, 1].set_xlabel("beta")
axes[0, 1].set_ylabel("gamma")
axes[0, 1].set_title("B. Peak response surface")
fig.colorbar(image, ax=axes[0, 1], label="Peak I")

# Panel 3: local sensitivity
colours = [
    "seagreen" if value >= 0 else "indianred"
    for value in sensitivity_values
]

axes[1, 0].bar(
    sensitivity_labels,
    sensitivity_values,
    color=colours
)
axes[1, 0].axhline(0, color="black", linewidth=0.8)
axes[1, 0].set_title("C. Local normalized sensitivity")
axes[1, 0].set_ylabel("Sensitivity")
axes[1, 0].tick_params(axis="x", rotation=25)
axes[1, 0].grid(axis="y", alpha=0.25)

# Panel 4: extra intervention-timing scenario
axes[1, 1].plot(
    intervention_days,
    intervention_peaks,
    "o-",
    label="Peak infectious"
)
axes[1, 1].plot(
    intervention_days,
    intervention_final_sizes,
    "s--",
    label="Final epidemic size"
)
axes[1, 1].set_title("D. Sensitivity to intervention timing")
axes[1, 1].set_xlabel("Intervention day")
axes[1, 1].set_ylabel("People")
axes[1, 1].grid(alpha=0.25)
axes[1, 1].legend()

fig.tight_layout()

Output

Run the code to see the result.

New programming tools

Code	Purpose
`np.zeros((rows, columns))`	Creates a two-dimensional numerical array.
`enumerate(values)`	Returns each item together with its index.
`elif`	Tests another condition when the preceding if condition was false.
`zip(a, b)`	Pairs corresponding entries from two collections.
`[a if condition else b for value in values]`	A list comprehension that creates one result for every value, choosing between two alternatives.
`np.round(array, 2)`	Rounds displayed values to two decimal places without changing the original array.
`imshow(matrix)`	Displays a two-dimensional array as a colour image.
`fig.colorbar(...)`	Adds a scale explaining how colours correspond to numerical values.

Extra scenario: sensitivity to intervention timing

Suppose transmission falls from \(\beta=0.30\) to \(\beta=0.12\) on an intervention day \(T_L\). We now treat \(T_L\) as the input being varied.

\[T_L\in\{10,20,30,40,50,60,70\}\quad\text{days}.\]

This scenario is useful because timing can have a nonlinear effect. Delaying an early intervention may strongly increase the peak, while delaying an intervention until after the uncontrolled peak may make little further difference to that peak.

Input	Outputs	Question
Intervention day \(T_L\)	Peak infectious population and final epidemic size	How rapidly is benefit lost as intervention is delayed?

Interpret the results carefully

Positive sensitivity

The outcome increases locally as the parameter increases. Transmission commonly has a positive effect on peak size near an epidemic baseline.

Negative sensitivity

The outcome decreases locally as the parameter increases. Faster recovery commonly reduces infectious burden.

The magnitude ranks local responsiveness only for the chosen model, baseline, outcome and perturbation. A small local sensitivity does not prove that a parameter is unimportant everywhere.

Sensitivity is not uncertainty

Sensitivity analysis	Uncertainty analysis
Asks how outputs respond when inputs change.	Asks how uncertain inputs produce uncertainty in outputs.
Can use deliberately chosen parameter changes.	Requires assumptions or evidence about plausible parameter distributions or ranges.
Does not assign probabilities to parameter values.	May propagate probability distributions through the model.

Neither sensitivity nor scenario comparison establishes that changing a real-world factor will cause exactly the modelled outcome. Conclusions remain conditional on the model structure and parameter interpretation.

Try changing the perturbation size

First change change_fraction=0.05 to 0.01 and then to 0.20. Next extend intervention_days to include days 80 and 100.

Compare the sensitivity estimates.
Identify which estimates remain stable.
Explain why a 20% change is less “local” than a 1% change.
Discuss whether one derivative-like number can summarise a strongly nonlinear response.