Fact-checked by Grok 2 weeks ago

Branching process

A branching process is a mathematical model in probability theory that describes the evolution of a population over discrete generations, where each individual independently produces a random number of offspring according to a fixed probability distribution, and then dies, with the process continuing indefinitely from the offspring.^[1] The simplest form, known as the Galton–Watson process, begins with a single individual (Z_0 = 1) and tracks the population size Z_n at generation n, where the offspring counts are independent and identically distributed.^[2] The concept originated in the mid-19th century with Irénée-Jules Bienaymé's 1845 work on family name extinction probabilities, though it gained prominence through Francis Galton's 1873 query on surname persistence in British peerages, solved collaboratively with Henry William Watson in 1874.^[3] Early developments included partial results on extinction, with the first complete proof of the criticality theorem—distinguishing subcritical (mean offspring μ < 1, certain extinction), critical (μ = 1, certain extinction), and supercritical (μ > 1, positive survival probability)—provided by J.F. Steffensen in 1930.^[3] Key analytical tools include probability generating functions, where the extinction probability η satisfies η = f(η) (with f the offspring generating function) and is the smallest non-negative solution to this equation, equaling 1 if μ ≤ 1 and less than 1 otherwise in the supercritical case.^[1] The expected population size grows as E[Z_n] = μ^n, highlighting exponential growth potential when μ > 1.^[2] Branching processes have broad applications across disciplines, modeling phenomena involving reproduction and proliferation. In biology and genetics, they simulate allele frequencies, mutation propagation, and cell division dynamics.^[2] In epidemiology, they approximate early-stage disease outbreaks, such as the spread of infectious agents where each case generates a random number of secondary infections, aiding in estimating the basic reproduction number R_0 and outbreak extinction risks.^[4] Nuclear physics employs them for neutron chain reactions in fission, determining critical mass thresholds for sustained reactions.^[1] Extensions include multi-type processes for heterogeneous populations (e.g., by age or type) and continuous-time variants for age-dependent reproduction, influencing fields like ecology, queueing theory, and computer science algorithms.^[5]

Introduction

Definition and Motivation

A branching process is a stochastic model that describes the evolution of a population over discrete generations, where each individual independently produces a random number of offspring according to a fixed probability distribution, and the process is structured as a Markov chain with the state space consisting of non-negative integers representing population sizes.^[6] The population size at generation n, denoted Z_n, is the sum of the offspring produced by each individual in the previous generation Z_{n-1}, with the process typically initialized at Z_0 = 1 for a single founding individual.^[7] This formulation assumes non-overlapping generations, meaning individuals reproduce once and then die, and that offspring production is independent across individuals.^[8] The motivation for branching processes stems from their ability to capture the random, independent nature of reproduction in systems where lineage continuation or proliferation is of interest, providing a framework to analyze long-term population behaviors such as growth, stability, or decline.^[6] Originating in the context of modeling the persistence of family surnames, as posed by Francis Galton in 1873, the model addresses questions like the probability that a lineage eventually dies out.^[8] In biology, branching processes simulate population dynamics, including the spread of genetic traits or species evolution under random reproduction.^[7] They also apply to nuclear chain reactions, where neutrons independently generate further neutrons analogous to offspring, helping predict reaction sustainability.^[9] Additionally, in cell biology, the model represents proliferation in contexts like cancerous growths, estimating the likelihood of a cell colony extinguishing before overtaking tissue.^[10] A key outcome in branching processes, particularly when the expected number of offspring exceeds one, is the potential for ultimate extinction despite initial growth, underscoring their utility in risk assessment for replicating systems.^[6]

Historical Development

The origins of branching process theory trace back to the mid-19th century, when French statistician Irénée-Jules Bienaymé first mathematically analyzed the phenomenon of family extinction based on varying family sizes, providing an early probabilistic framework for population persistence.^[3] Bienaymé's work, published in 1845, anticipated key results on extinction probabilities but remained relatively obscure until later rediscoveries.^[11] In 1873, British polymath Francis Galton posed a famous query in the Educational Times regarding the extinction of surnames among the aristocracy, motivated by concerns over the persistence of family lineages.^[3] This problem, solved in collaboration with mathematician H. W. Watson the following year, introduced the discrete-time model now known as the Galton-Watson process, focusing on the probability that a lineage dies out over generations.^[3] Their 1874 paper provided an initial solution using probability generating functions to address extinction, laying the groundwork for subsequent developments.^[11] Early 20th-century advances built on these foundations, particularly in genetics. Ronald A. Fisher applied branching models in 1930 to study the survival of genetic traits and progeny lines, adapting the framework to biological inheritance without directly referencing Galton and Watson.^[3] Danish statistician J. F. Steffensen delivered the first complete proof of the criticality theorem for extinction probabilities in 1930, resolving lingering gaps in the theory.^[3] The theory experienced a significant resurgence in the 1940s and 1950s, driven by applications in nuclear physics, where physicist Leo Szilard had independently reinvented branching models in the late 1930s to describe neutron multiplication in fission chain reactions.^[1] This period saw intensified research, exemplified by William Feller's comprehensive treatment in his 1950 book An Introduction to Probability Theory and Its Applications, which formalized branching processes within modern probability theory and highlighted the utility of generating functions in extinction analysis.^[11] Studies like the 1953 RAND Corporation report further extended the models to neutron transport, marking branching processes' transition from demographic curiosities to tools in physical sciences.^[12]

