A Spectral Koopman Approximation Framework for Stochastic Reaction Networks

TL;DR

This work tackles the challenge of analyzing stochastic reaction networks by introducing SKA, a spectral Stochastic Koopman Approximation. By exploiting the compactness of the Koopman operator, SKA separates dynamics into a state-independent decay-mode structure and small state-dependent coefficients, enabling continuous-time predictions for moments, event probabilities, and spectral metrics across all initial states. The framework provides computable a posteriori error bounds, robust parameter-sensitivity and cross-spectral-density estimations, and substantial computational speedups over trajectory-based baselines. Through biologically relevant networks and AIF-controlled systems, SKA demonstrates accurate, scalable analysis and opens avenues for initial-state inference and integration with learning-based methods.

Abstract

Stochastic reaction networks (SRNs) are a general class of continuous-time Markov jump processes used to model a wide range of systems, including biochemical dynamics in single cells, ecological and epidemiological populations, and queueing or communication networks. Yet analyzing their dynamics remains challenging because these processes are high-dimensional and their transient behavior can vary substantially across different initial molecular or population states. Here we introduce a spectral framework for the stochastic Koopman operator that provides a tractable, low-dimensional representation of SRN dynamics over continuous time, together with computable error estimates. By exploiting the compactness of the Koopman operator, we recover dominant spectral modes directly from simulated or experimental data, enabling efficient prediction of moments, event probabilities, and other summary statistics across all initial states. We further derive continuous-time parameter sensitivities and cross-spectral densities, offering new tools for probing noise structure and frequency-domain behavior. We demonstrate the approach on biologically relevant systems, including synthetic intracellular feedback controllers, stochastic oscillators, and inference of initial-state distributions from high-temporal-resolution flow cytometry. Together, these results establish spectral Koopman analysis as a powerful and general framework for studying stochastic dynamical systems across the biological, ecological, and computational sciences.

Figures (7)

• Figure 1: Analysis of the self-regulatory gene-expression network. (A) Schematic of the network. (B) Optimal cost $C_J^*$ of the convex optimization problem \ref{['defn_convex_optimization']} as a function of the number of decay modes $J$ (vertical axis on a logarithmic scale). (C) Identified decay modes on the complex plane, all with positive real parts and appearing in conjugate pairs. (D) Koopman trajectories for all observable functions in $\mathcal{F}$ (i.e., the first two moments \ref{['first_two_moments']}) generated by SKA, compared with SSA for an arbitrary initial state $x=(x_1,x_2) = (5, 10)$. (E) Parameter sensitivity estimates for selected observable functions with respect to a subset of model parameters, compared with CFD (same initial state as in panel D). (F) Cross-spectral density estimates for selected observable functions over the interval $[0,T]$ for two choices of $T$, compared with SSA–DFT (same initial state as in panel D). The full evolution of SKA-estimated CSDs from $T=0$ to $T=30$ seconds is provided in Supplementary Movie 1. When $f_1=f_2$, $\textnormal{CSD}_{f_1,f_2}(\omega,x,T)$ reduces to the real-valued power spectral density $\textnormal{PSD}_{f_1}(\omega,x,T)$; when $f_1 \neq f_2$, the magnitude and argument of $\textnormal{CSD}_{f_1,f_2}(\omega,x,T)$ are plotted separately. In each of the plots, the solid curve represents the mean while the shaded region represents the symmetric one standard deviation interval around the mean.
• Figure 2: Inference of the initial state distribution for the self-regulatory gene-expression network. (A) First two moments of the protein copy numbers computed from simulated empirical data, with each CTMC initialized from the product Poisson distribution \ref{['prod_poisson_distr']} with the true parameter vector $\hat{\eta} = (15, 50)$. (B) Marginal distributions of the inferred parameter vectors obtained from the SSA- and SKA-based approaches, each initialized from 100 randomly sampled parameter vectors. (C) Mean and standard deviation of the computational times required by the two approaches over 100 optimization runs.
• Figure 3: Analysis of the repressilator network. (A) Schematic of the network. Panels (B–F) correspond to the same analyses described in Figure \ref{['main_fig_self_reg_gene_ex']}, but here computed with the arbitrary initial state $x = (x_1,x_2,x_3) = (7, 10, 2)$. Corresponding to panel (F), the full evolution of SKA-estimated CSDs from $T=0$ to $T=30$ seconds is provided in Supplementary Movie 2.
• Figure 4: Effect of cooperativity on the repressilator network. We consider three repressilator networks with all Hill coefficients in the repression mechanism set to $H=1.0$, $1.5$, and $2.0$. (A) Decay modes of the three networks. (B) $\mathcal{L}_2(\hat{\pi})$ norm of the PSD of the $\mathbf{X}_1$ copy-number trajectory for the three networks, shown for two terminal times $T$, where $\hat{\pi}$ denotes the approximate stationary distribution. The full temporal dynamics from $T=0$ to $30$ seconds are provided in Supplementary Movie 5. (C) Probability distribution of the peak frequency ($\omega_0$), assuming the initial state is distributed according to $\hat{\pi}$, plotted at the two terminal times $T$. The full temporal dynamics from $T=0$ to $30$ seconds are provided in Supplementary Movie 6.
• Figure 5: Analysis of the closed-loop network consisting of constitutive gene-expression with the reference-based Antithetic Integral Feedback (rAIF) controller. (A) Schematic of the network. Panels (B–F) correspond to the same analyses shown in Figure \ref{['main_fig_self_reg_gene_ex']}, here evaluated for the arbitrary initial state $x = (x_1,x_2,z_1, z_2) = (5, 10, 11, 2)$. Corresponding to panel (F), the full evolution of SKA-estimated CSDs from $T=0$ to $T=30$ seconds is provided in Supplementary Movie 3.

Theorems & Definitions (11)

• Theorem 3.1: Frequency Domain Characterization of the Exact Finite Representation of the Koopman operator
• Remark 3.2
• Remark 3.3
• Remark 3.4
• Proposition A2.1
• proof
• proof : Proof of Theorem \ref{['thm_1']}
• Lemma A4.1
• proof
• Lemma A4.2