Search

Multi-scale systems often exhibit stochastic and deterministic dynamics. Capturing these aspects in a compartmental model is challenging. Notably, low occupancy compartments exhibit stochastic dynamics and high occupancy compartments exhibit deterministic dynamics. Failing to account for stochasticity in small populations can produce 'atto-foxes', e.g. in the Lotka-Volterra ordinary differential equation (ODE) model. This limitation becomes problematic when studying extinction of species or the clearance of infection, but it can be resolved by using discrete stochastic models e.g. continuous time Markov chains (CTMCs). Unfortunately, simulating CTMCs is infeasible for most realistic populations. We develop a novel mathematical framework, to couple continuous ODEs and discrete CTMCs: 'Jump-Switch-Flow' (JSF). In this framework populations can reach extinct states ("absorbing states"), thereby resolving atto-fox-type problems. JSF has the desired behaviours of exact CTMC simulation, but is substantially faster than existing alternatives. JSF's utility for simulation-based inference, particularly multi-scale problems, is demonstrated by several case-studies. In a simulation study, we demonstrate how JSF can enable a more nuanced analysis of the efficacy of public health interventions. We also carry out a novel analysis of longitudinal within-host data from SARS-CoV-2 infections to quantify the timing of viral clearance. JSF offers a novel approach to compartmental model development and simulation.

Phylogenetic trees describe the relationships between species in the evolutionary process, and provide information about the rates of diversification. To understand the mechanisms behind macroevolution, we consider a class of multitype branching processes called Markovian binary trees (MBTs). MBTs allow for trait-based variation in diversification rates, and provide a flexible and realistic probabilistic model for phylogenetic trees. We develop an approximate Bayesian computation (ABC) scheme to infer the rates of MBT parameters by exploiting the information in the shapes of phylogenetic trees. We evaluate the accuracy of this inference method using simulation studies, and find that our method is able to detect variation in the diversification rates, with accuracy comparable to, and generally better than, likelihood-based methods. In an application to a real-life phylogeny of squamata, we reinforce conclusions drawn from earlier studies, in particular supporting the existence of ovi-/viviparity transitions in both directions. Our method demonstrates the potential for more complex models of evolution to be employed in phylogenetic inference, in conjunction with likelihood-free schemes., Comment: 43 pages (17 pages for the main text)

Population-size dependent branching processes (PSDBP) and controlled branching processes (CBP) are two classes of branching processes widely used to model biological populations that exhibit logistic growth. In this paper we develop connections between the two, with the ultimate goal of determining when a population is more appropriately modelled with a PSDBP or a CBP. In particular, we state conditions for the existence of equivalent PSDBPs and CBPs, we then consider the subclass of CBPs with deterministic control functions (DCBPs), stating a necessary and sufficient condition for DCBP-PSDBP equivalence. Finally, we derive an upper bound on the total variation distance between non-equivalent DCBPs and PSDBPs with matching first and second moments and equal initial population size, and show that under certain conditions this bound tends to zero as the initial population size becomes large.

We derive the first conditionally consistent estimators for a class of parametric Markov population models with logistic growth, which are suitable for modelling endangered populations in restricted habitats with a carrying capacity. We focus on discrete-time parametric population-size-dependent branching processes, for which we propose a new class of weighted least-squares estimators based on a single trajectory of population size counts. We establish the consistency and asymptotic normality of our estimators, conditional on non-extinction up to time $n$, as $n\to\infty$. Since Markov population models with a carrying capacity become extinct almost surely under general conditions, our proofs rely on arguments distinct from those in the existing literature. Our results are motivated by conservation biology, where endangered populations are often studied precisely because they are still alive, leading to an observation bias. Through simulated examples, we show that our conditionally consistent estimators generally reduce this bias for key quantities such as a habitat's carrying capacity. We apply our methodology to estimate the carrying capacity of the Chatham Island black robin, a population reduced to a single breeding female in the 1970's, which has since recovered but has yet to reach the island's carrying capacity.