Discrete-Time Branching Processes

Galton-Watson Process Formulation

The Galton-Watson process is the foundational model for discrete-time, single-type branching processes, originally motivated by the problem of surname extinction in families.^[13] It describes the evolution of a population where each individual independently produces a random number of offspring and then dies, with generations succeeding non-overlappingly. The process is named after Francis Galton and Henry Watson, who introduced the concept in their 1875 analysis of family lineage persistence.^[13] Formally, let Z_n denote the population size at generation n, for n = 0, 1, 2, \dots. The process begins with Z_0 = 1 individual, though the formulation extends naturally to arbitrary non-negative integer initial conditions Z_0. Each individual in generation n independently produces a random number of offspring X_i, where the X_i are identically distributed non-negative integer-valued random variables. The next generation size is then given recursively by

Z_{n+1} = \sum_{i=1}^{Z_n} X_i,

with the understanding that the sum is zero if Z_n = 0. The offspring production across individuals and generations is assumed to be independent. The distribution of the offspring random variable X_i is characterized by its probability mass function p_k = P(X_i = k) for k = 0, 1, 2, \dots, where \sum_{k=0}^\infty p_k = 1 and p_0 > 0 to allow for extinction. Common examples include the geometric distribution, where p_k = (1-p)^k p for $0 < p < 1, modeling cases with constant probability of continued reproduction; the Poisson distribution, p_k = e^{-\lambda} \lambda^k / k! for \lambda > 0, often used for rare events like mutations; and binary fission, where p_0 = p_2 = 1/2, representing simple splitting into zero or two offspring. Key moments of the offspring distribution include the mean m = E[X_i] = \sum_{k=0}^\infty k p_k, which determines the average growth rate per generation. For analytical tractability in studying properties like convergence, the variance \sigma^2 = E[X_i^2] - m^2 is often assumed finite. Probability generating functions provide a powerful tool for analyzing the distribution of Z_n.

Probability Generating Functions

In the Galton-Watson branching process, the offspring distribution \{p_k\}_{k=0}^\infty, where p_k = \mathbb{P}(X = k) and X denotes the number of offspring produced by a single individual, is conveniently analyzed using its probability generating function (PGF), defined as

f(s) = \mathbb{E}[s^X] = \sum_{k=0}^\infty p_k s^k

for |s| \leq 1. This function encapsulates the entire distribution in a compact form, facilitating the computation of probabilities and moments through algebraic manipulations and series expansions.^[12]^[14] The PGF f(s) possesses several key properties that underpin its utility. Normalization of the probabilities ensures f(1) = 1, while the mean offspring number m = \mathbb{E}[X] equals the first derivative evaluated at the boundary, m = f'(1). The second derivative f''(s) = \sum_{k=2}^\infty k(k-1) p_k s^{k-2} is nonnegative for s \in [0,1], implying that f(s) is convex on this interval. Moments beyond the mean can also be extracted; for instance, the variance of the offspring distribution is given by

\operatorname{Var}(X) = f''(1) + m - m^2.

These properties enable efficient derivation of distributional characteristics without direct summation over the probabilities.^[12]^[14] To study the population size Z_n at generation n, starting from Z_0 = 1, the PGF f_n(s) = \mathbb{E}[s^{Z_n}] is obtained iteratively through functional composition: f_n(s) = f(f_{n-1}(s)), with initial conditions f_0(s) = s and f_1(s) = f(s). This recursive structure reflects the independent branching of each individual in the previous generation, allowing the distribution of Z_n to be computed by repeated application of f, which is particularly advantageous for analytical or numerical evaluation of tail probabilities and higher moments. As a motivational example, the limiting extinction probability \eta = \lim_{n \to \infty} \mathbb{P}(Z_n = 0) satisfies the functional equation \eta = f(\eta), solvable via fixed-point iteration of the PGF starting from s = 0.^[12]^[14] The PGF framework extends naturally to growth analysis; for example, differentiating f_n(s) and evaluating at s = 1 yields the mean population size \mathbb{E}[Z_n] = m^n, highlighting exponential growth when m > 1.^[14]

Fundamental Properties

Mean Population Size and Growth

The mean population size in a Galton-Watson branching process is governed by the mean number of offspring per individual, denoted m = f'(1), where f(s) is the probability generating function of the offspring distribution. The conditional expectation satisfies \mathbb{E}[Z_{n+1} \mid Z_n] = m Z_n. Taking expectations yields the recursion \mathbb{E}[Z_n] = m \mathbb{E}[Z_{n-1}], which solves to \mathbb{E}[Z_n] = m^n assuming an initial population of Z_0 = 1.^[1] The asymptotic behavior of the expected population size classifies the process into three regimes based on m. In the subcritical case (m < 1), \mathbb{E}[Z_n] = m^n decays exponentially to 0, reflecting inevitable population decline on average.^[15] In the critical case (m = 1), \mathbb{E}[Z_n] = 1 remains constant for all n, but the survival probability \mathbb{P}(Z_n > 0) \to 0 as n \to \infty, highlighting that the stable mean masks the stochastic tendency toward extinction.^[15] In the supercritical case (m > 1), \mathbb{E}[Z_n] = m^n grows exponentially to infinity, indicating sustained population expansion on average.^[1] In the supercritical regime, assuming the offspring distribution has finite variance \sigma^2 = f''(1) + m - m^2 < \infty, the population size exhibits exponential growth with fluctuations captured by the variance. The exact variance is given by

