Joint Neural Posterior Estimation and Model Selection for Chemical Reaction Networks

J. Arturo Esquivel; Jorge A. Campos-Gonzalez-Angulo; Sergio Pablo-García; Alán Aspuru-Guzik; Vianey Leos-Barajas

stat.MLOct 13, 2026

Joint Neural Posterior Estimation and Model Selection for Chemical Reaction Networks

Authors: J. Arturo Esquivel, Jorge A. Campos-Gonzalez-Angulo, Sergio Pablo-García, Alán Aspuru-Guzik, Vianey Leos-Barajas

Organizations: 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 · Queensland Quantum Decarbonisation Alliance, St Lucia, Queensland, Australia · The Matter Lab, University of Toronto, Toronto, Ontario, Canada · Institute of Micro and Nanotechnology · Departments of Chemistry, Computer Science, Materials Science & Engineering, Chemical Engineering & Applied Chemistry, and Institute of Medical Science, University of Toronto, Toronto, Ontario, Canada · Acceleration Consortium, Toronto, Ontario, Canada · Vector Institute for Artificial Intelligence, Toronto, Ontario, Canada · Canadian Institute for Advanced Research (CIFAR), Toronto, Ontario, Canada · NVIDIA, Toronto, Ontario, Canada · School of the Environment, University of Toronto, Toronto, Ontario, Canada

Abstract

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.

Figures & tables

Panel of Figure 1 — Figure 1: High-level architecture of the joint neural inference pipeline. Each time-series observation ( $4$ in total) goes through a parallel, shared-weight CNN to extract kinetic features. The architecture bifurcates to evaluate model evidence and condition $\mathcal{M}$ independent model-specific conditional Invertible Neural Networks to estimate parameter posteriors.

Observed Species	Time Points	Global Acc.	$M_{1}$ F1	$M_{2}$ F1	$M_{3}$ F1
Substrate & Product	20	94%	0.94	0.98	0.96
Substrate & Product	10	94%	0.93	0.97	0.93
Product Only	20	93%	0.91	0.97	0.92
Product Only	10	92%	0.89	0.95	0.90

Table 1: Model Selection Performance of Our Deep Learning Pipeline. We evaluate it across varying observational conditions: the impact of hidden chemical species (Product Only) and reduced temporal resolution (10 vs. 20 points). Global accuracy and model-specific F1-scores are reported for the

3,000

independent test trajectories.

Figure 2: Reliability diagram evaluating the calibration of the posterior probabilities across the 3,000 test trajectories. Predictions are grouped into 10 uniformly spaced probability bins. The dashed diagonal represents perfect calibration, where the estimated probability exactly matches empirical accuracy. All three mechanisms closely track the identity line, achieving low Expected Calibration Errors.

Parameter	$M_{1}$	$M_{2}$	$M_{3}$
$k_{1}$	$96.2\%$	$96.2\%$	$96.2\%$
$k_{-1}$	$96.4\%$	$95.2\%$	$96.0\%$
$k_{2}$	$98.1\%$	$97.3\%$	$95.1\%$
$k_{-2}$	$97.0\%$	$94.8\%$	$96.8\%$
$k_{3}$	–	$93.0\%$	$95.2\%$
$k_{-3}$	–	$96.2\%$	$96.0\%$

Table 2: Empirical Coverage of the

95\%

Posterior Credible Intervals. Results for the kinetic parameters and observation noise (

\sigma

) under the

20

-point (Substrate + Product) setting.

Figure 3: Posterior Predictive Check (PPC) for experimental protodeboronation data across three varying initial catalyst loadings. The $300$ predictive trajectories are sampled from the joint posterior, drawing from the candidate mechanisms proportionally to their estimated model probabilities.

Appendix figures & tables20 assets

Supplementary material from the paper’s appendix.

Appendix

Summary Network: CNN-Transformer
Component / Hyperparameter	Specification
1D-CNN filters (3 layers, kernel size 3)	$64\rightarrow 128\rightarrow 256$
Pooling Layer	Global Average Pooling (across time steps, masked for dropout points)
Dense layer ( $d_{model}$ )	$256$
Self-Attention heads	$4$
Feed Forward Network	$512\rightarrow 256$

Appendix

Table 3: Architecture and Training Hyperparameters for our Neural Inference Framework.

Panel of Figure 4 — Figure 4: Detailed architecture of the summary network $T_{\phi}$ . Independent experimental observations are processed through a shared 1D-CNN (indicated by dashed connections) to extract temporal signatures. The sequences are globally average pooled across time to ensure length invariance. The independent features are contextualized using a Set Attention Block (SAB), before being globally pooled across experiments and passed through dense layers to produce the unified embedding $\mathbf{e}$ .

