Quasi-Stationary Approximation Methods

• Quasi-stationary approximation is a method to compute the QSD of stochastic processes with absorbing states, modeling metastable equilibria conditioned on non-absorption.
• It encompasses techniques such as Fleming–Viot particle systems, stochastic approximation, and returned-process methods that provide rigorous error bounds and convergence analysis.
• These approximations have practical applications in population dynamics, epidemiology, network science, and statistical physics, offering insights into metastability and extinction phenomena.

A quasi-stationary approximation refers to any principled methodology for approximating or numerically evaluating the quasi-stationary distribution (QSD) of a Markovian or stochastic process with one or more absorbing states. These distributions encapsulate the limiting behavior of the process conditioned on non-absorption, rigorously describing the apparent long-term “metastable” equilibrium before eventual extinction or absorption. Quasi-stationary approximation has become central in population dynamics, statistical physics, epidemiology, network science, and computational probability, where direct computation or simulation of the QSD is often computationally intractable or analytically inaccessible. The field encompasses spectral, stochastic, particle, and mean-field techniques, enabled by advances in large deviations theory, stochastic approximation, and interacting particle systems.

1. Definition and Characterization of Quasi-Stationary Distributions

Given a Markov process $Z(t)$ on a countable or general state space $\Lambda\cup\{0\}$ with $0$ absorbing and irreducibility on $\Lambda$, the quasi-stationary distribution $\nu$ on $\Lambda$ is defined as the unique law such that, for all $t\geq 0$ and $A\subset\Lambda$,

$$\mathbb{P}_\nu (Z(t) \in A \mid \tau>t) = \nu(A), \quad \tau = \inf\{ t : Z(t)=0 \}.$$

This is equivalent to invariance under the conditioned semigroup $T_t$: $T_t\nu = \nu \qquad \text{for all } t\geq 0,$ where

$$T_t\mu(x) := \frac{\sum_{z\in\Lambda} \mu(z) P_t(z,x)}{1 - \sum_{z\in\Lambda} \mu(z) P_t(z,0)}, \quad P_t = e^{tQ}.$$

Differentiating this yields the nonlinear fixed-point (nonlinear eigenvalue) equation: $$0 = \sum_{y \in \Lambda} q(y,x)\nu(y) + \left( \sum_{y \in \Lambda} q(y,0)\nu(y) \right) \nu(x),\qquad x\in\Lambda,$$ which can also be viewed as a left-eigenvector equation for the substochastic generator restricted to $\Lambda$ with normalization. Typically, the decay rate $\theta_\nu = -\sum_{y\in\Lambda} q(y,0)\nu(y)>0$ so that $\mathbb{P}_\nu(\tau>t) = e^{-\theta_\nu t}$ (Groisman et al., 2012).

Existence and uniqueness depend on recurrence and mixing properties. There may be no QSD, a unique QSD (e.g., under a uniform Doeblin condition), or (as for birth-death and branching processes) a continuum of QSDs corresponding to minimal, maximal, and intermediate solutions (Barbour et al., 2011, Barbour et al., 2010).

2. Principle Methods and Algorithms for Quasi-Stationary Approximation

Quasi-stationary approximation comprises a suite of methodologies, including particle system methods, stochastic approximation, spectral/returned-process approaches, and functional fixed-point schemes.

2.1 Fleming–Viot Particle Systems

The Fleming–Viot particle system approximates the conditioned process by evolving $N$ independent copies of the underlying process; any particle reaching the absorbing state is instantaneously resurrected at the position of a uniformly chosen survivor. The empirical measure of the particles approximates the law conditioned on non-absorption. On finite time horizons, the error is $O(N^{-1/2})$ in expectation and vanishes as $N\to\infty$ (Groisman et al., 2012). On the stationary time scale, for processes with a unique QSD, the stationary empirical law of the FV system converges in total variation to the QSD as $N\to\infty$. For processes with infinitely many QSDs (e.g. subcritical branching), the FV system selects the minimal QSD (Groisman et al., 2012, Villemonais, 2014).

Recent extensions handle general Markov processes via normalized Feynman-Kac particle systems, with explicit uniform-in-time bias, variance, and exponential concentration bounds in $L^p$ and probability (Journel et al., 20 Dec 2024, Journel et al., 2019). For diffusions or soft-killing models, discretized particle approximations yield rigorous error bounds in Wasserstein distance, separated into contributions from simulation time $t$, sample size $N$, and time-step $\gamma$ (Journel et al., 2019).

2.2 Stochastic Approximation and Self-Interacting Chains

Particle-free algorithms based on self-interacting or reinforced random walks iteratively reinforce occupation or transition measures, approximating the QSD by a stochastic approximation algorithm. For finite-state Markov processes, occupation or evolution measures are updated recursively from a non-linear kernel incorporating “resurrected” moves, yielding almost sure convergence (and a central limit theorem) of the empirical measure to the QSD at polynomial rates (Benaïm et al., 2014). These methods extend to general compact spaces and even diffusions via continuous-time reinforced processes (Benaim et al., 2016, Benaïm et al., 2019). The algorithmic structure is typically:

• At each killing/absorption, resample the process at a state distributed according to the measure of previously visited states.
• Update the occupation/transition measure with a fixed or decreasing step-size.

Non-failable variants, such as the recycling resampling scheme, achieve robust uniform-in-time error bounds and prevent breakdown due to repeated simultaneous absorptions (Oçafrain et al., 2016).

2.3 Fixed-Point and Regeneration-Type Approximations

