PaperlayerView exact revision on arXiv

Source-linked AI summary

A scalable computational framework for establishing long-term behavior of stochastic reaction networks

Ankit Gupta, Corentin Briat, Mustafa Khammash

arXiv:1304.5404v3q-bio.MNcs.SYeess.SYmath.OCmath.PR

TL;DR

Long-term stability of stochastic reaction networks is difficult to establish directly, particularly when the state space is infinite and simulation or deterministic approximations do not settle the question. The paper develops sufficient theoretical conditions and computational methods for assessing ergodicity and moment behavior. It demonstrates these conditions on reaction networks and discusses a bimolecular-network condition and optimization complexity that bound their scope.

  • Problem

    For infinite-state stochastic reaction networks, direct stability analysis is difficult because the Chemical Master Equation is generally intractable and finite simulations cannot establish long-term behavior.

    Source: Introduction · p. 3

  • Method

    The paper uses sufficient linear Foster-Lyapunov conditions and computational optimization to assess ergodicity and moment properties of stochastic reaction networks.

    Source: Introduction · p. 4, Results for stochastic bimolecular… · p. 14

  • Results

    The authors demonstrate that their sufficient conditions establish stability properties for several reaction networks from the literature.

    Source: Introduction · p. 4

  • Takeaways & Limitations

    Moment bounds and convergence results can help assess stochastic-model suitability and the assumptions underlying moment-closure methods.

    Source: Introduction · p. 4

    For the bimolecular result, a non-empty Nq is required, equivalent to a conservation relation across all bimolecular reactions.

    Source: Results for stochastic bimolecular… · p. 14

Abstract

from arXiv · show

Reaction networks are systems in which the populations of a finite number of species evolve through predefined interactions. Such networks are found as modeling tools in many biological disciplines such as biochemistry, ecology, epidemiology, immunology, systems biology and synthetic biology. It is now well-established that, for small population sizes, stochastic models for biochemical reaction networks are necessary to capture randomness in the interactions. The tools for analyzing such models, however, still lag far behind their deterministic counterparts. In this paper, we bridge this gap by developing a constructive framework for examining the long-term behavior and stability properties of the reaction dynamics in a stochastic setting. In particular, we address the problems of determining ergodicity of the reaction dynamics, which is analogous to having a globally attracting fixed point for deterministic dynamics. We also examine when the statistical moments of the underlying process remain bounded with time and when they converge to their steady state values. The framework we develop relies on a blend of ideas from probability theory, linear algebra and optimization theory. We demonstrate that the stability properties of a wide class of biological networks can be assessed from our sufficient theoretical conditions that can be recast as efficient and scalable linear programs, well-known for their tractability. It is notably shown that the computational complexity is often linear in the number of species. We illustrate the validity, the efficiency and the wide applicability of our results on several reaction networks arising in biochemistry, systems biology, epidemiology and ecology. The biological implications of the results as well as an example of a non-ergodic biological network are also discussed.

Author Summary

Small populations make stochastic modeling necessary to capture random interactions in biological networks. The paper develops a framework for assessing stochastic reaction networks’ long-term behavior and stability.

Author Summary

  • Small populations make random effects significant, so stochastic models are needed to represent biological reaction networks.Such models are widely used, including to study heterogeneity in clonal cell populations.Source: Author Summary · p. 2
  • The paper addresses the lack of analytical tools for assessing long-term behavior and stability in stochastic reaction networks.Researchers have often relied on simulation-based approaches, which the authors describe as unsatisfactory.Source: Author Summary · p. 2
  • Its framework combines probability, linear algebra, and optimization, and is illustrated with examples from several biological areas.The paper also discusses biological implications of its results.Source: Author Summary · p. 2

Introduction

The introduction motivates direct analysis of stochastic reaction networks because their long-term behavior can differ from deterministic models, while infinite-state CME analysis and finite simulations have limitations. The paper presents sufficient theoretical conditions for stability and moment properties, including ergodicity.

