6.5 Birth-death processes

Definition of birth-death processes

A birth-death process is a continuous-time Markov chain where the state can only change by $+1$ (a "birth") or $-1$ (a "death") at any given transition. The rates at which these transitions happen are governed by state-dependent birth rates $\lambda_i$ and death rates $\mu_i$. These processes show up constantly in population dynamics, queueing theory, and epidemiology because many real systems evolve through incremental additions and removals.

Discrete-time vs continuous-time

Birth-death models exist in both discrete and continuous time. In discrete time, transitions happen at fixed intervals (every day, every generation, etc.). In continuous time, transitions can occur at any moment, and the waiting time in each state follows an exponential distribution. This unit focuses on the continuous-time version, where the exponential holding times give the process its memoryless property.

State space and transitions

The state space is typically $\{0, 1, 2, \ldots\}$, representing a count of something: individuals in a population, customers in a queue, infected people in an epidemic. The defining structural constraint is that transitions only occur between adjacent states:

• From state $i$ to state $i+1$ at rate $\lambda_i$ (birth)
• From state $i$ to state $i-1$ at rate $\mu_i$ (death)

No jumps of size 2 or larger are allowed. This gives the generator matrix $Q$ a tridiagonal structure, which is what makes many birth-death calculations tractable.

Birth and death rates

• Birth rates $\lambda_i$ govern how quickly the system moves upward from state $i$. Think of new customers arriving or new infections occurring.
• Death rates $\mu_i$ govern how quickly the system moves downward. Think of service completions or recoveries.

These rates can be constant (the same for every state) or state-dependent (varying with $i$). State-dependent rates are common: for example, in a population model, the total birth rate might scale linearly with population size, so $\lambda_i = i\lambda$.

Transition probabilities

Transition probabilities $P_{ij}(t) = P(X(t+s) = j \mid X(s) = i)$ describe the likelihood of being in state $j$ after time $t$, given the process started in state $i$. Computing these is central to understanding both the short-term dynamics and long-run behavior of the chain.

Chapman-Kolmogorov equations

The Chapman-Kolmogorov equations express the idea that you can break any time interval into two pieces and sum over all intermediate states:

$P_{ij}(t+s) = \sum_k P_{ik}(t) \, P_{kj}(s)$

This is just the law of total probability applied to the Markov property. In matrix form, it says $P(t+s) = P(t)P(s)$, which means the transition probability matrices form a semigroup. These equations are the foundation for deriving the Kolmogorov forward and backward differential equations, which give you the actual ODEs for computing $P_{ij}(t)$.

Transition rate matrix (generator matrix)

For continuous-time chains, the key object is the infinitesimal generator matrix $Q$, not a one-step transition probability matrix. For a birth-death process, $Q$ has the tridiagonal form:

• Off-diagonal entries: $q_{i,i+1} = \lambda_i$ and $q_{i,i-1} = \mu_i$
• Diagonal entries: $q_{ii} = -(\lambda_i + \mu_i)$
• All other entries are zero

The embedded discrete-time chain has transition probabilities $p_{i,i+1} = \lambda_i / (\lambda_i + \mu_i)$ and $p_{i,i-1} = \mu_i / (\lambda_i + \mu_i)$, which describe where the chain goes next (ignoring when).

Stationary distribution

The stationary distribution $\pi$ is the long-run probability of being in each state. For a continuous-time chain, you find it by solving:

$\pi Q = 0, \quad \sum_i \pi_i = 1$

For birth-death processes, the tridiagonal structure of $Q$ means you can use detailed balance (also called the local balance equations):

$\lambda_i \, \pi_i = \mu_{i+1} \, \pi_{i+1}, \quad i = 0, 1, 2, \ldots$

This gives the explicit solution:

$\pi_n = \pi_0 \prod_{i=0}^{n-1} \frac{\lambda_i}{\mu_{i+1}}$

where $\pi_0$ is determined by the normalization condition. A stationary distribution exists if and only if the chain is irreducible and the sum $\sum_{n=0}^{\infty} \prod_{i=0}^{n-1} \frac{\lambda_i}{\mu_{i+1}}$ converges (positive recurrence).

Poisson processes

The Poisson process is the simplest birth-death process: set all birth rates to a constant $\lambda$ and all death rates to zero. It models pure arrivals with no departures.

Relationship to birth-death processes

