Programming a DTMC SIS model

The SIS model describes an infection that does not give lasting immunity. Recovery returns an infectious person to the susceptible compartment, so the same person may later become infected again.

The new biological idea: reinfection

Susceptible\(S\)

⇄

Infectious\(I\)

\[ S\xrightarrow{\text{infection}}I, \qquad I\xrightarrow{\text{recovery}}S. \]

Unlike the SIR model, SIS has no removed compartment. Recovery does not create permanent immunity. SIS models can be useful for infections for which immunity is absent, short-lived or ignored at the chosen modelling scale.

Closed population and a one-dimensional state

Assume a closed population of constant size \(N\):

\[ S_n+I_n=N. \]

Knowing \(I_n\) automatically determines the susceptible count:

\[ S_n=N-I_n. \]

Therefore the DTMC can use the single integer state \(X_n=I_n\) with state space:

\[ \mathcal S=\{0,1,2,\ldots,N\}. \]

This reduction is not an approximation. It follows exactly from population conservation.

Possible movements from state (i)

Suppose \(I_n=i\). Under the one-event-per-step assumption, the next infectious count can be only:

Infection: upward movement

\[i\longrightarrow i+1.\]

One susceptible person becomes infectious.

Recovery: downward movement

\[i\longrightarrow i-1.\]

One infectious person becomes susceptible again.

No change: stay

\[i\longrightarrow i.\]

No infection or recovery occurs in the step.

A chain with only upward, downward and stay movements between neighbouring integer states is commonly called a birth–death DTMC. Here “birth” means an increase in the infectious count, not a biological birth.

Derive the SIS transition probabilities

At state \(i\), the susceptible count is \(N-i\). The infection and recovery event rates are:

\[ \lambda_i=\beta\frac{(N-i)i}{N}, \qquad \mu_i=\gamma i. \]

For sufficiently short fixed step \(\Delta t\), the DTMC probabilities are:

\[ \begin{aligned} p_i&=\Pr(I_{n+1}=i+1\mid I_n=i) =\beta\frac{(N-i)i}{N}\Delta t,\\ q_i&=\Pr(I_{n+1}=i-1\mid I_n=i) =\gamma i\Delta t,\\ r_i&=\Pr(I_{n+1}=i\mid I_n=i) =1-p_i-q_i. \end{aligned} \]

Symbol	Meaning	State change
\(p_i\)	Infection probability at state \(i\)	\(+1\)
\(q_i\)	Recovery probability at state \(i\)	\(-1\)
\(r_i\)	No-change probability at state \(i\)	\(0\)

Why the probabilities change with the state

This DTMC is time-homogeneous when \(\beta\), \(\gamma\) and \(\Delta t\) remain fixed. However, its transition probabilities still depend on the current infectious count \(i\).

When \(i\) is small, there are many susceptible people but few infectious people.
At intermediate values, both groups are present, so infection probability can be larger.
When \(i\) is near \(N\), few susceptible people remain, so infection probability falls.
Recovery probability increases linearly with \(i\).

Boundary states

Extinction: (i=0)

\[ p_0=0,\quad q_0=0,\quad r_0=1. \]

With no infectious people, neither infection nor recovery can occur. State 0 is absorbing.

Everyone infectious: (i=N)

\[ p_N=0,\quad q_N=\gamma N\Delta t. \]

There is nobody susceptible to infect, but recovery can move the process down to \(N-1\). State \(N\) is not absorbing when \(\gamma>0\).

Check one fixed time step over the whole state space

A simulation may visit many values of \(i\). Therefore it is not enough to check the probabilities only at the initial state. For every possible state:

\[ p_i+q_i\le1, \qquad i=0,1,\ldots,N. \]

The code calculates this sum for every integer state and rejects \(\Delta t\) if any no-change probability is negative.

Satisfying the inequality makes the probabilities valid. It does not by itself prove that the one-event approximation is accurate. A smaller time step makes multiple unrepresented events within one step less likely.