Introduction

  • For infinite state spaces, solving the Chemical Master Equation is generally difficult, and truncation methods face a curse of dimensionality.The CME has one equation per state; numerical approaches work well only for relatively small systems, although bounded-rank tensor representations can scale linearly in species number.Source: Introduction · p. 2
  • Finite Monte Carlo simulations cannot establish long-term stability directly, while large-population SDE approximations may fail for low-copy species and diverge in long-term behavior.SSA simulation time also depends linearly on the number of reactions during the simulated period.Source: Introduction · p. 3
  • The framework uses linear Foster-Lyapunov functions to establish ergodicity, moment bounds and convergence, and attractive sets for moments.Its sufficient conditions also support structural conclusions for classes of networks and are demonstrated on networks from the literature.Source: Introduction · p. 4
  • Stochastic stability cannot generally be inferred from deterministic stability: examples include ergodic stochastic toggle-switch and repressilator models, and a non-ergodic model with a stable deterministic fixed point.The paper also describes stochastic trajectories escaping a deterministic region of attraction and becoming unbounded.Source: Introduction · p. 4
  • The framework checks whether process moments remain bounded and converge, and provides bounds uniform in time and bounds on steady-state moments.These results can help assess model suitability and support biological controller design; they can also test moment-closure assumptions, including existence of a moment-generating function for all times.Source: Introduction · pp. 4–5
  • The deterministic analogues DC1–DC3 describe bounded trajectories, an eventual common attractor set, and convergence to a globally attracting fixed point.The authors note that dynamical-systems methods can verify these conditions without simulation.Source: Introduction · p. 5
  • For an infinite state space, explicit or numerical CME solutions are generally nearly impossible, motivating analysis that avoids solving the distribution equations directly.The Markov process evolves by discrete reaction events, with states given by molecular-count vectors.Source: Introduction · pp. 5–6
  • The paper frames stochastic stability through SC1–SC3, from compactness of distributions to a globally attracting stationary distribution.SC3 corresponds to ergodicity; then long-run trajectory time proportions can reveal stationary probabilities, with applications to single-cell and population measurements.Source: Introduction · pp. 6–7

Results

The framework uses a weighted norm satisfying the drift-diffusivity condition to establish bounds and long-term properties of stochastic reaction networks. Under stated conditions, it gives moment bounds, exponential ergodicity, moment convergence, and light-tailedness results.

Preliminaries

  • Condition DD bounds drift and diffusivity in a weighted norm; above c1/c2, the norm process experiences negative drift.The weight vector v is chosen so that the norm process has negative drift at large values.Source: Preliminaries · pp. 7–9
  • The weight vector v can be determined for a large class of networks by solving suitably constructed optimization problems.The selected norm is used to measure state size when verifying long-term behavior conditions.Source: Preliminaries · pp. 7–8

General Results

  • For each positive integer r < rmax, Condition DD yields uniform and asymptotic upper bounds on the process's r-th moment.The constants in these bounds can be explicitly computed using a recursive relationship.Source: General Results · pp. 9–10
  • When Condition DD holds with c5 = 0, the process is uniformly light-tailed from each initial state and the stationary distribution is light-tailed.In this case, rmax is infinite.Source: General Results · p. 10
  • Under Condition DD, moments of integer order r < rmax−1 converge to finite stationary moments; with irreducibility and DD1, the process is exponentially ergodic.The distribution converges exponentially to a unique stationary distribution under irreducibility and DD1.Source: General Results · p. 10
  • For networks satisfying Condition DD, the framework's results apply broadly, including monomolecular and bimolecular mass-action networks and networks with Hill-function propensities.The framework describes itself as applicable to any network satisfying Condition DD.Source: General Results · p. 10

Methods

The framework uses Condition DD to establish long-term properties of stochastic reaction networks, then turns sufficient conditions into tractable optimization problems. For unimolecular networks, a positive-vector test characterizes ergodicity and supports scalable checks, including under some rate uncertainties.

Methods

  • Condition DD is the framework’s condition for establishing uniform-in-time moment bounds and global convergence to equilibrium, respectively, through Theorems 2 and 5.The Methods section develops a theoretical and computational framework for checking this condition.Source: Methods · p. 11

Results for stochastic unimolecular reaction networks

  • The unimolecular results cover general reaction types and can determine whether a unique stationary distribution exists without solving the Chemical Master Equation.The approach also informs the analysis of bimolecular networks.Source: Results for stochastic unimolecular… · pp. 11–12
  • For irreducible unimolecular networks, Hurwitz stability of A is equivalent to the existence of v > 0 with Av < 0; either condition gives exponential ergodicity and bounded, globally converging moments.The positive-vector condition can be checked with a linear program.Source: Results for stochastic unimolecular… · p. 12
  • If moments need only be bounded and convergent, A may have zero eigenvalues, and their limits may depend on the initial conditions.Source: Results for stochastic unimolecular… · p. 12
  • With uncertain rates but known network structure, an upper-bound matrix A+ can support a robust ergodicity test when it exists and the process is irreducible for every allowed parameter value.Under these conditions, the process is robustly exponentially ergodic, with bounded and globally converging moments; checking the family reduces to checking A+.Source: Results for stochastic unimolecular… · p. 12
  • The first-order moment-bound optimization uses d + 1 variables and 2d + 1 constraints, so its size scales linearly with species count and is independent of reaction count.A bisection over z solves a linear program at each iteration.Source: Results for stochastic unimolecular… · p. 13

