Epidemic threshold

Final size equation

Definition

We define the final size $z$ as the proportion of cumulative infected individuals in the entire population for $t \to \infty$ , which is

z = \frac{R_\infty}{N}

where $R_\infty$ is the number of recovered individual at $t \to \infty$ and $N$ is the size of the entire population.

The final size $z$ follows the relationship below:

z = 1 - \exp (- R_0 z)

where $R_0$ is the basic reproduction number. This equation is known as the final size equation.

Derivation

Let us derive the final equation. We first start with the set of SIR differential equations.

\begin{equation*} \left\{ \begin{aligned} \frac{dS(t)}{dt} &= -\beta \frac{I(t)}{N}S(t) \\ \frac{dI(t)}{dt} &= \beta \frac{I(t)}{N}S(t) - \mu I(t) \\ \frac{dR(t)}{dt} &= \mu I(t) \end{aligned} \right. \end{equation*}

Notice that $I$ appears in both $\frac{dS(t)}{dt}$ and $\frac{dR(t)}{dt}$ . By eliminating $I$ , we form a differential relation between $R$ and $S$ . Dividing $\frac{dR(t)}{dt}$ by $\frac{dS(t)}{dt}$ , we get

\begin{aligned} \frac{dR(t)}{dS(t)} &= \frac{\frac{dR(t)}{dt}}{\frac{dS(t)}{dt}} \\ &= \frac{\mu I(t)}{-\beta \frac{I(t)}{N}S(t)} \\ &= -\frac{\mu N}{\beta S(t)}. \end{aligned}

Since the basic reproduction number is $R_0 = \frac{\beta}{\mu}$ ,

\frac{dR(t)}{dS(t)} = -\frac{N}{R_0 S(t)},

Rearranging the equation, we get

dR(t)= -\frac{N}{R_0} \frac{dS(t)}{S}.

By integrating both sides with t from $t=0$ to $t \to \infty$ ,

\int_{R(0)}^{R(\infty)} dR(u) = -\frac{N}{R_0} \int_{S(0)}^{S(\infty)} \frac{dS(u)}{S}

Note that following relationship holds:

At $t=0$ , $R(0)=0$ and $S(0)=S_0$
At $t \to \infty$ , $R(\infty)=R_\infty$ and $S(\infty)=S_\infty$

We get

\begin{aligned} R_\infty &= -\frac{N}{R_0} \bigg[ \ln S \bigg] _{S_0}^{S_\infty} \\ &= -\frac{N}{R_0} \left( \ln S_\infty - \ln S_0 \right) \\ &= -\frac{N}{R_0} \ln \frac{S_\infty}{S_0} \\ &= \frac{N}{R_0} \ln \frac{S_0}{S_\infty} \end{aligned}

Since the total population is conserved, we have:

N = S_\infty + I(\infty) + R_\infty

At the end of the epidemic, $I(\infty) \approx 0$ , so:

N = S_\infty + R_\infty \\ S_\infty = N - R_\infty

Substitute this into the previous equation:

R_\infty = \frac{N}{R_0} \ln \frac{S_0}{N - R_\infty}

Assuming that almost the entire population is susceptible at $t=0$ , i.e., $S_0 \approx N$ ,

R_\infty = \frac{N}{R_0} \ln \frac{N}{N - R_\infty}

Dividing through by $N$ gives:

\frac{R_\infty}{N} = \frac{1}{R_0} \ln \frac{1}{1 - \frac{R_\infty}{N}}

Now, let $z = \frac{R_\infty}{N}$ which is the fraction of the cumulative population that were infected at $t \to \infty$ . The equation becomes

\begin{aligned} z &= \frac{1}{R_0} \ln \frac{1}{1 - z} \\ &= \frac{1}{R_0} \bigg( \ln 1 - \ln (1-z) \bigg) \\ &= - \frac{1}{R_0} \ln (1-z) \\ \end{aligned}

Multiplying by $R_0$ and taking exponent on both sides, we get

\begin{aligned} R_0 z &= - \ln (1-z) \\ \exp( - R_0 z ) &= 1-z \end{aligned}

Rewriting the equation, we finally get the final size equation for the SIR model as:

z = 1 - \exp( - R_0 z ) \tag{1}

Now, let $s = \frac{s_\infty}{N}$ which is the fraction of the susceptible population at $t \to \infty$ , where $s = 1 - z$ holds. The final size equation can also be written as:

1-s = 1 - \exp( - R_0 (1-s) )

s = \exp( - R_0 (1-s) ) \tag{2}

Threshold in theory

Now, let's solve the final size equation for $z$ , both analytically and numerically.

# Define functions to compute Z numerically and analytically

from scipy.optimize import root_scalar
from scipy.special import lambertw

def final_size_equation(z, R0):
    """
    Final size equation in the SIR model.
    This function represents the equation z = 1 - exp(-R0 * z).
    """
    # LHS - RHS = 0
    return z - (1 - np.exp(-R0 * z)) 

def calc_Z_numerically(R0):
    """
    Compute z numerically using the root-finding method.
    This function uses the Brent's method to find the root of the final size equation.
    Since the equation has multiple solutions for z in R≥1, we set the range for the solution to [1e-6, 1.0] to find the non-trivial solution and avoid trivial solutions (z=0).

    Parameters:
    R0 (float or numpy.ndarray): The basic reproduction number(s).

    Returns:
    float: The computed value of Z for the given R0.
    """
    if R0 <= 1:
        return 0.0
    result = root_scalar(
        final_size_equation,
        args=(R0,),
        bracket=[1e-6, 1.0], # Set range for solution. 
        method='brentq'
    )
    return result.root if result.converged else np.nan

def calc_Z_analytically(R0):
    """
    Compute Z analytically using Lambert W function.
    For R0 > 0, Z = 1 + (1/R0) * W(-R0 * exp(-R0)), taking the principal branch.

    Parameters:
    R0 (float or numpy.ndarray): The basic reproduction number(s).

    Returns:
    float or numpy.ndarray: The corresponding value(s) of Z.
    """
    if R0 <= 0:
        return 0.0
    return 1 + (1 / R0) * lambertw(-R0 * np.exp(-R0)).real

# Range for R0
R0_values = np.linspace(0.0, 3.0, 200)

# Calculate Z values for both numerical and analytical methods
Z_values_numerical = [calc_Z_numerically(R0) for R0 in R0_values] # Z values calculated numerically for each R0
Z_values_analytical = [calc_Z_analytically(R0) for R0 in R0_values] # Z values calculated analytically for each R0

The figure below compares the the numerical and analytical solutions of $z$ value. We observe that the final size $z$ is zero for $R_0<1$ , and experiences a phase transision at $R_0=1$ , showing the threshold property.

In other words, $R_0 ≥ 1$ has to hold for an outbreak to take off.

# Compare the results from numerical and analytical methods
fig, ax = plt.subplots(1,1, figsize=(4,3))

ax.plot(
    R0_values,
    Z_values_analytical,
    color=colors[2],
    label='Numerical'
)
ax.plot(
    R0_values,
    Z_values_analytical,
    ':',
    color='black',
    label='Analytical'
)

ax.set_xlabel(r'$R_0$')
ax.set_ylabel(r'$z$')
ax.set_title(r'Epidemic threshold', loc='left')

ax.legend()

plt.close()

show_fig_svg(fig)

Now, recall that we had to set a range of solution for $z$ in the calc_Z_numerically(R0) function above. This is because there is a trivial solution (z=0) besides the non-trivial solution that we want to solve for. Here, we'll make an animation plotting two functions below to observe the solutions (point of intersection) as we change the value of $R_0$ .

\begin{equation*} \left\{ \begin{aligned} y &= z \\ y &= 1 - \exp(-R_0 z) \end{aligned} \right. \end{equation*}

# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/
colors_cbfp = ["#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"]

import matplotlib.animation as animation

Z_values_ani = np.linspace(-0.1, 1.05, 500) # Z-value for animation plot

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

# Left: Intersection of two functions
line1, = ax1.plot([], [], label=r"$y = z$", color=colors_cbfp[5])
line2, = ax1.plot([], [], label=r"$y = 1 - \exp(- R_0 z)$", color=colors_cbfp[2])
point_trivial, = ax1.plot([], [], 'o', color=colors_cbfp[1], label="Trivial ($z=0$)")
point_nontrivial, = ax1.plot([], [], 'o', color=colors_cbfp[7], label=r"Nontrivial ($\hat{z}$)")
title = ax1.text(0, 1.02, "", transform=ax1.transAxes, ha="left")
ax1.set_xlim(-0.05, 1.05)
ax1.set_ylim(-0.05, 1.05)
ax1.set_xlabel("$z$")
ax1.set_ylabel("$y$")
ax1.legend(fontsize=8)

# Right: Z vs R0
ax2.plot(R0_values, Z_values_numerical, '-', color="black")
point_current, = ax2.plot([], [], 'o', color=colors_cbfp[7], label=r"$\hat{z}$ (current)")
ax2.set_xlim(-0.05, 3.05)
ax2.set_ylim(-0.05, 1.05)
ax2.set_xlabel(r"$R_0$")
ax2.set_ylabel(r"$z$")
ax2.legend(fontsize=8, loc='upper left')

def init_finalsize_animation():
    line1.set_data([], [])
    line2.set_data([], [])
    point_trivial.set_data([], [])
    point_nontrivial.set_data([], [])
    point_current.set_data([], [])
    title.set_text("")
    return line1, line2, point_trivial, point_nontrivial, point_current, title

def animate_finalsize(i):
    R0 = R0_values[i]
    y1 =Z_values_ani 
    y2 = 1 - np.exp(-R0 * Z_values_ani)
    Z_star = Z_values_numerical[i]

    line1.set_data(Z_values_ani, y1)
    line2.set_data(Z_values_ani, y2)
    point_trivial.set_data([0], [0])
    point_nontrivial.set_data([Z_star], [Z_star])
    point_current.set_data([R0], [Z_star])
    title.set_text(f"$R_0$ = {R0:.2f}")

    return line1, line2, point_trivial, point_nontrivial, point_current, title

ani = animation.FuncAnimation(
    fig, animate_finalsize, frames=len(R0_values), init_func=init_finalsize_animation,
    blit=True, interval=100, repeat=False
)

plt.close(fig)

# We use this custom function to have the animation embedded as mp4. This allows the extraction of videos when we build the document.
show_animation_video_base64(ani)

# Use this when you don't need to convert to Markdown
# HTML(ani.to_jshtml())

Threshold in simulation

beta_list = np.arange(0,0.4,0.005)

data_final_size_determ = []

for beta in beta_list:
    _, _, R_dict = SIR_model(
        S_initial=params["S_initial"],
        I_initial=params["I_initial"],
        R_initial=params["R_initial"],
        n_steps=params["n_steps"],
        beta=beta,
        gamma=params["gamma"]
    )
    final_size = list(R_dict.values())[-1]
    data_final_size_determ.append(
        {
            'final_size': final_size,
            'beta': beta,
            'gamma': params["gamma"]
        }
    )
    # final_size_dict_determ[beta] = final_size / params["N"]

df_final_size_determ = pd.DataFrame(data_final_size_determ)
df_final_size_determ['final_size'] = df_final_size_determ['final_size'] / params["N"]

df_final_size_determ['R0'] = df_final_size_determ.apply(
    lambda d:
        calc_R0_homogeneous(beta = d['beta'], gamma = d['gamma']),
    axis=1
)

# df_final_size_determ['R0'] = df_final_size_determ['beta'] / df_final_size_determ['gamma']

df_final_size_determ

	final_size	beta	gamma	R0
0	0.000400	0.000	0.1	0.00
1	0.000421	0.005	0.1	0.05
2	0.000444	0.010	0.1	0.10
3	0.000471	0.015	0.1	0.15
4	0.000500	0.020	0.1	0.20
...	...	...	...	...
75	0.979013	0.375	0.1	3.75
76	0.980273	0.380	0.1	3.80
77	0.981459	0.385	0.1	3.85
78	0.982573	0.390	0.1	3.90
79	0.983622	0.395	0.1	3.95

80 rows × 4 columns

data_final_size_stoc = []

# N = 50000
for beta in tqdm(beta_list):
    for it in range(params["n_stoc_sim"]):
        _, _, R_dict = SIR_model_stochastic(
            S_initial=params["S_initial"],
            I_initial=params["I_initial"],
            R_initial=params["R_initial"],
            n_steps=params["n_steps"],
            beta=beta,
            gamma=params["gamma"]
        )
        final_size = list(R_dict.values())[-1]
        data_final_size_stoc.append(
            {
                'iter': it,
                'final_size': final_size,
                'beta': beta,
                'gamma': params["gamma"]
            }
        )
    # final_size_dict_stoc[beta] = final_size / params["N"]

0%|          | 0/80 [00:00<?, ?it/s]

100%|██████████| 80/80 [00:01<00:00, 48.82it/s]

df_final_size_stoc = pd.DataFrame(data_final_size_stoc)
df_final_size_stoc['final_size'] = df_final_size_stoc['final_size'] / params["N"]

df_final_size_stoc['R0'] = df_final_size_stoc.apply(
    lambda d:
        calc_R0_homogeneous(beta = d['beta'], gamma = d['gamma']),
    axis=1
)

summary_df_final_size_stoc = df_final_size_stoc.groupby('R0')['final_size'].agg(['mean', q25, q75])

Comparison

fig, ax = plt.subplots(1,1, figsize=(4,3))

# Analytical solution when n-> \infty
ax.plot(R0_values, Z_values_analytical, color=colors[0], label=r'Analytical $( n \rightarrow \infty)$', zorder=-1)

# Stochastic
ax.fill_between(
    summary_df_final_size_stoc.index,
    summary_df_final_size_stoc['q25'],
    summary_df_final_size_stoc['q75'],
    color=colors[2],
    alpha=0.2,
)
ax.plot(
    summary_df_final_size_stoc.index,
    summary_df_final_size_stoc['mean'],
    color=colors[2],
    label="Stochastic",
)

# Deterministic
ax.plot(
    df_final_size_determ['R0'],
    df_final_size_determ['final_size'],
    color="#2e2e2e",
    label="Deterministic",
    linestyle=(0,(1,2))
)

# R0=1 line
ax.axvline(x=1, color='gray', linestyle='--', zorder=-10)


ax.set_title("Epidemic threshold", loc='left')
ax.set_xlabel("$R_0$")
ax.set_ylabel(r"$z$")
ax.legend()

ax.set_xlim(-0.1,3.1)

fig.tight_layout()
# plt.savefig('fig/fig_name.pdf', bbox_inches='tight')
plt.close()

Final size equation​

Definition​

Derivation​

Threshold in theory​

Threshold in simulation​

Comparison​