\operatorname{Var}(Z_n) = \sigma^2 m^{n-1} \frac{m^n - 1}{m - 1}.

For large n, this asymptotes to \frac{\sigma^2}{m(m-1)} m^{2n}. Equivalently, the normalized process satisfies \operatorname{Var}(Z_n / m^n) \to \frac{\sigma^2}{m(m-1)}, corresponding to the variance of the almost sure limit W of Z_n / m^n, where \mathbb{E}[W] = 1. Under conditioning on non-extinction, the normalized size converges in distribution to the law of W given W > 0, preserving the scaling of growth and variability relative to the mean.^[16]^[15]

Extinction Probability and Criticality

The extinction probability \eta of a Galton-Watson branching process, starting from a single individual, is defined as \eta = \lim_{n \to \infty} P(Z_n = 0), where Z_n denotes the population size at generation n. This probability equals the smallest non-negative solution to the fixed-point equation \eta = f(\eta), where f(s) is the probability generating function of the offspring distribution. If the mean number of offspring m \leq 1, then \eta = 1; if m > 1 and p_0 + p_1 < 1 (where p_0 and p_1 are the probabilities of zero and one offspring, respectively), then $0 \leq \eta < 1.^[14] Branching processes are classified into three regimes based on the value of m: subcritical (m < 1), in which extinction occurs with probability 1 almost surely; critical (m = 1), in which extinction also occurs with probability 1 unless the process is trivial (p_1 = 1); and supercritical (m > 1), in which the extinction probability \eta < 1 and the survival probability $1 - \eta > 0. In the supercritical case, conditioned on non-extinction, the population grows exponentially with rate m. These classifications determine the long-term fate of the population, with certain extinction in subcritical and critical regimes due to the drift toward zero population size.^[14]^[12] To illustrate, consider a simple binary offspring distribution where each individual produces zero offspring with probability p_0 = 1/4 or two offspring with probability p_2 = 3/4, yielding mean m = 1.5 > 1. The pgf is f(s) = 1/4 + (3/4)s^2, and solving \eta = 1/4 + (3/4)\eta^2 gives the quadratic equation \eta^2 - (4/3)\eta + 1/3 = 0, with roots $1 and $1/3; thus, the extinction probability is the smaller root \eta = 1/3. For a geometric offspring distribution with p_k = (1-p) p^k for k = 0,1,2,\dots and mean m = p/(1-p), the pgf is f(s) = (1-p)/(1 - p s); if m > 1 (i.e., p > 1/2), then \eta = (1-p)/p < 1. These examples highlight how explicit solutions to \eta = f(\eta) reveal survival chances in supercritical settings.^[14] Near criticality, additional asymptotics describe the decay of survival probability. In the critical case (m=1) with finite positive variance \sigma^2 = f''(1) < \infty, the probability of survival up to generation n satisfies P(Z_n > 0) \sim 2/(\sigma^2 n) as n \to \infty. This slow polynomial decay underscores the precarious balance at criticality, where extinction is inevitable but the process lingers with small positive probability for large n.^[14]

Continuous-Time Branching Processes

Poissonian Offspring Formulation

In the Poissonian offspring formulation of continuous-time branching processes, each individual in the population lives for an exponentially distributed lifetime with rate parameter \lambda > 0, after which it produces a random number k of offspring according to a probability mass function p_k for k = 0, 1, 2, \dots, and then dies.^[14] This setup assumes that the offspring distribution has finite mean m = \sum_k k p_k and that reproduction events occur independently for each individual, modeling scenarios such as nuclear chain reactions or biological populations where generations overlap in time.^[12] The process begins with Z(0) = 1 individual, and Z(t) denotes the population size at time t \geq 0, evolving without discrete generational boundaries. The dynamics are captured through the probability generating function (PGF) g(t, s) = \mathbb{E}[s^{Z(t)} \mid Z(0) = 1] for s \in [0, 1], where the offspring PGF is f(s) = \sum_{k=0}^\infty p_k s^k. This PGF satisfies the backward Kolmogorov differential equation

\frac{\partial g}{\partial t}(t, s) = \lambda \left( f(s) - g(t, s) \right),

with initial condition g(0, s) = s.^[14] The solution to this equation provides the distribution of Z(t), and for a starting population of n individuals, the PGF is [g(t, s)]^n due to the branching property. This formulation generalizes the discrete-time Galton-Watson process to continuous time via a limit as the time step approaches zero.^[12] This Poissonian model is a special case of the more general Bellman-Harris process, where lifetimes follow an arbitrary distribution; here, the exponential lifetime ensures the Markov property, making the process memoryless and analyzable as a continuous-time Markov chain.^[14] Specifically, Z(t) is a Markov chain on the non-negative integers with state 0 absorbing, where from state n \geq 1, the transition rate to state n + k - 1 is n \lambda p_k for each k \geq 0, reflecting that any of the n individuals reproduces at rate \lambda and is replaced by k offspring.^[12] The infinitesimal generator of the chain encodes these rates, facilitating derivations of extinction probabilities and growth behaviors when m > 1.^[14]

Yule Process and Birth-Death Interpretation

