Simulating birth-death processes

The birth-death process $\{N(t)\}$ with birth and death rates $\lambda$ and $\mu$, respectively, and initial population size $N(0)=n_0$ is characterized by the following differential-difference equations for the probability mass functions $p_n(t)=Pr(N(t)=n)$:

$$\begin{array}{rcl}\frac{d{p}_0(t)}{dt}& =& \mu p_1(t),\\ \text{(Eq:BDn)}& \frac{d{p}_n(t)}{dt}& =& \lambda (n-1)p_{n-1}(t)-(\lambda +\mu )n{p}_n(t)+\mu (n+1)p_{n+1}(t),(n\ge 1).\end{array}$$

with the initial condition

$$p_n(0)={\delta }_{n,n_0}.$$

See Combining birth and death processes for the detail of how to analytically solve the above differential-difference equations. Here, we show how to simulate this process numerically as we did for the birth process and death process; See also:

• Simulating birth processes
• Simulating death processes

First, we discretize the time steps with some small interval $\delta t$: $N(0),N(\delta t),N(2\delta t),N(3\delta t),\cdots ,N(k\delta t),\cdots$. 

Next, we exploit the Markov property of the birth-death process: Every birth or death is the first birth or death since the last one. That is, given $N(\tau )=n_{\tau }$ at time $t=\tau$, we may regard this as the "initial" condition for the next step: $N(\tau +\delta t)=n_{\tau }$ or $n_{\tau }+1$ or $n_{\tau }-1$. 

Unlike the simple birth or death process, we have three mutually disjoint events. That is, we assume that during the short period $\delta t$, one of the following occurs: one birth $N(\tau +\delta t)=n_{\tau }+1$, one death $N(\tau +\delta t)=n_{\tau }-1$, or nothing $N(\tau +\delta t)=n_{\tau }$. Thus, we have the following differential-difference equations with $p_n(\tau )={\delta }_{n,n_{\tau }}$

$$\begin{array}{rcl}\frac{d{p}_{n_{\tau }-1}(t)}{dt}& =& \mu n_{\tau }{p}_{n_{\tau }}(t),\\ \frac{d{p}_{n_{\tau }}(t)}{dt}& =& -(\lambda +\mu )n_{\tau }{p}_{n_{\tau }}(t),\\ \frac{d{p}_{n_{\tau }+1}(t)}{dt}& =& \lambda n_{\tau }{p}_{n_{\tau }}(t).\end{array}$$

Solving these, we have

$$\begin{array}{rcl}{p}_{n_{\tau }-1}(\tau +\delta t)& =& \frac{\mu }{\lambda +\mu }(1-e^{-(\lambda +\mu )n_{\tau }\delta t})\\ p_{n_{\tau }}(\tau +\delta t)& =& e^{-(\lambda +\mu )n_{\tau }\delta t},\\ p_{n_{\tau }+1}(\tau +\delta t)& =& \frac{\lambda }{\lambda +\mu }(1-e^{-(\lambda +\mu )n_{\tau }\delta t}).\end{array}$$

Note that $p_{n_{\tau }-1}(\tau +\delta t)+p_{n_{\tau }}(\tau +\delta t)+p_{n_{\tau }+1}(\tau +\delta t)=1.$ 

Based on this, we have the following algorithm for simulating the birth-death process:

1. Set $N(0)=n_0$; set $k=0$.
2. Generate a uniformly distributed random number $r\in [0,1)$.
1. If $r<\frac{\mu }{\lambda +\mu }(1-e^{-(\lambda +\mu )N(k\delta t)\delta t})$, then $N((k+1)\delta t)=N(k\delta t)-1$;
2. Otherwise, if $r<1-e^{-(\lambda +\mu )N(k\delta t)\delta t}$, then $N((k+1)\delta t)=N(k\delta t)+1$;
3. Otherwise, $N((k+1)\delta t)=N(k\delta t)$.
3. Update $k:=k+1$, go to Step 2, and repeat.

A sample path obtained from such simulations is shown below:

A sample path of a birth-death process with birth rate $\lambda =0.1$ and death rate $\mu =0.09$.