Multiscale Modelling of Birth-Death Processes

TL;DR

This work tackles the challenge of multiscale stochastic–deterministic modelling in biology by developing a principled method to choose the Jump-Switch-Flow (JSF) threshold $\Omega$ that bounds extinction-probability error $P_{ext}$. By reducing near-extinction dynamics to a time-inhomogeneous birth–death process and deriving backward Riccati equations for extinction, the authors decompose error into early and late excursions and provide both a rigorous bound and a practical late-term heuristic $\Delta_\Omega(T) \lesssim [q_Y(s_*)]^{\Omega+1}$. They validate the approach on a stochastic Lotka–Volterra model, showing the heuristic reliably upper-bounds empirical errors and enabling selection of the smallest $\Omega$ to meet a target tolerance. The framework offers a principled, efficient pathway for principled multiscale simulation in stochastic biology, with potential extensions to adaptive thresholds and other rare-event observables.

Abstract

Many biological systems exhibit multiscale dynamics, where some species occur in high copy numbers while others remain rare. This heterogeneity necessitates hybrid modelling approaches: deterministic models are computationally efficient but inaccurate for low-count species, while fully stochastic simulations are accurate but prohibitively expensive. Hybrid methods like the Jump-Switch-Flow (JSF) algorithm address this by simulating low-count species stochastically and high-count species deterministically. However, selecting regime-switching thresholds to control errors for specific observables remains an open challenge. We develop a principled framework for threshold selection targeting extinction probability. We formalise JSF as a piecewise-deterministic Markov process and derive backward equations for extinction under exact and hybrid dynamics. Near extinction boundaries, complex nonlinear dynamics reduce to tractable time-inhomogeneous linear birth-death processes. This structure yields a rigorous error decomposition based on early and late excursions. Isolating the dominant error term motivates a fast, actionable heuristic. We demonstrate via Monte Carlo studies on a stochastic Lotka-Volterra model that our heuristic reliably upper-bounds empirical errors in extinction probability. This enables users to select the smallest threshold that satisfies a target error tolerance. This work paves the way for principled, efficient multiscale modelling and simulation in stochastic biological systems.

Figures (5)

• Figure 1: Illustration of the critical time $t_c$ and point of no return $s_*$. (a) The net growth rate $r(t) = \lambda(t) - \mu(t)$ transitions from negative (death-dominated) to positive (birth-dominated) at $t_c$. The point of no return $s_*$ marks the latest time an up-crossing can occur and still be expected to return below $\Omega$ before $t_c$. (b) Illustrative trajectories: early excursions (green) are corrected by negative drift, while late excursions (red) persist due to favourable growth conditions.
• Figure 2: Comparison of full stochastic and hybrid simulations for the predator–prey model with decoupled predator decay. Parameters: $\alpha = 1.10$, $\beta = 0.05$, $\gamma = 0.4$, initial conditions $(x_{1,0}, x_{2,0}) = (10, 50)$. Top: 1000 exact SSA trajectories (solid black: mean). Bottom: 1000 JSF trajectories.
• Figure 3: As Fig. \ref{['fig:example_system_1']}, with increased prey birth rate $\alpha=1.60$. The higher birth rate lowers extinction probability.
• Figure 4: Validation of the practical heuristic bound for the Lotka-Volterra system. Blue circles show empirical error $\Delta_\Omega(T)$ between exact SSA and JSF hybrid simulations (2000 replications each), with Wilson score confidence intervals. Red dashed line shows the late-only heuristic bound $[q_Y(s_*)]^{\Omega+1}$. The simplified bound provides a reliable upper envelope, demonstrating the effectiveness of the practical approximation. The small uptick in the empirical $\Delta$ at $\Omega=35$ is a statistical fluctuation. Parameters: $\alpha = 1.10$, $\beta = 0.04$, $\gamma = 0.1$, $(x_{1,0}, x_{2,0}) = (10, 50)$, $T = 7$.
• Figure 5: Parameter space exploration for $\Omega=10$ across prey birth rate $\alpha$ and predation rate $\beta$. (a) Empirical extinction probability gap $\Delta_{\text{emp}}$ from 2000 simulations per grid point. (b) Deterministic theoretical bound values showing monotonic dependence on $\alpha$ but near-invariance to $\beta$. (c) Margin (bound slack) with green indicating positive margins; the black dashed contour at margin=0 indicates the bound-violation boundary. (d) Statistical significance (z-score = margin/error) with green indicating the bound holds and red indicating violations. Of the 36 parameter combinations, 35 show positive margins; the single marginal violation (z-score $= -0.10\sigma$) is not statistically significant, validating the theoretical framework across diverse dynamical regimes from high-extinction to high-survival scenarios.