The Gillespie algorithm draws a waiting time with total event rate , then chooses infection with probability and recovery otherwise. Each event conserves , and is absorbing in this closed stochastic SIR model.
Choose a horizon of 100 days, which is longer than all extinction times in the runs reported here. This independent implementation also avoids executing a sampled event beyond the horizon:
import numpy as np
import matplotlib.pyplot as plt
EPIDEMIC_RNG = np.random.default_rng(2345)
def sir_sim(t_max,X0,beta,gamma,rng=None):
if rng is None:rng=EPIDEMIC_RNG
state=np.array(X0,dtype=int);population=state.sum();time=0.
times=[time];states=[state.copy()]
while time<t_max and state[1]>0:
infection=beta*state[0]*state[1]/population;recovery=gamma*state[1]
rate=infection+recovery
if rate==0:
times.append(float(t_max));states.append(state.copy());break
next_time=time+rng.exponential(1/rate)
if next_time>t_max:
times.append(float(t_max));states.append(state.copy());break
if rng.random()<infection/rate:state+=np.array([-1,1,0])
else:state+=np.array([0,-1,1])
time=next_time;times.append(time);states.append(state.copy())
return np.array(times),np.array(states)
rng = np.random.default_rng(2345)
trajectories = [sir_sim(100, [99,1,0], 1, .25, rng) for _ in range(10)]
for times, states in trajectories:
plt.step(times, states[:,1], where="post", alpha=.8)
plt.xlabel("Time (days)")
plt.ylabel("Infected individuals")
plt.show()A step plot represents the piecewise constant infected population between events. The trajectories show both substantial outbreaks and early epidemic extinction despite the large basic reproduction number .
Ten original stochastic SIR trajectories, with frequency-dependent transmission and one initial infection
. Details. 
One simulation contributes its largest observed infected population and its last recorded time. Repeating this operation produces one row per realization:
def sir_sim_stats(n_sim,t_max,X0,beta,gamma,rng=None):
result=np.empty((n_sim,2))
for run in range(n_sim):
t,s=sir_sim(t_max,X0,beta,gamma,rng)
result[run]=s[:,1].max(),t[-1]
return resultThe maximum must be taken over the infected column,
states[:,1], rather than across all compartments. With the implementation in part (a), the final time is the epidemic extinction time if extinction occurs first, and otherwise the specified horizon; an outbreak still active at that horizon has a censored extinction time.Continue the random-number stream used for the ten trajectories, then generate 1000 further independent realizations with the same model parameters. Plot each statistic in its own histogram:
stats = sir_sim_stats(1000, 100, [99,1,0], 1, .25, rng)
fig, axes = plt.subplots(1, 2, figsize=(9, 4), layout="constrained")
axes[0].hist(stats[:,0], bins=np.arange(.5,65.6,2))
axes[1].hist(stats[:,1], bins=np.arange(0,61,2))
axes[0].set(xlabel="Peak infected population", ylabel="Number of outbreaks")
axes[1].set(xlabel="Final recorded time (days)", ylabel="Number of outbreaks")
plt.show()Both histograms have a group of early failures separated from a group of large outbreaks. For this reproducible run the mean peak is 33.73 individuals, and the mean final time is 22.165 days; the largest final time is 56.916 days, so all 1000 runs reached epidemic extinction before the 100-day horizon.
Peak size and extinction time for 1000 original stochastic SIR realizations, with operational failure cutoffs
. Details. 
“Failed to take off” needs an operational definition; the two separated groups in the histograms suggest a low peak or an early end. For example, classify a run as a failure if its peak is below ten infected individuals, or alternatively if its extinction time is below ten days:
failed_by_peak = stats[:,0] < 10
failed_by_time = stats[:,1] < 10
print(failed_by_peak.mean(), failed_by_time.mean())Both definitions give in this run. Here they also classify the same 255 realizations, although agreement is not guaranteed for other simulations or cutoffs. These are empirical failure fractions for the stated rules, not exact model constants.
The early branching process approximation provides a useful check: while susceptible depletion is negligible, gives eventual extinction probability about from one initial infection. Allowing for the initial susceptible fraction gives a rate ratio ; both approximations are compatible with this simulation, but neither defines the finite-population cutoff exactly.
Articles by others on the same topic
There are currently no matching articles.