Figure 5: Convergence diagnostics for SMC-ABC inference. The solid lines (left axis) track the shrinkage of the discrepancy threshold (

\epsilon

) across generations for candidate mechanisms M1, M2, and M3. The dashed lines (right axis) denote the corresponding acceptance rates required to maintain the target population of

N=3,000

particles per generation. As the posterior approximation refines and the

\epsilon

tolerance tightens, the acceptance rate exhibits a rapid decline, reflecting the increasing computational cost of proposing valid kinetic parameters from the transition kernel.

Figure 6: PPC for the four experimental observations, where SMC-ABC was implemented. Black lines denote the true observed synthetic data, with solid and dashed lines representing the substrate (S) and product (P) concentrations, respectively. Faint colored lines depict the 100 inferred trajectories generated by drawing from the network’s estimated posterior.

		Predicted Model			Performance Metrics
Setup	Truth	$M_{1}$	$M_{2}$	$M_{3}$	Precision	Recall	$F_{1}$ -score
(1) S & P, $T=20$	$M_{1}$	954	9	37	0.93	0.95	0.94
	$M_{2}$	18	976	6	0.99	0.98	0.98
	$M_{3}$	51	5	944	0.96	0.94	0.95
(2) S & P, $T=10$	$M_{1}$	936	17	47	0.92	0.94	0.93
	$M_{2}$	19	965	16	0.98	0.96	0.97

Appendix

Table 4: Classification Performance. Results for the three-model setting (

\mathcal{M}=\{M_{1},M_{2},M_{3}\}

) across all four observational setups. The left pane displays the confusion matrix across the holdout test set, and the right pane details the corresponding performance metrics.

	(1) S & P, $T=20$			(2) S & P, $T=10$			(3) P only, $T=20$			(4) P only, $T=10$
Param	$M_{1}$	$M_{2}$	$M_{3}$	$M_{1}$	$M_{2}$	$M_{3}$	$M_{1}$	$M_{2}$	$M_{3}$	$M_{1}$	$M_{2}$	$M_{3}$
$k_{1}$	96.2	96.2	96.2	95.0	96.2	96.0	97.7	94.9	98.0	95.9	92.6	97.4
$k_{-1}$	96.4	95.2	96.0	97.0	94.9	96.4	95.8	94.4	94.9	96.5	94.8	94.0
$k_{2}$	98.1	97.3	95.1	97.5	94.3	95.4	97.8	95.2	93.9	96.5	94.6	93.9
$k_{-2}$	97.0	94.8	96.8	97.7	94.8	97.3	97.0	95.4	94.3	96.6	94.6	94.8
$k_{3}$	–	93.0	95.2	–	89.5	93.1	–	97.7	94.0	–	96.5	93.4

Appendix

Table 6: Empirical Coverage for Three Models. Coverage (%) for the

95\%

credible intervals for the kinetic parameters and observation noise across all four observational setups. Values close to

95\%

indicate well-calibrated marginal posteriors.

Figure 8: Estimated posterior distributions, using

2,000

samples, for all parameters from a test case where the true mechanism was

M_{3}

. The true value (in the log-space) is marked with black dashed line.

Figure 9: Inference calibration across all testing cases from

M_{1}

