Data-driven Effective Modeling of Stochastic Chemical Reaction Networks

Yuan Chen; Weize Mao; Dongbin Xiu

math.NAAug 26, 2026

Data-driven Effective Modeling of Stochastic Chemical Reaction Networks

Authors: Yuan Chen, Weize Mao, Dongbin Xiu

Organizations: Department of Mathematics, The Ohio State University, Columbus, OH 43210, USA.

Abstract

The Stochastic Simulation Algorithm (SSA), widely considered an exact algorithm for stochastic chemical reaction networks, suffers from high computational cost. In this work, we propose a data-driven effective model that operates on a user-defined coarse time step independent of the underlying microscopic reaction-event scale. This is accomplished by directly approximating the finite-time transition kernel of the continuous-time Markov chain induced by SSA, using a generative machine learning model trained on short bursts of SSA simulation data. The trained model constructs a stochastic propagator that recursively generates statistically consistent trajectories at the constant coarse time step, with significantly reduced computational cost. In this paper, we employ conditional normalizing flow as the stochastic propagator. A comprehensive set of numerical examples is presented to demonstrate the accuracy and efficiency of the proposed method.

Figures & tables

Figure 1: The DNN model structure for the proposed normalizing flow sFML method ( 12 ).

Example	$T$	sFML Steps	SSA		Adaptive Tau-leaping
Example	$T$	sFML Steps	Steps	Ratio	Steps	Ratio
Transfer	$10$	$1.0\times 10^{2}$	$1.9\times 10^{2}$	$1.93\times$	$1.8\times 10^{2}$	$1.75\times$
LV (slow)	$20$	$2.0\times 10^{2}$	$7.2\times 10^{3}$	$36\times$	$2.5\times 10^{3}$	$12.3\times$
LV (fast)	$20$	$2.0\times 10^{3}$	$6.0\times 10^{5}$	$300\times$	$3.2\times 10^{4}$	$15.8\times$
Brusselator	$15$	$1.5\times 10^{3}$	$1.6\times 10^{6}$	$1070\times$	$2.1\times 10^{5}$	$143\times$
Autocatalysis	$10$	$1.0\times 10^{3}$	$3.2\times 10^{5}$	$317\times$	$7.6\times 10^{3}$	$7.65\times$

Table 1: Comparison of simulation steps. SSA and tau-leaping steps are estimated as averages over 10,000 simulated trajectories and rounded to two significant figures. sFML steps equal

T/\Delta

exactly. The ratios for SSA and adaptive tau-leaping are defined as their respective steps divided by the sFML steps.

Example	Min MMD	Median MMD	Max MMD
Transfer	$1.89\times 10^{-2}$	$5.59\times 10^{-2}$	$1.05\times 10^{-1}$
LV (slow)	$1.72\times 10^{-2}$	$4.46\times 10^{-2}$	$9.60\times 10^{-2}$
LV (fast)	$1.55\times 10^{-2}$	$5.07\times 10^{-2}$	$8.72\times 10^{-2}$
Brusselator	$8.09\times 10^{-2}$	$2.29\times 10^{-1}$	$2.83\times 10^{-1}$
Autocatalysis	$1.81\times 10^{-2}$	$4.18\times 10^{-2}$	$9.35\times 10^{-2}$
Oregonator	$5.28\times 10^{-2}$	$1.53\times 10^{-1}$	$3.09\times 10^{-1}$

Table 2: One-step transition-kernel comparison. For each example, MMD is computed at 30 conditioning states sampled from one SSA test trajectory. The table reports the minimum, median, and maximum MMD over the selected states.

Example	$T$	$\Delta$	Method	$E_{\mu}$	$E_{\sigma}$
Transfer	$10$	$0.1$	sFML	$1.78\times 10^{-3}$	$7.92\times 10^{-2}$
Transfer	$10$	$0.1$	Tau-leaping	$7.64\times 10^{-3}$	$4.81\times 10^{-2}$
LV (slow)	$20$	$0.1$	sFML	$2.16\times 10^{-2}$	$3.07\times 10^{-2}$
LV (slow)	$20$	$0.1$	Tau-leaping	$1.70\times 10^{-1}$	$3.04\times 10^{0}$
LV (fast)	$20$	$0.01$	sFML	$8.31\times 10^{-3}$	$5.43\times 10^{-2}$
LV (fast)	$20$	$0.01$	Tau-leaping	unstable	unstable

Table 3: Global trajectory comparison with fixed-step tau-leaping. Tau-leaping uses the same coarse step size