Interactive Python laboratory

The program verifies all states, simulates one SIS trajectory, records event counts and creates two graphs: the realised epidemic path and the transition probabilities across the state space.

Interactive PythonDTMC SIS epidemic

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def sis_probabilities(i, N, beta, gamma, dt):
    susceptible = N - i

p_infection = beta * susceptible * i / N * dt
    p_recovery = gamma * i * dt
    p_no_change = 1 - p_infection - p_recovery

return p_infection, p_recovery, p_no_change

def one_sis_transition(i, N, beta, gamma, dt, rng):
    p_infection, p_recovery, p_no_change = sis_probabilities(
        i, N, beta, gamma, dt
    )

if p_no_change < -1e-12:
        raise ValueError("dt gives invalid transition probabilities.")

u = rng.random()

if u < p_infection:
        return i + 1, "infection"
    elif u < p_infection + p_recovery:
        return i - 1, "recovery"
    else:
        return i, "no change"

# Population, parameters and initial state
N = 200
beta = 0.25
gamma = 0.20
dt = 0.01
initial_I = 5
max_days = 300
seed = 7

if N <= 0:
    raise ValueError("N must be positive.")

if initial_I < 0 or initial_I > N:
    raise ValueError("initial_I must be between 0 and N.")

if beta < 0 or gamma < 0 or dt <= 0:
    raise ValueError("Rates must be non-negative and dt positive.")

# Check dt at every possible infectious-count state
state_grid = np.arange(N + 1)
infection_probabilities = (
    beta * (N - state_grid) * state_grid / N * dt
)
recovery_probabilities = gamma * state_grid * dt
total_event_probabilities = (
    infection_probabilities + recovery_probabilities
)

maximum_total_probability = total_event_probabilities.max()
worst_state = int(total_event_probabilities.argmax())

if maximum_total_probability > 1 + 1e-12:
    raise ValueError(
        "dt is invalid somewhere in the SIS state space."
    )

print("Whole-state-space time-step check")
print(f"Maximum infection + recovery probability: "
      f"{maximum_total_probability:.4f}")
print(f"This maximum occurs at I = {worst_state}")
print()

# Store and simulate one trajectory
rng = np.random.default_rng(seed)
I = initial_I
times = [0.0]
I_values = [I]
events = ["initial state"]

max_steps = int(max_days / dt)

for step in range(max_steps):
    I, event = one_sis_transition(
        I, N, beta, gamma, dt, rng
    )

times.append((step + 1) * dt)
    I_values.append(I)
    events.append(event)

if I == 0:
        break

trajectory = pd.DataFrame({
    "time": times,
    "S": N - np.array(I_values),
    "I": I_values,
    "event": events
})

# Validate the stored trajectory
if not ((trajectory["S"] + trajectory["I"]) == N).all():
    raise ValueError("Population conservation failed.")

if (trajectory[["S", "I"]] < 0).any().any():
    raise ValueError("A compartment became negative.")

final_time = trajectory["time"].iloc[-1]
final_I = trajectory["I"].iloc[-1]
peak_I = trajectory["I"].max()
peak_time = trajectory.loc[
    trajectory["I"].idxmax(), "time"
]

if final_I == 0:
    outcome = "infection became extinct"
else:
    outcome = "infection remained at the time horizon"

event_counts = (
    trajectory["event"]
    .value_counts()
    .reindex(["infection", "recovery", "no change"], fill_value=0)
)

print("SIS trajectory summary")
print(f"Seed: {seed}")
print(f"Basic ratio beta/gamma: {beta / gamma:.3f}")
print(f"Final time: {final_time:.2f} days")
print(f"Outcome: {outcome}")
print(f"Peak infectious count: {peak_I}")
print(f"First time at peak: {peak_time:.2f} days")
print()
print("Event counts")
print(event_counts.to_string())
print()
print("First 8 states")
print(trajectory.head(8).to_string(index=False))
print()
print("Last 8 states")
print(trajectory.tail(8).to_string(index=False))

