Pure birth process and Yule Process

Introduction

You need to be familiar with continuous Markov Chain, Poisson process, especially about the notations from John’s notes.

Motivation

In this vignette, we want to model the number of patients at the initial stage of the pandemic, when unluckily, no one is paying attention to COVID or wearing a mask to prevent the spread. Intuitively, the spread of COVID will be closely related to the number of current infected patients, who have the ability to infect more people (i.e. give birth to more patients). As more and more people get infected, the time interval before the ‘birth’ of next patient will become shorter and shorter, and thus a traditional Poisson process, which assumes a constant birth rate at any time, fails to model this scenario. This leads us to Pure-birth process and Yule process.

Definition and key properties

Pure-birth process is a special case of continuous Markov-Chain process. It only has one type of state transition: ‘birth’, which increases the current state variable by one. The birth rate at each state is only dependent on the current state. It’s a generalization of Poisson process, which assumes a constant birth rate at any state.

\(X(t)\) is known as birth process with rate \(\lambda_n, n \in N\) if:

\[ \forall t, h >0, X(t+h) \geq X(t) \\ P\{X(t+h) - X(t) = 1 | X(t) = n \} = \lambda_{n}h + o(h) \\ P\{X(t+h) - X(t) = 0 | X(t) = n \} = 1- \lambda_n h + o(h) \\ P\{X(t+h) - X(t) > 1 | X(t) = n \} = o(h) \\ \]

Accordingly, the transition matrix \(Q\) of a pure-birth process can be expressed as \[ \begin{pmatrix} -\lambda_0 & \lambda_0 & 0 & 0 & \cdots \\ 0 & -\lambda_1 & \lambda_1 & 0 & \cdots \\ 0 & 0 & -\lambda_2 & \lambda_2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{pmatrix} \]

Clearly, pure birth process is equal to Poisson process on the condition that \(\lambda_i = \lambda (i=0, 1, \cdots)\).

Defining a Yule process

Yule process is another special case of pure birth process, which assumes that at any time, each individual in the population gives birth to an offspring at a constant rate \(\lambda\), independently from the rest of the population. Using the previous notation, the birth rate at state \(n\) is \(\lambda_n=\lambda n\). That is to say the birth rate of the population is proportional to the number of individuals with coefficient \(\lambda\). This characteristic makes it a suitable model for modelling the COVID patients in the initial stage.

Differential equations

The system equations for Yule process can be derived as:

\[ P_n(t+h) = P_{n-1}(t) ((n-1) \lambda h + o(h)) + P_{n}(t) (1-n \lambda h + o(h)) \]

which if we subtract \(P_n(t)\) from each side and further take limit as \(h\) goes to 0, we can get the differential equations: \[ \frac{d}{dt}P_n(t) = (n-1) \lambda P_{n-1}(t) - n \lambda P_{n}(t) \]

Typically, Yule process doesn’t start from state 0 (when no individual is in the population and birth rate is 0 based on the assumption). Under initial condition \(N(0)=n_0\), the solution to the differential equations is:

\[ P_n(t) = P(N(t)=n) = \begin{pmatrix} n-1 \\ n_0-1 \\ \end{pmatrix} e^{-\lambda n_0t} (1-e^{-\lambda t})^{n-n_0}, \]

where \(n=n_0, n_0+1, \cdots\). If you’re familiar with negative binomial distribution, you’ll surprisingly find \(N(t)\) follows a negative binomial distribution with \(k=n_0\), \(p=e^{-\lambda t}\). When \(n_0=1\), it degenerates to a geometric distribution with \(p=e^{-\lambda t}\).

Simulation of Yule process

Birth times \(T_n\) and inter-birth intervals \(T_{n+1}-T_{n}\)

Suppose \(T_n\) (\(n = n_0+1, \cdots\)) is the birth times of \(n\)-th individual in a Yule process, where \(n_0\) is the initial state of the population at time \(t=0\). Then \(T_{n+1}-T_{n}\) is the inter-birth interval during which \(n\)-many patients spread COVID to the next patient. This interval would follow an exponential distribution with rate \(\lambda_n = n \lambda\) using the according birth rate of the population. That is to say \(T_{n+1}-T_{n} \sim Exp(n\lambda)\).

Another way to interpret the event is that if we now look at \(n\) individuals after time \(T_n\) separately, each of them has a birth rate \(\lambda\), and thus, the birth time of the next individual \(T_{n+1}\) is going be the minimum birth times given by these \(n\) patients, namely \(T_{n+1}-T_{n}=\min\{Z_i : Z_1, \cdots, Z_n \overset{i.i.d}\sim Exp( \lambda ) \}\). Next, we prove the equivalence of these two explanations.