The Yule process, introduced by G. Udny Yule in 1925, serves as a foundational model for continuous-time branching processes emphasizing unconstrained population growth through pure births. In this process, each individual independently produces one offspring at an exponential rate \lambda > 0, with no deaths occurring, effectively modeling binary fission where the parent continues to exist alongside the new individual. Starting from a single ancestor, the population size Z(t) at time t follows a geometric distribution:

P(Z(t) = k \mid Z(0) = 1) = e^{-\lambda t} (1 - e^{-\lambda t})^{k-1}, \quad k = 1, 2, \dots

This distribution arises from the memoryless property of exponential inter-birth times, leading to an exponentially growing lineage tree. The expected population size grows deterministically as E[Z(t)] = e^{\lambda t}, reflecting the supercritical nature of the model where growth accelerates without bound almost surely.^[17]^[18] The Yule process exemplifies how general continuous-time branching processes align with linear birth-death frameworks, where demographic events are Poissonian and independent across individuals. In a linear birth-death process, the birth rate for a population of size k is \lambda_k = k \beta and the death rate is \mu_k = k \delta, with \beta, \delta \geq 0 denoting per capita rates; this structure captures branching by treating births as offspring production and deaths as lineage termination. The Yule process emerges as the special case \delta = 0 and \beta = \lambda, reducing to pure birth dynamics. Such interpretations facilitate analysis of population trajectories as Markov chains, with transitions driven by individual-level events.^[19]^[20] Key properties of linear birth-death processes derive from solving the Kolmogorov forward equations, yielding the mean population size m(t) = E[Z(t) \mid Z(0) = 1] = e^{(\beta - \delta)t}, which indicates subcritical decay (\beta < \delta), critical stasis (\beta = \delta), or supercritical explosion (\beta > \delta) based on the net growth rate \beta - \delta. For extinction probabilities, starting from one individual, the ultimate extinction probability \eta satisfies \eta = 1 if \beta \leq \delta, and \eta = \delta / \beta if \beta > \delta > 0; in the pure birth Yule case (\delta = 0), \eta = 0, ensuring persistence with probability 1. These results stem from the fixed-point equation of the probability generating function and highlight the process's alignment with discrete-time branching criticality, adapted to continuous parameters.^[19]^[20]

Multitype Branching Processes

Vector Formulation and Matrices

In multitype branching processes, populations consist of d distinct types of individuals, labeled $1 through d. The state of the population at discrete generation n is described by the nonnegative integer vector \mathbf{Z}_n = (Z_n^{(1)}, \dots, Z_n^{(d)})^\top, where Z_n^{(i)} denotes the number of individuals of type i. Each individual reproduces independently, and the process evolves across generations according to the offspring distributions specific to each type. This formulation generalizes the single-type Galton-Watson process, which corresponds to the special case d=1.^[14] The reproduction of an individual of type i is modeled by a random vector \mathbf{Y}_i = (Y_i^{(1)}, \dots, Y_i^{(d)})^\top, where Y_i^{(j)} is the number of offspring of type j produced by that parent. These offspring vectors for individuals of the same type are independent and identically distributed. The probability generating function (PGF) for the offspring distribution of type i is defined as

f_i(s_1, \dots, s_d) = \mathbb{E}\left[ \prod_{j=1}^d s_j^{Y_i^{(j)}} \right],

where \mathbf{s} = (s_1, \dots, s_d) with |s_j| \leq 1 for all j. The functions f_i satisfy f_i(1, \dots, 1) = 1 and are convex, analytic in the unit polydisc.^[14] The recursion for the population vector is given by

\mathbf{Z}_{n+1} = \sum_{i=1}^d \sum_{k=1}^{Z_n^{(i)}} \mathbf{Y}_{i,k},

where \mathbf{Y}_{i,k} are i.i.d. copies of \mathbf{Y}_i for each type i and individual k. This expresses \mathbf{Z}_{n+1} as the vector sum of the offspring contributions from all individuals in generation n. Conditional on \mathbf{Z}_n, the components of \mathbf{Z}_{n+1} are independent across types, with the offspring from type i individuals determining the increments to each type j.^[14] The expected behavior is captured by the mean matrix M = (m_{ij})_{i,j=1}^d, where m_{ij} = \mathbb{E}[Y_i^{(j)}] is the expected number of type j offspring from a type i parent. These entries are finite under standard moment assumptions and given by

m_{ij} = \left. \frac{\partial f_i(s_1, \dots, s_d)}{\partial s_j} \right|_{s_1 = \dots = s_d = 1}.

The conditional expectation satisfies \mathbb{E}[\mathbf{Z}_{n+1} \mid \mathbf{Z}_n] = M \mathbf{Z}_n, so the unconditional mean evolves as \mathbb{E}[\mathbf{Z}_n] = M^n \mathbb{E}[\mathbf{Z}_0]. The matrix M is nonnegative, and its spectral properties, such as the dominant eigenvalue, determine long-term growth rates.^[14] The probability generating function for the population vector at generation n, starting from an arbitrary \mathbf{Z}_0, is F_n(\mathbf{s}) = \mathbb{E}\left[ \prod_{j=1}^d s_j^{Z_n^{(j)}} \right]. For a single ancestor of type i, this reduces to the n-fold functional iterate F_n(\mathbf{s}) = f_i^{(n)}(\mathbf{s}). In vector notation, define the offspring PGF vector \mathbf{f}(\mathbf{s}) = (f_1(\mathbf{s}), \dots, f_d(\mathbf{s}))^\top. The iterates are obtained componentwise via \mathbf{s}^{(0)} = \mathbf{s} and \mathbf{s}^{(n)} = \mathbf{f}(\mathbf{s}^{(n-1)}) for n \geq 1, so F_n(\mathbf{s}) involves the n-th iterate \mathbf{s}^{(n)}. This composition structure facilitates analysis of distributional properties across generations.^[14] Key assumptions for the model include finite first moments (\mathbb{E}[|\mathbf{Y}_i|] < \infty for all i) to ensure the mean matrix is well-defined. Additionally, irreducibility of the types is typically required, meaning the directed graph with vertices as types and edges if P(Y_i^{(j)} > 0) > 0 for some i,j is strongly connected; equivalently, some power of M has all positive entries (positive regularity). These conditions ensure meaningful interactions among types and avoid trivial decompositions into independent subprocesses.^[14]