. Left Panel: Posterior mean estimates vs. the true values. Right Panel: Rank histograms for all parameters (normalized to

(0,1)

). The dashed black line indicated the level for uniform density (

1

Figure 10: Inference calibration across all testing cases from

M_{2}

. Left Panel: Posterior mean estimates vs. the true values. Right Panel: Rank histograms for all parameters (normalized to

(0,1)

). The dashed black line indicated the level for uniform density (

1

Figure 11: Inference calibration across all testing cases from

M_{3}

. Left Panel: Posterior mean estimates vs. the true values. Right Panel: Rank histograms for all parameters (normalized to

(0,1)

). The dashed black line indicated the level for uniform density (

1

Figure 12: PPC for a test trajectory generated by mechanism $M_{3}$ under the high-resolution setup (1). Black dots denote the $20$ observations for both species across four distinct initial conditions. The solid lines represent $100$ posterior predictive trajectories sampled from the network’s inferred joint posterior.

Figure 13: PPC for a test trajectory generated by mechanism $M_{1}$ under the low-resolution setup (4). Black dots denote the $10$ observations for the product across four distinct initial conditions. The gray dots indicate the corresponding substrate, which was not observed by the model. The solid lines represent $100$ posterior predictive trajectories sampled from the network’s inferred joint posterior.

	Predicted Model				Performance Metrics
Truth	$M_{1}$	$M_{2}$	$M_{3}$	$M_{5}$	Precision	Recall	$F_{1}$ -score
$M_{1}$	789	6	19	186	0.75	0.79	0.77
$M_{2}$	17	977	6	0	0.98	0.98	0.98
$M_{3}$	54	9	928	9	0.93	0.93	0.93
$M_{5}$	190	4	43	763	0.80	0.76	0.78

Appendix

Table 7: Classification Performance for Four Models. Results for the expanded four-model setting (

\mathcal{M}=\{M_{1},M_{2},M_{3},M_{5}\}

) under the main observational setup (both species observed,

T=20

). The left pane displays the confusion matrix across the holdout test set, and the right pane details the corresponding performance metrics.

Figure 15: PPC for a test trajectory generated by mechanism $M_{1}$ under the high-resolution setup (1) with four candidate models. Black dots denote the $20$ discrete observations for both species across four distinct initial conditions. The solid lines represent $100$ posterior predictive trajectories sampled from the network’s inferred joint posterior.

Figure 16: Inference calibration across all testing cases from $M_{5}$ . Left Panel: Posterior mean estimates vs. the true values. Right Panel: Rank histograms for all parameters (normalized to $(0,1)$ ). The dashed black line indicated the level for uniform density ( $1$ ).

	$M_{14}$ (3.71%)		$M_{15}$ (92.19%)		$M_{16}$ (4.10%)
Parameter	Mean	Median	Mean	Median	Mean	Median
$\log k_{1}$	$3.9572$	$2.2728$	$3.5822$	$1.0382$	$4.2962$	$2.6632$
$\log k_{-1}$	$-5.6564$	$-6.7500$	$-6.1921$	$-7.2830$	$-5.8420$	$-6.8342$
$\log k_{2}$	$-1.8568$	$-3.9131$	$-1.8872$	$-2.5498$	$-2.0608$	$-4.0818$
$\log k_{-2}$	$-0.9392$	$-1.4539$	$-0.2301$	$-1.2468$	$-1.0966$	$-1.8808$
$\log k_{-3}$	$-2.1533$	$-4.2273$	$9.9645$	$10.5861$	$2.2611$	$0.2867$

Appendix

Table 9: Parameter Estimates. Posterior parameter estimates for the protodeboronation structural hypotheses.

Figure 17: Estimated posterior distributions for

M_{15}

, which was estimated to be the most likely mechanism from our candidate models. We use

2,000

samples, for all parameters estimated from real protodeboronation data. The estimated posterior mean (in the log-space) is marked with black dashed line.

Figure 18: Bivariate and marginal posterior distributions (

2,000

samples) for the most probable protodeboronation mechanism. The diagonal panels display the marginal posterior distributions. The lower-triangular panels feature bivariate KDE contour maps, illustrating the pairwise uncertainty between the kinetic constants.

Explore similar work

Jun 26, 2026

Reduction of Probabilistic Chemical Reaction Networks

Programming adaptive behaviors at the cellular level is a long-standing goal that raises the question of how probabilistic computation can be implemented in biochemical systems. Chemical reaction networks (CRNs) provide such a substrate and have been shown to realize probabilistic models, including hidden Markov models and factor graphs, with dynamics reproducing Bayesian inference and belief propagation. However, encoding these algorithms typically requires prohibitively large reaction networks, and classical CRN reduction techniques do not directly apply. By recovering the factor graph structure encoded in Napp--Adams-compiled CRNs, we transport recent factor-graph reduction results to their chemical implementations, obtaining significantly smaller CRNs while preserving the belief-propagation fixed points on surviving variables.

Mauricio Montes, Gregoire Sergeant-Perthuis

COSAM, Auburn University, Auburn, AL, USA · CQSB, Sorbonne Universite, Paris, France

Probabilistic Graphical Models Reduced-Order Modeling

Ask PDF

May 27, 2025

Causal Posterior Estimation

We present Causal Posterior Estimation (CPE), a novel method for Bayesian inference in simulator models, where evaluating the likelihood function is intractable or computationally expensive, but generating outputs given parameter values is straightforward. CPE approximates the posterior distribution using flow matching while directly incorporating the conditional dependence structure induced by the model's graphical representation into the neural network architecture. Across extensive experiments, we demonstrate that hard-coding these conditional dependencies into the network, rather than requiring them to be learned from data, enables CPE to achieve highly accurate posterior inference that matches or outperforms state-of-the-art baselines.

Simon Dirmeier, Antonietta Mira

Swiss Data Science Center, ETH Zurich, Switzerland · Università della Svizzera italiana, Switzerland · Insubria University, Italy

Bayesian Inference Simulation-Based Inference Neural Posterior Estimation

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