Birth/birth-death processes and their computable transition probabilities with biological applications

Abstract

Birth-death processes track the size of a univariate population, but many biological systems involve interaction between populations, necessitating models for two or more populations simultaneously. A lack of efficient methods for evaluating finite-time transition probabilities of bivariate processes, however, has restricted statistical inference in these models. Researchers rely on computationally expensive methods such as matrix exponentiation or Monte Carlo approximation, restricting likelihood-based inference to small systems, or indirect methods such as approximate Bayesian computation. In this paper, we introduce the birth/birth-death process, a tractable bivariate extension of the birth-death process, where rates are allowed to be nonlinear. We develop an efficient algorithm to calculate its transition probabilities using a continued fraction representation of their Laplace transforms. Next, we identify several exemplary models arising in molecular epidemiology, macro-parasite evolution, and infectious disease modeling that fall within this class, and demonstrate advantages of our proposed method over existing approaches to inference in these models. Notably, the ubiquitous stochastic susceptible-infectious-removed (SIR) model falls within this class, and we emphasize that computable transition probabilities newly enable direct inference of parameters in the SIR model. We also propose a very fast method for approximating the transition probabilities under the SIR model via a novel branching process simplification, and compare it to the continued fraction representation method with application to the 17th century plague in Eyam. Although the two methods produce similar maximum a posteriori estimates, the branching process approximation fails to capture the correlation structure in the joint posterior distribution.

Appendices

A Continued fractions

In this section, we give some basic definitions and properties related to continued fractions.

Definition A1

A continued fraction $\phi _0$ is a scalar quantity expressed in

$$\begin{aligned} \phi _0 = \frac{x_1}{ y_1 + \frac{x_2}{ y_2 + \frac{x_3}{ y_3 + \cdots }~,}} \end{aligned}$$

(A.1)

where $\{ x_i \}_{i=1}^\infty$ and $\{ y_i \}_{i=1}^\infty$ are infinite sequences of complex numbers.

Definition A2

The $n^{\mathrm{th}}$ convergent of $\phi _0$ is

$$\begin{aligned} \frac{X_n}{Y_n} = \frac{x_1}{ y_1 + \frac{x_2}{ y_2 + \frac{x_3}{ y_3 + \cdots + \frac{x_n}{y_n }~.}}} \end{aligned}$$

(A.2)

Definition A3

We define the corresponding sequence $\{\phi _n\}_{n=0}^\infty$ of a continued fraction (A.1) by the following recurrence formulae

$$\begin{aligned} \begin{aligned}&\phi _1 = x_1 - y_1\phi _0,~\text{ and } \\&\phi _n = x_n \phi _{n-2} - y_n \phi _{n-1}~\text{ for }~n \ge 2. \end{aligned} \end{aligned}$$

(A.3)

Murphy and O’Donohoe (1975) provided the following sufficient condition for the convergence of (A.1):

Lemma A1

Assume that there exists N such that $\inf _{n>N}|Y_n| > 0$ and $\lim _{n \rightarrow \infty } \phi _n = 0$. Then, the continued fraction (A.1) is convergent. Moreover,

$$\begin{aligned} \phi _n = \prod _{i=1}^{n}{x_i} \frac{x_{n+1}}{ Y_{n+1} + \frac{x_{n+2} Y_n}{ y_{n+2} + \frac{x_{n+3}}{ y_{n+3} + \frac{x_{n+4}}{ y_{n+4} + \cdots }~.}}} \end{aligned}$$

(A.4)

Now, if we consider a more general recurrence formulae

$$\begin{aligned} \begin{aligned}&\phi ^{(m)}_1 = - y_1\phi ^{(m)}_0 + k_1 1_{\{m = 0\}} \\&\phi ^{(m)}_n = x_n \phi ^{(m)}_{n-2} - y_n \phi ^{(m)}_{n-1} + k_{m+1} 1_{\{m = n - 1\}}~\text {for}~n \ge 2, \end{aligned} \end{aligned}$$

(A.5)

then under the assumption of Lemma A.4, we have the following lemma:

Lemma A2

The solution for (A.5) is