A Poisson process with rate $\lambda$ is a birth-death chain on $\{0, 1, 2, \ldots\}$ with $\lambda_i = \lambda$ for all $i$ and $\mu_i = 0$ for all $i$. The state only increases, so it counts the cumulative number of events. The number of events in an interval of length $t$ follows a Poisson distribution with mean $\lambda t$.

Poisson process properties

Three properties characterize a Poisson process:

• Memoryless property: The time until the next event doesn't depend on how long you've already waited. This comes directly from the exponential distribution of inter-arrival times.
• Stationary increments: The distribution of the number of events in any interval depends only on the interval's length, not its location in time.
• Independent increments: Event counts in non-overlapping intervals are independent random variables.

Exponential inter-arrival times

The time between consecutive events follows an exponential distribution with rate $\lambda$:

$f(t) = \lambda e^{-\lambda t}, \quad t \geq 0$

• Expected inter-arrival time: $E(T) = \frac{1}{\lambda}$
• Variance: $\text{Var}(T) = \frac{1}{\lambda^2}$

The memoryless property of the exponential distribution ($P(T > t + s \mid T > t) = P(T > s)$) is what makes the Poisson process Markovian.

Applications of birth-death processes

Population growth models

Birth-death processes model population dynamics where each individual independently gives birth at rate $\lambda$ and dies at rate $\mu$. With state-dependent rates $\lambda_i = i\lambda$ and $\mu_i = i\mu$, you get a linear birth-death process. If $\lambda > \mu$, the population tends to grow exponentially; if $\lambda < \mu$, it tends toward extinction. More sophisticated models add carrying capacity through nonlinear rate functions, similar in spirit to the logistic growth model.

Queueing systems

Queueing theory is one of the richest application areas. The classic M/M/1 queue has Poisson arrivals at rate $\lambda$ and exponential service at rate $\mu$, giving a birth-death process with constant rates. The stationary distribution exists when the traffic intensity $\rho = \lambda / \mu < 1$, and the expected number of customers in the system is $\frac{\rho}{1 - \rho}$. The M/M/c queue extends this to $c$ servers, where the death rate becomes $\min(i, c) \cdot \mu$ in state $i$.

Epidemiology and disease spread

In the stochastic SIS model, infected individuals recover at rate $\mu$ and infect susceptible individuals at a rate proportional to the number of infected. This creates a birth-death process on the number of infected individuals. The stochastic SIR model is similar but with permanent recovery. Birth-death formulations let you compute quantities like the probability of a major outbreak and the expected time to disease extinction, which deterministic ODE models cannot provide.

Limiting behavior

The limiting behavior describes what happens to the process as $t \to \infty$. Does it settle into an equilibrium? Does it drift to infinity? Does it get absorbed?

Transient vs recurrent states

• A state is recurrent if the process returns to it with probability 1. It's positive recurrent if the expected return time is finite.
• A state is transient if there's a positive probability of never returning.

For an irreducible birth-death process on $\{0, 1, 2, \ldots\}$, all states share the same classification. The chain is positive recurrent (and has a stationary distribution) if and only if $\sum_{n=1}^{\infty} \prod_{j=1}^{n} \frac{\lambda_{j-1}}{\mu_j} < \infty$.

Absorbing states and absorption probabilities

If state 0 has $\lambda_0 > 0$ but $\mu_0 = 0$, it's not absorbing. But if $\lambda_0 = 0$, then state 0 is absorbing: once the process reaches 0, it stays there forever. Absorption probabilities from state $i$ can be computed by solving a system of linear equations derived from first-step analysis on the embedded chain, or equivalently using the fundamental matrix of the transient portion of the chain.

Long-run distribution

For an irreducible, positive recurrent birth-death process, the long-run distribution equals the stationary distribution found by solving $\pi Q = 0$ with $\sum_i \pi_i = 1$. The detailed balance solution gives:

$\pi_n = \pi_0 \prod_{i=0}^{n-1} \frac{\lambda_i}{\mu_{i+1}}$

If the chain is not positive recurrent (the sum diverges), no stationary distribution exists, and the process is either null recurrent or transient.

First passage times

The first passage time $T_{ij}$ is the time it takes the process to reach state $j$ for the first time, starting from state $i$. These random variables capture the timescales of the system: how long until a queue empties, how long until a population reaches a threshold, etc.

Definition and properties

• $T_{ij}$ is always non-negative.
• By the Markov property, the distribution of $T_{ij}$ depends only on $i$ and $j$, not on the prior history.
• For birth-death processes, reaching state $j$ from state $i$ requires passing through every intermediate state (since only $\pm 1$ transitions are allowed). This means $T_{ij}$ can be decomposed as a sum of independent passage times through consecutive states.