Extinction Probabilities and Perron-Frobenius Analysis

In multitype branching processes, the extinction probabilities are captured by the vector \eta = (\eta_1, \dots, \eta_d)^\top \in [0,1]^d, where \eta_j denotes the probability of eventual extinction of the entire population starting from a single individual of type j. This vector is the smallest (componentwise) non-negative solution to the fixed-point equations

\eta_j = f_j(\eta_1, \dots, \eta_d), \quad j = 1, \dots, d,

with f_j(s_1, \dots, s_d) being the probability generating function (PGF) for the offspring distribution of a type-j individual, defined as f_j(\mathbf{s}) = \sum_{\mathbf{k} \in \mathbb{N}_0^d} p_{j,\mathbf{k}} \prod_{i=1}^d s_i^{k_i}, where p_{j,\mathbf{k}} is the probability of producing \mathbf{k} = (k_1, \dots, k_d)^\top offspring of each type.^[14] The solution \eta = \mathbf{1} (the all-ones vector) always satisfies the equations, but the minimal solution determines the true extinction probabilities, obtained iteratively as the limit \eta^{(n)} \to \eta where \eta^{(0)} = \mathbf{0} and \eta^{(n+1)}_j = f_j(\eta^{(n)}). The process experiences global extinction with probability 1 if at least one component \eta_i = 1, which occurs in subcritical and critical regimes; otherwise, the probabilities of long-term survival starting from type i are $1 - \eta_i > 0. This formulation extends the single-type case (where d=1) by solving a system of multivariate equations rather than a scalar one.^[14] The mean behavior of the process is governed by the offspring mean matrix M = (m_{ij}), where m_{ij} is the expected number of type-j offspring produced by a type-i individual, so the j-th row of M is the gradient \nabla f_j(\mathbf{1}). The long-term population growth rate is dictated by the spectral radius \rho(M): the process is subcritical if \rho(M) < 1 (extinction probability \eta = \mathbf{1}), critical if \rho(M) = 1 (also \eta = \mathbf{1} under mild conditions), and supercritical if \rho(M) > 1 (where \eta < \mathbf{1} componentwise, provided there is positive probability of producing no offspring, i.e., f_j(\mathbf{0}) > 0 for some j).^[14] In the supercritical case, the expected population size grows asymptotically as \rho^n times a random limit proportional to the positive left and right Perron eigenvectors. When M is a positive irreducible nonnegative matrix, the Perron-Frobenius theorem ensures that \rho(M) is a simple positive real eigenvalue (the spectral radius), with strictly positive right eigenvector v > \mathbf{0} (normalized so v^\top \mathbf{1} = 1) and positive left eigenvector u > \mathbf{0} (with u^\top v = 1), dominating all other eigenvalues in modulus. The extinction vector satisfies \eta < \mathbf{1} if and only if \rho(M) > 1 and the boundary conditions hold (e.g., the probability of immediate extinction is positive for at least one type); otherwise, \eta = \mathbf{1}. This spectral analysis classifies the process's criticality and links extinction to the principal eigenvalue, providing a matrix-theoretic foundation for multivariate growth and survival.^[14] For illustration, consider a two-type (d=2) branching process modeling mutual production between types (e.g., two interacting species or particle types like electrons and photons in cascades), where type 1 has PGF f_1(s,t) and type 2 has f_2(s,t). The extinction vector (\eta_1, \eta_2)^\top solves the system \begin{align*} \eta_1 &= f_1(\eta_1, \eta_2), \ \eta_2 &= f_2(\eta_1, \eta_2), \end{align*} with \eta = \mathbf{1} always a solution, but the minimal solution is found iteratively or numerically if the PGFs are not analytically invertible. The mean matrix M then determines if \eta < \mathbf{1} via \rho(M) > 1, as in models of interacting populations where cross-production ensures irreducibility.^[14]

Extensions and Variations

Size-Dependent Branching Processes

In size-dependent branching processes, also known as density-dependent branching processes, the offspring distribution of each individual depends on the current population size Z_n = n, allowing the model to incorporate regulatory mechanisms such as reduced fertility in crowded populations, which contrasts with the fixed offspring distribution in the standard Galton-Watson process. This dependence breaks the independence assumption across generations while maintaining that, conditional on the current size, individuals reproduce independently. Such models are particularly useful in ecology and population biology to simulate logistic-like growth where reproduction rates decline with density. The process is formally defined by the recursion

Z_{n+1} = \sum_{i=1}^{Z_n} X_i(n),

