Tracking large chemical reaction networks and rare events by neural networks

Jiayu Weng; Xinyi Zhu; Jing Liu; Linyuan Lü; Pan Zhang; Ying Tang

q-bio.MNDec 11, 2025

Tracking large chemical reaction networks and rare events by neural networks

Authors: Jiayu Weng, Xinyi Zhu, Jing Liu, Linyuan Lü, Pan Zhang, Ying Tang

Organizations: 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 · School of Physical Science and Technology, Beijing University of Posts and Telecommunications, Beijing 102206, China · Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China · School of Cyber Science and Technology, University of Science and Technology of China, Hefei 230027, China · School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China · Institute of Fundamental and Frontier Sciences, University of Electronic Science and Technology of China, Chengdu 611731, China · Non-classical Information Science Basic Discipline Research Center of Sichuan Province, University of Electronic Science and Technology of China, Chengdu 611731, China

Abstract

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.

Figures & tables

Figure 1: The present work efficiently tracks large reaction networks and rare events. (a) The chemical reaction system follows the CME, whose exponentially large state space makes the exact probability distribution intractable. (b) The variational autoregressive network factorizes the joint distribution into conditional probabilities. We use NADE-based architectures to improve training efficiency over sequential RNNs. (c) Our advances include both reducing computational cost to deal with large networks and quantifying rare-event statistics. (d) The VAN is trained by minimizing the KL-divergence between time steps. Moving from stochastic gradient descent (SGD), we now use natural gradient (NG), which improves convergence by rescaling updates with the inverse Fisher information matrix and the time-dependent variational principle (TDVP) with the projected temporal evolution. In the genetic toggle switch model [ 24 , 25 ] , NG and TDVP achieve significant speedup over SGD. (e) To capture rare events, we augment the training procedure with enhanced-sampling strategies, including mixture sampling with configurations drawn from a uniform distribution

P_{U}

, diffusive sampling based on a kernel

\mathcal{K}

, and

\alpha

sampling using an exponent-based modification. In the Schlögl model, these sampling schemes markedly reduce the Hellinger distance to the baseline from expensive Gillespie simulations.

Example	Optimize by	Species	Reactions	$M$	$T_{step}$	$\delta t$	Depth	Width	Batch	Epochs per step	$\eta$	Time (hr)
Toggle switch	SGD	4	8	80	8001	0.005	1	16	2000	100	0.005	3.234
	NG	4	8	80	8001	0.005	1	16	2000	5	0.5	0.806
	TDVP	4	8	80	8001	0.005	1	16	2000	1	1.0	0.185
MAPK cascade	NG	16	35	10	$10^{6}$	0.01	1	8	2000	5	0.5	127.584
1D Schlögl (8 Sites)	NG	8	46	85	$1.25\times 10^{5}$	$8\times 10^{-6}$	1	16	5000	5	0.8	79.433
2D Schlögl (2x4 Sites)	NG	8	52	85	$1.25\times 10^{5}$	$8\times 10^{-6}$	1	16	5000	5	0.8	102.82

Table 1: The models solved by NNCME-2 and the corresponding computational cost under the chosen hyperparameters. The time-step length

\delta t

is measured in units of the inverse reaction rates, and the physical time is

t=\delta t\,T_{\mathrm{step}}

. The upper count limit is

M

. Depth and width refer to the NADE architecture, and

\eta

is the learning rate for optimization algorithms. NG denotes natural gradient and TDVP denotes the time-dependent variational principle. All computations were performed on a single Tesla V100 GPU.

Figure 2: NNCME-2 achieves at least a 5-fold acceleration while maintaining high accuracy, as demonstrated in the genetic toggle switch. (a) A schematic of the chemical reactions. (b) The time evolution of the average counts for the genes and proteins, comparing the VAN (dots), Gillespie simulation (solid line), FSP results from SSIT (dashed), and ACME results (dash-dotted). (c) Comparison of the standard deviations obtained from VAN, FSP (SSIT), and ACME against Gillespie at time points $t=0,2,...,38,40$ . The Gillespie simulation has $1\times 10^{4}$ trajectories. (d) Joint distribution of the two proteins from VAN at time points $t=5,40$ , showing four stable states in the long-time limit. The initial distribution assigns equal probability to the bound and unbound promoter states for each gene, with zero initial protein counts. The results here are from the NADE architecture with natural gradient optimization. Parameters are $s_{x}=s_{y}=50$ , $d_{x}=d_{y}=1$ , $b_{x}=b_{y}=10^{-4}$ , and $u_{x}=u_{y}=0.1$ , with hyperparameters in Table I. (e) Training performance of different optimizers on the NADE architecture. (f) Comparison of computational cost. NADE reduces overall training time compared to RNN, and both NG and TDVP significantly accelerate training.