$$\begin{aligned} \phi ^{(m)}_n = {\left\{ \begin{array}{ll} \frac{(-1)^{m-n}k_{m+1}}{\prod _{i=1}^{m+1}{x_i}} Y_n \phi _m, &{} \text{ if } n \le m \\ \frac{k_{m+1}}{\prod _{i=1}^{m+1}{x_i}} Y_m \phi _n, &{} \text{ if } n \ge m. \end{array}\right. } \end{aligned}$$

(A.6)

B Modified Lentz method

Modified Lentz method (Lentz 1976; Thompson and Barnett 1986) is an efficient algorithm to finitely approximate the infinite expression of the continued fraction $\phi _0$ in (A.1) to within a prescribed error tolerance. Let $\phi _0^{(n)}$ be the $n^{{ \mathrm \tiny th}}$ convergence of $\phi _0$, that is $\phi _0^{(n)} = X_n/Y_n.$ The main idea of Lentz’s algorithm lies in using the ratios

$$\begin{aligned} A_n = \frac{X_n}{X_{n-1}} ~~ \text{ and } ~~ B_n = \frac{Y_{n-1}}{Y_n} \end{aligned}$$

(B.1)

to stabilize the computation of $\phi _0^{(n)}$. We can calculate $A_n$, $B_n$, and $\phi _0^{(n)}$ recursively as follows:

$$\begin{aligned} \begin{aligned} A_n&= y_n + \frac{x_n}{A_{n-1}} \\ B_n&= \frac{1}{y_n + x_n B_{n-1}} \\ \phi _0^{(n)}&= \phi _0^{(n-1)} A_n B_n. \end{aligned} \end{aligned}$$

(B.2)

If $\phi _0^{(n)}$ converges to $\phi _0$, then Craviotto et al. (1993) show that

$$\begin{aligned} \left| \phi _0^{(n)} - \phi _0 \right| \le \frac{|Y_n/Y_{n-1}|}{\mathcal{I} [Y_n/Y_{n-1}]} \left| \phi _0^{(n)} - \phi _0^{(n-1)} \right| = \frac{|1/B_n|}{\mathcal{I} [1/B_n]} \left| \phi _0^{(n)} - \phi _0^{(n-1)} \right| , \end{aligned}$$

(B.3)

where $\mathcal{I} [Y_n/Y_{n-1}]$ is the imaginary part of $Y_n/Y_{n-1}$ and is assumed to be non-zero. Hence, the Lentz’s algorithm terminates when

$$\begin{aligned} \frac{|1/B_n|}{\mathcal{I} [1/B_n]} \left| \phi _0^{(n)} - \phi _0^{(n-1)} \right| \end{aligned}$$

(B.4)

is small enough. However, $A_n$ and $B_n$ can equal zero themselves and cause problem. Hence, Thompson and Barnett (1986) propose a modification for Lentz’s algorithm by setting $A_n$ and $B_n$ to a very small number, such as $10^{-16}$, whenever they equal zero. In practice, the algorithm often terminates after small number of iterations. However, in some rare cases where the numerical computation is unstable, it might take too long before the algorithm terminates. So, we set a predefined maximum number of iterations H as a fallback for these cases.

C Convergence results of increasing the truncation level

Let $f_{ab}^{(B)}(s)$ be the output of the approximation scheme (19) in Theorem 2. In this section, we prove that $f_{ab}^{(B)}(s)$ converges to $f_{ab}(s)$ as B goes to infinity. To do so, let us consider a truncated birth/birth-death process $\mathbf {X}^{(B)}(t) = (X^{(B)}_1(t), X^{(B)}_2(t) )$ at truncation level B such that it executes the same process as $\mathbf {X}(t)$ on the state $\{ a_0, a_0 + 1, a_0 + 2, \ldots \} \times \{ 0, 1, 2, \ldots , B\}$ except that $\lambda ^{(2)}_{aB} = 0$. Define $P_{ab}^{a_0 b_0, (B)}(t)$ be the transition probabilities of $\mathbf {X}^{(B)}(t)$ and $T_B$ be the hitting time at which $X_2(t)$ first reach state $B+1$. For any set $S \subset {\mathbb {N}}^2$, we have

$$\begin{aligned} \Pr (\mathbf {X}(t) \in S)&= \Pr (\mathbf {X}(t) \in S |T_B> t) \Pr (T>t) {+} \Pr (\mathbf {X}(t) \in S|T_B \le t) \Pr (T_B \le t) \\&= \Pr (\mathbf {X}^{(B)}(t) \in S) \Pr (T_B > t) + \Pr (\mathbf {X}(t) \in S ~|~ T_B \le t) \Pr (T_B \le t) \\&= \Pr (\mathbf {X}^{(B)}(t) \in S) + [\Pr (\mathbf {X}(t) \in S ~|~ T_B \le t)\\&\quad - \Pr (\mathbf {X}^{(B)}(t) \in S)] \Pr (T_B \le t) \end{aligned}$$

Therefore $|\Pr (\mathbf {X}(t) \in S) - \Pr (\mathbf {X}^{(B)}(t) \in S)| \le \Pr (T_B \le t)$. Note that $f_{ab}^{(B)}(s)$ is the Laplace transform of $P_{ab}^{a_0 b_0, (B)}(t)$. Hence

$$\begin{aligned} |f_{ab}^{(B)}(s)- & {} f_{ab}(s)| \le \int _{0}^\infty {|P_{ab}^{a_0 b_0, (B)}(t) - P_{ab}^{a_0 b_0}(t)| e^{-st}dt}\nonumber \\\le & {} \quad \int _{0}^\infty {Pr(T_B \le t) e^{-st}dt} \end{aligned}$$

(C.1)

By Dominated convergence theorem and the fact that $\lim _{B \rightarrow \infty } Pr(T_B \le t) = 0$, we deduce that $\lim _{B \rightarrow \infty } f_{ab}^{(B)}(s) = f_{ab}(s)$.

D Branching SIR approximation

Here we derive and solve the Kolmogorov backward equations of the two-type branching process necessary for evaluating the probability generating functions (PGFs) whose coefficients yield transition probabilities.

1.1 D.1 Deriving the PGF

Our two-type branching process is represented by a vector $(X_1(t), X_2(t))$ that denotes the numbers of particles of two types at time t. Let the quantities $a_1(k,l)$ denote the rates of producing k type 1 particles and l type 2 particles, starting with one type 1 particle, and $a_2(k,l)$ be analogously defined but beginning with one type 2 particle. Given a two-type branching process defined by instantaneous rates $a_i(k,l)$, denote the following pseudo-generating functions for $i = 1,2$ as

$$\begin{aligned} u_i(s_1,s_2) = \sum _k \sum _l a_i(k,l)s_1^k s_2^l . \end{aligned}$$

(D.1)

We may expand the probability generating functions in the following form:

$$\begin{aligned} \phi _{10}(t, s_1, s_2)&= E (s_1^{X_1(t)} s_2^{X_2(t)} | X_1(0) = 1, X_2(0) = 0) \nonumber \\&= \sum _{k=0}^\infty \sum _{l=0}^\infty P_{1,0}^{kl} (t) s_1^k s_2^l \nonumber \\&= \sum _{k=0}^\infty \sum _{l=0}^\infty ( \mathbf {1}_{k=1, l = 0} + a_1(k,l) t + o(t) ) \nonumber s_1^k s_2^l \\&= s_1 + u_1(s_1, s_2) t + o(t) . \end{aligned}$$

(D.2)

We have an analogous expression for $\phi _{01}(t, s_1, s_2)$ beginning with one particle of type 2 instead of type 1. For short, we will write $\phi _{10} := \phi _1, \phi _{01} := \phi _2$. Thus, we have the following relation between the functions $\phi$ and u:

$$\begin{aligned} \frac{d \phi _1}{dt} (t, s_1, s_2) |_{t=0}&= u_1(s_1, s_2) \text { and}\nonumber \\ \frac{d \phi _2}{dt} (t, s_1, s_2) |_{t=0}&= u_2(s_1, s_2) . \end{aligned}$$

(D.3)

To derive the backwards and forward equations, Chapman–Kolmogorov arguments yield the symmetric relations

$$\begin{aligned} \phi _1(t+h, s_1, s_2)&= \phi _1(t, \phi _1(h, s_1, s_2), \phi _2(h, s_1, s_2)) \nonumber \\&= \phi _1(h, \phi _1(t, s_1, s_2), \phi _2(t, s_1, s_2)) . \end{aligned}$$

(D.4)

First, we derive the backward equations by expanding around t and applying (D.3):

$$\begin{aligned} \phi _1(t+h, s_1, s_2)&= \phi _1(t, s_1, s_2) + \frac{d \phi _1}{dh}(t+h, s_1, s_2) | _{h=0} h + o(h) \nonumber \\&= \phi _1(t, s_1, s_2) + \frac{d \phi _1}{dh}(h, \phi _1(t, s_1, s_2), \phi _2(t, s_1, s_2) |_{h=0} h + o(h) \nonumber \\&= \phi _1(t,s_1,s_2) + u_1( \phi _1(t,s_1, s_2) , \phi _2(t, s_1, s_2) h + o(h) ) . \end{aligned}$$

(D.5)

Since an analogous argument applies for $\phi _2$, we arrive at the system

$$\begin{aligned} \frac{d}{dt} \phi _1(t, s_1, s_2)&= u_1( \phi _1(t, s_1, s_2), \phi _2(t, s_1, s_2) ) \text { and} \nonumber \\ \frac{d}{dt} \phi _2(t, s_1, s_2)&= u_2( \phi _1(t, s_1, s_2), \phi _2(t, s_1, s_2) ) , \end{aligned}$$

(D.6)

with initial conditions $\phi _1(0, s_1, s_2) = s_1, \phi _2(0, s_1, s_2) = s_2$.

Recall in our SIR approximation, we use the initial population $X_2(0)$ as a constant that scales the instantaneous rates over any time interval $[t_0, t_1)$. The only nonzero rates specifying this proposed model, in the notation above, are

$$\begin{aligned} a_1(0,1) = \beta X_2(0), \quad \quad a_1(1,0) = -\beta X_2(0), \quad \quad a_2(0,1) = -\alpha , \quad \quad a_2(0,0) = \alpha . \nonumber \\ \end{aligned}$$

(D.7)

For simplicity, call $X_2(0) := I_0$, the constant representing the infected population at the beginning of the time interval. Thus, the corresponding pseudo-generating functions have a simple form:

$$\begin{aligned} u_1(s_1, s_2) = \beta I_0 s_2 - \beta I_0 s_1 \text { and} \nonumber \\ u_2(s_1, s_2) = \alpha - \alpha s_2 = \alpha (1 - s_2) . \end{aligned}$$

(D.8)

Plugging into the backward equations, we obtain

$$\begin{aligned} \frac{d}{dt} \phi _1(t,s_1,s_2)&= \beta I_0 \big ( \phi _2(t, s_1, s_2) - \phi _1(t, s_1, s_2) \big ) \text { and} \nonumber \\ \frac{d}{dt} \phi _2(t,s_1, s_2)&= \alpha - \alpha \phi _2(t,s_1,s_2). \end{aligned}$$

(D.9)

The $\phi _2$ differential equation corresponds to a pure death process and is immediately solvable; suppressing the arguments of $\phi _2$ for notational convenience, we obtain

$$\begin{aligned} \frac{d}{dt} \phi _2&= \alpha - \alpha \phi _2 \nonumber \\ \frac{d}{dt} \phi _2( \frac{1}{1 - \phi _2})&= \alpha \nonumber \\ \ln (1 - \phi _2)&= -\alpha t + C \nonumber \\ \phi _2&= 1 - \exp (-\alpha t + C) . \end{aligned}$$

(D.10)

Plugging in $\phi _2(0, s_1, s_2) = s_2$, we obtain $C = \ln (1-s_2)$, and we arrive at

$$\begin{aligned} \phi _2(t, s_1, s_2) = 1 + (s_2 - 1)\exp (-\alpha t) \end{aligned}$$

(D.11)

Substituting this solution into the first differential equation and applying the integrating factor method provides

$$\begin{aligned} \phi _1 e^{\beta I_0 t}&= \int \beta I_0 e^{\beta I_0 t} (1 + \frac{s_2-1}{e^{\alpha t}}) \, dt = e^{\beta I_0 t} + \beta I_0 (s_2 -1) \int e^{(\beta I_0 - \alpha )t} \, dt \nonumber \\&= e^{\beta I_0 t} + \beta I_0 (s_2-1) \frac{e^{(\beta I_0 - \alpha ) t}}{\beta I_0 - \alpha } + C . \end{aligned}$$

(D.12)

Plugging in the initial condition $\phi _1(0,s_1,s_2) = s_1$ and rearranging yields

$$\begin{aligned} \phi _1 = 1 + \frac{ \beta I_0 (s_2 - 1)}{\beta I_0 - \alpha } e^{-\alpha t} + e^{-\beta I_0 t} ( s_1 - 1 - \frac{\beta I_0 (s_2 - 1) }{\beta I_0 - \alpha } ) . \end{aligned}$$

(D.13)

1.2 D.2 Transition probability expressions

Transition probabilities are related to the PGF via repeated partial differentiation; note that

$$\begin{aligned} P_{kl}^{mn}(t)&= \frac{1}{k!}\frac{1}{l!}\frac{\partial ^k}{\partial s_1^k} \frac{\partial ^l}{\partial s_2^l} \phi _{mn}(t, s_1, s_2) \bigg |_{s_1=s_2=0}\nonumber \\&= \frac{1}{k!}\frac{1}{l!}\frac{\partial ^k}{\partial s_1^k} \frac{\partial ^l}{\partial s_2^l} \phi _1^m(t, s_1, s_2) \phi _2^n(t, s_1, s_2) \bigg |_{s_1=s_2=0} \nonumber \\&= \frac{\partial ^l}{\partial s_2^l} \sum _{i=0}^k {k \atopwithdelims ()i} \frac{ \partial ^{k-i}}{\partial s_1^ {k-i} } \phi _1^m(t, s_1, s_2) \frac{ \partial ^i}{\partial s_1^i} \phi _2^n(t, s_1, s_2) \bigg |_{s_1=s_2=0} . \end{aligned}$$

(D.14)

This expression is generally unwieldy, but notice $\frac{ \partial ^i}{\partial s_1^i} \phi _2^n(t, s_1, s_2) \bigg |_{s_1 = 0} = 0 \text { for all } i > 0$ in our model. Remarkably, this allows us to further simplify and ultimately arrive at closed-form expressions. Continuing, we see

$$\begin{aligned} P_{kl}^{mn}(t)&= \frac{\partial ^l}{\partial s_2^l} \left[ { k \atopwithdelims ()0 } \phi _2^n(t, s_1, s_2) \frac{ \partial ^k}{\partial s_1^k} \phi _1^m(t, s_1, s_2) \right] \bigg |_{s_1=s_2=0} \nonumber \\&= \frac{\partial ^l}{\partial s_2^l} \bigg \{ \phi _2^n(t, s_1, s_2) \cdot \frac{m!}{(m-k)!} e^{-k \beta I_0 t} \bigg [ 1 + \frac{ \beta I_0 (s_2 -1)}{\beta I_0 - \alpha } e^{-\alpha t} \nonumber \\&- e^{-\beta I_0 t} \bigg ( 1 + \frac{ \beta I_0 (s_2 - 1)}{ \beta I_0 - \alpha } \bigg ) \bigg ] ^{m-k} \bigg \} \bigg |_{s_1=s_2=0} \nonumber \\&:= \frac{\partial ^l}{\partial s_2^l} \left[ \phi _2^n(t, s_1, s_2) \cdot h(t, s_1, s_2) \right] \bigg |_{s_1=s_2=0} \nonumber \\&= \sum _{i=0}^l {l \atopwithdelims ()i} \frac{ \partial ^{l-i}}{\partial s_2^ {l-i}} h(t, s_1, s_2) \frac{\partial ^i}{\partial s_2^i} \phi _2^n(t, s_1, s_2) \nonumber \\&:= \sum _{i=0}^l {l \atopwithdelims ()i} A(l-i) B(i) . \end{aligned}$$

(D.15)

From here, it is straightforward to take partial derivatives of $h(t, s_1, s_2)$ and our closed-form expression of $\phi _2^n(t, s_1, s_2)$ to arrive at Conditions (41) and (42). A heatmap visualization of the difference between transition probabilities under the branching approximation and those computed using the continued fraction method for the SIR model is included below (Fig. 8).

Fig. 8

Heatmap visualizations of transition probabilities near the region of support across methods for $t=0.5, 1$. We see that the branching approximation is noticeably different from the Monte Carlo ground truth when we increase t to 1, while the continued fraction approach remains accurate