Results for stochastic bimolecular reaction networks

For bimolecular networks, the paper extends unimolecular tests using conservation relations or additional matrix conditions. The resulting criteria yield ergodicity and moment conclusions, with computational checks ranging from linear programs to semidefinite programs.

Results for stochastic bimolecular reaction networks

  • When the bimolecular reactions admit a non-empty Nq, any v in Nq with Av < 0 implies ergodicity and bounded, globally converging moments.For v in Nq, the quadratic term vanishes and the resulting expression matches the unimolecular case.Source: Results for stochastic bimolecular… · p. 14
  • Nq must be non-empty for this criterion; this is equivalent to a conserved positive linear combination of species populations across every bimolecular reaction.A necessary condition is that Sq is not full-row rank.Source: Results for stochastic bimolecular… · p. 14
  • When the Nq criterion is unavailable, either Av < 0 with M(v) ≤ 0 or a positive v with negative-definite M(v) implies ergodicity and bounded, globally converging moments up to the stated order.The moment order is ⌊1 + 2c2/c5⌋−2.Source: Results for stochastic bimolecular… · p. 14
  • The first alternative is checkable by linear programming, while the second requires a semidefinite program.Source: Results for stochastic bimolecular… · p. 14
  • The optimization problem for the conservation-constrained criterion scales linearly with the number of species, allowing checks on large networks.The paper identifies it as the computational translation of Proposition 10.Source: Results for stochastic bimolecular… · pp. 14–15
  • The alternative optimization involving M(v) can require quadratically many constraints in the number of species in the worst case, though its linear-programming structure keeps it tractable.Source: Results for stochastic bimolecular… · p. 15

Qualitative differences between deterministic and stochastic dynamics

The examples show that stochastic reaction networks can have long-term behavior sharply different from their deterministic counterparts, so deterministic stability does not in general establish stochastic ergodicity or moment convergence.

Jumping potential wells

  • In the potential-well example, stochastic drift becomes positive for x ≥ 1 + x+, and the drift-negative region exceeds the deterministic basin because the bimolecular propensity differs between models.The largest root is reported as x+ ≃ 26.194 > β.Source: Jumping potential wells · pp. 15–16
  • From the same initial state, all 1000 stochastic simulations were stopped before 100 seconds as the reaction rate led to unbounded trajectories, despite deterministic convergence.The simulations stopped when the bimolecular propensity exceeded the rate corresponding to 15,000 molecules.Source: Jumping potential wells · p. 16, Figure 1 · p. 16

Globally stable deterministic dynamics does not imply moments stability

  • Despite a globally exponentially stable deterministic equilibrium, network (26) is not ergodic and all its moments grow unboundedly with time.For initial conditions with X1(0) − X2(0) = α > 0, E[X1(t) − X2(t)] remains α for all t ≥ 0.Source: Globally stable deterministic dynamics… · p. 17

Finding an attractive compact set for the first-order moments

For a bimolecular network, the authors use optimization to bound first-order moments within an attractive compact set, and show how that set can also assess moment-closure predictions.

Finding an attractive compact set for the first-order moments

  • The set is computed for the first-order moment of ⟨v, X(t)⟩ using optimization problems (17) or (20), addressing bimolecular networks whose steady-state moments lack analytical expressions.The stated obstacle is the moment closure problem.Source: Finding an attractive compact set for… · p. 17
  • For the irreducible dimerization network, the optimization yields exponential ergodicity, with all moments bounded and converging.The reported bound satisfies 1/c*_2 = 5, matching the theoretical value k/min{γ1, γ2}.Source: Finding an attractive compact set for… · pp. 17–18
  • The computed attractive compact set contains the unique steady-state first-order moments, which converge close to its boundary across initial conditions.Figure 3 averages over 5000 cells.Source: Finding an attractive compact set for… · p. 18, Figure 3 · p. 19
  • Comparing a moment-closure equilibrium with the attractive set can test consistency, but the set gives only a lower bound on closure error when that equilibrium lies outside it.The discrepancy between closed-system and stochastic equilibrium moments is expected because the closure neglects variance.Source: Finding an attractive compact set for… · pp. 18–19