We introduce a class of branching processes in which the reproduction or lifetime distribution at a given time depends on the total cumulative number of individuals who have been born in the population until that time. We focus on a continuous-time version of these processes, called total-progeny-dependent birth-and-death processes, and study some of their properties through the analysis of their fluid (deterministic) approximation. These properties include the maximum population size, the total progeny size at extinction, the time to reach the maximum population size, and the time until extinction. As the fluid approach does not allow us to approximate the time until extinction directly, we propose several methods to complement this approach. We also use the fluid approach to study the behaviour of the processes as we increase the magnitude of the individual birth rate.

Birth-and-death processes (BDPs) form a class of continuous-time Markov chains that are particularly suited to describing the changes in the size of a population over time. Population-size-dependent BDPs (PSDBDPs) allow the rate at which a population grows to depend on the current population size. The main purpose of our new Python package BirDePy is to provide easy-to-use functions that allow the parameters of discretely-observed PSDBDPs to be estimated. The package can also be used to estimate parameters of continuously-observed PSDBDPs, simulate sample paths, approximate transition probabilities, and generate forecasts. We describe in detail several methods which have been incorporated into BirDePy to achieve each of these tasks. The usage and effectiveness of the package is demonstrated through a variety of examples of PSDBDPs, as well as case studies involving annual population count data of two endangered bird species.

Stochastic models for performance analysis, optimization and control of queues hinge on a multitude of alternatives for input point processes. In case of bursty traffic, one very popular model is the \textit{Markov Modulated Poisson Process} (MMPP), however it is not the only option. Here, we introduce an alternative that we call \textit{Markovian transition counting process} (MTCP). The latter is a point process counting the number of transitions of a finite continuous-time Markov chain. For a given MTCP one can establish an MMPP with the same first and second moments of counts. In this paper, we show the other direction by establishing a duality in terms of first and second moments of counts between MTCPs and a rich class of MMPPs which we refer to as slow MMPPs (modulation is slower than the events). Such a duality confirms the applicability of the MTCP as an alternative to the MMPP which is superior when it comes to moment matching and finding the important measures of the inter-event process. We illustrate the use of such equivalence in a simple queueing example, showing that the MTCP is a comparable and competitive model for performance analysis.

We consider Galton-Watson branching processes with countable typeset $\mathcal{X}$. We study the vectors ${\bf q}(A)=(q_x(A))_{x\in\mathcal{X}}$ recording the conditional probabilities of extinction in subsets of types $A\subseteq \mathcal{X}$, given that the type of the initial individual is $x$. We first investigate the location of the vectors ${\bf q}(A)$ in the set of fixed points of the progeny generating vector and prove that $q_x(\{x\})$ is larger than or equal to the $x$th entry of any fixed point, whenever it is different from 1. Next, we present equivalent conditions for $q_x(A)< q_x (B)$ for any initial type $x$ and $A,B\subseteq \mathcal{X}$. Finally, we develop a general framework to characterise all \emph{distinct} extinction probability vectors, and thereby to determine whether there are finitely many, countably many, or uncountably many distinct vectors. We illustrate our results with examples, and conclude with open questions., Comment: 33 pages, 7 figures

We consider population-size-dependent branching processes (PSDBPs) which eventually become extinct with probability one. For these processes, we derive maximum likelihood estimators for the mean number of offspring born to individuals when the current population size is $z\geq 1$. As is standard in branching process theory, an asymptotic analysis of the estimators requires us to condition on non-extinction up to a finite generation $n$ and let $n\to\infty$; however, because the processes become extinct with probability one, we are able to demonstrate that our estimators do not satisfy the classical consistency property ($C$-consistency). This leads us to define the concept of $Q$-consistency, and we prove that our estimators are $Q$-consistent and asymptotically normal. To investigate the circumstances in which a $C$-consistent estimator is preferable to a $Q$-consistent estimator, we then provide two $C$-consistent estimators for subcritical Galton-Watson branching processes. Our results rely on a combination of linear operator theory, coupling arguments, and martingale methods.

