Total event rate and exponential waiting time

Several biological events may be possible from the current state. Their rates combine to determine how long the CTMC waits before the next event of any type.

Separate “when?” from “which?”

When does the next event occur?

The total event rate determines the random waiting time.

This lesson

Which event occurs?

Relative event rates determine whether the jump is infection, recovery or another event.

Lesson 4

Keeping these questions separate makes the Gillespie algorithm easier to understand and program.

Add competing event rates

At SIS state \(I=i\), infection and recovery are both possible:

Infection rate\(b(i)\)

Recovery rate\(d(i)\)

Total rate\(a(i)=b(i)+d(i)\)

The total event rate is the sum of the rates of all events currently possible. It measures the instantaneous rate at which the process leaves its present state.

Worked current state

For \(N=1000\), \(i=10\), \(\beta=0.3\) and \(\gamma=0.1\):

\[ b(10)=2.97,\qquad d(10)=1.00, \] \[ a(10)=2.97+1.00=3.97\text{ events/day}. \]

A larger total rate produces a shorter typical waiting time. A smaller total rate produces a longer typical waiting time.

The exponential waiting-time model

If the process is currently in a state with total rate \(a>0\), the waiting time \(\tau\) to the next event is exponentially distributed:

\[ \tau\sim\operatorname{Exponential}(a). \]

The rate \(a\) remains fixed only while the state remains unchanged. When an event occurs, the state changes, the rates are recalculated, and a new waiting time is generated.

Three descriptions of the same distribution

Survival function\[\Pr(\tau>t)=e^{-at}.\]

Probability that no event has occurred by time \(t\).

Cumulative distribution\[\Pr(\tau\le t)=1-e^{-at}.\]

Probability that the next event occurs by time \(t\).

Density\[f(t)=ae^{-at},\quad t\ge0.\]

Describes how waiting-time probability is distributed over continuous time.

The density can exceed 1 for some parameter values; it is not a point probability. Probabilities are obtained from areas under the density curve.

Mean and median waiting times

The mean waiting time is:

\[ \mathbb E[\tau]=\frac{1}{a}. \]

At \(a=3.97\) events/day:

\[ \mathbb E[\tau] =\frac{1}{3.97} \approx0.2519\text{ day} \approx6.05\text{ hours}. \]

The median solves \(\Pr(\tau\le t)=0.5\):

\[ \operatorname{median}(\tau) =\frac{\log 2}{a}. \]

The mean exceeds the median because the exponential distribution has a long right tail: most waits are relatively short, but occasional long waits raise the average.

Connection with the previous lesson

The exact no-event probability over length \(h\) is:

\[ \Pr(\tau>h)=e^{-ah}. \]

For very small \(h\), the exponential expansion gives:

\[ e^{-ah}=1-ah+O(h^2). \]

This recovers the first-order no-event approximation \(1-ah\) from Lesson 2. The exponential distribution provides the exact waiting-time description without advancing through artificial small fixed steps.

Generate a waiting time from a uniform number

Let \(U\sim\operatorname{Uniform}(0,1)\). Inverse-transform sampling gives:

\[ \tau=-\frac{\log U}{a}. \]

Because \(0<U\le1\), \(\log U\le0\) and therefore \(\tau\ge0\).

Generate one uniform value.
Take its natural logarithm.
Change the sign.
Divide by the current total rate.

Why use `1 - rng.random()`?

NumPy's rng.random() generates values in \([0,1)\). The logarithm of zero is not finite. Therefore the code uses:

u = 1.0 - rng.random()

This gives \(U\in(0,1]\) and guarantees that np.log(u) is defined. The transformation has the same uniform distribution apart from irrelevant endpoint conventions.

The memoryless property

\[ \Pr(\tau>s+t\mid\tau>s) = \Pr(\tau>t). \]

Conditional on no event having occurred during the first \(s\) time units, the distribution of the additional waiting time is the same as it was initially.

Memorylessness does not mean the process has no state. It means that while the state and its rates remain unchanged, elapsed waiting time alone does not alter the future waiting-time distribution.

Competing-clock intuition

Infection clock

Imagine an exponential clock with rate \(b(i)\).

Recovery clock

Imagine another exponential clock with rate \(d(i)\).

The first clock to ring gives the next-event time. The minimum of these independent exponential clocks is exponential with rate \(b(i)+d(i)\).

This intuition explains why rates add. Lesson 4 uses the same competing rates to determine which clock wins.

What if the total rate is zero?

\[ a=0\quad\Longrightarrow\quad\tau=\infty. \]

No event is possible from the current state. The process remains there indefinitely, so the code must not divide by zero or try to generate a finite waiting time.

For the closed SIS model, this occurs at \(I=0\), the absorbing extinction state.

Interactive Python laboratory

The code generates 10,000 waiting times at the current state and compares their empirical mean, histogram and survival curve with the exponential theory. It does not choose event types.

