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.
Similar content being viewed by others
References
Abate J, Whitt W (1992) The Fourier-series method for inverting transforms of probability distributions. Queueing Syst 10(1–2):5–87
Andersson H, Britton T (2000) Stochastic epidemic models and their statistical analysis, vol 4. Springer, New York
Andrieu C, Doucet A, Holenstein R (2010) Particle Markov chain Monte Carlo methods. J R Stat Soc Ser B (Stat Methodol) 72(3):269–342
Blum MG, Tran VC (2010) HIV with contact tracing: a case study in approximate Bayesian computation. Biostatistics 11(4):644–660
Brauer F (2008) Compartmental models in epidemiology. Mathematical epidemiology. Springer, Berlin, pp 19–79
Britton T (2010) Stochastic epidemic models: a survey. Math Biosci 225(1):24–35
Cauchemez S, Ferguson N (2008) Likelihood-based estimation of continuous-time epidemic models from time-series data: application to measles transmission in London. J R Soc Interface 5(25):885–897
Correia-Gomes C, Economou T, Bailey T, Brazdil P, Alban L, Niza-Ribeiro J (2014) Transmission parameters estimated for salmonella typhimurium in swine using susceptible-infectious-resistant models and a bayesian approach. BMC Vet Res 10(1):101
Craviotto C, Jones WB, Thron W (1993) A survey of truncation error analysis for Padé and continued fraction approximants. Acta Appl Math 33(2–3):211–272
Crawford FW, Suchard MA (2012) Transition probabilities for general birth-death processes with applications in ecology, genetics, and evolution. J Math Biol 65(3):553–580
Crawford FW, Minin VN, Suchard MA (2014) Estimation for general birth-death processes. J Am Stat Assoc 109(506):730–747
Crawford FW, Weiss RE, Suchard MA (2015) Sex, lies, and self-reported counts: Bayesian mixture models for longitudinal heaped count data via birth-death processes. Ann Appl Stat 9:572–596
Crawford FW, Stutz TC, Lange K (2016) Coupling bounds for approximating birth-death processes by truncation. Stat Prob Lett 109:30–38
Csilléry K, Blum MG, Gaggiotti OE, François O (2010) Approximate Bayesian computation (ABC) in practice. Trends Ecol Evol 25(7):410–418
Doss CR, Suchard MA, Holmes I, Kato-Maeda M, Minin VN (2013) Fitting birth-death processes to panel data with applications to bacterial DNA fingerprinting. Ann Appl Stat 7(4):2315–2335
Drovandi CC, Pettitt AN (2011) Estimation of parameters for macroparasite population evolution using approximate Bayesian computation. Biometrics 67(1):225–233
Dukic V, Lopes HF, Polson NG (2012) Tracking epidemics with Google flu trends data and a state-space SEIR model. J Am Stat Assoc 107(500):1410–1426
Earn DJ (2008) A light introduction to modelling recurrent epidemics. Mathematical epidemiology. Springer, Berlin, pp 3–17
Ephraim Y, Mark BL (2012) Bivariate Markov processes and their estimation. Found Trends Signal Process 6(1):1–95
Feller W (1968) An introduction to probability theory and its applications, vol 1. Wiley, Hoboken
Finkenstädt B, Grenfell B (2000) Time series modelling of childhood diseases: a dynamical systems approach. J R Stat Soc Ser C (Appl Stat) 49(2):187–205
Golightly A, Wilkinson DJ (2005) Bayesian inference for stochastic kinetic models using a diffusion approximation. Biometrics 61(3):781–788
Griffiths D (1972) A bivariate birth-death process which approximates to the spread of a disease involving a vector. J Appl Probab 9(1):65–75
Hitchcock S (1986) Extinction probabilities in predator-prey models. J Appl Probab 23(1):1–13
Iglehart DL (1964) Multivariate competition processes. Ann Math Stat 35(1):350–361
Ionides E, Bretó C, King A (2006) Inference for nonlinear dynamical systems. Proc Natl Acad Sci USA 103(49):18,438–18,443
Ionides EL, Nguyen D, Atchadé Y, Stoev S, King AA (2015) Inference for dynamic and latent variable models via iterated, perturbed Bayes maps. Proce Natl Acad Sci USA 112(3):719–724
Jahnke T, Huisinga W (2007) Solving the chemical master equation for monomolecular reaction systems analytically. J Math Biol 54(1):1–26
Karev GP, Berezovskaya FS, Koonin EV (2005) Modeling genome evolution with a diffusion approximation of a birth-and-death process. Bioinformatics 21(Suppl 3):iii12–iii19
Keeling M, Ross J (2008) On methods for studying stochastic disease dynamics. J R Soc Interface 5(19):171–181
Kermack W, McKendrick A (1927) A contribution to the mathematical theory of epidemics. Proc R Soc Lond Ser A 115(772):700–721
Lentz WJ (1976) Generating Bessel functions in Mie scattering calculations using continued fractions. Appl Opt 15(3):668–671
Levin D (1973) Development of non-linear transformations for improving convergence of sequences. Int J Comput Math 3(1–4):371–388
Martin AD, Quinn KM, Park JH (2011) MCMCpack: Markov chain Monte Carlo in R. J Stat Softw 42(9):22. http://www.jstatsoft.org/v42/i09/
McKendrick A (1926) Applications of mathematics to medical problems. Proce Edinb Math Soc 44:98–130
McKinley T, Cook AR, Deardon R (2009) Inference in epidemic models without likelihoods. Int J Biostat 5(1):1557–4679
Moler C, Loan C (2003) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Rev 45:3–49
Murphy J, O’Donohoe M (1975) Some properties of continued fractions with applications in Markov processes. IMA J Appl Math 16(1):57–71
Novozhilov AS, Karev GP, Koonin EV (2006) Biological applications of the theory of birth-and-death processes. Brief Bioinform 7(1):70–85
Owen J, Wilkinson DJ, Gillespie CS (2015) Scalable inference for Markov processes with intractable likelihoods. Stat Comput 25(1):145–156
Rabier CE, Ta T, Ané C (2014) Detecting and locating whole genome duplications on a phylogeny: a probabilistic approach. Mol Biol Evol 31(3):750–762
Raggett G (1982) A stochastic model of the Eyam plague. J Appl Stat 9(2):212–225
Renshaw E (2011) Stochastic population processes: analysis, approximations, simulations. Oxford University Press, Oxford
Reuter GEH (1957) Denumerable Markov processes and the associated contraction semigroups on l. Acta Math 97(1):1–46
Reuter GEH (1961) Competition processes. In: Proc. 4th Berkeley Symp. Math. Statist. Prob, vol 2, pp 421–430
Riley S, Donnelly CA, Ferguson NM (2003) Robust parameter estimation techniques for stochastic within-host macroparasite models. J Theor Biol 225(4):419–430
Robert CP, Cornuet JM, Marin JM, Pillai NS (2011) Lack of confidence in approximate Bayesian computation model choice. Proc Natl Acad Sci 108(37):15,112–15,117
Rosenberg NA, Tsolaki AG, Tanaka MM (2003) Estimating change rates of genetic markers using serial samples: applications to the transposon IS6110 in Mycobacterium tuberculosis. Theor Popul Biol 63(4):347–363
Schranz HW, Yap VB, Easteal S, Knight R, Huttley GA (2008) Pathological rate matrices: from primates to pathogens. BMC Bioinform 9(1):550
Sidje RB (1998) Expokit: a software package for computing matrix exponentials. ACM Trans Math Softw (TOMS) 24(1):130–156
Sunnåker M, Busetto AG, Numminen E, Corander J, Foll M, Dessimoz C (2013) Approximate Bayesian computation. PLoS Comput Biol 9(1):e1002,803
Thompson I, Barnett A (1986) Coulomb and Bessel functions of complex arguments and order. J Comput Phys 64(2):490–509
van den Eshof J, Hochbruck M (2006) Preconditioning lanczos approximations to the matrix exponential. SIAM J Sci Comput 27(4):1438–1457
Xu J, Guttorp P, Kato-Maeda M, Minin VN (2015) Likelihood-based inference for discretely observed birth-death-shift processes, with applications to evolution of mobile genetic elements. Biometrics 71(4):1009–1021
Author information
Authors and Affiliations
Corresponding author
Additional information
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.
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
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
Definition A3
We define the corresponding sequence \(\{\phi _n\}_{n=0}^\infty \) of a continued fraction (A.1) by the following recurrence formulae
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,
Now, if we consider a more general recurrence formulae
then under the assumption of Lemma A.4, we have the following lemma:
Lemma A2
The solution for (A.5) is
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
to stabilize the computation of \(\phi _0^{(n)}\). We can calculate \(A_n\), \(B_n\), and \(\phi _0^{(n)}\) recursively as follows:
If \(\phi _0^{(n)}\) converges to \(\phi _0\), then Craviotto et al. (1993) show that
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
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
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
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
We may expand the probability generating functions in the following form:
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:
To derive the backwards and forward equations, Chapman–Kolmogorov arguments yield the symmetric relations
First, we derive the backward equations by expanding around t and applying (D.3):
Since an analogous argument applies for \(\phi _2\), we arrive at the system
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
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:
Plugging into the backward equations, we obtain
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
Plugging in \(\phi _2(0, s_1, s_2) = s_2\), we obtain \(C = \ln (1-s_2)\), and we arrive at
Substituting this solution into the first differential equation and applying the integrating factor method provides
Plugging in the initial condition \(\phi _1(0,s_1,s_2) = s_1\) and rearranging yields
1.2 D.2 Transition probability expressions
Transition probabilities are related to the PGF via repeated partial differentiation; note that
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
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).
Rights and permissions
About this article
Cite this article
Ho, L.S.T., Xu, J., Crawford, F.W. et al. Birth/birth-death processes and their computable transition probabilities with biological applications. J. Math. Biol. 76, 911–944 (2018). https://doi.org/10.1007/s00285-017-1160-3
Received:
Revised:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00285-017-1160-3