We present a new algorithm for computing the quasi-stationary distribution of subcritical Galton--Watson branching processes. This algorithm is based on a particular discretization of a well-known functional equation that characterizes the quasi-stationary distribution of these processes. We provide a theoretical analysis of the approximate low-rank structure that stems from this discretization, and we extend the procedure to multitype branching processes. We use numerical examples to demonstrate that our algorithm is both more accurate and more efficient than other approaches.

We consider a class of multitype Galton-Watson branching processes with a countably infinite type set $\mathcal{X}_d$ whose mean progeny matrices have a block lower Hessenberg form. For these processes, the probability $\boldsymbol{q}(A)$ of extinction in subsets of types $A\subseteq \mathcal{X}_d$ may differ from the global extinction probability $\boldsymbol{q}$ and the partial extinction probability $\tilde{\boldsymbol{q}}$. After deriving partial and global extinction criteria, we develop conditions for $\boldsymbol{q}<\boldsymbol{q}(A)<\tilde{\boldsymbol{q}}$. We then present an iterative method to compute the vector $\boldsymbol{q}(A)$ for any set $A$. Finally, we investigate the location of the vectors $\boldsymbol{q}(A)$ in the set of fixed points of the progeny generating vector.

Birth-and-death processes are widely used to model the development of biological populations. Although they are relatively simple models, their parameters can be challenging to estimate, because the likelihood can become numerically unstable when data arise from the most common sampling schemes, such as annual population censuses. Simple estimators may be based on an embedded Galton-Watson process, but this presupposes that the observation times are equi-spaced. We estimate the birth, death, and growth rates of a linear birth-and-death process whose population size is periodically observed via an embedded Galton-Watson process, and by maximizing a saddlepoint approximation to the likelihood. We show that a Gaussian approximation to the saddlepoint-based likelihood connects the two approaches, we establish consistency and asymptotic normality of quasi-likelihood estimators, compare our estimators on some numerical examples, and apply our results to census data for two endangered bird populations and the H1N1 influenza pandemic.

We consider a class of branching processes with countably many types which we refer to as Lower Hessenberg branching processes. These are multitype Galton-Watson processes with typeset $\mathcal{X}=\{0,1,2,\dots\}$, in which individuals of type $i$ may give birth to offspring of type $j\leq i+1$ only. For this class of processes, we study the set $S$ of fixed points of the progeny generating function. In particular, we highlight the existence of a continuum of fixed points whose minimum is the global extinction probability vector $\boldsymbol{q}$ and whose maximum is the partial extinction probability vector $\boldsymbol{\tilde{q}}$. In the case where $\boldsymbol{\tilde{q}}=\boldsymbol{1}$, we derive a global extinction criterion which holds under second moment conditions, and when $\boldsymbol{\tilde{q}}<\boldsymbol{1}$ we develop necessary and sufficient conditions for $\boldsymbol{q}=\boldsymbol{\tilde{q}}$.

We consider a class of branching processes called Markovian binary trees, in which the individuals lifetime and reproduction epochs are modeled using a transient Markovian arrival process (TMAP). We estimate the parameters of the TMAP based on population data containing information on age-specific fertility and mortality rates. Depending on the degree of detail of the available data, a weighted non-linear regression method or a maximum likelihood method is applied. We discuss the optimal choice of the number of phases in the TMAP, and we provide confidence intervals for the model outputs. The results are illustrated using real data on human and bird populations.

Deriving the time-dependent expected reward function associated with a continuous-time Markov chain involves the computation of its transient deviation matrix. In this paper we focus on the special case of a finite quasi-birth-and-death (QBD) process, motivated by the desire to compute the expected revenue lost in a MAP/PH/1/C queue. We use two different approaches in this context. The first is based on the solution of a finite system of matrix difference equations; it provides an expression for the blocks of the expected reward vector, the deviation matrix, and the mean first passage time matrix. The second approach, based on some results in the perturbation theory of Markov chains, leads to a recursive method to compute the full deviation matrix of a finite QBD process. We compare the two approaches using some numerical examples.

