Programming one CTMC event

We now combine the total rate, exponential waiting time, conditional event selection and biological state change to produce exactly one complete continuous-time transition.

The unique purpose of this lesson

1. Calculate ratesUse the current state only.

2. Generate waiting timeDetermine when the jump occurs.

3. Select eventDetermine which jump occurs.

4. Update stateApply one integer compartment change.

This page completes this mechanism once. It does not repeat the mechanism into a trajectory; that loop is the unique purpose of Lesson 6.

Scenario

At current time \(t=12\) days, use:

\[ (S,E,I,R)=(990,5,5,0), \qquad N=1000, \]

\[ \beta=0.3,\qquad \sigma=0.2,\qquad \gamma=0.1\text{ per day}. \]

The current event rates are:

\[ a_1=1.485,\qquad a_2=1.000,\qquad a_3=0.500, \]

for infection, progression and recovery, respectively. The total rate is \(a_0=2.985\) events/day.

Use two independent random numbers

Waiting-time draw

\[\tau=-\frac{\log U_1}{a_0}.\]

\(U_1\) determines when the event occurs.

Event-selection draw

\[U_2\sim\operatorname{Uniform}(0,1).\]

\(U_2\) determines which event occurs.

Use separate independent draws. Reusing the same uniform number would create an artificial relationship between waiting time and event type that is not part of the CTMC model.

Advance continuous time

The current state remains unchanged throughout the waiting period. At its end, the selected event occurs instantaneously:

\[ t_{\mathrm{next}}=t+\tau. \]

Unlike a DTMC, this increment is not a fixed \(\Delta t\). Every CTMC event generally produces a different positive waiting time.

Select the event from rate ratios

The cumulative conditional boundaries are:

\[ c_1=\frac{a_1}{a_0}, \qquad c_2=\frac{a_1+a_2}{a_0}. \]

If \(U_2<c_1\), select infection.
If \(c_1\le U_2<c_2\), select progression.
Otherwise, select recovery.

There is no no-change choice because waiting has already been represented by \(\tau\).

Apply exactly one state update

Infection

\[(S,E,I,R)\to(S-1,E+1,I,R).\]

Progression

\[(S,E,I,R)\to(S,E-1,I+1,R).\]

Recovery

\[(S,E,I,R)\to(S,E,I-1,R+1).\]

The event rate controls selection, but the selected event controls the integer-valued compartment changes.

Do not update while calculating rates

All rates must be calculated from the same current state. Only after both random decisions have been made should the program form the next state.

Using next_S, next_E, next_I and next_R keeps the current and next states separate and avoids mixing values from different times.

What if no event is possible?

If \(a_0=0\), the current state is absorbing:

there is no finite waiting time;
no event type can be normalised or selected; and
the state remains unchanged indefinitely.

A robust function should detect this case before calculating \(-\log(U_1)/a_0\).

Validation after the transition

Non-negative countsNo compartment may fall below zero.

Population conservationThe next counts must still sum to \(N\).

Increasing timeFor \(a_0>0\), the generated event time must not precede the current time.

Interactive Python laboratory

Click Run code to perform one complete event. The result shows rates, two random numbers, waiting time, selected event, event time and before/after state.

Interactive PythonOne complete SEIR CTMC event

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

def one_seir_ctmc_event(
    time, S, E, I, R, beta, sigma, gamma, rng
):
    N = S + E + I + R

# All rates use the same current state
    infection_rate = beta * S * I / N
    progression_rate = sigma * E
    recovery_rate = gamma * I

rates = np.array([
        infection_rate,
        progression_rate,
        recovery_rate
    ])

if (rates < 0).any():
        raise ValueError("Event rates cannot be negative.")

total_rate = rates.sum()

if total_rate == 0:
        return {
            "absorbing": True,
            "time": time,
            "waiting_time": float("inf"),
            "event": "no possible event",
            "rates": rates,
            "probabilities": np.zeros(3),
            "u_time": None,
            "u_event": None,
            "next_state": (S, E, I, R)
        }

probabilities = rates / total_rate

# Independent random draw 1: waiting time
    u_time = 1.0 - rng.random()
    waiting_time = -np.log(u_time) / total_rate
    next_time = time + waiting_time

# Independent random draw 2: event type
    u_event = rng.random()
    infection_end = probabilities[0]
    progression_end = probabilities[0] + probabilities[1]

if u_event < infection_end:
        event = "infection"
        next_S, next_E, next_I, next_R = (
            S - 1, E + 1, I, R
        )
    elif u_event < progression_end:
        event = "progression"
        next_S, next_E, next_I, next_R = (
            S, E - 1, I + 1, R
        )
    else:
        event = "recovery"
        next_S, next_E, next_I, next_R = (
            S, E, I - 1, R + 1
        )