Figure 3: NNCME-2 handles large reaction networks efficiently. (a) A schematic of a representative large reaction network (MAPK cascade, 16 species). Bidirectional arrows indicate reversible reactions, and the detailed reactions and rates are in Appendix G.2 . (b) The left panel displays the time evolution of the mean counts of the principal species (MKP3, K, and Kpp), comparing Gillespie simulations (lines), VAN (circles), and ACME (squares). The two panels on the right compare the means and standard deviations obtained from VAN and ACME with Gillespie at time points

t=0,1000,\ldots,10000

. (c) Marginal distributions of MKP3, K, and Kpp at

t=2000

and

t=10000

, obtained from the three methods. The inset reports the Hellinger distance

D_{\mathrm{HD}}^{G,V}

between the Gillespie and VAN distributions. (d) Joint distributions of K and Kpp from the VAN at

t=0,2000,

and

10000

. The system evolves from the initial distribution

P(\text{K}=3,\text{MKP3}=1,\text{others}=0)=1

to the bimodal state. The Gillespie results are based on

10^{4}

trajectories. The hyperparameters of the VAN are in Table I.

Figure 4: NNCME-2 captures rare events in the one-dimensional spatially extended Schlögl model (6 sites). (a) A schematic of the reaction-diffusion system, where each site hosts a well-mixed Schlögl model with fixed species $A$ and $B$ , and stochastic species $X_{i}$ diffusing between neighboring sites with rate $d$ . (b) Comparison between the vanilla VAN (top) and the one with enhanced sampling (bottom). The results obtained from the VAN (dots) closely match those from Gillespie simulations (lines) for the time evolution of the mean counts at sites $X_{1}$ , $X_{2}$ , $X_{4}$ , and $X_{6}$ (left), standard deviations at time points $t=0,0.1,...,1.0$ (middle), and the marginal distribution of $X_{6}$ at $t=1$ with the Hellinger distance (right). The inset highlights the low-probability regime ( $X\in[0,20]$ , probability $<5\times 10^{-4}$ ). (c) Joint distributions of $X_{5}$ and $X_{6}$ from the VAN with enhanced sampling, with probabilities shown by color in a logarithmic scale. White boxes indicate the transition of the probability peak. The initial distribution is a delta peak at $X=50$ . The Gillespie results are based on $5\times 10^{5}$ trajectories. The hyperparameters of the VAN are in Table I. Mixture ES is used in this example, with exploratory samples comprising 15% of each batch.

Figure 5: The scaling of rare probabilities and computational time for the one-dimensional Schlögl model. (a) Probability mass

P_{\mathcal{B}}

\mathcal{B}=\{\mathbf{X}\mid 0\leq X_{i}\leq 15,\ \forall i\}

across lattice size

L_{\text{lattice}}

, from vanilla Gillespie simulations (red), NNCME-2 (blue), and Gillespie (SSIT [ 82 ] ) with an optimized SSA implementation (orange). (b) Total computation time (log scale), with a linear-scale view shown in the inset. Both Gillespie implementations used

1\times 10^{6}

trajectories. The original Gillespie simulations were run on an Intel Xeon Gold 6248 CPU node, SSIT calculations were performed in MATLAB on a single-processor MacBook Pro, and NNCME-2 computations were performed on an NVIDIA Tesla V100 GPU node. Enhanced sampling is applied to all lattice sizes using the mixture ES scheme, with 15% exploratory samples per training batch.

Figure 6: Absolute error of the rare-event probability

P_{B}

for the one-dimensional Schlögl model with