where the X_i(n) are independent and identically distributed random variables with probability mass function P(X_i(n) = k) = p_k(n) for k = 0, 1, 2, \dots, and the distribution \{p_k(n)\} is parameterized by the current population size n. The mean number of offspring per individual is \mu(n) = \sum_{k=0}^\infty k p_k(n), so the conditional expected population size is E[Z_{n+1} \mid Z_n = n] = n \mu(n). For many applications, the variance of X_i(n) is also considered, often assumed to stabilize or grow mildly with n to ensure well-behaved limits. Key properties revolve around the behavior of \mu(n) as n \to \infty. In subcritical regimes where \mu(n) \to m < 1, the population tends to extinction almost surely. When \mu(n) \to 1, the process may exhibit critical-like behavior, with extinction probability 1 if \mu(n) approaches 1 from above sufficiently fast (e.g., \mu(n) - 1 \sim n^{- \gamma} with \gamma > 1); conversely, if \mu(n) > 1 and approaches 1 from above at a suitable rate (e.g., \mu(n) - 1 \sim n^{-\alpha} with $0 < \alpha < 1), the normalized population Z_n / \prod_{k=1}^{n-1} \mu(k) converges almost surely to a positive random limit. For stability around an equilibrium, models often impose \mu(1) = 1 and \mu'(1) < 0 when density is normalized by a carrying capacity K, ensuring the population lingers near K in large-K limits via convergence to an Ornstein-Uhlenbeck process for fluctuations. A representative example is the logistic branching process, a model that combines branching reproduction with density-dependent mortality (e.g., quadratic death rates), leading to regulated growth and convergence to an Ornstein-Uhlenbeck process in large populations.^[21] Another instance is the Cannings model, which enforces exchangeable offspring partitions dependent on a fixed population size N, effectively introducing size-dependent reproduction through the constraint that total offspring sums to N, as studied in population genetics for analyzing genetic drift. Branching processes with environmental stochasticity tied to size further extend this by letting the distribution \{p_k(n)\} incorporate random environmental factors scaling with n, such as variable resource availability.

Branching Processes with Migration or Competition

Branching processes with migration incorporate external arrivals of individuals, modifying the standard Galton-Watson dynamics by adding an independent immigration component at each generation. In this extension, the population size at generation n+1, denoted Z_{n+1}, is given by the sum of the offspring from the current population Z_n and an independent immigration random variable B_{n+1}, where the B_i are i.i.d. with some distribution. This model captures scenarios such as population growth influenced by ongoing recruitment from outside the reproducing group, common in ecological or demographic studies.^[22] The mean population size in such processes satisfies a modified recurrence relation. Let m_n = \mathbb{E}[Z_n] denote the expected size at generation n, \mu the mean offspring number per individual, and \nu = \mathbb{E}[B_{n+1}] the mean immigration rate (assumed constant for simplicity). Then, m_{n+1} = \mu m_n + \nu. For subcritical cases (\mu < 1), this linear equation implies convergence to a stationary mean m_\infty = \nu / (1 - \mu), ensuring long-term persistence despite individual extinction risk.^[23] Equilibrium analysis further reveals a stationary distribution for the process, often compound-geometric in form, which describes the limiting population size under immigration. Competition introduces density-dependent effects that limit growth, reflecting resource constraints or intra- and inter-individual rivalry. A common formulation adjusts fertility or survival probabilities based on current population size, yielding a logistic-like regulation. This density-dependent branching process stabilizes population fluctuations around an equilibrium, contrasting pure branching models that grow unboundedly or extinct.^[24] In multitype settings, competition extends to interactions across types, where offspring production rates for type i depend on the total population vector \mathbf{Z}_n. For example, logistic competition modifies the mean matrix to include terms penalizing total density, leading to bounded growth and coexistence analysis via eigenvalue perturbations of the base reproduction matrix.^[25] Properties include altered extinction probabilities, often lower than in non-competitive multitype processes due to density feedback, and stationary distributions supported on compact sets reflecting competitive equilibria.^[25] An illustrative example is the branching process with absorption, where population size is constrained by an absorbing barrier at zero (extinction state), but immigration prevents total collapse, or with reflecting barriers that redirect excess individuals. These boundary conditions alter recurrence relations for means and probabilities, modeling regulated environments like finite habitats.^[26]

Computational Methods

Simulation Techniques