# Figure 1: one realised SIS path
plt.figure(figsize=(9, 4.8))
plt.step(trajectory["time"], trajectory["I"],
         where="post", label="Infectious", color="#dc2626")
plt.step(trajectory["time"], trajectory["S"],
         where="post", label="Susceptible", color="#2563eb",
         alpha=0.75)
plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("One stochastic SIS trajectory")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: probabilities at every state
no_change_probabilities = 1 - total_event_probabilities

plt.figure(figsize=(9, 4.8))
plt.plot(state_grid, infection_probabilities,
         label="Infection: i to i + 1", color="#dc2626")
plt.plot(state_grid, recovery_probabilities,
         label="Recovery: i to i - 1", color="#2563eb")
plt.plot(state_grid, no_change_probabilities,
         label="No change: i to i", color="#16a34a")
plt.xlabel("Current infectious count i")
plt.ylabel("One-step probability")
plt.title("SIS transition probabilities depend on the state")
plt.ylim(0, 1.02)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Read the two graphs differently

Graph	Question answered	What it represents
Trajectory graph	What happened in this simulation?	One randomly realised path through time.
Probability-by-state graph	How does the model behave from each possible state?	The model's transition rule, not one realised path.

Do not interpret the probability curves as compartment trajectories. Their horizontal axis is the current state \(i\), not time.

What does \(\beta/\gamma\) tell us?

For this simple SIS model, the ratio:

\[ \mathcal R_0=\frac{\beta}{\gamma} \]

compares transmission with recovery near the disease-free state. If it exceeds 1, infection has an initial tendency to grow in the corresponding deterministic model.

\(\mathcal R_0>1\) does not guarantee growth in every stochastic realisation. With only a few infectious people, early recoveries can still send the DTMC to the absorbing state \(I=0\).

Why event counts are not counts of different people

In an SIS trajectory, one person can be infected, recover and later become infected again. Therefore the total number of infection events may exceed the population size.

This differs from a closed SIR model, where a person enters the removed compartment after recovery and cannot be infected again. SIS has no conventional “final epidemic size” measured by a permanently removed compartment.

Understand the model-specific code

Code	Meaning
`susceptible = N - i`	Reconstructs \(S\) from the one-dimensional infectious state.
`return i + 1`	An infection moves the birth–death chain one state upward.
`return i - 1`	A recovery returns one person to susceptibility and moves the chain downward.
`np.arange(N + 1)`	Creates every infectious-count state from 0 through \(N\), including both boundaries.
`.argmax()`	Returns the position where the total event probability is largest.
`N - np.array(I_values)`	Calculates the complete susceptible trajectory from the stored infectious counts.

Try these experiments

Change seed = 7. Compare early extinction with longer persistence.
Set initial_I = 1. Explain why early chance events become especially influential.
Use beta = 0.15 and gamma = 0.20, so \(\beta/\gamma<1\). Compare the trajectory.
Increase dt gradually. Find when the whole-state-space validity check rejects it.
Set initial_I = N. Confirm from the probability graph that infection probability is zero at the right boundary.

Biological interpretation

The SIS model permits repeated circulation because recovery replenishes the susceptible group. Nevertheless, a finite closed stochastic population eventually can reach \(I=0\), after which the infection cannot return without importation.

A trajectory that remains infectious at 300 days demonstrates persistence during the observation window. It does not prove permanent persistence. Similarly, an early extinction path does not prove that every introduction will fail.

What this lesson has added

SIS recovery returns people to susceptibility.
Population conservation reduces the state to the single count \(I_n\).
The process moves between neighbouring states \(i-1,i,i+1\).
State 0 is absorbing, while state \(N\) allows recovery.
A fixed time step must be valid across the whole state space.
One person may contribute multiple infection events through reinfection.

Lesson 8 turns to the SIR model, where recovery is irreversible and the removed compartment makes final epidemic size meaningful.