We consider the extinction events of Galton-Watson processes with countably infinitely many types. In particular, we construct truncated and augmented Galton-Watson processes with finite but increasing sets of types. A pathwise approach is then used to show that, under some sufficient conditions, the corresponding sequence of extinction probability vectors converges to the global extinction probability vector of the Galton-Watson processes with countably infinitely many types. This gives rise to a number of iterative methods for the computation of the global extinction probability vector., Comment: 35 pages, 11 figures; References 4 and 6 updated

We perform an analytical sensitivity analysis for a model of a continuous-time branching process evolving on a fixed network. This allows us to determine the relative importance of the model parameters to the growth of the population on the network. We then apply our results to the early stages of an influenza-like epidemic spreading among a set of cities connected by air routes in the United States. We also consider vaccination and analyze the sensitivity of the total size of the epidemic with respect to the fraction of vaccinated people. Our analysis shows that the epidemic growth is more sensitive with respect to transmission rates within cities than travel rates between cities. More generally, we highlight the fact that branching processes offer a powerful stochastic modeling tool with analytical formulas for sensitivity which are easy to use in practice., Comment: 17 pages (30 with SI), Journal of Complex Networks, Feb 2015

In this paper, we use a finite-state continuous-time Markov chain with one absorbing state to model an individual's lifetime. Under this model, the time of death follows a phase-type distribution, and the transient states of the Markov chain are known as phases. We then attempt to provide an answer to the simple question "What is the conditional age distribution of the individual, given its current phase"? We show that the answer depends on how we interpret the question, and in particular, on the phase observation scheme under consideration. We then apply our results to the computation of the age pyramid for the endangered Chatham Island black robin Petroica traversi during years of intensive conservation efforts in 1980-1989.

We consider multitype Markovian branching processes evolving in a Markovian random environment. To determine whether or not the branching process becomes extinct almost surely is akin to computing the maximal Lyapunov exponent of a sequence of random matrices, which is a notoriously difficult problem. We define dual processes and we construct bounds for the Lyapunov exponent. The bounds are obtained by adding or by removing information: to add information results in a lower bound, to remove information results in an upper bound and we show that to add more information gives smaller lower bounds. We give a few illustrative examples and we observe that the upper bound is generally more accurate than the lower bound.

We consider queueing output processes of some elementary queueing models such as the M/M/1/K queue and the M/G/1 queue. An important performance measure for these counting processes is their variance curve, indicating the variance of the number of served customers over a time interval. Recent work has revealed some non-trivial properties dealing with the asymptotic rate at which the variance curve grows. In this paper we add to the results by finding explicit expressions for the second order approximation of the variance curve, namely the y-intercept of the linear asymptote. For M/M/1/K queues our results are based on the Drazin inverse of the generator. It turns out that by viewing output processes as MAPs (Markovian Arrival Processes) and considering the Drazin inverse, one can obtain explicit expressions for the y-intercept, together with some further insight regarding the BRAVO effect (Balancing Reduces Asymptotic Variance of Outputs). For M/G/1 queues our results are based on a classic transform of D.J. Daley. In this case we represent the y-intercept of the variance curve in terms of the first three moments of the service time distribution. A further performance measure that we are able to calculate for both models, is the asymptotic covariance between the queue length and the number of arrivals or departures. In addition we shed light on a classic conjecture of Daley, dealing with characterization of stationary M/M/1 queues within the class of stationary M/G/1 queues, based on the variance curve.

We present two iterative methods for computing the global and partial extinction probability vectors for Galton-Watson processes with countably infinitely many types. The probabilistic interpretation of these methods involves truncated Galton-Watson processes with finite sets of types and modified progeny generating functions. In addition, we discuss the connection of the convergence norm of the mean progeny matrix with extinction criteria. Finally, we give a sufficient condition for a population to become extinct almost surely even though its population size explodes on the average, which is impossible in a branching process with finitely many types. We conclude with some numerical illustrations for our algorithmic methods.

Matrix-form Poisson probability distributions were recently introduced as one matrix generalization of Panjer distributions. We show in this paper that under the constraint that their representation is to be nonnegative, they have a physical interpretation as extensions of PH distributions, and we name this restricted family Phase-type Poisson. We use our physical interpretation to construct an EM algorithm-based estimation procedure.