Simulation of branching processes is essential for empirical investigation, particularly when analytical solutions are intractable or for validating theoretical predictions through Monte Carlo methods. In discrete-time branching processes, such as the Galton-Watson model, simulation begins by initializing the population size Z_0 = 1. For each subsequent generation n, the number of offspring produced by each of the Z_n individuals is independently sampled from the specified offspring probability distribution p_k, where k is the number of offspring. The population size for the next generation is then Z_{n+1} = \sum_{i=1}^{Z_n} X_i, with X_i denoting the offspring of the i-th individual; this process repeats until either extinction (Z_n = 0) or a predetermined maximum number of generations is reached to prevent infinite runs in supercritical cases. Sampling can employ direct methods for simple distributions or inverse cumulative distribution function techniques for more complex ones, enabling the generation of multiple independent realizations to estimate quantities like extinction probabilities or growth trajectories.^[27]^[28] For continuous-time branching processes, the Gillespie algorithm provides an exact stochastic simulation method by modeling the process as a continuous-time Markov chain with event-driven updates. Each individual waits an exponentially distributed time until reproduction or death, with rates determined by the offspring and survival distributions; the algorithm selects the next event time across all individuals using the minimum of their waiting times and updates the population accordingly upon occurrence. This direct method ensures statistical exactness for the underlying master equation without time discretization approximations, making it suitable for birth-death interpretations like the Yule process. Variations, such as the first-reaction or next-reaction methods, optimize efficiency for large populations by avoiding redundant rate recalculations.^[29] In multitype branching processes, simulation extends the single-type algorithms by tracking populations of distinct types, where each individual of type i produces a vector of offspring across all types according to a joint probability distribution or mean matrix M. Discrete-time versions sample these vectors independently for each individual, aggregating counts per type to form the next generation's composition; continuous-time implementations apply the Gillespie algorithm with type-specific event rates for transitions between types. This approach captures interactions, such as in ecological or epidemiological models with multiple species or states.^[30] Estimating extinction probabilities via simulation, especially in supercritical regimes where the probability is small and rare, introduces bias from finite run lengths or incomplete conditioning; correction techniques include importance sampling, which biases the offspring distribution toward fewer progeny to increase extinction events, followed by likelihood reweighting, or conditioning simulations on paths leading to extinction. These methods reduce variance and improve accuracy for rare-event probabilities compared to naive Monte Carlo.^[31]^[28] Software implementations facilitate these simulations: the R package 'branching' provides functions for single- and multi-type discrete-time processes, including extinction probability estimation, while Python users can leverage SciPy's random number generators for custom offspring sampling in both discrete and continuous settings. Simulations often visualize outcomes as family trees to illustrate lineage branching and extinction paths.

Numerical Solution of Extinction Equations