\[\begin{align} P(\min\{Z_i : Z_i \sim Exp(\lambda)\}\leq z) &= 1-P(\min\{Z_i : Z_i \sim {} Exp(\lambda)\} > z) \\ &= 1-P(Z_1 > z, Z_2 > z, \cdots, Z_n > z) \\ &= 1-(P(Z_i > z))^n \\ &= 1-e^{-n\lambda z} \end{align}\]

Simulating Yule process

Next, we’ll simulate the infected times (or birth times) of Yule process using the previous result \(T_{n+1}-T_n \sim Exp(n \lambda)\). The following Yule_sim function simulates the process of infection before time T_end from the initial state n0.

# ' @param n0 initial state
# ' @param lambda individual birth rate 
# ' @param T_end 
Yule_sim <- function(lambda, n0, T_end){
  Z = c()
  n = c()
  i <- 1
  Z[i] <- rexp(1, ((n0+i-1)*lambda))
  n[i] <- n0 + i
  while (sum(Z) < T_end){
    i = i+1
    Z[i] <- rexp(1, ((n0+i-1)*lambda))
    n[i] <- n0 + i
  }
  return(list(Z=Z, n=n))
}

Example: The trajectory of Yule process

Here we use Yule_sim function to repeat Yule process sim=100 times: in each time, the infection rate for a patient is lambda=0.2, and the initial number of patients is n0=10. We stop simulating after T_end=10, and plot the trajectory.

set.seed(1234)
par(mfrow=c(1,1))
sim <- 100
res <- Yule_sim(lambda=0.2, n0=10, T_end =10)
plot(cumsum(res$Z), res$n, type='l', col=1, ylim=c(0,100), xlab='Time', ylab='population')
for (i in 2:sim){
  res <- Yule_sim(lambda=0.2, n0=10, T_end =10)
  lines(cumsum(res$Z), res$n, type='l', col=i)
}

Example: The influence of individual brith rate (mutation of the virus)

Here we explore the influence of individual infection rate caused by mutation of the virus on \(N(t)\). We choose a set of quite similar \(\lambda=0.1, 0.15, 0.2, 0.25\), and set \(T\_end=2\).

set.seed(123)
lambda0 <- c(0.1, 0.15, 0.2, 0.25)
sim <- 1000
par(mfrow=c(2,2))
for (i in 1:length(lambda0)){
  rec <- rep(0, sim)
  for (j in 1:sim){
    res <- Yule_sim(lambda=lambda0[i], n0=5, T_end = 2)
    rec[j] <- rev(res$n)[1]
  }
  hist(rec, freq=TRUE, xlab='population', main=paste0('lambda=',lambda0[i]), breaks=10)
}

You may conclude from the plot that the individual infection rate doesn’t seem to have that much influence on \(N(t)\). Hold on! Let’s take a look at the plot after \(T\_end=10\).

set.seed(123)
lambda0 <- c(0.1, 0.15, 0.2, 0.25)
sim <- 1000
par(mfrow=c(2,2))
for (i in 1:length(lambda0)){
  rec <- rep(0, sim)
  for (j in 1:sim){
    res <- Yule_sim(lambda=lambda0[i], n0=5, T_end =10)
    rec[j] <- rev(res$n)[1]
  }
  hist(rec, freq=TRUE, xlab='population', main=paste0('lambda=',lambda0[i]), breaks=10)
}

There seems to be an obvious difference among the four plots now, especially with respect to the x axis. What about after \(T\_end=20\)?

set.seed(123)
lambda0 <- c(0.1, 0.15, 0.2, 0.25)
sim <- 1000
par(mfrow=c(2,2))
for (i in 1:length(lambda0)){
  rec <- rep(0, sim)
  for (j in 1:sim){
    res <- Yule_sim(lambda=lambda0[i], n0=5, T_end =20)
    rec[j] <- rev(res$n)[1]
  }
  hist(rec, freq=TRUE, xlab='population', main=paste0('lambda=',lambda0[i]), breaks=10)
}

The difference becomes even larger! With time going on, even a small difference of the individual infection rate caused by the mutation can result in a huge difference.

Example: Testify the distribution of \(N(t)\)

Finally, let’s testify that the distribution of \(N(t)\) follows a negative binomial distribution with parameter \(k=n_0\) and \(p=e^{-\lambda t}\). Here we set lambda=0.2, n0=10, and T_end=20.

sim <- 10000
lambda = 0.2
n0 = 5
T_end = 20
p = exp(-T_end * lambda)
rec <- rep(0, sim)
par(mfrow=c(1,1))
for (j in 1:sim){
  res <- Yule_sim(lambda, n0, T_end)
  rec[j] <- rev(res$n)[1]
}
x <- seq(min(rec), max(rec), 1)
hist(rec, probability=TRUE, xlab='population', breaks=50, main='N(t)')
lines(x,dnbinom(x, size=n0, prob=p), col=2)

The simulation results fits pretty well with the theoretical result!