L=8

as a function of computation time. Results from the vanilla Gillespie implementation, Gillespie (SSIT), and NNCME-2 are compared with the extrapolated reference

P_{B}^{\mathrm{ref}}=4.75\times 10^{-6}

. The Gillespie and SSIT curves with shaded regions show the variation over 20 resampled estimates. The NNCME-2 points correspond to different training settings. The inset shows the same results on linear axes. The computational platforms follow those described for Fig. 5 .

Figure 7: NNCME-2 captures rare events in the two-dimensional spatially extended Schlögl model ( $2\times 2$ grid). (a) A schematic of the 2d well-mixed Schlögl model. The parameters are the same as Fig. 4 . (b) Comparison between the vanilla VAN (top) and the enhanced VAN (bottom), using mixture enhanced sampling with 15% exploratory samples per training batch. The left and middle panels show the time evolution of the mean counts and the corresponding standard deviations of sites $X_{1}$ - $X_{4}$ , respectively. The right panels display the marginal distribution of $X_{3}$ at $t=1$ , with the Hellinger distance $D_{\mathrm{HD}}$ between VAN and Gillespie results indicated. (c) Joint distributions of $X_{3}$ and $X_{4}$ from the enhanced VAN at $t=1$ , shown on a logarithmic probability scale.

Figure 8: The scaling of rare probabilities and computational time for the two-dimensional spatially extended Schlögl model. (a) Probability mass

P_{\mathcal{B}}

obtained using the vanilla Gillespie implementation, NNCME-2, and Gillespie (SSIT). Both Gillespie implementations use

1\times 10^{6}

trajectories. (b) The total computation times of the three methods on a logarithmic scale, with the inset showing the same data on a linear scale. The computational platforms follow those described for Fig. 5 . All NNCME-2 results use mixture enhanced sampling with 15% exploratory samples per training batch.

Appendix figures & tables15 assets

Supplementary material from the paper’s appendix.

Appendix

Symbol	Description
$\eta$	Learning rate used in optimization.
$\delta t$	Discrete time step of CME propagation.
$T_{\mathrm{step}}$	Number of discrete time steps in the evolution.
$t$	Physical time, $t=\delta t\,T_{\mathrm{step}}$ .
$N_{b}$	Number of samples drawn per training step.
$N_{p}$	Total number of trainable parameters.

Appendix

Table 2: Mathematical symbols used throughout the manuscript.

Figure 9: Effect of network width on the learning dynamics and accuracy of the toggle-switch system. (a) Evolution of the average loss during training of the initial distribution (

t=0

) for different network widths. (b) Evolution of the Hellinger distance over time between the VAN-generated distribution and the Gillespie reference for different network widths. In these experiments, the VAN model employs NADE architecture (

\text{depth}=1

and

N_{b}=2000

). Training utilizes NG and reverse KL, with

\text{epoch}=50

\eta

= 0.5 at

t=0

and

\text{epoch}=5

\eta

= 0.5 for

t>0

Figure 10: Evaluation of KL-divergence loss functions in learning the dynamics of the toggle switch system. (a) Mean of training loss at the initial time step

t=0

. (b) The accuracy of the VAN-generated distribution over time, quantified by the Hellinger distance from the benchmark

10^{4}

Gillespie trajectories. (c) Joint probability heatmap of

P_{x}

and

P_{y}

for the benchmark Gillespie simulation at

t=1

(left) and

t=40

(right). (d) The time evolution of the average count for the genes and proteins, from the VAN (lines) and the Gillespie simulation (dotted line) (left). Standard deviations comparisons (middle). The heatmaps depict the joint distributions of

P_{x}

and

P_{y}

at time points

t=1

and

t=40

(right). In these experiments, the VAN model employs NADE architecture (

\text{width}=16

\text{depth}=1

and

N_{b}=2000

). Training utilizes SGD, with

\text{epoch}=2000

\eta

= 0.005 at

t=0

and

\text{epoch}=100

\eta

= 0.005 for

t>0