Feedback loop

The feedback-loop example models a dimer acting on gene expression through a bounded nonnegative feedback function and establishes stability for that broad class of feedback rules.

Feedback loop

  • For any positive rate parameters and any bounded nonnegative feedback function, the dimerization feedback loop is ergodic, with all moments bounded and globally converging.The dimer S3 acts back on gene expression through f(·).Source: Feedback loop · p. 20

Stochastic switch

The stochastic-switch example establishes ergodicity and global moment convergence for arbitrary bounded nonnegative feedback functions and positive rate parameters.

Stochastic switch

  • For any positive rate parameters and bounded nonnegative functions f1(·) and f2(·), the stochastic switch is ergodic, with all moments bounded and globally converging.The network represents a stochastic switch described by.Source: Stochastic switch · pp. 20–21

Repressilator

The section considers the stochastic N-gene repressilator and establishes ergodicity and convergence of all moments for every positive choice of its rate parameters.

  • The model is an N-gene repressilator whose reaction network is specified in the section.Source: Repressilator · p. 21
  • Each gene has a corresponding mRNA and protein species.Source: Repressilator · p. 22
  • For any positive rate parameters, the stochastic N-gene repressilator is ergodic, and all moments are bounded and globally converging.Source: Repressilator · p. 22

Stochastic SIR model

The SIR and circadian-clock examples are shown to be ergodic with bounded, globally converging moments for positive rate parameters. For the circadian clock, sample averages converge to stationary values that equal the asymptotic time averages.

Stochastic SIR model

  • For any positive rate parameters, the SIR model is ergodic, and all moments are bounded and globally converging.Source: Stochastic SIR model · p. 22

Circadian clock

  • For any positive rate parameters, the 9-species, 18-reaction circadian clock model is ergodic, and all moments are bounded and globally converging.Source: Circadian clock · p. 22
  • The circadian-clock protein sample averages converge to stationary values that, by ergodicity, are globally attracting and equal the asymptotic time averages.For A, R, and C, the reported steady-state averages are 222.1797, 534.8853, and 549.7195, respectively.Source: Circadian clock · p. 22

p53 model

The section applies the framework to p53, Lotka–Volterra, and Schlögl models, establishing ergodicity under the stated conditions. For the Schlögl model, moment stability beyond first-order convergence remains unresolved by the stated condition.

p53 model

  • The p53 model defines S1 as p53 molecules, S2 as precursor Mdm2 molecules, and S3 as Mdm2 molecules, with nonlinear feedback on p53 degradation.Source: p53 model · p. 22
  • For any positive rate parameters, the oscillatory p53 model is ergodic, and all moments are bounded and globally converging.Source: p53 model · p. 22

Lotka-Volterra model

  • The stochastic Lotka–Volterra network is ergodic with bounded, globally converging moments under either the stated positive-definiteness condition or the copositivity condition plus β_i − δ_i < 0 for all i.Source: Lotka-Volterra model · p. 23

Schl¨ogl model

  • For positive rate parameters and positive XA and XB, the Schlögl model is exponentially ergodic.Source: Schl¨ogl model · p. 24
  • For the Schlögl model, the stated condition does not establish moment stability beyond first-order convergence because the model contains a cubic term.Source: Schl¨ogl model · p. 24

Discussion

The discussion connects ergodicity and moment bounds to experimental interpretation and model assessment, while illustrating that a biological network can be non-ergodic for some parameters. The framework also has potential applications to synthetic-biology design, whose specifics are outside this article’s scope.

  • The framework is presented as applicable beyond mass-action kinetics, including kinetics described by Hill functions.Source: Discussion · p. 24
  • Ergodicity can make single-cell observations over time and population snapshots yield the same stationary-distribution information.Source: Discussion · p. 25
  • For the carcinogenesis model, k1 > γ2 gives unbounded species trajectories, illustrating non-ergodicity for that parameter choice.Source: Discussion · p. 26
  • The analysis may inform synthetic biological circuit design when desired behavior must hold despite unknown or difficult-to-control initial conditions.Source: Discussion · p. 26
  • Moment bounds can assess whether higher-order moments remain bounded, while uniform light-tailedness under certain conditions ensures the moment-generating function exists for all time.Source: Discussion · pp. 26–27
  • The specifics of using these techniques to design circuits that steer moments to a target steady-state value are outside this article’s scope.Source: Discussion · p. 27
Loading 1304.5404v3…