A Sensitivity Analysis Methodology for Rule-Based Stochastic Chemical Systems

TL;DR

This work develops a gradient-based sensitivity analysis framework for stochastic chemical systems, where intrinsic noise shapes observable outcomes. It employs forward finite-difference gradient estimation across a parameter subset, with adaptive replication guided by a gradient-uncertainty angular range $\omega$ to balance precision and cost. The method produces local sensitivity vector fields and global sensitivity coefficients for arbitrary user-defined observables, and validates the approach on Michaelis–Menten kinetics and a rule-based formose chemistry model implemented in MØD. By exploring parameter space and adaptively controlling simulation effort, the approach offers a robust, interpretable, and computationally efficient tool for understanding robustness and parameter importance in complex stochastic chemical networks.

Abstract

In this study, we introduce a sensitivity analysis methodology for stochastic systems in chemistry, where dynamics are often governed by random processes. Our approach is based on gradient estimation via finite differences, averaging simulation outcomes, and analyzing variability under intrinsic noise. We characterize gradient uncertainty as an angular range within which all plausible gradient directions are expected to lie. A key feature of our approach is that this uncertainty measure adaptively guides the number of simulations performed for each nominal-perturbation pair of points in order to minimize unnecessary computations while maintaining robustness. Systematically exploring a range of parameter values across the parameter space, rather than focusing on a single value, allows us to identify not only sensitive parameters but also regions of parameter space associated with different levels of sensitivity. These results are visualized through vector field plots to offer an intuitive representation of local sensitivity across parameter space. Additionally, global sensitivity coefficients over sampled points in the parameter space are computed to capture overall trends. Flexibility regarding the choice of output observable measures is another key feature of our method: while traditional sensitivity analyses often focus on species concentrations, our framework allows for the definition of a large range of problem-specific observables. This makes it broadly applicable in diverse chemical and biochemical scenarios. We demonstrate our approach on two systems: classical Michaelis-Menten kinetics and a rule-based model of the formose reaction, using the cheminformatics software MØD for Gillespie-based stochastic simulations.

Figures (12)

• Figure 1: Data Collection Workflow. For each nominal–perturbation pair, we first run a batch of $N$ simulations, evaluate the gradient uncertainty angle $\omega$ (see Figure \ref{['fig:angle_uncertainty']}), and, if necessary, refine the estimates by running additional batches of $N$ simulations for both points. Because the stopping criterion depends on the pair, the same number of extra simulations is assigned to each point. This strategy could be further optimized by: i) distributing extra simulations according to each point’s individual uncertainty, and ii) reusing nominal-point simulations for multiple perturbation points.
• Figure 2: Angular range of gradient uncertainty $\omega$. \ref{['fig:angle_a']} The individual confidence intervals (CIs) for both the nominal point and its perturbation in one dimension. \ref{['fig:angle_b']} The extreme possible gradients. \ref{['fig:angle_d']} The angular range within which the possible gradients exist.
• Figure 3: Input parameter points in the $2$-simplex.
• Figure 4: Stochastic simulation of the evolution (simulation time units) of Michaelis-Menten dynamics.
• Figure 5: Sensitivity analysis of the Michaelis-Menten Kinetics. \ref{['fig:fig1']} Gradient vectors projected onto the 2-simplex plane. Each point represents a sampled parameter set; the arrow at each point shows the local sensitivity (gradient vector) of the observable (simulation time to consume $S$) with respect to all three rate constants. The color of each point and vector indicates the value of the observable. Points with very small or zero gradient magnitude appear without a visible arrow. \ref{['fig:fig2']} Each sampled point colored according to the magnitude of the full, non-projected gradient vector; no arrows are shown in this panel. This highlights regions of high and low sensitivity. \ref{['fig:product']} Gradient vectors projected onto the three basic 2D planes of the 3D coordinate system, with colors indicating the observable as in panel (a). From these plots, we observe that sensitivity is largest with respect to the catalysis rate constant, especially when it is close to zero, whereas binding and unbinding rate constants contribute little to the variation in the observable.