Critical Community Size (CCS)¶

import numpy as np
import matplotlib.pyplot as plt

from laser.core import PropertySet
from laser.core.utils import grid
from laser.generic.utils import ValuesMap
from laser.cohorts import Model, Campaign, Intervention
from laser.cohorts.vitaldynamics import NonDiseaseMortality, ConstantPopBirths
import laser.cohorts.SIR as SIR


class Importation(Intervention):
    """Add a fixed count of infectious individuals to targeted nodes each scheduled tick."""

    @property
    def states(self):
        return []

    @property
    def properties(self):
        return []

    def execute(self, tick, who, where, params, notes):
        count   = int(params.get("count", 1))
        target  = where if where is not None else list(range(len(self.model.scenario)))
        I_idx   = self.model.states.get_state_index("I")
        if I_idx is None:
            return
        for node in target:
            self.model.states[tick + 1][I_idx, node] += count


Campaign.register(Importation)

def persistence_fraction(N, nticks, r_mortality, beta, gamma, n_reps):
    """Estimate the fraction of endemic time with I > 0 across n_reps replicates.

    Runs n_reps independent SIR simulations for a single patch of size N with
    vital dynamics and annual importation of 1 case.  After a 1-year burn-in,
    measures the fraction of ticks with at least one infectious individual.

    Args:
        N: Population size.
        nticks: Simulation length in ticks (days).
        r_mortality: Daily non-disease mortality rate.
        beta: Transmission rate (per day).
        gamma: Recovery rate (per day).
        n_reps: Number of independent replicates.

    Returns:
        Mean persistence fraction across replicates (0 = always extinct, 1 = always present).
    """
    fracs         = []
    I0            = max(1, N // 1000)
    import_period = 365
    when_ticks    = list(range(0, nticks, import_period))

    for _ in range(n_reps):
        scenario = grid(M=1, N=1)
        scenario["S"] = N - I0
        scenario["I"] = I0
        scenario["R"] = 0

        params = PropertySet({"nticks": nticks})
        model  = Model(scenario, params)

        schedule = [{"who": "*", "what": "Importation", "when": when_ticks,
                     "where": [0], "parameters": {"count": 1}, "notes": ""}]
        campaign = Campaign(model, schedule)

        betas        = ValuesMap.from_scalar(beta,  nticks, 1)
        r_recoveries = ValuesMap.from_scalar(gamma, nticks, 1)

        model.components = [
            SIR.Susceptible(model),
            SIR.Infectious(model, r_recovery=r_recoveries),
            SIR.Recovered(model),
            NonDiseaseMortality(model, r_mortality=r_mortality),
            ConstantPopBirths(model),
            campaign,
            SIR.Transmission(model, beta=betas),
        ]
        model.run()

        I_ts = model.states.I[:, 0]

        # Discard first year as burn-in
        burn = 365
        fracs.append((I_ts[burn:] > 0).mean())

    return np.mean(fracs)

# Find the approximate CCS as the population size where persistence crosses 0.5
cross_idx = np.where(persistence >= 0.5)[0]
CCS_obs   = pop_sizes[cross_idx[0]] if len(cross_idx) > 0 else None

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(pop_sizes, persistence, "o-", color="steelblue", markersize=5,
            label="Persistence fraction")
ax.axhline(0.5, color="firebrick", linestyle="--", linewidth=1,
           label="Persistence = 0.5 (CCS threshold)")

if CCS_obs is not None:
    ax.axvline(CCS_obs, color="orange", linestyle=":", linewidth=1.5,
               label=f"Estimated CCS ≈ {CCS_obs:,}")

ax.set_xlabel("Population size N")
ax.set_ylabel("Fraction of time with I > 0  (after 1-year burn-in)")
ax.set_title(f"Critical Community Size  (R₀={beta/gamma:.0f}, {n_replicates} replicates)")
ax.set_ylim(-0.05, 1.05)
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

if CCS_obs:
    print(f"Observed CCS (persistence ≥ 0.5): N ≈ {CCS_obs:,}")

# Rough theoretical CCS comparison
# One widely cited approximation (Keeling & Grenfell 1997) is:
#   CCS ~ N* such that the inter-epidemic susceptible replenishment
#   barely refills the epidemic threshold between outbreaks.
# For measles-like parameters (R0 ~ 5, CDR ~ 20/1000) the empirical CCS is ~300,000.
# A simple heuristic: CCS ~ (birth rate) / (extinction probability per wave)
# We can compute the endemic I* and see how quickly it drives to zero.

R0     = beta / gamma
mu     = r_mortality          # birth ≈ death in stable population
I_star = mu * (R0 - 1) / beta  # endemic equilibrium (fraction of N)

print(f"R0           = {R0:.1f}")
print(f"Theoretical endemic I*/N  = {I_star:.5f}")
print()
print("For a population of N, the expected endemic infectious count is:")
for N_ex in [10_000, 50_000, 100_000, 300_000, 500_000, 1_000_000]:
    print(f"  N={N_ex:>8,}  →  I* ≈ {I_star * N_ex:.1f} individuals")

print()
print("When I* < ~1, stochastic extinction becomes likely between waves.")
N_extinction = 1.0 / I_star
print(f"I*/N · N = 1  at  N = {N_extinction:,.0f}  (rough heuristic CCS)")
if CCS_obs:
    print(f"Observed CCS from simulation: ~{CCS_obs:,}")