Loss Type	Loss
ForwardKL	$\mathcal{L}_{F}=\mathbb{E}_{s\sim\mathbb{T}\hat{P}^{\theta_{t}}_{t}}\big[\ln\mathbb{T}\hat{P}^{\theta_{t}}_{t}(s)-\ln\hat{P}^{\theta_{t+\delta t}}_{t+\delta t}(s)\big]$
ReverseKL	$\mathcal{L}_{R}=\mathbb{E}_{s\sim\hat{P}^{\theta_{t+\delta t}}_{t+\delta t}}\big[\ln\hat{P}^{\theta_{t+\delta t}}_{t+\delta t}(s)-\ln\mathbb{T}\hat{P}^{\theta_{t}}_{t}(s)\big]$
ReverseKL2	$\mathcal{L}_{R_{2}}=\mathbb{E}_{s\sim\mathbb{T}\hat{P}^{\theta_{t}}_{t}}\!\left[\frac{\hat{P}^{\theta_{t+\delta t}}_{t+\delta t}(s)}{\mathbb{T}\hat{P}^{\theta_{t}}_{t}(s)}\big(\ln\hat{P}^{\theta_{t+\delta t}}_{t+\delta t}(s)-\ln\mathbb{T}\hat{P}^{\theta_{t}}_{t}(s)\big)\right]$

Appendix

Table 3: Summary of loss functions.

Figure 11: Validation of the Hellinger distance as a representative evaluation metric for the toggle-switch system. The VAN is trained using different loss functions based on variants of the KL divergence, and model performance is evaluated using multiple distance measures, including Total Variation (TV), Kullback–Leibler (KL), Jensen–Shannon (JS), and Hellinger distance (HD). The reference distribution is based on

10^{4}

trajectories generated via the Gillespie algorithm. In these experiments, the VAN model employs NADE architecture (

\text{width}=16,\text{depth}=1

and

N_{b}=2000

). Training utilizes SGD, with

\text{epoch}=2000

\eta

= 0.005 at

t=0

and

\text{epoch}=100

\eta

= 0.005 for

t>0

Figure 12: Effect of the learning rate

\eta

on NG. (a) Mean of loss at

t=0

for different values of

\eta

. (b) Standard deviation of the loss at

t=0

for different values of

\eta

. The red dashed line marks a threshold of

1\times 10^{-2}

, and the colored vertical lines indicate the first epoch at which the threshold is reached for the corresponding

\eta

. (c) Mean of loss at

t=5

. (d) The accuracy of the VAN-generated distribution for various

\eta

over time, quantified by the Hellinger distance from the

10^{4}

Gillespie trajectories. In these experiments, the VAN model employs NADE architecture (

\text{width}=8,\text{depth}=1

) and is optimized using the reverse KL-divergence. Training utilizes NG, with

\text{epoch}=50

for

\eta\in\{0.05,0.1\}

and

\text{epoch}=5

for

\eta=0.5

, and

N_{b}=2000

Figure 13: Comparison of the training efficiency and convergence behavior of different optimization algorithms for the VAN in the toggle-switch system. (a) Means of training loss at the initial time step $t=0$ for SGD (green) and NG (blue) with 1000 epochs. (b) The accuracy of the VAN-generated distribution using three optimization methods over time, quantified by the Hellinger distance from the benchmark $10^{4}$ Gillespie trajectories. (c) Joint probability heatmap of $P_{x}$ and $P_{y}$ for the benchmark Gillespie simulation at $t=1$ (left) and $t=40$ (right). (d) The left panel shows the time evolution of the mean counts of the genes and proteins obtained from the VAN (lines) and from Gillespie simulations (dotted lines). The middle panel presents a comparison of the corresponding standard deviations. The right panel displays heatmaps of the joint distributions of $P_{x}$ and $P_{y}$ at time points $t=1$ and $t=40$ . All experiments were trained with reverse KL-divergence, the optimizer parameters: learning rate $\eta$ = 0.005, epoch = 100 at $t>0$ , epoch = 2000 at $t=0$ for SGD; $\eta$ = 0.5, epoch = 5 for $t>0$ , epoch = 50 for $t=0$ for NG, in NADE architecture (width $16$ , $N_{b}=2000$ ).

Figure 14: Comparative analysis of training loss for NG and SGD. (a) Means of training loss at the initial time step