Returned-process or regeneration-based approximations replace absorption with an immediate return to a typical or reference state or distribution (often centered at or near a metastable equilibrium). The stationary law of the returned process serves as a tractable proxy for the QSD. When explicit return-time bounds (mean and hitting probability) are satisfied, the total variation distance between the QSD and any returned-process stationary law is sharply bounded by explicit expressions in these parameters, often exponentially small in system size for large, metastable systems (Barbour et al., 2010, Barbour et al., 2011). Truncation to a finite “core” and return-to-core operators enable low-dimensional linear algebra solutions.

Self-interacting chains and McKean-Vlasov dynamics, conditioned on non-extinction or non-absorption, can be analyzed as fixed-points of measure-valued evolution equations (Assadeck et al., 26 Jun 2025). Recursion and backward error-analysis further enable systematic control over bias and error.

3. Approximation Schemes in Practice: Examples

The practical implementation of quasi-stationary approximations has been illustrated in broad application domains:

• SIS/epidemic models: Particle, ODE-based, and mean-field quasi-stationary approximations have been intensively studied for the logistic SIS model and network-based epidemics (Nåsell, 2022, Overton et al., 2022, Cantwell et al., 15 Sep 2025). The three-stage Ovaskainen approximation $q^{(OV3)}$ for the QSD in logistic SIS yields exponentially small uniform error in system size $N$ (i.e., $O(e^{-I_1(R_0)N})$ for $R_0>1$), outperforming the auxiliary process $p^{(0)}$ (Nåsell, 2022).
• Birth–death processes: Minimal QSDs and their approximation via Fleming–Viot and Lyapunov techniques control the domain of attraction for a variety of initial conditions (Villemonais, 2014, Barbour et al., 2010).
• Diffusions and McKean–Vlasov SDEs: Euler schemes, weighted occupation measures, or reinforced processes deliver consistent QSD approximations in bounded or compact sets, with convergence in the weak topology and uniform error bounds (Benaïm et al., 2019, Wang et al., 2018, Assadeck et al., 26 Jun 2025, Benaim et al., 2016).
• Chemical reaction networks and mass-action kinetics: Sensitivity and scaling of QSDs against diffusion approximations can be estimated via pathwise Poisson–diffusion couplings, giving $O(\log V/V)$ 1-Wasserstein error in system size $V$ (Li et al., 2022).
• Bayesian computation: Regeneration-based quasi-stationary Monte Carlo schemes sample from posteriors via target QSD simulation in killed diffusions (Wang et al., 2018).

4. Theoretical Properties: Existence, Uniqueness, and Error

Existence and uniqueness are controlled by irreducibility, recurrence, and mixing-type hypotheses:

• Sufficient conditions: uniform lower bounds on return probabilities to a “core,” bounded expected return times, and Doeblin or small-set minorization ensure unique QSD (Barbour et al., 2010, Barbour et al., 2011, Benaïm et al., 2014).
• Error bounds: For FV-type schemes, $O(N^{-1/2})$ convergence rates for empirical measures are standard; for stochastic approximation, polynomial or exponential rates apply depending on eigenvalue spectral gaps and step-size schedules (Benaïm et al., 2014, Journel et al., 20 Dec 2024).
• Particle approximations of QSDs for diffusions on compact domains admit explicit uniform-in-time bounds (Wasserstein or total variation) that decouple simulation-time, sample size, and discretization step (Journel et al., 2019, Journel et al., 20 Dec 2024).

When QSD is not unique, returned (quasi-equilibrium) distributions can still serve as effective approximations with error controlled by renewal and coupling analytic tools (Barbour et al., 2011). For models exhibiting metastability, the time-window over which the process stays near a QSD is both explicit and exponentially long in the relevant system-size parameter.

5. Comparison of Approaches and Computational Complexity

Method	Core Idea	Error Bounds
Fleming–Viot Particle	$N$ clones with resurrection	$O(N^{-1/2})$
Stochastic Approximation	Reinforced/self-interacting chain	$O(n^{-\theta})$, CLT
Returned-Process Stationary	Absorption replaced by return	Explicit TV bounds
Euler/occupation (diffusion)	Occupation of killed Euler schemes	Weak convergence; $O(h^{1/2})$

Particle schemes typically scale linearly in the sample size $N$ and can reach high accuracy but require synchronization and storage overhead. Stochastic approximation and returned-process methods scale with the number of iterations and are thus suited to very high-dimensional or continuous-state problems, with the trade-off of statistical noise and slower convergence.

6. Applications and Model-Specific Implications

Quasi-stationary approximation underpins modeling in:

• Ecology (population persistence times, extinction)
• Epidemiology (endemic levels, threshold analysis)
• Network science (epidemic and voter models, persistent activity)
• Physics (metastability in spin and transport models)
• Bayesian computation (Rare-event likelihoods, posterior simulation in killed diffusions)

In SIS models on networks, quasi-stationary approximations using ODE or pair-closure frameworks recapitulate both endemic finite positive prevalence and accurate epidemic thresholds inaccessible to traditional mean-field approximations (Overton et al., 2022, Cantwell et al., 15 Sep 2025). In mass-action and chemical kinetics, rigorous error estimation on QSDs informs the validity of diffusion or fluid approximations (Li et al., 2022).

7. Future Directions and Open Problems

Ongoing research seeks:

• Uniform-in-time and pathwise control for more general classes of Markov processes (including non-compact or non-uniformly elliptic models)
• Multi-scale and rare-event asymptotics for QSDs under singular perturbation or heterogeneous dynamics
• Algorithmic advances in particle approximations to prevent failure in degenerate regimes and accelerate convergence
• Tight complexity and performance guarantees for large-scale or real-time quasi-stationary simulation
• Rigorous sensitivity analysis for QSDs under model uncertainty or parameter perturbation

The field remains a nexus of probabilistic, analytical, and computational advances, with direct impact on both the theoretical understanding and numerical implementation of metastable phenomena in stochastic dynamical systems.