In single-type branching processes, the extinction probability \eta satisfies the fixed-point equation \eta = f(\eta), where f(s) is the probability generating function of the offspring distribution. When a closed-form solution is unavailable, numerical computation often begins with the iterative scheme starting from \eta_0 = 0 and updating \eta_{k+1} = f(\eta_k); this monotone increasing sequence converges to \eta from below provided the mean offspring number m \leq 1 or m > 1 with \eta < 1. To accelerate convergence, particularly when the derivative f'(s) is known or computable, Newton's method can be applied: \eta_{k+1} = \eta_k - \frac{f(\eta_k) - \eta_k}{f'(\eta_k)}, which exhibits quadratic convergence near the root under standard smoothness conditions on f. For multitype branching processes, the extinction probability vector \boldsymbol{\eta} solves \boldsymbol{\eta} = \mathbf{f}(\boldsymbol{\eta}), where \mathbf{f} is the vector of probability generating functions. Componentwise fixed-point iteration \boldsymbol{\eta}_{k+1} = \mathbf{f}(\boldsymbol{\eta}_k), initialized at the zero vector, converges monotonically to \boldsymbol{\eta} in the subcritical case (spectral radius of the mean matrix m < 1), leveraging the contraction mapping theorem on the unit cube. In more general settings, including supercritical regimes with finite types, polynomial-time fixed-point solvers such as value iteration or specialized algorithms exploiting the nonnegative matrix structure ensure efficient computation. When the process is near criticality (m \approx 1), direct iteration may converge slowly, prompting approximation methods for the ultimate extinction probability or finite-time survival P(Z_n > 0). Saddlepoint approximations provide asymptotic expansions for tail probabilities, integrating the cumulant generating function to yield accurate estimates for large n. Diffusion approximations model the rescaled population size as a Feller process, capturing the quasi-stationary behavior and yielding P(Z_n > 0) \approx 2 / (\sigma^2 n) for critical cases with variance \sigma^2 > 0, where the probability generating function admits the series expansion f(s) = s + \frac{\sigma^2 (1-s)^2}{2} + o((1-s)^2) near s=1. In continuous-time branching processes, the probability generating function g(t,s) evolves according to the backward Kolmogorov equation \frac{\partial g}{\partial t} = Q(g), where Q is the infinitesimal generator derived from the reproduction rates; the ultimate extinction probability is \lim_{t \to \infty} g(t,0). Numerical solutions employ integration schemes like Runge-Kutta methods to solve this partial differential equation over a discretized grid, allowing evaluation of the long-time limit. Convergence analysis for these methods reveals linear rates for fixed-point iteration, with the constant governed by \sup |f'(\eta)| in single-type cases or the spectral radius in multitype settings, while Newton's method achieves quadratic convergence, reducing error by squaring it at each step near the solution. These numerical approaches can be validated against Monte Carlo simulations from the underlying stochastic process for empirical verification.

References

[1]
[PDF] BRANCHING PROCESSES Galton-Watson processes were ...
Galton-Watson processes were introduced by Francis Galton in 1889 as a simple mathemat- ical model for the propagation of family names. They were reinvented by ...
[2]
[PDF] Branching Processes and Their Applications - Lancaster University
Branching processes are mathematical representations of population growth over generations, where members reproduce a random number of offspring.
[3]
[PDF] THREE PAPERS ON THE HISTORY OF BRANCHING PROCESSES
ABSTRACf. The first complete proof of the criticality theorem for Bienayme-Galton-Watson branching processes was published in Danish in 1930 by J. F ...
[4]
Branching Processes: Their Role in Epidemiology - PubMed Central
Branching processes are stochastic individual-based processes leading consequently to a bottom-up approach.
[5]
[PDF] Branching processes
Mar 6, 2021 · Branching processes arise in stochastic processes on trees and graphs, like Galton-Watson processes, which model population growth.
[6]
[PDF] Essentials of Stochastic Processes - Duke Mathematics Department
Consider a branching process as defined in Example 7.2, in which each ... Lawler, G.F. (2006) Introduction to Stochastic Processes. Second edition ...
[7]
[PDF] Stochastic Processes Amir Dembo (revised by Kevin Ross) April 12 ...
Apr 12, 2021 · . Stochastic processes are typically defined only up to a set of zero probability, and ... Definition 4.6.1 (Branching Process). The ...
[8]
[PDF] COURSE NOTES STATS 325 Stochastic Processes
• formal definition of stochastic processes. 1.1 Revision: Sample spaces and ... Definition: A branching process is defined as follows. • Single ...
[9]
[PDF] Coalescence in Bellman-Harris and multi-type branching processes
because of the analogy between the growth of families and nuclear chain reactions and also partly because of the increased general interest in applications of ...Missing: names proliferation
[10]
[PDF] MULTI-TYPE BRANCHING PROCESSES WITH TIME-DEPENDENT ...
May 27, 2018 · Branching processes are central in the study of evolution of various populations such as bacteria, cancer cells, carriers of a particular form ...Missing: proliferation | Show results with:proliferation
[11]
(PDF) Branching Processes: A Personal Historical Perspective
It gives a—personally biased—sketch of the development of branching processes, from the mid nineteenth century to 2010, emphasizing relations to bioscience and ...
[12]
[PDF] The Theory of Branching Process - RAND
A branching process is a mathematical representation of a population's development where members reproduce and die, subject to chance.
[13]
None
Nothing is retrieved...<|control11|><|separator|>
[14]
Branching Processes | SpringerLink
The purpose of this book is to give a unified treatment of the limit theory of branching processes ... Ney. Krishna B. Athreya. University of Wisconsin ...
[15]
[PDF] Chapter 6: Branching Processes: The Theory of Reproduction
A branching process starts with one individual, each living one unit, then producing Y offspring (0, 1, 2,...) and dying. The number of offspring is random.Missing: mathematics | Show results with:mathematics
[16]
[PDF] Branching Processes I: Supercritical growth and population structure
is a martingale and converges to a random variable W a.s. as n → ∞. The characterization of the limit W in the supercritical case, m > 1, under minimal.
[17]
Reading Yule in light of the history and present of macroevolution
Feb 20, 2025 · Yule's 1925 paper introducing the branching model that bears his name was a landmark contribution to the biodiversity sciences.
[18]
23. Continuous-Time Branching Chains - Random Services
The Yule Process. Next we consider the pure birth branching chain in which each particle, at the end of its life, is replaced by 2 new particles. Equivalently ...
[19]
[PDF] THE LINEAR BIRTH–DEATH PROCESS - Columbia University
The linear birth-death process models population growth/extinction, where individuals birth at rate λ and die at rate μ, independently.
[20]
On the Generalized "Birth-and-Death" Process - Project Euclid
In this paper, I shall give the complete solution of the equations governing the generalised birth-and-death process in which the birth and death rates λ(t) λ ...
[21]
https://doi.org/10.1214/105051605000000098
[22]
https://www.jstor.org/stable/3211835
[23]
Branching Processes with Immigration - jstor
O < ao, ao + a1, bo <1. The first two conditions allow us to use the results that have been developed for the ordinary branching process (b.p.), see ...
[24]
Estimation of the Means in the Branching Process with Immigration
Let {Xn} { X n } be the branching process with immigration and let m m and λ λ be the means of the offspring and immigration distributions, respectively.
[25]
The branching process with logistic growth - Project Euclid
In order to model random density-dependence in population dynamics, we construct the random analogue of the well-known logistic process in the branching ...
[26]
On multitype Branching Processes with Logistic Competition
Aug 27, 2025 · In this paper we take a completely different approach and consider multi-type continuous-state branching process, now allowing for up to a ...
[27]
A condition for the extinction of a branching process with ... - PubMed
A branching process with an absorbing lower barrier is considered. This is a Galton-Watson process with the condition that at any generation the number of ...Missing: reflecting | Show results with:reflecting
[28]
Discrete-time branching processes and their applications to ...
Aug 2, 2024 · Applications of branching processes include the model of cell reproduction, the propagation of neutrons in a nuclear reactor, or the spread of ...
[29]
Branching Processes - Cambridge University Press & Assessment
Branching Processes: Variation, Growth, and Extinction of Populations. Search within full text. Access. Patsy Haccou, Rijksuniversiteit Leiden, The Netherlands.
[30]
[PDF] The Gillespie Algorithm: Simulating Stochasticity
The Gillespie algorithm is a generic useful algorithm that's great for simulating Markov Process where events happen at exponential rate.Missing: continuous- | Show results with:continuous-
[31]
Multitype branching process method for modeling complex ...
Mar 17, 2022 · The MTBP formulation of the dynamics allowed us to define the mean matrix, M , which describes the average behavior of the system. Using M , a ...
[32]
Unbiased Simulation of Rare Events in Continuous Time
Nov 2, 2021 · This provides unbiased estimations of the probability in question. We discuss the practical feasibility of the algorithm with reference to ...