\tau=\Delta

as sFML. Errors are estimated from 10,000 trajectories.

Panel of Figure 2 — Figure 2: Mean and standard deviation (STD) of Example 4.1 with initial condition $(83,26,69)$ : upper left: $X_{1}$ , upper right: $X_{2}$ , lower left: $X_{3}$ . Lower right: sample trajectories of Example 4.1 with initial condition $(83,26,69)$ .

Figure 3: Test sample trajectories of Example 4.2 with initial condition

(100,100)

and parameter

c_{1}=1.0

c_{2}=0.005

and

c_{3}=0.6

up to

T=20

Figure 4: Mean and standard deviation (STD) of Example 4.2 with initial condition

(100,100)

and parameters

c_{1}=1.0

c_{2}=0.005

and

c_{3}=0.6

up to

T=20

: left: number of

X_{1}

, right: number of

X_{2}

Panel of Figure 5 — Figure 5: Sample trajectories produced by trained model of Example 4.2 with initial condition $(1000,1000)$ and parameter $c_{1}=10$ , $c_{2}=0.01$ and $c_{3}=10$ up to $T=20$ . Upper left, upper right, lower left: 3 independent simulations, lower right: phase plot of upper left.

Panel of Figure 6 — Figure 6: Sample trajectories produced by SSA of Example 4.2 with initial condition $(1000,1000)$ and parameter $c_{1}=10$ , $c_{2}=0.01$ and $c_{3}=10$ up to $T=20$ . Upper left, upper right, lower left: 3 independent simulations, lower right: phase plot of upper left.

Panel of Figure 7 — Figure 7: Comparison of one-step conditional distribution of Example 4.2 with condition $\mathbf{x}=(540,734)$ and parameter $c_{1}=10$ , $c_{2}=0.01$ and $c_{3}=10$ . Upper left: ground truth of joint distribution of $(X_{1},X_{2})$ , upper right: learned joint distribution of $(X_{1},X_{2})$ , lower left: comparison of the marginal distribution for $X_{1}$ , lower right: comparison of the marginal distribution for $X_{2}$ .

Panel of Figure 8 — Figure 8: Comparison of sample trajectories of Example 4.3 with initial condition $(1000,2000)$ and parameter $c_{1}\overline{Y}_{1}=5000$ , $c_{2}\overline{Y}_{2}=50$ , $c_{3}=5\times 10^{-5}$ and $c_{4}=5$ up to $T=15$ . Upper left: a sample trajectory in state space produced by SSA, upper right: a sample trajectory in state space produced by learned model, lower left: phase plot of upper left, lower right: phase plot of upper right.

Panel of Figure 9 — Figure 9: Comparison of one-step conditional distribution of Example 4.3 with initial condition $(1959,971)$ and parameter $c_{1}\overline{Y}_{1}=5000$ , $c_{2}\overline{Y}_{2}=50$ , $c_{3}=5\times 10^{-5}$ and $c_{4}=5$ . Upper left: ground truth of joint distribution of $(X_{1},X_{2})$ , upper right: learned joint distribution of $(X_{1},X_{2})$ , lower left: comparison of the marginal distribution for $X_{1}$ , lower right: comparison of the marginal distribution for $X_{2}$ .

Figure 10: Comparison of Fourier modes for sFML generated trajectory and ground truth of Example 4.3 with initial condition

(1000,2000)

and parameter

c_{1}\overline{Y}_{1}=5000

c_{2}\overline{Y}_{2}=50

c_{3}=5\times 10^{-5}

and

c_{4}=5

up to

T=15

. Left: Fourier modes of SSA trajectory (upper left of Figure 8 ), right: Fourier modes of model trajectory (upper right of Figure 8 ).

Panel of Figure 11 — Figure 11: Comparison of sample trajectories of Example 4.4 with initial condition $(3000,2000,2000)$ and parameter $c_{1}=0.002$ , $c_{2}=0.002$ and $c_{3}=0.002$ up to $T=10$ . Upper left, upper right: 2 independent sample trajectories produced by SSA, lower left, lower right: 2 independent sample trajectories produced by the trained model.

Figure 12: Comparison of Fourier modes for sFML generated trajectory and ground truth of Example 4.4 with initial condition

(3000,2000,2000)

and parameter

c_{1}=0.002

c_{2}=0.002

and

c_{3}=0.002

up to

T=10

. Left: Fourier modes of SSA trajectory (upper right of Figure 11 ), right: Fourier modes of model trajectory (lower right of Figure 11 ).

