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.

Keywords

This work was partially supported by the National Institutes of Health (R01 HG006139, R01 AI107034, and U54 GM111274) and the National Science Foundation (IIS 1251151, DMS 1264153, DMS 1606177). We thank Christopher Drovandi, Edwin Michael, and David Denham for access to the Brugia pahangi count data.

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

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

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

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:

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

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).

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