Global parameter identification of stochastic reaction networks from single trajectories

We consider the problem of inferring the unknown parameters of a stochastic biochemical network model from a single measured time-course of the concentration of some of the involved species. Such measurements are available, e.g., from live-cell fluorescence microscopy in image-based systems biology. In addition, fluctuation time-courses from, e.g., fluorescence correlation spectroscopy provide additional information about the system dynamics that can be used to more robustly infer parameters than when considering only mean concentrations. Estimating model parameters from a single experimental trajectory enables single-cell measurements and quantification of cell–cell variability. We propose a novel combination of an adaptive Monte Carlo sampler, called Gaussian Adaptation, and efficient exact stochastic simulation algorithms that allows parameter identification from single stochastic trajectories. We benchmark the proposed method on a linear and a non-linear reaction network at steady state and during transient phases. In addition, we demonstrate that the present method also provides an ellipsoidal volume estimate of the viable part of parameter space and is able to estimate the physical volume of the compartment in which the observed reactions take place.

💡 Research Summary

This paper addresses the challenging problem of inferring kinetic parameters and the reaction compartment volume of stochastic biochemical networks when only a single noisy time‑course of a subset of species is available. Such data arise naturally from live‑cell fluorescence microscopy and fluorescence correlation spectroscopy (FCS), which provide both mean concentration trajectories and fluctuation spectra at single‑cell resolution. Traditional parameter estimation based solely on mean trajectories suffers from non‑identifiability, especially in low‑copy‑number regimes where stochastic fluctuations carry essential information.

The authors propose a novel framework that couples an adaptive Monte‑Carlo sampler—Gaussian Adaptation (GaA)—with exact stochastic simulation algorithms based on the partial‑propensity formulation of Gillespie’s SSA. GaA operates on a multivariate Gaussian proposal distribution N(m, Σ) whose mean vector m and covariance matrix Σ are iteratively updated to maintain a predefined acceptance rate while maximizing the entropy of the search distribution. This makes GaA robust to noisy, non‑differentiable objective functions, a key requirement when the forward model is stochastic and the data consist of a single trajectory.

The forward model is the full Chemical Master Equation (CME). Rather than solving the CME analytically, the authors generate synthetic trajectories by exact SSA, ensuring that each sampled parameter set yields a statistically exact realization of the stochastic dynamics. To evaluate the match between simulated and experimental data, they introduce a composite distance metric that combines the mean‑squared error of the concentration time‑course with a term quantifying the discrepancy between the power spectral densities of the experimental and simulated fluctuations. By incorporating both mean dynamics and fluctuation spectra, the metric exploits all available information, dramatically reducing the viable parameter region.

The algorithm proceeds as follows: (1) initialize m and Σ; (2) draw a candidate parameter vector θ from N(m, Σ); (3) simulate a single stochastic trajectory using partial‑propensity SSA; (4) compute the composite distance f(θ); (5) accept θ if f(θ) falls below a dynamically adjusted threshold that enforces the target acceptance rate; (6) update m and Σ based on accepted samples; repeat steps 2‑6 for many iterations. Accepted samples not only provide point estimates of the kinetic rates and volume Ω but also define an ellipsoidal confidence region via the final covariance matrix, allowing quantitative assessment of parameter uncertainty.

The methodology is benchmarked on two models: (i) a linear chain of five reversible reactions, representing a monostable network with analytically tractable dynamics; and (ii) a nonlinear colloidal aggregation model, which exhibits more complex stochastic behavior. In both cases, the algorithm accurately recovers the true kinetic parameters and the reaction volume, even when only a single noisy trajectory is used. Inclusion of the fluctuation‑based term in the distance metric is shown to be crucial for the nonlinear case, where mean‑only fitting leads to ambiguous solutions. Moreover, the ellipsoidal volume estimate reliably captures the region of parameter space that yields acceptable fits, providing a useful diagnostic for model identifiability.

Key contributions of the work are:

1. Demonstration that single‑trajectory data, when combined with fluctuation information, suffice for global stochastic parameter identification.
2. Introduction of Gaussian Adaptation as a robust, gradient‑free optimizer and approximate Bayesian sampler tailored to noisy, black‑box objective functions.
3. Integration of partial‑propensity exact SSA to keep computational costs tractable while preserving exact stochastic dynamics.
4. Provision of quantitative uncertainty estimates (ellipsoidal confidence regions) and simultaneous inference of the physical compartment volume Ω, which is often experimentally inaccessible.

The authors discuss extensions such as hierarchical Bayesian inference across multiple cells, application to multistable or oscillatory systems, and online variants of GaA for real‑time parameter tracking. Overall, the paper delivers a practical, theoretically sound tool for quantitative systems biology at the single‑cell level, opening new avenues for model calibration, experimental design, and the study of cell‑to‑cell variability.