Interactive PythonExponential CTMC waiting times

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

# Current SIS state and rates
N = 1000
i = 10
beta = 0.30
gamma = 0.10

infection_rate = beta * (N - i) * i / N
recovery_rate = gamma * i
total_rate = infection_rate + recovery_rate

if total_rate < 0:
    raise ValueError("Total rate cannot be negative.")

if total_rate == 0:
    raise ValueError(
        "No event is possible: the waiting time is infinite."
    )

theoretical_mean = 1 / total_rate
theoretical_median = np.log(2) / total_rate

# One waiting time by inverse-transform sampling
seed = 7
rng = np.random.default_rng(seed)

u = 1.0 - rng.random()
one_waiting_time = -np.log(u) / total_rate

print("Current rates")
print(f"Infection rate: {infection_rate:.3f} events/day")
print(f"Recovery rate: {recovery_rate:.3f} events/day")
print(f"Total rate: {total_rate:.3f} events/day")
print()
print(f"Uniform value U: {u:.6f}")
print(f"Generated waiting time: {one_waiting_time:.6f} days")
print(f"Generated waiting time: "
      f"{one_waiting_time * 24:.3f} hours")
print()

# Many waiting times from the same current state
number_of_draws = 10000
u_values = 1.0 - rng.random(number_of_draws)
waiting_times = -np.log(u_values) / total_rate

summary_table = pd.DataFrame({
    "quantity": [
        "Theoretical mean (days)",
        "Empirical mean (days)",
        "Theoretical median (days)",
        "Empirical median (days)",
        "Minimum simulated wait (days)",
        "Maximum simulated wait (days)"
    ],
    "value": [
        theoretical_mean,
        waiting_times.mean(),
        theoretical_median,
        np.median(waiting_times),
        waiting_times.min(),
        waiting_times.max()
    ]
})

print("Waiting-time summary")
print(summary_table.round(6).to_string(index=False))
print()
print("First 10 waiting times (days)")
print(np.round(waiting_times[:10], 6))

# Figure 1: histogram and theoretical density
t_grid = np.linspace(
    0,
    np.percentile(waiting_times, 99.5),
    300
)
density = total_rate * np.exp(-total_rate * t_grid)

plt.figure(figsize=(9, 4.9))
plt.hist(
    waiting_times,
    bins=45,
    density=True,
    color="#9fc8b1",
    edgecolor="white",
    label="Simulated waiting-time density"
)
plt.plot(
    t_grid,
    density,
    color="#174c34",
    linewidth=2.5,
    label="Theoretical exponential density"
)
plt.xlabel("Waiting time to next event (days)")
plt.ylabel("Probability density")
plt.title("Exponential waiting times at the current state")
plt.legend()
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

# Figure 2: empirical and theoretical survival
survival_times = np.linspace(0, 1.5, 180)
empirical_survival = np.array([
    (waiting_times > t).mean()
    for t in survival_times
])
theoretical_survival = np.exp(
    -total_rate * survival_times
)

plt.figure(figsize=(9, 4.9))
plt.step(
    survival_times,
    empirical_survival,
    where="post",
    color="#2563eb",
    label="Empirical Pr(wait > t)"
)
plt.plot(
    survival_times,
    theoretical_survival,
    color="#dc2626",
    linestyle="--",
    linewidth=2.2,
    label="Theoretical exp(-a t)"
)
plt.xlabel("Time t (days)")
plt.ylabel("Probability no event has occurred")
plt.title("Exponential survival function")
plt.ylim(0, 1.03)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the new code

Code	Meaning
`np.log(u)`	Calculates the natural logarithm required by inverse-transform sampling.
`-np.log(u) / total_rate`	Transforms one uniform value into an exponential waiting time.
`rng.random(10000)`	Generates 10,000 uniform values in one NumPy array.
`density=True`	Scales histogram area to 1 so it can be compared with a probability-density curve.
`(waiting_times > t).mean()`	Calculates the proportion of simulated waits exceeding \(t\), an empirical survival probability.
`np.percentile(..., 99.5)`	Limits the plotted range so a few extreme waits do not compress the main histogram.

Biological interpretation

At the current state, the CTMC waits a random amount of continuous time before either infection or recovery occurs. The mean wait is about 6.05 hours, but individual waits may be much shorter or considerably longer.

The exponential model determines only when the next event occurs. It does not yet say whether that event is infection or recovery. That conditional selection is the unique purpose of Lesson 4.

What this lesson has—and has not—done

You can now:

add competing transition rates to obtain the total leaving rate;
interpret exponential survival, cumulative probability and density;
calculate mean and median waiting times;
generate waiting times by inverse transformation;
explain memorylessness and the zero-total-rate case; and
connect the exact survival function to the short-interval approximation.

Not yet: no event type has been selected and no state has been updated. Lesson 4 converts competing rates into conditional event probabilities.