Figure 13: Comparison of sample trajectories of Example 4.5 in state space with initial condition

(500,1000,2000)

and parameter

c_{1}\overline{Y}_{1}=2

c_{2}=0.1

c_{3}\overline{Y}_{2}=104

c_{4}=0.016

and

c_{5}\overline{Y}_{3}=26

up to

T=6

. Left: a sample trajectory produced by SSA, right: a sample trajectory produced by the trained model.

Panel of Figure 14 — Figure 14: Comparison of sample trajectories of Example 4.5 in phase space with initial condition $(500,1000,2000)$ and parameter $c_{1}\overline{Y}_{1}=2$ , $c_{2}=0.1$ , $c_{3}\overline{Y}_{2}=104$ , $c_{4}=0.016$ and $c_{5}\overline{Y}_{3}=26$ up to $T=6$ . Upper: sample trajectory produced by SSA, lower: sample trajectory produced by trained model.

Figure 15: Comparison of Fourier modes for two sample trajectories of Example 4.5 with initial condition

(500,1000,2000)

and parameter

c_{1}\overline{Y}_{1}=2

c_{2}=0.1

c_{3}\overline{Y}_{2}=104

c_{4}=0.016

and

c_{5}\overline{Y}_{3}=26

up to

T=6

. Left: Fourier modes of SSA trajectory (left of Figure 13 ), right: Fourier modes of model trajectory (right of Figure 13 ).

Explore similar work

Tracking large chemical reaction networks and rare events by neural networks
Dec 11, 2025Jiayu Weng, Xinyi Zhu, Jing Liu +3Dynamical Systems Neural Surrogate Modeling
Institute of Data Science, University of Hong Kong, Hong Kong · Department of Systems Science, Faculty of Arts and Sciences, Beijing Normal University, Zhuhai 519087, China · School of Physics, University of Electronic Science and Technology of China, Chengdu 611731, China +6
Exact Discrete Stochastic Simulation with Deep-Learning-Scale Gradient Optimization
Feb 23, 2026Jose M. G. Vilar, Leonor SaizDifferentiable Physics
Biofisika Institute (CSIC, UPV/EHU), University of the Basque Country (UPV/EHU), P.O. Box 644, 48080 Bilbao, Spain · IKERBASQUE, Basque Foundation for Science, 48011 Bilbao, Spain · Department of Biomedical Engineering, University of California, 451 E. Health Sciences Drive, Davis, CA 95616, USA
Conditional Score-Based Modeling of Effective Langevin Dynamics
Apr 27, 2026Ludovico T. GiorginiLangevin Dynamics Model Calibration
Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Dec 11, 2025J/K move · Enter open · S save

Tracking large chemical reaction networks and rare events by neural networks

Jiayu Weng, Xinyi Zhu, Jing Liu, Linyuan Lü+2

Institute of Data Science, University of Hong Kong, Hong Kong · Department of Systems Science, Faculty of Arts and Sciences, Beijing Normal University, Zhuhai 519087, China · School of Physics, University of Electronic Science and Technology of China, Chengdu 611731, China+6

Dynamical Systems Neural Surrogate Modeling Markov Models

Open paper AskPDF arXiv

Chemical reaction networks are widely used to model stochastic dynamics in biophysical processes. Their exact probabilistic description via the chemical master equation becomes intractable for high-dimensional systems due to the exponential growth of the state space. This difficulty is further exacerbated when metastable transitions and rare events dominate the dynamics. Here, we present an advanced neural-network approach that effectively tackles large reaction networks and rare-event dynamics within a unified framework. Building on autoregressive representations of the full joint probability distribution, we achieve efficient training through second-order optimization methods, including natural gradient descent and the time-dependent variational principle. We leverage enhanced-sampling strategies that enable accurate resolution of low-probability regions governing rare transitions. We demonstrate the method on a multistable genetic toggle switch and a mitogen-activated protein kinase (MAPK) cascade consisting of 16 species and 35 reactions with dense coupling. We further quantify rare events in spatially extended Schlögl models in both one- and two-dimensional lattices, reaching system sizes with extremely large state spaces up to

3\times10^{15}

. These results highlight the potential of neural modeling as a flexible and scalable framework for studying high-dimensional chemical reaction networks and rare events.

Abstract

Explore similar work

Tracking large chemical reaction networks and rare events by neural networks

Exact Discrete Stochastic Simulation with Deep-Learning-Scale Gradient Optimization

Conditional Score-Based Modeling of Effective Langevin Dynamics