Worked Birth-Death Simulations with Gillespie's Algorithm

From a Birth–Death Model to Simulated Cell Paths

A deterministic population model describes a smooth change in cell number. A stochastic birth–death model instead represents individual cell events: one cell divides, or one cell dies, and the population changes by one. Even when the event rates are fixed, repeated simulations can produce different event times and different population sizes. This variation is relevant when cell counts are small or when a particular simulated trajectory matters.

Here, “finite population” means the initial count is finite and the process tracks integer cell counts. The basic birth–death model below has no fixed upper bound on the population; it is observed over a finite time horizon. Its rates are proportional to the current number of cells. The model assumes events occur one at a time, with rates that remain constant between events. It does not include delayed effects or changes in rates between events.

The rates \(\lambda\) and \(\mu\) have units of inverse time, such as day\(^{-1}\), and the propensities have units of events per unit time. Birth and death events are conditionally independent given the current state. At \(n=0\), both propensities are zero: without an immigration process, zero cells is an absorbing state. For \(n>0\), the next event is simulated from the current propensities, then the propensities are recalculated for the new state.

Why the Direct Method Generates These Events

The earlier material on Gillespie’s direct method establishes how exponential waiting times and reaction propensities determine the next event. For this birth–death application, if the current state is \(n>0\), the waiting time to a birth is exponential with rate \(a_B(n)\), and the waiting time to a death is exponential with rate \(a_D(n)\). The first of these competing event times determines what happens next.

The probability that neither event has occurred by elapsed time \(t\) is the product of their survival probabilities. Thus the waiting time to either event has an exponential distribution with rate \(a_0(n)\). The density that the next event occurs at time \(t\) and is a birth is \(a_B(n)\exp[-a_0(n)t]\). Integrating this density over all \(t\geq 0\) gives the probability that the next event is a birth.

Consequently, with independent draws \(U_1,U_2\sim\operatorname{Uniform}(0,1)\), the waiting time is \(\tau=-\log(U_1)/a_0(n)\). Select a birth if \(U_2a_0(n)<a_B(n)\); otherwise select a death. The event time is advanced by \(\tau\), the cell count is changed by one, and the calculation repeats. These formulas require finite propensities that stay constant while waiting for the next event.

Representing and Summarizing a Trajectory

Between events, \(N(t)\) is constant. A useful record therefore includes the event time, the event type, the count immediately before the event, and the count immediately after it. The trajectory is usually taken to be right-continuous: at an event time, its value is the count after that event. For example, the preceding simulation has a death at day 0.4581, a birth at day 1.5311, and a death at day 1.7865. Recording only the final count would discard the timing and order of those events.

A time-weighted average population over \([0,T]\) is obtained by multiplying each population level by the time spent at that level, summing, and dividing by \(T\). This differs from the unweighted average of the counts observed at event times, because event times need not be evenly spaced.

A simulation can also be saved at selected observation times, such as \(t=0.5\) and \(t=1.0\), by reading the current state of the event-driven path at each time. If no event occurs between two observation times, the recorded population is unchanged. These saved values are useful for comparing trajectories, but they do not replace the event-time record when the timing of events is important.

Replicate Variability and Its Interpretation

A replicate is a separate simulation using the same model, parameter values, initial condition, and time horizon, but new random draws. Conditional on those shared inputs, differences among replicates are demographic stochasticity: the randomness in the birth and death events. They are not, by themselves, evidence that the rate parameters differ across cells or subjects.

With more replicates, a histogram or empirical quantiles can show the range and shape of simulated counts at a selected time. Individual trajectories can also reveal whether variability is concentrated early or late, or whether paths frequently reach zero. Any such summary is conditional on the specified model, parameter values, initial state, and simulation horizon.

Assumptions, Diagnostics, and Alternatives

The direct simulation is exact for the specified continuous-time Markov birth–death model when its assumptions hold: the state is an integer count, events change the count by one, the stated propensities are correct, and rates do not change during a waiting interval. Exactness means that the event-driven algorithm samples the model’s distribution; it does not mean that the model is a complete description of real cell biology. On a finite horizon, the linear birth–death process starting from a finite count has a finite count and finitely many events with probability one.

Check implementation and output before interpreting a simulation. The following checks are particularly useful:

• Every event time must be later than the previous one, and no event after the chosen horizon should be applied.
• Each event must change the count by exactly \(+1\) or \(-1\). A death must not be selected when \(n=0\), and the count must never become negative.
• Recompute \(a_B=\lambda n\), \(a_D=\mu n\), and \(a_0=a_B+a_D\) after every event. Using rates from the prior state gives incorrect waiting times and event probabilities.
• Check units: a propensity in day\(^{-1}\) produces a waiting time in days. A rate intended per cell must be multiplied by the current cell count.
• Compare replicated runs with different random draws. Identical paths can indicate that the random stream is being reused; unusually large changes can indicate a unit or rate error.

A common mistake is to simulate every event but keep the population fixed when calculating the next waiting time. For example, after a death from \(n=4\) to \(n=3\), the total propensity in the worked model changes from \(2.00\) to \(1.50\) day\(^{-1}\). Reusing \(2.00\) would make the next waiting time too short on average and distort the trajectory. The correction is to update the state first, then calculate all propensities from that state.

The exact direct method is event-driven and preserves integer counts, making it a natural choice when event timing or small populations matter. A deterministic ordinary differential equation instead gives a smooth trajectory and does not represent individual random event histories. Time-step approximations can reduce computational cost in some settings, but require care about step size and possible invalid count changes. The appropriate alternative depends on whether the goal is individual event timing, a smooth population trend, or computational efficiency.