Expected first passage times

For a birth-death process, the expected first passage time from state $i$ to state $j$ satisfies:

$E(T_{ij}) = \frac{1}{\lambda_i + \mu_i} + \frac{\lambda_i}{\lambda_i + \mu_i} E(T_{i+1,j}) + \frac{\mu_i}{\lambda_i + \mu_i} E(T_{i-1,j})$

for $i \neq j$, with boundary condition $E(T_{jj}) = 0$ (you're already there). The $\frac{1}{\lambda_i + \mu_i}$ term is the expected holding time in state $i$, and the remaining terms account for which direction the chain jumps.

Because transitions are only to adjacent states, the expected time to go from state $i$ to state $i+1$ (one step up) can be computed directly, and longer passages are sums of these one-step times.

Variance of first passage times

The variance $\text{Var}(T_{ij})$ measures the spread around the expected passage time. It can be computed by setting up a second system of equations for $E(T_{ij}^2)$ using the same first-step conditioning approach, then applying $\text{Var}(T_{ij}) = E(T_{ij}^2) - [E(T_{ij})]^2$. High variance indicates the passage time is unpredictable, which matters in applications like service guarantees in queueing.

Extinction probabilities

Extinction means the process reaches state 0 and stays there (when state 0 is absorbing). This is especially relevant in population models: will a small population die out, or will it survive?

Conditions for extinction

For a birth-death process on $\{0, 1, 2, \ldots\}$ with state 0 absorbing, extinction from state $i$ is certain (probability 1) if and only if:

$\sum_{n=1}^{\infty} \prod_{j=1}^{n} \frac{\mu_j}{\lambda_j} = \infty$

Intuitively, if death rates dominate birth rates sufficiently, the process will always drift down to 0. If birth rates are strong enough that this sum converges, there's a positive probability of escaping to infinity (survival).

Calculation of extinction probabilities

The extinction probability $\rho_i$ from state $i$ satisfies the recursion from first-step analysis on the embedded chain:

$\rho_i = \frac{\mu_i}{\lambda_i + \mu_i} \, \rho_{i-1} + \frac{\lambda_i}{\lambda_i + \mu_i} \, \rho_{i+1}$

with boundary condition $\rho_0 = 1$. This is a second-order linear difference equation. The solution is:

$\rho_i = \begin{cases} 1 & \text{if } \sum_{n=1}^{\infty} \prod_{j=1}^{n} \frac{\mu_j}{\lambda_j} = \infty \\ \frac{\sum_{n=i}^{\infty} \prod_{j=1}^{n} \frac{\mu_j}{\lambda_j}}{1 + \sum_{n=1}^{\infty} \prod_{j=1}^{n} \frac{\mu_j}{\lambda_j}} & \text{otherwise} \end{cases}$

For the simple case where $\lambda_i = i\lambda$ and $\mu_i = i\mu$ (linear birth-death), the extinction probability from state $i$ is $(\mu/\lambda)^i$ when $\lambda > \mu$, and 1 when $\lambda \leq \mu$.

Time to extinction

The expected time to extinction $E(T_0 \mid X(0) = i)$ can be computed using the first passage time framework with state 0 as the target. For the linear birth-death process starting from state 1 with $\lambda < \mu$, the expected extinction time is $\frac{1}{\mu - \lambda} \ln\!\left(\frac{\mu}{\mu - \lambda}\right)$. In general, these calculations involve solving the appropriate system of equations using the generator matrix.

Birth-death processes with immigration

Adding immigration means individuals enter the system from an external source at some rate, independent of the current state. This prevents the population from staying at zero and often guarantees the existence of a stationary distribution.

Definition and properties

The transition rates become:

• State $i \to i+1$: rate $\lambda_i + \nu$, where $\nu$ is the immigration rate
• State $i \to i-1$: rate $\mu_i$ (for $i > 0$)

The immigration rate $\nu$ acts as a constant upward push regardless of the current state. Because of this, state 0 is never absorbing (even if $\lambda_0 = 0$, the rate out of state 0 is at least $\nu > 0$).

Stationary distribution with immigration

The stationary distribution is found from $\pi Q = 0$ with $\sum_i \pi_i = 1$, using the modified rates. Detailed balance gives:

$(\lambda_i + \nu) \, \pi_i = \mu_{i+1} \, \pi_{i+1}$

so $\pi_n = \pi_0 \prod_{i=0}^{n-1} \frac{\lambda_i + \nu}{\mu_{i+1}}$. The immigration term makes it more likely that this sum converges (and a stationary distribution exists), since the numerator grows but the death rates in the denominator can still dominate.

For example, the M/M/1 queue is a birth-death process with constant arrival rate $\lambda$ (which you can think of as pure immigration with $\lambda_i = 0$ and $\nu = \lambda$) and service rate $\mu$. Its stationary distribution is geometric: $\pi_n = (1 - \rho)\rho^n$ where $\rho = \lambda/\mu < 1$.

Applications in queueing theory

• M/M/1 queue with immigration: Models a single server where customers arrive from outside (Poisson arrivals at rate $\lambda$) and are served at rate $\mu$. The immigration framework applies when arrivals are purely external.
• M/M/$\infty$ queue: Every arriving customer gets their own server immediately. With arrival rate $\lambda$ and per-server service rate $\mu$, the stationary distribution is Poisson with mean $\lambda/\mu$. This is a birth-death process with $\lambda_i = \lambda$ (immigration) and $\mu_i = i\mu$.

Branching processes

Branching processes model populations where each individual independently produces a random number of offspring. They're closely related to birth-death processes but allow for multiple offspring per reproduction event.

Relationship to birth-death processes

In a branching process, each individual in generation $n$ independently produces a random number of offspring according to some distribution $\{p_k\}_{k=0}^{\infty}$. A birth-death process is the special case where the offspring distribution is concentrated on $\{0, 1\}$: each individual either dies (0 offspring) or is replaced by two (itself plus one new, effectively 1 net birth). The branching process framework is more general because it allows jumps larger than $+1$ in a single generation.

Generating functions

The probability generating function (PGF) of the offspring distribution is:

$f(s) = \sum_{k=0}^{\infty} p_k s^k$

This is the main analytical tool. If $Z_n$ is the population size at generation $n$, its PGF satisfies the recursion:

$F_{n+1}(s) = f(F_n(s)), \quad F_0(s) = s$

From the PGF you can extract moments: $f'(1) = E[\text{offspring}] = m$ (the mean number of offspring per individual), and $f''(1) + f'(1) - [f'(1)]^2 = \text{Var}(\text{offspring})$.

Extinction probabilities in branching processes

The extinction probability $q$ (starting from a single individual) is the smallest non-negative solution to:

$f(s) = s$

The classification depends on the mean offspring number $m = f'(1)$:

• Subcritical ($m < 1$): Extinction is certain ($q = 1$). The population shrinks on average each generation.
• Critical ($m = 1$): Extinction is still certain ($q = 1$), but it takes longer on average.
• Supercritical ($m > 1$): There's a positive probability of survival ($q < 1$). The extinction probability $q$ is the unique root of $f(s) = s$ in $[0, 1)$.

Simulation of birth-death processes

When analytical solutions are hard to obtain (complex state-dependent rates, finite populations with boundaries, etc.), simulation provides a way to estimate stationary distributions, passage times, and extinction probabilities numerically.

Monte Carlo methods

Monte Carlo simulation for birth-death processes means generating many independent sample paths and using them to estimate quantities of interest. For example, to estimate the extinction probability from state 5, you simulate thousands of trajectories starting at state 5 and compute the fraction that reach state 0. The accuracy improves with more samples, scaling as $1/\sqrt{N}$ where $N$ is the number of trajectories.

Gillespie algorithm

The Gillespie algorithm (also called the stochastic simulation algorithm) is the standard method for simulating continuous-time Markov chains exactly. For a birth-death process, it works as follows:

1. Initialize: Set the current state $X = i$ and time $t = 0$.

2. Compute the total rate: $r = \lambda_i + \mu_i$ (the total rate of leaving state $i$).

3. Generate the holding time: Draw $\tau$ from an exponential distribution with rate $r$. Update $t \leftarrow t + \tau$.

4. Determine the event type: With probability $\lambda_i / r$, a birth occurs ($X \leftarrow X + 1$). With probability $\mu_i / r$, a death occurs ($X \leftarrow X - 1$).

5. Repeat from step 2 until a stopping condition is met (e.g., $t$ exceeds a time limit, or the process reaches an absorbing state).

This algorithm is exact: it produces sample paths with the correct statistical properties, unlike approximate time-discretization schemes. Each step requires only two random number draws (one for the waiting time, one for the event type), making it computationally efficient.