t=0

. (b) Means of training loss at an early training stage

t=5

. (c) Standard deviations of training loss at the initial time step

t=0

, the horizontal dotted lines indicate convergence threshold (

5\times 10^{-2}

). (d) Standard deviations of training loss at an early training stage

t=5

, the convergence threshold is set as

10^{-2}

. Red arrows quantify the epoch reduction ratio. All experiments were trained with reverse KL-divergence, the optimizer parameters: learning rate

\eta

= 0.005, epoch = 100 at

t>0

, epoch = 2000 at

t=0

for SGD;

\eta

= 0.5, epoch = 5 for

t>0

, epoch = 50 for

t=0

for NG, in NADE architecture (width

16

N_{b}=2000

Algorithm	Parameter update
Stochastic Gradient Descent	$\Delta\theta=-\eta\dfrac{\partial\mathcal{L}}{\partial\theta}$
Natural Gradient	$\Delta\theta=-\eta S(\theta_{t})^{-1}\dfrac{\partial\mathcal{L}}{\partial\theta}$
Time-Dependent Variational Principle	$\dot{\theta}=S(\theta_{t})^{-1}F$

Appendix

Table 4: Comparison of optimization algorithms and their update mechanisms. SGD follows the negative gradient direction scaled by the learning rate

\eta

; NG incorporates the FIM

S

for preconditioning; TDVP employs a force term

F

derived from the system’s Hamiltonian. All methods iteratively minimize the reverse KL-divergence loss

\mathcal{L}

Figure 15: Comparison of probability distributions for the genetic toggle-switch model obtained from Gillespie simulations, VAN, SSIT-FSP, and ACME. (a) Distribution of the promoter state $G_{x}$ at $t=5$ . (b) Marginal distribution of protein $P_{x}$ at $t=5$ . (c) Joint distributions of $P_{x}$ and $P_{y}$ at $t=1$ and $t=40$ , shown on a logarithmic probability scale.

Figure 16: Joint distributions of K and Kpp for the MAPK cascade at $t=0$ , $2000$ , and $10000$ , obtained from Gillespie simulations, ACME, and VAN. The distribution evolves from the delta-initialized state to a bimodal long-time distribution, with close agreement among the three methods.

Figure 17: Effect of hyperparameter

k

on diffusive ES for the Schlögl (2 sites) model. (a) The time evolution of the average count for

X_{1}

with different

k

from the VAN (dot) and the Gillespie simulation (gray line). (b) The accuracy of the VAN-generated distribution over time, quantified by the Hellinger distance from the benchmark

10^{5}

Gillespie trajectories. (c) Marginal distribution of species

X_{1}

from the VAN at

t=1

for different

k

values (colored) and compared with the Gillespie simulation (gray). The learning parameters for the system are fixed at

T_{\text{step}}=1\times 10^{5}

and

\delta t=1\times 10^{-5}

. All experiments are performed using the NADE architecture with

\text{width}=16

and

\text{depth}=1

. The Gillespie simulation has

10^{5}

trajectories. In these experiments, the VAN model employs NADE architecture (

\text{width}=16

\text{depth}=1

, and

N_{b}=2000

). Training utilizes reverse KL, with

\text{epoch}=50

\eta

= 0.5 at

t=0

using NG and

\text{epoch}=1

\eta

=1 for

t>0

using TDVP.

Figure 18: Effect of hyperparameter

\alpha

\alpha

ES for the Schlögl (2 sites) model. Smaller values of

\alpha

correspond to stronger exploration in the sampling process. The layout of the figure and the hyperparameters of VAN are the same as those of Fig. 17 .

Figure 19: Results of four ES strategies- Vanilla , Mixture ES, Diffusive ES and

\alpha

ES used in training VAN on the Schlögl (2 sites) system. (a) Time evolution of the means of species

X_{1}

estimated by each method. (b) Comparison of the means of

X_{1}

between VAN-based estimates and Gillespie simulation. (c) Comparison of the standard deviation of

X_{1}

between VAN-based estimates and Gillespie simulation. (d) Marginal distributions of species

X_{1}

and

X_{2}

t=1

(first and second rows), and joint probability distribution heatmaps of

X_{1}

and

X_{2}

t=1

(third row). The Gillespie simulation has

10^{5}

trajectories. In these experiments, the VAN model employs NADE architecture (

\text{width}=16

\text{depth}=1

, and

N_{b}=2000

). Training utilizes reverse KL, with

\text{epoch}=50

\eta

= 0.5 at

t=0

using NG and

\text{epoch}=1

\eta

=1 for

t>0

using TDVP.

Figure 20: Comparison of computational efficiency and precision between PyTorch and JAX frameworks, framework implementations are distinguished by color (red: PyTorch, blue: JAX) and line style (solid: TDVP, dashed: NG algorithm). Experiments were conducted using the NADE net (

\text{width}=16,\text{depth}=1

). Performance metrics were evaluated through reverse KL-divergence measurements. Computational efficiency was evaluated through Hellinger distance measurements relative to

10^{4}

Gillespie simulated trajectories. All experiments were executed on an NVIDIA Tesla V100 GPU under identical computational constraints.

Explore similar work

Aug 26, 2026

Data-driven Effective Modeling of Stochastic Chemical Reaction Networks

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.

Yuan Chen, Weize Mao, Dongbin Xiu

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

Surrogate Modeling

Ask PDF

Mar 12, 2026

Chemical Reaction Networks Learn Better than Spiking Neural Networks

We mathematically prove that chemical reaction networks without hidden layers can solve tasks for which spiking neural networks require hidden layers. Our proof uses the deterministic mass-action kinetics formulation of chemical reaction networks. Specifically, we prove that a certain reaction network without hidden layers can learn a classification task previously proved to be achievable by a spiking neural network with hidden layers. We provide analytical regret bounds for the global behavior of the network and analyze its asymptotic behavior and Vapnik-Chervonenkis dimension. In a numerical experiment, we confirm the learning capacity of the proposed chemical reaction network for classifying handwritten digits in pixel images, and we show that it solves the task more accurately and efficiently than a spiking neural network with hidden layers. This provides a motivation for machine learning in chemical computers and a mathematical explanation for how biological cells might exhibit more efficient learning behavior within biochemical reaction networks than neuronal networks.

Sophie Jaffard, Ivo F. Sbalzarini

Max Planck Institute of Molecular Cell Biology and Genetics Center for Systems Biology Dresden Dresden, Germany · Dresden University of Technology, Faculty of Computer Science Max Planck Institute of Molecular Cell Biology and Genetics Center for Systems Biology Dresden Dresden, Germany

Classification Dynamical Systems Spiking Neural Networks

Ask PDF

Oct 13, 2026

Joint Neural Posterior Estimation and Model Selection for Chemical Reaction Networks

Simulation-based Bayesian inference for Chemical Reaction Networks (CRNs) must simultaneously solve two computationally demanding problems: discriminating between overlapping reaction pathways and estimating reaction rate constants. Classical approaches, such as Sequential Monte Carlo ABC, rely on handcrafted summary statistics and scale poorly across large parameter spaces. To overcome this, we introduce a neural inference architecture. Our pipeline utilizes shared 1D-CNNs and a permutation-invariant Set Attention mechanism to learn low-dimensional summary statistics from variable-length, multi-experiment time-series data. These statistics simultaneously condition a model selection network to compute posterior model probabilities and model-specific conditional Invertible Neural Networks (cINNs) that estimate parameter posterior distributions. We demonstrate the practical utility of our approach on a real protodeboronation example, illustrating how it handles structural ambiguity and extracts kinetic parameters from competing catalyst deactivation pathways.

J. Arturo Esquivel, Jorge A. Campos-Gonzalez-Angulo, Sergio Pablo-García +2

Department of Statistical Sciences, University of Toronto, Toronto, Ontario, Canada · Schools of Mathematics & Physics and Chemistry & Molecular Biosciences, The University of Queensland, St Lucia, Queensland, Australia · Australian Research Council Centre of Excellence in Quantum Biotechnology, St Lucia, Queensland, Australia +9

Bayesian Inference Neural Posterior Estimation Simulation-Based Inference Dynamical Systems

Ask PDF