# Validate the next state
    if min(next_S, next_E, next_I, next_R) < 0:
        raise ValueError("A compartment became negative.")

if next_S + next_E + next_I + next_R != N:
        raise ValueError("Population conservation failed.")

if next_time < time:
        raise ValueError("Event time moved backwards.")

return {
        "absorbing": False,
        "time": next_time,
        "waiting_time": waiting_time,
        "event": event,
        "rates": rates,
        "probabilities": probabilities,
        "u_time": u_time,
        "u_event": u_event,
        "next_state": (next_S, next_E, next_I, next_R)
    }

# Current time, state and parameters
time = 12.0
S, E, I, R = 990, 5, 5, 0
beta = 0.30
sigma = 0.20
gamma = 0.10

seed = 7
rng = np.random.default_rng(seed)

result = one_seir_ctmc_event(
    time, S, E, I, R, beta, sigma, gamma, rng
)

event_names = ["infection", "progression", "recovery"]

rate_table = pd.DataFrame({
    "event": event_names,
    "rate (events/day)": result["rates"],
    "conditional probability": result["probabilities"]
})

print("Current-event calculations")
print(rate_table.round(6).to_string(index=False))
print()
print(f"Waiting-time uniform U1: {result['u_time']:.6f}")
print(f"Event-choice uniform U2: {result['u_event']:.6f}")
print(f"Waiting time: {result['waiting_time']:.6f} days")
print(f"Current time: {time:.6f} days")
print(f"Next event time: {result['time']:.6f} days")
print(f"Selected event: {result['event']}")
print()

next_S, next_E, next_I, next_R = result["next_state"]

state_table = pd.DataFrame({
    "compartment": ["Susceptible", "Exposed", "Infectious", "Removed", "Total"],
    "before": [S, E, I, R, S + E + I + R],
    "after": [
        next_S, next_E, next_I, next_R,
        next_S + next_E + next_I + next_R
    ],
    "change": [
        next_S - S,
        next_E - E,
        next_I - I,
        next_R - R,
        (next_S + next_E + next_I + next_R)
        - (S + E + I + R)
    ]
})

print("State update")
print(state_table.to_string(index=False))

labels = ["Susceptible", "Exposed", "Infectious", "Removed"]
changes = [
    next_S - S,
    next_E - E,
    next_I - I,
    next_R - R
]
colours = [
    "#dc2626" if value < 0
    else "#39815e" if value > 0
    else "#9ca3af"
    for value in changes
]

plt.figure(figsize=(8.5, 4.8))
bars = plt.bar(labels, changes, color=colours)
plt.axhline(0, color="black", linewidth=1)
plt.bar_label(bars, fmt="%+g")
plt.ylabel("Change caused by one event")
plt.ylim(-1.4, 1.4)
plt.title(
    f"One CTMC state update: {result['event']} "
    f"at t={result['time']:.3f} days"
)
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the combined code

Code	Purpose
`rates.sum()`	Calculates the current total rate.
`u_time = 1.0 - rng.random()`	Generates a positive uniform value for the logarithm.
`u_event = rng.random()`	Uses a separate random draw for event selection.
`next_time = time + waiting_time`	Advances continuous time by the generated waiting duration.
`next_S, next_E, next_I, next_R = ...`	Assigns the complete next state in the selected branch.
`float("inf")`	Represents indefinite waiting when no transition rate is positive.

Read the result correctly

The generated waiting time tells how long the process remains in the current state. At the reported next-event time, the selected event changes the state by one allowed transition.

Running the code with another seed can change both the waiting time and event type. Neither result alone is an average prediction; each is one valid random CTMC event.

Try these experiments

Change the seed and observe that the two random draws and resulting event change.
Set \(E=0\) while keeping \(I>0\). Confirm that progression probability becomes zero.
Use \((S,E,I,R)=(990,10,0,0)\). Progression remains possible even when infection and recovery rates are zero.
Use \((S,E,I,R)=(990,0,0,10)\). Confirm that the function reports an absorbing state rather than dividing by zero.
Temporarily reuse u_time as u_event, then restore the independent draw and explain why coupling them is inappropriate.

What this lesson has—and has not—done

You can now:

calculate all rates from one current state;
generate one exponential waiting time;
select one event using an independent uniform draw;
advance continuous time and update the state;
handle absorbing states; and
validate population counts and time direction.

Not yet: the mechanism has been executed only once. Lesson 6 repeats it, recalculates rates after every jump and stores a complete Gillespie trajectory.