Synergizing Physically Constrained MCMC and Chemical-Informed Gaussian Processes for Reaction Network Discovery

# Synergizing Physically Constrained MCMC and Chemical-Informed Gaussian Processes for Reaction Network Discovery
Source: [https://arxiv.org/html/2606.23757](https://arxiv.org/html/2606.23757)
###### Abstract

Extracting interpretable governing equations from sparse, noisy chemical time\-series data remains difficult because discrete reaction topology and continuous kinetic parameters are tightly coupled\. We present PC\-MCMC\-CIGP, a reproducible gray\-box workflow that combines spike\-and\-slab topology sampling, hard conservation and thermodynamic screening, and a Chemical\-Informed Gaussian Process \(CIGP\) residual model for parameter calibration and experimental design\. The methodological contribution is not a new MCMC or GP family in isolation; rather, it is the integration of these components into a physically constrained workflow with explicit uncertainty\-aware acquisition choices\. On theH2\+Br2\\mathrm\{H\_\{2\}\+Br\_\{2\}\}benchmark, the constrained sampler distinguishes elementary radical pathways from deceptive phenomenological fits in our experiments\. On styrene epoxidation, the CIGP optimization loop improves final yield by12\.5%12\.5\\%over the reported GP\-BO baseline\. A new 10\-seed acquisition study shows that EI, GWU, PC\-EI, uncertainty sampling, discrepancy hunting, and random search have different trade\-offs: PC\-EI substantially reduces low\-yield BO suggestions, while EI\-style criteria give the strongest final\-yield performance\.

reaction network discovery, Bayesian inference, Gaussian processes, experimental design

## 1Introduction

The discovery of governing equations from noisy time\-series data remains a fundamental challenge in scientific machine learning\(Bruntonet al\.,[2016](https://arxiv.org/html/2606.23757#bib.bib1); Schmidt and Lipson,[2009](https://arxiv.org/html/2606.23757#bib.bib3); Championet al\.,[2019](https://arxiv.org/html/2606.23757#bib.bib4)\)\. Although deep learning paradigms approximate continuous dynamics effectively, interpretability is often sacrificed for black\-box interpolation\(Cranmer,[2023](https://arxiv.org/html/2606.23757#bib.bib5)\)\. In chemical kinetics, the primary objective is*structure identification*, defined as the reconstruction of the discrete reaction topology governing system evolution to recover a physically interpretable causal graph\(Burnhamet al\.,[2008](https://arxiv.org/html/2606.23757#bib.bib11); Jianget al\.,[2022](https://arxiv.org/html/2606.23757#bib.bib15)\)\.

Mechanism discovery is mathematically formulated as an ill\-posed inverse problem characterized by the strong coupling between discrete structures and continuous kinetic constants\(Walter and Pronzato,[1997](https://arxiv.org/html/2606.23757#bib.bib19); Davidescu and Jørgensen,[2008](https://arxiv.org/html/2606.23757#bib.bib16); Audolyet al\.,[2002](https://arxiv.org/html/2606.23757#bib.bib18)\)\. Conventional Bayesian inference is frequently hindered by this interdependence\. Standard parameter estimation techniques lack explicit physical encoding\(Williams and Rasmussen,[2006](https://arxiv.org/html/2606.23757#bib.bib31)\), whereas structure sampling is often impeded by vast combinatorial spaces and insufficient physical constraints, leading to mathematically valid yet physically prohibited solutions\(Encisoet al\.,[2021](https://arxiv.org/html/2606.23757#bib.bib51)\)\.

To address these limitations, a unified physically constrained Bayesian framework is proposed to synergize discrete sampling on physical manifolds with chemically informed continuous modeling\(Brubakeret al\.,[2012](https://arxiv.org/html/2606.23757#bib.bib44)\)\. The problem is treated as the cooperative solution of structural inference and parameter estimation\. A multi\-physically constrained Spike\-and\-Slab Markov Chain Monte Carlo \(MCMC\) sampler is introduced for structural identification\. Constraints regarding mass conservation, electron conservation, and thermodynamic potential loops are encoded as boundary conditions to prune the posterior distribution and eliminate spurious mechanisms\.

For parameter estimation, a*Chemical\-Informed Gaussian Process*\(CIGP\) is utilized to construct a hybrid gray\-box model\. The mechanism\-based ordinary differential equation is embedded as the prior mean function to capture dominant kinetic trends, while structural discrepancies are modeled via non\-parametric kernels\(Raissiet al\.,[2017](https://arxiv.org/html/2606.23757#bib.bib33); Chang and Zeng,[2023](https://arxiv.org/html/2606.23757#bib.bib38); Maet al\.,[2020](https://arxiv.org/html/2606.23757#bib.bib34)\)\. Physical parameters and kernel hyperparameters are jointly optimized to ensure robust estimation\(Crosset al\.,[2024](https://arxiv.org/html/2606.23757#bib.bib37); Daltonet al\.,[2024](https://arxiv.org/html/2606.23757#bib.bib53)\)\. Additionally, physics\-aware acquisition functions are developed to facilitate automatic experimental design\. Figure effig:framework summarizes the two\-stage workflow: PC\-MCMC first identifies a physically admissible reaction topology, and CIGP then uses the identified ODE model as a GP mean function during calibration and active learning\.

The main contributions of this work are summarized as follows:

1. 1\.A fixed\-dimensional spike\-and\-slab MCMC workflow is formulated for candidate reaction selection, with mass, charge, and detailed\-balance constraints used as hard admissibility checks\.
2. 2\.A Chemical\-Informed Gaussian Process calibration stage is used as a gray\-box residual model, preserving the mechanistic ODE as the prior mean while modeling systematic discrepancy\.
3. 3\.A reproducible active\-learning benchmark is provided, including a 10\-seed comparison of PC\-EI, EI, GWU, discrepancy hunting, uncertainty sampling, and random search\.

![Refer to caption](https://arxiv.org/html/2606.23757v1/fig1.jpg)Figure 1:Overview of the proposed PC\-MCMC\-CIGP framework\. The upper panel illustrates the physically constrained structure discovery stage, where candidate reaction networks are sampled using a Spike\-and\-Slab MCMC scheme subject to mass conservation and thermodynamic constraints\. The lower panel depicts the Chemical\-Informed Gaussian Process \(CIGP\), which embeds the discovered mechanistic ODE model as the GP mean function while modeling structural discrepancies non\-parametrically\. The unified framework enables simultaneous mechanism identification, parameter inference, and physics\-aware active learning\.
## 2Related Work

### 2\.1Phenomenological versus Mechanistic Discovery

Data\-driven discovery of kinetic laws is commonly approached through symbolic regression or sparse identification techniques\(Bruntonet al\.,[2016](https://arxiv.org/html/2606.23757#bib.bib1); Rudyet al\.,[2017](https://arxiv.org/html/2606.23757#bib.bib2); Schmidt and Lipson,[2009](https://arxiv.org/html/2606.23757#bib.bib3); Burlacuet al\.,[2020](https://arxiv.org/html/2606.23757#bib.bib6)\)\. Algorithms such as genetic programming and sparse regression can recover accurate macroscopic expressions from observational data, including fractional\-order kinetics that arise in radical chain reactions\(Cranmer,[2023](https://arxiv.org/html/2606.23757#bib.bib5); Otte,[2014](https://arxiv.org/html/2606.23757#bib.bib52)\)\. Their limitation for this work is that an accurate expression is not necessarily a reaction mechanism: a half\-order rate law may fitHBr\\mathrm\{HBr\}formation while leaving the elementary radical topology and hidden intermediates unidentified\.

### 2\.2Reaction\-Network and Neural\-ODE Models

Reaction network reconstruction has also been formulated as deterministic superstructure selection using MINLP, flux\-balance\-style constraints, or hypergraph search\(Langary and Nikoloski,[2019](https://arxiv.org/html/2606.23757#bib.bib13); Searsonet al\.,[2007](https://arxiv.org/html/2606.23757#bib.bib14); Willis and von Stosch,[2016](https://arxiv.org/html/2606.23757#bib.bib12); Palet al\.,[2025](https://arxiv.org/html/2606.23757#bib.bib29); Bonvin and Rippin,[1990](https://arxiv.org/html/2606.23757#bib.bib20)\)\. These methods can encode stoichiometry explicitly, but they often provide limited posterior uncertainty over alternative topologies\. Recent neural ODE and chemical\-reaction neural\-network models improve flexible dynamics learning for kinetic systems\(Kimet al\.,[2021](https://arxiv.org/html/2606.23757#bib.bib54); Owoyele and Pal,[2021](https://arxiv.org/html/2606.23757#bib.bib55); Chang and Zeng,[2023](https://arxiv.org/html/2606.23757#bib.bib38)\)\. They are complementary to our objective: rather than learning an unconstrained neural vector field, PC\-MCMC\-CIGP assumes a candidate elementary\-reaction set and estimates which reactions are active, so the output remains an interpretable subset of reaction steps\.

### 2\.3Probabilistic Modeling and Experimental Design

Gaussian Processes \(GPs\) are widely used as probabilistic surrogates for calibration and Bayesian optimization\(Williams and Rasmussen,[2006](https://arxiv.org/html/2606.23757#bib.bib31); Kocijan,[2016](https://arxiv.org/html/2606.23757#bib.bib35); Snoeket al\.,[2012](https://arxiv.org/html/2606.23757#bib.bib49); Shieldset al\.,[2021](https://arxiv.org/html/2606.23757#bib.bib30)\)\. Standard GP\-BO is effective for yield maximization, but it usually treats the experiment as a black\-box response surface\. Physics\-informed and hybrid GP models incorporate mechanistic structure into the mean or kernel\(Kennedy and O’Hagan,[2001](https://arxiv.org/html/2606.23757#bib.bib32); Raissiet al\.,[2017](https://arxiv.org/html/2606.23757#bib.bib33); Maet al\.,[2020](https://arxiv.org/html/2606.23757#bib.bib34); Crosset al\.,[2024](https://arxiv.org/html/2606.23757#bib.bib37)\)\. Our CIGP component follows this gray\-box tradition\. The novelty claimed here is therefore deliberately scoped: we provide a unified, physically constrained reaction\-discovery and active\-learning workflow, not a fundamentally new GP inference algorithm\.

## 3The PC\-MCMC\-CIGP Framework

### 3\.1Problem Formulation

We consider a chemical system withNsN\_\{s\}species and a pre\-enumerated candidate set ofNrN\_\{r\}elementary reactionsℛ=\{Rj\}j=1Nr\\mathcal\{R\}=\\\{R\_\{j\}\\\}\_\{j=1\}^\{N\_\{r\}\}\. This candidate set is assumed known; the unknowns are the binary activity vector𝜸∈\{0,1\}Nr\\boldsymbol\{\\gamma\}\\in\\\{0,1\\\}^\{N\_\{r\}\}, the active reaction graph𝒢​\(𝜸\)\\mathcal\{G\}\(\\boldsymbol\{\\gamma\}\), and kinetic parameters𝐤\\mathbf\{k\}\. Let𝐜​\(t\)∈ℝ≥0Ns\\mathbf\{c\}\(t\)\\in\\mathbb\{R\}\_\{\\geq 0\}^\{N\_\{s\}\}denote concentrations and let𝐒​\(𝒢\)∈ℤNs×Nr\\mathbf\{S\}\(\\mathcal\{G\}\)\\in\\mathbb\{Z\}^\{N\_\{s\}\\times N\_\{r\}\}be the stoichiometric matrix whosejj\-th column is the net stoichiometric vector of reactionRjR\_\{j\}and is masked byγj\\gamma\_\{j\}during simulation\. The mass\-action dynamics are

d​𝐜d​t=𝐒​\(𝒢\)​𝐫​\(𝐜​\(t\);𝐤,𝜸\),\\frac\{d\\mathbf\{c\}\}\{dt\}=\\mathbf\{S\}\(\\mathcal\{G\}\)\\,\\mathbf\{r\}\\bigl\(\\mathbf\{c\}\(t\);\\mathbf\{k\},\\boldsymbol\{\\gamma\}\\bigr\),\(1\)whererj​\(𝐜;𝐤,𝜸\)=γj​kj​∏iciνi​j−r\_\{j\}\(\\mathbf\{c\};\\mathbf\{k\},\\boldsymbol\{\\gamma\}\)=\\gamma\_\{j\}k\_\{j\}\\prod\_\{i\}c\_\{i\}^\{\\nu\_\{ij\}^\{\-\}\}for reactant stoichiometric ordersνi​j−\\nu\_\{ij\}^\{\-\}\. Reaction rates follow either Arrhenius parameterizationk=A​exp⁡\(−Ea/R​T\)k=A\\exp\(\-E\_\{a\}/RT\)or the detailed\-balance parameterization described below\(Horn and Jackson,[1972](https://arxiv.org/html/2606.23757#bib.bib21); Vlad and Ross,[2009](https://arxiv.org/html/2606.23757#bib.bib25)\)\.

The mechanism\-discovery data are sparse noisy trajectories𝒟=\{\(tm,𝐲m\)\}m=1M\\mathcal\{D\}=\\\{\(t\_\{m\},\\mathbf\{y\}\_\{m\}\)\\\}\_\{m=1\}^\{M\}, where𝐲m\\mathbf\{y\}\_\{m\}contains observed components of𝐜​\(tm\)\\mathbf\{c\}\(t\_\{m\}\)under a Gaussian observation model\. The active\-learning stage uses a distinct notation:𝐮∈𝒰\\mathbf\{u\}\\in\\mathcal\{U\}denotes controllable experimental conditions such as initial concentrations, temperature, and residence time, andy​\(𝐮\)y\(\\mathbf\{u\}\)denotes the scalar objective such as product yield\. This separates state variables𝐜​\(t\)\\mathbf\{c\}\(t\)from design variables𝐮\\mathbf\{u\}\.

### 3\.2Physically Constrained Topology Search

To mitigate the combinatorial explosion inherent in reaction network discovery, we formulate a Bayesian variable selection problem on a physically admissible manifold\(Mitchell and Beauchamp,[1988](https://arxiv.org/html/2606.23757#bib.bib45); George and McCulloch,[1993](https://arxiv.org/html/2606.23757#bib.bib46)\)\. The existence of thejj\-th candidate reaction is controlled by a binary latent variableγj∈\{0,1\}\\gamma\_\{j\}\\in\\\{0,1\\\}\.

Effective Rate Formulation\.To ensure rigorous structural sparsity consistent with the governing kinetics, we define the effective rate constant as a masked variable:

kje​f​f=γj⋅exp⁡\(θj\),k\_\{j\}^\{eff\}=\\gamma\_\{j\}\\cdot\\exp\(\\theta\_\{j\}\),\(2\)whereθj\\theta\_\{j\}is the auxiliary logarithmic kinetic parameter\. This multiplicative formulation ensures that when a reaction is deactivated \(γj=0\\gamma\_\{j\}=0\), its contribution to the differential equations \(Eq\.[1](https://arxiv.org/html/2606.23757#S3.E1)\) vanishes exactly \(kje​f​f≡0k\_\{j\}^\{eff\}\\equiv 0\), regardless of the value ofθj\\theta\_\{j\}\.

Prior Specification\.To induce strict sparsity, a Spike\-and\-Slab prior is imposed on the auxiliary parameterθj\\theta\_\{j\}to regularize the search space:

p​\(θj∣γj\)=\(1−γj\)​𝒩​\(θj;0,ϵ2\)\+γj​𝒩​\(θj;μθ,σθ2\),p\(\\theta\_\{j\}\\mid\\gamma\_\{j\}\)=\(1\-\\gamma\_\{j\}\)\\mathcal\{N\}\(\\theta\_\{j\};0,\\epsilon^\{2\}\)\+\\gamma\_\{j\}\\mathcal\{N\}\(\\theta\_\{j\};\\mu\_\{\\theta\},\\sigma\_\{\\theta\}^\{2\}\),\(3\)where the spike component \(approximated by a narrow Gaussian withϵ→0\\epsilon\\to 0\) constrains the auxiliary parameters of inactive pathways for identifiability, while the slab component captures the parameter uncertainty of active reactions\(Ishwaran and Rao,[2005](https://arxiv.org/html/2606.23757#bib.bib47); Piironen and Vehtari,[2017](https://arxiv.org/html/2606.23757#bib.bib10)\)\. Posterior inference over𝜸\\boldsymbol\{\\gamma\}and𝜽\\boldsymbol\{\\theta\}is performed using a Metropolis–Hastings Markov Chain Monte Carlo \(MCMC\) sampler\(Metropoliset al\.,[1953](https://arxiv.org/html/2606.23757#bib.bib41); Peskun,[1973](https://arxiv.org/html/2606.23757#bib.bib42); Girolami and Calderhead,[2011](https://arxiv.org/html/2606.23757#bib.bib43)\)\.

#### Physical Constraints Enforcement\.

To eliminate mathematically admissible but physically prohibited solutions, three classes of hard constraints are imposed via rejection sampling\.

Mass and electron conservation\.Let𝐀∈ℝNa×Ns\\mathbf\{A\}\\in\\mathbb\{R\}^\{N\_\{a\}\\times N\_\{s\}\}denote the atomic composition matrix\. A valid reactionjjmust satisfy

𝐀​𝐒⋅j=𝟎,\\mathbf\{A\}\\,\\mathbf\{S\}\_\{\\cdot j\}=\\mathbf\{0\},\(4\)ensuring strict conservation of atomic species and charge\.

Thermodynamic consistency\.For each detected forward–reverse pair, detailed balance constrains the rate ratio through latent dimensionless chemical potentials𝝁\\boldsymbol\{\\mu\}:

ln⁡kfwdkrev=−∑i=1Nsνi​j​μi,\\ln\\frac\{k\_\{\\mathrm\{fwd\}\}\}\{k\_\{\\mathrm\{rev\}\}\}=\-\\sum\_\{i=1\}^\{N\_\{s\}\}\\nu\_\{ij\}\\mu\_\{i\},\(5\)whereνi​j\\nu\_\{ij\}is the net stoichiometric coefficient of speciesiifor the forward reaction\. The potentials\{μi\}\\\{\\mu\_\{i\}\\\}are inferred jointly with kinetic parameters under broad box priors, which avoids requiring tabulated thermodynamic data for transient radicals\. Candidate moves that violate mass, charge, rate bounds, or detailed\-balance constraints are rejected before likelihood evaluation\. This rejection\-based implementation is simple but may reduce acceptance rates as the number of reversible pairs grows; we report this as a limitation rather than assuming scalability is automatic\.

### 3\.3Chemical\-Informed Gaussian Processes

To decouple parameter estimation from numerical integration error and systematic model mismatch, we adopt a hybrid gray\-box modeling strategy\(Kristensenet al\.,[2004](https://arxiv.org/html/2606.23757#bib.bib40)\)\. Let𝐮=\[𝐜0⊤,T,t\]⊤∈𝒰\\mathbf\{u\}=\[\\mathbf\{c\}\_\{0\}^\{\\top\},T,t\]^\{\\top\}\\in\\mathcal\{U\}denote the design vector comprising initial concentrations, temperature, and reaction time\. The observation model decomposes the system response into a mechanistic component, a structural discrepancy term, and observation noise:

𝐲​\(𝐮\)=fphys​\(𝐮;𝐖\)\+gerr​\(𝐮\)\+ϵ,\\mathbf\{y\}\(\\mathbf\{u\}\)=f\_\{\\mathrm\{phys\}\}\(\\mathbf\{u\};\\mathbf\{W\}\)\+g\_\{\\mathrm\{err\}\}\(\\mathbf\{u\}\)\+\\boldsymbol\{\\epsilon\},\(6\)wherefphysf\_\{\\mathrm\{phys\}\}denotes the solution of the ODE system defined by the inferred reaction topology,gerr​\(⋅\)g\_\{\\mathrm\{err\}\}\(\\cdot\)is a zero\-mean Gaussian Process with a radial basis function kernel, andϵ∼𝒩​\(𝟎,σn2​𝐈\)\\boldsymbol\{\\epsilon\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\sigma\_\{n\}^\{2\}\\mathbf\{I\}\)represents i\.i\.d\. measurement noise\.

#### Semi\-Parametric Inference\.

The physical model is embedded as the mean function of the Gaussian Process\. For a test design𝐮∗\\mathbf\{u\}\_\{\*\}, the predictive distribution is Gaussian with mean:

μpred=\\displaystyle\\mu\_\{\\mathrm\{pred\}\}=\{\}fphys​\(𝐮∗;𝐖\)\\displaystyle f\_\{\\mathrm\{phys\}\}\(\\mathbf\{u\}\_\{\*\};\\mathbf\{W\}\)\(7\)\+𝐤∗⊤​\(𝐊\+σn2​𝐈\)−1\\displaystyle\+\\mathbf\{k\}\_\{\*\}^\{\\top\}\(\\mathbf\{K\}\+\\sigma\_\{n\}^\{2\}\\mathbf\{I\}\)^\{\-1\}\(𝐲−fphys​\(𝐗;𝐖\)\)\.\\displaystyle\\quad\\bigl\(\\mathbf\{y\}\-f\_\{\\mathrm\{phys\}\}\(\\mathbf\{X\};\\mathbf\{W\}\)\\bigr\)\.To ensure the physical interpretability of the hybrid model, we strictly control the complexity of the non\-parametric component\. The physical parameters𝐖\\mathbf\{W\}and GP hyperparameters𝚯gp\\boldsymbol\{\\Theta\}\_\{\\mathrm\{gp\}\}are jointly optimized by solving a constrained maximization problem:

max𝐖,𝚯gplog⁡p​\(𝐲∣𝐗,𝐖,𝚯gp\)−λ​σf2\\displaystyle\\max\_\{\\mathbf\{W\},\\boldsymbol\{\\Theta\}\_\{\\mathrm\{gp\}\}\}\\quad\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{W\},\\boldsymbol\{\\Theta\}\_\{\\mathrm\{gp\}\}\)\-\\lambda\\sigma\_\{f\}^\{2\}\(8\)s\.t\.0<σf2≤ϵt​o​l,\\displaystyle\\text\{s\.t\.\}\\quad 0<\\sigma\_\{f\}^\{2\}\\leq\\epsilon\_\{tol\},whereσf2\\sigma\_\{f\}^\{2\}is the kernel signal variance\. This formulation combines a soft sparsity penalty \(λ\\lambda\) with hard safety bounds \(ϵt​o​l\\epsilon\_\{tol\}\)\(Hewinget al\.,[2019](https://arxiv.org/html/2606.23757#bib.bib39)\), compelling the Gaussian Process to capture only significant structural residuals that cannot be explained by the physical mechanism\.

![Refer to caption](https://arxiv.org/html/2606.23757v1/fig2.jpg)Figure 2:Geometric intuition of physics\-aware acquisition strategies\. The landscape illustrates the sampling objectives of different acquisition functionsα​\(𝐮\)\\alpha\(\\mathbf\{u\}\)over a noisy reaction yield curve\.\(1\) GWU \(Gradient\-Weighted Uncertainty\):Targets regions of maximal physical sensitivity \(teal arrow,‖∇f‖2\\\|\\nabla f\\\|\_\{2\}\) to maximize information gain about kinetic parameters, typically near the steepest gradient\.\(2\) EI \(Expected Improvement\):Targets the global optimum \(diamond marker\) for direct yield maximization based on the predictive mean\.\(3\) DH \(Discrepancy Hunter\):Identifies regions of significant structural model mismatch \(vertical lines,\|μerr\|\|\\mu\_\{\\mathrm\{err\}\}\|\) to correct systematic bias\.\(4\) PC\-EI \(Physically Constrained EI\):Restricts the search space to physically valid regions, where the infeasible domain \(grey shaded area\) is weighted by the soft feasibility gateσs​i​g​m​o​i​d\\sigma\_\{sigmoid\}\.![Refer to caption](https://arxiv.org/html/2606.23757v1/fig3.jpg)Figure 3:Evaluation of structural identifiability and predictive stability on theH2\+Br2\\mathrm\{H\_\{2\}\+Br\_\{2\}\}benchmark\.\(a\) Posterior Inclusion Probabilities \(PIP\):Derived from the PC\-MCMC sampler, showing high posterior support for the reported elementary steps and low support for spurious pathways\.\(b\) Trajectory Reconstruction:Predicted concentrations by the proposed PC\-MCMC framework align robustly with noisy experimental observations across all chemical species\.\(c\) PySR Parity Plot:Reaction rates predicted by symbolic regression exhibit high statistical goodness\-of\-fit, despite the model being purely phenomenological\.\(d\) SINDy Trajectory Prediction:Unconstrained sparse regression leads to immediate numerical divergence due to violation of non\-negativity constraints, illustrating instability\.\(e\) Comparative Summary:Overview of inferred model representations, physical validity, and outcome classification, highlighting the distinction between true mechanism discovery \(Ours\) and curve fitting \(Baselines\)\.

### 3\.4Active Learning Strategy

To guide autonomous experimental design, four physics\-aware acquisition functionsα​\(𝐮\)\\alpha\(\\mathbf\{u\}\)are constructed based on the CIGP posterior\.

#### Expected Improvement \(EI\)

is used as a baseline to balance exploration and exploitation in a black\-box optimization setting\(Joneset al\.,[1998](https://arxiv.org/html/2606.23757#bib.bib48); Shahriariet al\.,[2015](https://arxiv.org/html/2606.23757#bib.bib50)\)\.

#### Gradient\-Weighted Uncertainty \(GWU\)

prioritizes regions that are simultaneously uncertain and highly informative for physical parameters:

αGWU​\(𝐮\)=σpred​\(𝐮\)​‖∇𝐖fphys​\(𝐮;𝐖\)‖2\.\\alpha\_\{\\mathrm\{GWU\}\}\(\\mathbf\{u\}\)=\\sigma\_\{\\mathrm\{pred\}\}\(\\mathbf\{u\}\)\\,\\left\\\|\\nabla\_\{\\mathbf\{W\}\}f\_\{\\mathrm\{phys\}\}\(\\mathbf\{u\};\\mathbf\{W\}\)\\right\\\|\_\{2\}\.\(9\)

#### Discrepancy Hunter \(DH\)\.

To actively refine the model’s validity, this strategy targets regions where the structural discrepancy is significant or uncertain\. The acquisition score is defined as:

αD​H​\(𝐮\)=\|μe​r​r​\(𝐮\)\|\+β​σe​r​r​\(𝐮\),\\alpha\_\{DH\}\(\\mathbf\{u\}\)=\|\\mu\_\{err\}\(\\mathbf\{u\}\)\|\+\\beta\\sigma\_\{err\}\(\\mathbf\{u\}\),\(10\)whereμe​r​r\\mu\_\{err\}andσe​r​r\\sigma\_\{err\}are the posterior mean and standard deviation of the GP discrepancy termge​r​rg\_\{err\}\. We setβ=1\.0\\beta=1\.0to equally prioritize the exploration of high\-bias \(mean\) and high\-variance regions\.

#### Physically Constrained EI \(PC\-EI\)\.

To prevent the optimizer from exploiting unphysical regimes \(e\.g\., negative concentrations predicted by the unconstrained GP\), we modulate the standard Expected Improvement \(EI\) with a soft physical feasibility gate:

αP​C−E​I​\(𝐮\)=αE​I​\(𝐮\)⋅σs​i​g​m​o​i​d​\(γ​\(fp​h​y​s​\(𝐮\)−τ\)\),\\alpha\_\{PC\-EI\}\(\\mathbf\{u\}\)=\\alpha\_\{EI\}\(\\mathbf\{u\}\)\\cdot\\sigma\_\{sigmoid\}\\bigl\(\\gamma\(f\_\{phys\}\(\\mathbf\{u\}\)\-\\tau\)\\bigr\),\(11\)whereσs​i​g​m​o​i​d​\(z\)=\(1\+e−z\)−1\\sigma\_\{sigmoid\}\(z\)=\(1\+e^\{\-z\}\)^\{\-1\}acts as a differentiable mask\. The parametersγ\\gamma\(sharpness\) andτ\\tau\(threshold\) ensure the acquisition value vanishes smoothly as the physical model predictionfp​h​y​sf\_\{phys\}drops below a feasibility threshold, avoiding gradient discontinuities associated with hard indicators\.

## 4Experiments

### 4\.1Experimental Setup

#### Datasets\.

The proposed framework is evaluated on two benchmarks designed to assess structural identifiability and optimization efficiency\.

*Hydrogen–Bromine \(H2\+Br2\\mathrm\{H\_\{2\}\+Br\_\{2\}\}\)\.*This benchmark targets mechanistic structure discovery in a stiff radical chain system\. Isothermal simulations atT=600​KT=600\\,\\mathrm\{K\}were generated using Latin Hypercube Sampling \(LHS\) to produce high\-information trajectories spanning three orders of magnitude in concentration\. Thermal effects were intentionally excluded to isolate topological ambiguity from temperature\-induced variability\. In the released benchmark script, initialH2\\mathrm\{H\_\{2\}\}andBr2\\mathrm\{Br\_\{2\}\}concentrations are sampled by LHS over a bounded concentration range and integrated with a stiff ODE solver; the plotted proof\-of\-concept uses additive Gaussian concentration noise\. Gaussian noise is an idealization matching the GP likelihood\. Shot\-noise and heteroscedastic measurement models are not claimed to be covered and are listed as limitations for experimental deployment\.

*Styrene Epoxidation\.*This benchmark represents an industrial closed\-loop optimization task characterized by complex competitive side reactions, including catalyst self\-poisoning\. The objective is to maximize product yield under a strictly limited experimental budget, reflecting realistic laboratory constraints\.

#### Baselines\.

We compare against state\-of\-the\-art methods across structure discovery, surrogate modeling, and experimental design\. PySR and SINDy are used as representative symbolic and sparse regression approaches for mechanism discovery\. A standard Gaussian Process with an RBF kernel combined with Expected Improvement \(EI\) is the GP\-BO baseline in Fig\.[4](https://arxiv.org/html/2606.23757#S4.F4)\. For the proposed CIGP loop, Fig\.[4](https://arxiv.org/html/2606.23757#S4.F4)uses PC\-EI; Table[3](https://arxiv.org/html/2606.23757#S4.T3)reports the added comparison among all acquisition functions, random search, and uncertainty sampling\.

Table 1:Comparison of discovered kinetic parameters versus ground truth\.ParameterSymbolGround TruthDiscovered \(Ours\)Relative ErrorReaction OrdernSty,nPAAn\_\{\\mathrm\{Sty\}\},\\,n\_\{\\mathrm\{PAA\}\}1\.0,1\.01\.0,\\;1\.01\.0,1\.01\.0,\\;1\.00\.00%0\.00\\%Main Activation EnergyEa,mainE\_\{a,\\mathrm\{main\}\}55\.00​kJ/mol55\.00\\,\\mathrm\{kJ/mol\}57\.41​kJ/mol57\.41\\,\\mathrm\{kJ/mol\}4\.38%4\.38\\%Side Activation EnergyEa,sideE\_\{a,\\mathrm\{side\}\}85\.00​kJ/mol85\.00\\,\\mathrm\{kJ/mol\}90\.14​kJ/mol90\.14\\,\\mathrm\{kJ/mol\}6\.05%6\.05\\%Main Pre\-exponentialAmainA\_\{\\mathrm\{main\}\}1\.0×1061\.0\\times 10^\{6\}2\.66×1062\.66\\times 10^\{6\}Order CorrectSide Pre\-exponentialAsideA\_\{\\mathrm\{side\}\}1\.0×10101\.0\\times 10^\{10\}4\.15×10104\.15\\times 10^\{10\}Order CorrectTable 2:Quantitative ablation study on structural identifiability and physical consistency\.Model VariantPriorDistributionThermodynamicConstraintsStructureF1PhysicalValidityPrimary FailureModePC\-MCMC \(Ours\)Spike\-and\-SlabYes1\.0†1\.0^\{\\dagger\}100%100\\%Minor pathway pruningWithout SparsityGaussian \(ℒ2\\mathcal\{L\}\_\{2\}\)Yes0\.730\.73100%100\\%Dense topology \(overfitting\)Without PhysicsSpike\-and\-SlabNo0\.500\.500%0\\%Violated reversibility
†\\daggerRecovers the dominant elementary\-step set in the reported benchmark\.

Table 3:Acquisition\-function comparison for styrene epoxidation over 10 random seeds\. Final best yield is reported as mean±\\pmstandard deviation\. BO violations count suggestions below0\.2​M0\.2\\,\\mathrm\{M\}after the initial LHS design\.StrategyFinal best yieldBO violationsInterpretationEI0\.707±0\.0590\.707\\pm 0\.0591\.41\.4Strong final\-yield baseline for exploitation\.GWU0\.706±0\.0310\.706\\pm 0\.0311\.91\.9Stable final yield while emphasizing parameter\-sensitive regions\.PC\-EI0\.699±0\.0600\.699\\pm 0\.0600\.80\.8Lowest low\-yield violation rate among model\-based strategies\.Uncertainty0\.700±0\.0540\.700\\pm 0\.0548\.58\.5High exploration but many unsafe low\-yield suggestions\.DH0\.668±0\.0610\.668\\pm 0\.0614\.24\.2Useful for model\-error probing, weaker for direct yield maximization\.Random0\.652±0\.0540\.652\\pm 0\.0547\.67\.6Lowest final yield and frequent low\-yield experiments\.

### 4\.2Robust Mechanism Discovery

#### Structural Identifiability and Topological Disentanglement\.

The structural identifiability of the proposed framework is rigorously evaluated on theH2\+Br2\\mathrm\{H\_\{2\}\+Br\_\{2\}\}benchmark using a sparse dataset generated via Latin Hypercube Sampling\. Despite the inherent stiffness of the radical chain mechanism, decisive structural selectivity is achieved by the PC\-MCMC sampler\. As summarized inFig\. 3e, the posterior inclusion probabilities \(PIP\) for the governing elementary steps—including initiation, propagation, and termination—show high posterior support in the reported run\. Crucially, the direct molecular collision pathwayH2\+Br2→2​H​B​r\\mathrm\{H\_\{2\}\+Br\_\{2\}\\rightarrow 2HBr\}, which is mathematically plausible but kinetically prohibited, is successfully rejected \(PIP<0\.01\\mathrm\{PIP\}<0\.01\)\. These results confirm that the explicit enforcement of mass conservation \(𝐀⋅𝐒=0\\mathbf\{A\}\\cdot\\mathbf\{S\}=0\) and thermodynamic consistency effectively prunes the combinatorial search space, isolating the causal micro\-kinetic mechanism from phenomenological correlations\.

#### Numerical Instability of Unconstrained Sparse Regression\.

The necessity of physical constraints is further illustrated by the failure of the SINDy baseline\. As evidenced by the inferred equations, unconstrained sparse regression retrieves an ill\-posed model dominated by spurious fractional terms \(e\.g\.,x0​x10\.5x\_\{0\}x\_\{1\}^\{0\.5\}\) and negative coefficients, violating the non\-negativity of chemical concentrations\. As a consequence, severe numerical instability is observed: the ODE solver diverges after only two time steps \(2/202/20\) during trajectory reconstruction, rendering the learned dynamics physically invalid for extrapolation\.

#### Phenomenological Overfitting by Symbolic Regression\.

The distinction between curve fitting and true mechanism discovery is further highlighted by the PySR results\. Symbolic regression successfully retrieves the macroscopic rate law

d​\[HBr\]d​t≈0\.315​\[H2\]​\[Br2\]0\.5,\\frac\{d\[\\mathrm\{HBr\}\]\}\{dt\}\\approx 0\.315\\,\[\\mathrm\{H\_\{2\}\}\]\\,\[\\mathrm\{Br\_\{2\}\}\]^\{0\.5\},\(12\)with a low training mean squared error of1\.65×10−31\.65\\times 10^\{\-3\}\(complexity=6=6\)\. However, as summarized inFig\. 3e, this solution remains purely phenomenological\. PySR fails to uncover the underlying radical reaction topology involvingH⋅\\mathrm\{H\\\!\\cdot\}andBr⋅\\mathrm\{Br\\\!\\cdot\}intermediates that gives rise to the observed fractional\-order dependence\. In contrast, the proposed framework naturally reconstructs this behavior through the quasi\-equilibrium of bromine dissociation, yielding a physically interpretable derivation fully consistent with the law of mass action\.

![Refer to caption](https://arxiv.org/html/2606.23757v1/fig4.jpg)Figure 4:Comparative evaluation of optimization performance between the Standard GP baseline and the proposed CIGP framework\.\(a–b\) Exploration Trajectories:Three\-dimensional sampling paths in the design space are visualized over the ground\-truth yield landscape\. The Standard GP baseline\(a\)exhibits a scattered exploration pattern with frequent sampling in low\-yield regions, whereas CIGP\(b\)effectively constrains exploration to the high\-yield manifold via physical priors\.\(c\) Optimization Efficiency:Best\-so\-far yield trajectories are compared across iterations\. CIGP \(blue\) achieves significantly faster convergence, approaching the theoretical maximum \(green dashed line\) within five iterations, while the Standard GP baseline \(red\) converges more slowly\.\(d\) Process Stability and Safety:Observed yield per iteration is reported\. The Standard GP demonstrates substantial volatility and multiple violations of the efficiency threshold \(gray shaded region,<0\.2​M<0\.2\\,\\mathrm\{M\}\), whereas the PC\-EI BO phase of CIGP avoids low\-yield suggestions in this run\.

### 4\.3Efficient Closed\-Loop Optimization and Mechanism Identification

#### Optimization Efficiency and Convergence Analysis\.

We evaluate closed\-loop optimization performance of the proposed CIGP framework against a Standard Bayesian Optimization \(GP\-BO\) baseline on the styrene epoxidation benchmark\. The exploration trajectories shown in Fig\.[4](https://arxiv.org/html/2606.23757#S4.F4)\(a–b\) reveal fundamentally different search behaviors\. The GP\-BO baseline exhibits a scattered and unguided exploration pattern, with frequent redundant sampling in low\-yield regions\. In contrast, CIGP constrains exploration to the high\-yield manifold, exhibiting a gradient\-guided search behavior induced by physical priors\.

Quantitative convergence results in Fig\.[4](https://arxiv.org/html/2606.23757#S4.F4)\(c\) further highlight this disparity\. Using PC\-EI, CIGP enters the high\-yield regime \(\>0\.70​M\>0\.70\\,\\mathrm\{M\}\) within four BO iterations after the initial LHS design and reaches a maximum observed yield of0\.7788​M0\.7788\\,\\mathrm\{M\}by the fifth BO iteration\. By comparison, the GP\-BO baseline fails to exceed the0\.70​M0\.70\\,\\mathrm\{M\}threshold even after fifteen iterations, achieving a best observed yield of only0\.6925​M0\.6925\\,\\mathrm\{M\}\. Overall, the physics\-informed framework achieves a3\.75×3\.75\\timesacceleration in convergence and a12\.5%12\.5\\%improvement in final yield over the purely data\-driven baseline\.

#### Process Stability and Safety Assessment\.

Optimization robustness is assessed via per\-iteration observed yields \(Fig\.[4](https://arxiv.org/html/2606.23757#S4.F4)\(d\)\)\. The GP\-BO baseline incurs four severe violations, defined as yields below0\.2​M0\.2\\,\\mathrm\{M\}\(gray shaded region\), corresponding to experimentally infeasible operating conditions\. Such volatility implies substantial resource waste and elevated operational risk in industrial settings\. In contrast, no violations are observed during the PC\-EI BO phase; low\-yield points can still occur in the shared initial LHS design\. The resulting74%74\\%reduction in cumulative regret in this benchmark indicates that physical inductive biases effectively regularize exploration during the data\-sparse initialization phase, preventing excursions into physically invalid regions\.

#### Parameter Identification Accuracy\.

The superior optimization efficiency of CIGP is attributed to accurate identification of the underlying reaction kinetics\. Unlike black\-box surrogate models, micro\-kinetic parameters are inferred jointly within the CIGP optimization loop\. As reported in Table[1](https://arxiv.org/html/2606.23757#S4.T1), the reaction orders of styrene and PAA are correctly identified as unity, consistent with the theoretical bimolecular mechanism\. Furthermore, relative errors in the inferred activation energies are limited to4\.38%4\.38\\%for the main reaction and6\.05%6\.05\\%for side reactions\. These results demonstrate that rapid convergence is driven by a physically consistent reconstruction of the reaction energy landscape rather than numerical curve fitting alone\.

### 4\.4Ablation Study

#### Setup\.

An ablation study is conducted on theH2\+Br2\\mathrm\{H\_\{2\}\+Br\_\{2\}\}benchmark to quantify the contributions of sparsity\-inducing priors and thermodynamic constraints\. Two variants are evaluated: \(i\)*without Sparsity*, using a continuous Gaussian prior; and \(ii\)*without Physics*, disabling thermodynamic consistency checks\.

#### Effect of the Sparsity Prior\.

As reported in Table[2](https://arxiv.org/html/2606.23757#S4.T2), the full PC\-MCMC framework achieves a structural F1 score of 1\.0, corresponding to the recovery of all ground\-truth elementary steps\. Crucially, the reverse propagation stepHBr\+H⋅→Br⋅\+H2HBr\+H\\cdot\\rightarrow Br\\cdot\+H\_\{2\}, which is often challenging to identify due to thermodynamic constraints, is successfully retrieved with high confidence \(PIP\>0\.95\\mathrm\{PIP\}\>0\.95, see Fig\. 3a\)\. This demonstrates that the Spike\-and\-Slab prior, combined with thermodynamic regularization, effectively resolves the ambiguity between low\-flux reversible pathways and spurious noise, a capability that is lost when the sparsity prior is removed \(F1 drops to 0\.73\)\.

#### Role of Thermodynamic Constraints\.

The importance of physical regularization is further evidenced by the failure of the without Physics variant\. Relaxing thermodynamic constraints reduces the structural F1 score to0\.500\.50and eliminates physical validity of inferred parameters\. Specifically, microscopic reversibility is violated: the forward initiation stepBr2→2Br⋅\\mathrm\{Br\_\{2\}\\rightarrow 2Br\\cdot\}is selected \(PIP=1\.0\\mathrm\{PIP\}=1\.0\), while its reverse termination step is rejected \(PIP=0\.0\\mathrm\{PIP\}=0\.0\)\. To compensate for missing pathways, the rate constant of the rate\-determining stepBr⋅\+H2\\mathrm\{Br\\cdot\+H\_\{2\}\}is overestimated by three orders of magnitude \(kest≈9\.8×102k\_\{\\mathrm\{est\}\}\\approx 9\.8\\times 10^\{2\}vs\.ktrue=1\.0k\_\{\\mathrm\{true\}\}=1\.0\)\. These results demonstrate that thermodynamic regularization is indispensable for ensuring interpretability, numerical stability, and extrapolative validity of the discovered mechanisms \(see Appendix Fig\.[5](https://arxiv.org/html/2606.23757#A1.F5)\)\.

## 5Conclusion

We proposed PC\-MCMC\-CIGP, a unified framework that combines physically constrained Bayesian topology inference with Chemical\-Informed Gaussian Processes\. In our benchmarks, hard conservation and detailed\-balance checks improve interpretability of mechanism discovery, while gray\-box active learning improves styrene epoxidation optimization relative to the reported GP\-BO and random\-search baselines\. The added acquisition comparison clarifies that the acquisition functions have different purposes: EI\-style criteria are competitive for final yield, whereas PC\-EI reduces low\-yield experimental suggestions\.

Limitations and Future Directions\.The present evidence is based on synthetic and oracle\-based benchmarks rather than wet\-lab closed\-loop campaigns\. The PC\-MCMC stage relies on a pre\-enumerated candidate reaction set, and rejection sampling may become inefficient when many reversible pairs impose detailed\-balance constraints\. The observation model currently assumes Gaussian noise; Poisson shot noise, heteroscedastic sensor error, and asynchronous trajectory measurements require additional likelihood models\. Future work will address these scalability and realism gaps using adaptive proposals, variational approximations, and automated reaction\-template generation\.

## Impact Statement

This work is a step toward trustworthy and autonomous scientific discovery in the chemical sciences\. By synergizing physical laws with probabilistic machine learning, our PC\-MCMC\-CIGP framework addresses critical challenges in interpretability, safety, and efficiency that have hindered the adoption of black\-box AI in high\-stakes laboratory environments\.

Accelerating Green Chemistry\.In the styrene epoxidation benchmark, the proposed framework demonstrates faster process optimization and a 74% reduction in cumulative regret compared to the reported GP\-BO baseline\. In industrial applications, such efficiency gains directly translate to reduced consumption of raw materials and energy, facilitating the rapid development of sustainable chemical processes and supporting global efforts toward carbon neutrality\.

Safety\-Critical AI\.Unlike purely data\-driven approaches that may suggest hazardous operating conditions, our framework explicitly enforces thermodynamic and conservation constraints\. This ”safety\-by\-design” philosophy prevents the algorithm from exploring physically prohibited or dangerous regimes, making it suitable for direct integration into autonomous robotic laboratories where human supervision is limited\.

Democratization of Scientific AI\.We have demonstrated that rigorous mechanism discovery is feasible on standard CPUs without requiring massive GPU clusters \(see Appendix[B](https://arxiv.org/html/2606.23757#A2)\)\. This low computational barrier democratizes access to advanced AI tools, empowering experimentalists in resource\-constrained settings to leverage data\-driven insights for mechanistic research\.

Broader Applicability\.While validated on chemical kinetics, the core principle of embedding domain\-specific differential equations into Bayesian priors is transferable to other scientific domains, such as systems biology and climate modeling, where data is sparse but physical knowledge is abundant\.

## References

- S\. Audoly, G\. Bellu, L\. D’Angio, M\. P\. Saccomani, and C\. Cobelli \(2002\)Global identifiability of nonlinear models of biological systems\.IEEE Transactions on biomedical engineering48\(1\),pp\. 55–65\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p2.1)\.
- D\. Bonvin and D\. Rippin \(1990\)Target factor analysis for the identification of stoichiometric models\.Chemical Engineering Science45\(12\),pp\. 3417–3426\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- M\. Brubaker, M\. Salzmann, and R\. Urtasun \(2012\)A family of mcmc methods on implicitly defined manifolds\.InArtificial intelligence and statistics,pp\. 161–172\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p3.1)\.
- S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz \(2016\)Discovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the national academy of sciences113\(15\),pp\. 3932–3937\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p1.1),[§2\.1](https://arxiv.org/html/2606.23757#S2.SS1.p1.1)\.
- B\. Burlacu, G\. Kronberger, and M\. Kommenda \(2020\)Operon c\+\+ an efficient genetic programming framework for symbolic regression\.InProceedings of the 2020 genetic and evolutionary computation conference companion,pp\. 1562–1570\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.1](https://arxiv.org/html/2606.23757#S2.SS1.p1.1)\.
- S\. C\. Burnham, D\. P\. Searson, M\. J\. Willis, and A\. R\. Wright \(2008\)Inference of chemical reaction networks\.Chemical Engineering Science63\(4\),pp\. 862–873\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p1.1)\.
- K\. Champion, B\. Lusch, J\. N\. Kutz, and S\. L\. Brunton \(2019\)Data\-driven discovery of coordinates and governing equations\.Proceedings of the National Academy of Sciences116\(45\),pp\. 22445–22451\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p1.1)\.
- C\. Chang and T\. Zeng \(2023\)A hybrid data\-driven\-physics\-constrained gaussian process regression framework with deep kernel for uncertainty quantification\.Journal of Computational Physics486,pp\. 112129\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p4.1),[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- M\. Cranmer \(2023\)Interpretable machine learning for science with pysr and symbolicregression\. jl\.arXiv preprint arXiv:2305\.01582\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p1.1),[§2\.1](https://arxiv.org/html/2606.23757#S2.SS1.p1.1)\.
- E\. J\. Cross, T\. J\. Rogers, D\. J\. Pitchforth, S\. J\. Gibson, S\. Zhang, and M\. R\. Jones \(2024\)A spectrum of physics\-informed gaussian processes for regression in engineering\.Data\-Centric Engineering5,pp\. e8\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p4.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- S\. Da Veiga and A\. Marrel \(2012\)Gaussian process modeling with inequality constraints\.InAnnales de la Faculté des sciences de Toulouse: Mathématiques,Vol\.21,pp\. 529–555\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- D\. Dalton, A\. Lazarus, H\. Gao, and D\. Husmeier \(2024\)Boundary constrained gaussian processes for robust physics\-informed machine learning of linear partial differential equations\.Journal of Machine Learning Research25\(272\),pp\. 1–61\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p4.1)\.
- F\. P\. Davidescu and S\. B\. Jørgensen \(2008\)Structural parameter identifiability analysis for dynamic reaction networks\.Chemical Engineering Science63\(19\),pp\. 4754–4762\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p2.1)\.
- M\. Ederer and E\. D\. Gilles \(2007\)Thermodynamically feasible kinetic models of reaction networks\.Biophysical journal92\(6\),pp\. 1846–1857\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- G\. Enciso, R\. Erban, and J\. Kim \(2021\)Identifiability of stochastically modelled reaction networks\.European Journal of Applied Mathematics32\(5\),pp\. 865–887\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p2.1)\.
- M\. Feinberg \(1987\)Chemical reaction network structure and the stability of complex isothermal reactors—i\. the deficiency zero and deficiency one theorems\.Chemical engineering science42\(10\),pp\. 2229–2268\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- E\. I\. George and R\. E\. McCulloch \(1993\)Variable selection via gibbs sampling\.Journal of the American Statistical Association88\(423\),pp\. 881–889\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p1.2)\.
- M\. Girolami and B\. Calderhead \(2011\)Riemann manifold langevin and hamiltonian monte carlo methods\.Journal of the Royal Statistical Society Series B: Statistical Methodology73\(2\),pp\. 123–214\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p3.4)\.
- F\. Hase, L\. M\. Roch, C\. Kreisbeck, and A\. Aspuru\-Guzik \(2018\)Phoenics: a bayesian optimizer for chemistry\.ACS central science4\(9\),pp\. 1134–1145\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- L\. Hewing, J\. Kabzan, and M\. N\. Zeilinger \(2019\)Cautious model predictive control using gaussian process regression\.IEEE Transactions on Control Systems Technology28\(6\),pp\. 2736–2743\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.3](https://arxiv.org/html/2606.23757#S3.SS3.SSS0.Px1.p1.6)\.
- F\. Horn and R\. Jackson \(1972\)General mass action kinetics\.Archive for rational mechanics and analysis47\(2\),pp\. 81–116\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.1](https://arxiv.org/html/2606.23757#S3.SS1.p1.14)\.
- H\. Ishwaran and J\. S\. Rao \(2005\)Spike and slab variable selection: frequentist and bayesian strategies\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p3.4)\.
- G\. Jenkinson, X\. Zhong, and J\. Goutsias \(2010\)Thermodynamically consistent bayesian analysis of closed biochemical reaction systems\.BMC bioinformatics11\(1\),pp\. 547\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- R\. Jiang, P\. Singh, F\. Wrede, A\. Hellander, and L\. Petzold \(2022\)Identification of dynamic mass\-action biochemical reaction networks using sparse bayesian methods\.PLoS computational biology18\(1\),pp\. e1009830\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p1.1)\.
- D\. R\. Jones, M\. Schonlau, and W\. J\. Welch \(1998\)Efficient global optimization of expensive black\-box functions\.Journal of Global optimization13\(4\),pp\. 455–492\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.4](https://arxiv.org/html/2606.23757#S3.SS4.SSS0.Px1.p1.1)\.
- M\. C\. Kennedy and A\. O’Hagan \(2001\)Bayesian calibration of computer models\.Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)63\(3\),pp\. 425–464\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- S\. Kim, W\. Ji, S\. Deng, Y\. Ma, and C\. Rackauckas \(2021\)Stiff neural ordinary differential equations\.Chaos: An Interdisciplinary Journal of Nonlinear Science31\(9\),pp\. 093122\.Cited by:[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- J\. Kocijan \(2016\)Modelling and control of dynamic systems using gaussian process models\.Springer\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- N\. R\. Kristensen, H\. Madsen, and S\. B\. Jørgensen \(2004\)Parameter estimation in stochastic grey\-box models\.Automatica40\(2\),pp\. 225–237\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.3](https://arxiv.org/html/2606.23757#S3.SS3.p1.1)\.
- D\. Langary and Z\. Nikoloski \(2019\)Inference of chemical reaction networks based on concentration profiles using an optimization framework\.Chaos: An Interdisciplinary Journal of Nonlinear Science29\(11\)\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- T\. Lubitz, M\. Schulz, E\. Klipp, and W\. Liebermeister \(2010\)Parameter balancing in kinetic models of cell metabolism\.The Journal of Physical Chemistry B114\(49\),pp\. 16298–16303\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- T\. Ma, D\. A\. Barajas\-Solano, R\. Tipireddy, and A\. M\. Tartakovsky \(2020\)Physics\-informed gaussian process regression for probabilistic states estimation and forecasting in power grids\.arXiv preprint arXiv:2010\.04591\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p4.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- N\. Metropolis, A\. W\. Rosenbluth, M\. N\. Rosenbluth, A\. H\. Teller, and E\. Teller \(1953\)Equation of state calculations by fast computing machines\.The journal of chemical physics21\(6\),pp\. 1087–1092\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p3.4)\.
- T\. J\. Mitchell and J\. J\. Beauchamp \(1988\)Bayesian variable selection in linear regression\.Journal of the american statistical association83\(404\),pp\. 1023–1032\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p1.2)\.
- C\. Otte \(2014\)Interpretable semi\-parametric regression models with defined error bounds\.Neurocomputing143,pp\. 1–6\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.1](https://arxiv.org/html/2606.23757#S2.SS1.p1.1)\.
- O\. Owoyele and P\. Pal \(2021\)ChemNODE: a neural ordinary differential equations approach for chemical kinetics solvers\.Energy and AI5,pp\. 100082\.Cited by:[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- A\. Pal, R\. Fagerberg, J\. L\. Andersen, C\. Flamm, P\. Dittrich, and D\. Merkle \(2025\)Finding thermodynamically favorable pathways in chemical reaction networks using flows in hypergraphs and mixed\-integer linear programming\.Journal of Chemical Information and Modeling\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- T\. Park and G\. Casella \(2008\)The bayesian lasso\.Journal of the american statistical association103\(482\),pp\. 681–686\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- P\. H\. Peskun \(1973\)Optimum monte\-carlo sampling using markov chains\.Biometrika60\(3\),pp\. 607–612\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p3.4)\.
- J\. Piironen and A\. Vehtari \(2017\)Sparsity information and regularization in the horseshoe and other shrinkage priors\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.2](https://arxiv.org/html/2606.23757#S3.SS2.p3.4)\.
- M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis \(2017\)Machine learning of linear differential equations using gaussian processes\.Journal of Computational Physics348,pp\. 683–693\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p4.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- S\. H\. Rudy, S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz \(2017\)Data\-driven discovery of partial differential equations\.Science advances3\(4\),pp\. e1602614\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.1](https://arxiv.org/html/2606.23757#S2.SS1.p1.1)\.
- N\. Sani, J\. Lee, R\. Nabi, and I\. Shpitser \(2020\)A semiparametric approach to interpretable machine learning\.arXiv preprint arXiv:2006\.04732\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- M\. Schmidt and H\. Lipson \(2009\)Distilling free\-form natural laws from experimental data\.science324\(5923\),pp\. 81–85\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p1.1),[§2\.1](https://arxiv.org/html/2606.23757#S2.SS1.p1.1)\.
- D\. P\. Searson, M\. J\. Willis, S\. J\. Horne, and A\. R\. Wright \(2007\)Inference of chemical reaction networks using hybrid s\-system models\.Chemical Product and Process Modeling2\(1\)\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.
- B\. Shahriari, K\. Swersky, Z\. Wang, R\. P\. Adams, and N\. De Freitas \(2015\)Taking the human out of the loop: a review of bayesian optimization\.Proceedings of the IEEE104\(1\),pp\. 148–175\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.4](https://arxiv.org/html/2606.23757#S3.SS4.SSS0.Px1.p1.1)\.
- B\. J\. Shields, J\. Stevens, J\. Li, M\. Parasram, F\. Damani, J\. I\. M\. Alvarado, J\. M\. Janey, R\. P\. Adams, and A\. G\. Doyle \(2021\)Bayesian reaction optimization as a tool for chemical synthesis\.Nature590\(7844\),pp\. 89–96\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- J\. Snoek, H\. Larochelle, and R\. P\. Adams \(2012\)Practical bayesian optimization of machine learning algorithms\.Advances in neural information processing systems25\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- R\. Tibshirani \(1996\)Regression shrinkage and selection via the lasso\.Journal of the Royal Statistical Society Series B: Statistical Methodology58\(1\),pp\. 267–288\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- A\. Van der Schaft, S\. Rao, and B\. Jayawardhana \(2015\)Complex and detailed balancing of chemical reaction networks revisited\.Journal of Mathematical Chemistry53\(6\),pp\. 1445–1458\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- M\. O\. Vlad and J\. Ross \(2009\)Thermodynamically based constraints for rate coefficients of large biochemical networks\.Wiley Interdisciplinary Reviews: Systems Biology and Medicine1\(3\),pp\. 348–358\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§3\.1](https://arxiv.org/html/2606.23757#S3.SS1.p1.14)\.
- E\. Walter and L\. Pronzato \(1997\)Identification of parametric models: from experimental data\.\(No Title\)\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p2.1)\.
- D\. J\. Warne, R\. E\. Baker, and M\. J\. Simpson \(2019\)Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state\-of\-the\-art\.Journal of the Royal Society Interface16\(151\),pp\. 20180943\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1)\.
- C\. K\. Williams and C\. E\. Rasmussen \(2006\)Gaussian processes for machine learning\.Vol\.2,MIT press Cambridge, MA\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§1](https://arxiv.org/html/2606.23757#S1.p2.1),[§2\.3](https://arxiv.org/html/2606.23757#S2.SS3.p1.1)\.
- M\. J\. Willis and M\. von Stosch \(2016\)Inference of chemical reaction networks using mixed integer linear programming\.Computers & Chemical Engineering90,pp\. 31–43\.Cited by:[§B\.2](https://arxiv.org/html/2606.23757#A2.SS2.p3.1),[§2\.2](https://arxiv.org/html/2606.23757#S2.SS2.p1.1)\.

## Appendix APC\-MCMC Inference Algorithm

The detailed procedure for the Physically Constrained MCMC \(PC\-MCMC\) sampler is outlined in Algorithm[1](https://arxiv.org/html/2606.23757#alg1)\. This algorithm ensures that all sampled reaction networks strictly adhere to mass conservation, electron conservation, and thermodynamic consistency \(detailed balance\) while efficiently exploring the combinatorial topology space using a Spike\-and\-Slab prior\.

Fixed\-Dimension Parameter Space\.Contrary to trans\-dimensional approaches \(e\.g\., Reversible Jump MCMC\), our Spike\-and\-Slab formulation maintains a fixed parameter dimension corresponding to the total number of candidate reactions\|ℛ\|\|\\mathcal\{R\}\|\. The auxiliary parameters𝜽\\boldsymbol\{\\theta\}remain instantiated even for inactive reactions \(γj=0\\gamma\_\{j\}=0\), allowing standard Metropolis\-Hastings updates without dimension\-matching Jacobians\.

Structure Move \(Topology Sampling\)\.The proposal distributionq​\(𝜸′∣𝜸\(t−1\)\)q\(\\boldsymbol\{\\gamma\}^\{\\prime\}\\mid\\boldsymbol\{\\gamma\}^\{\(t\-1\)\}\)consists of selecting a reaction indexjjuniformly at random and flipping its status \(γj′=1−γj\(t−1\)\\gamma^\{\\prime\}\_\{j\}=1\-\\gamma\_\{j\}^\{\(t\-1\)\}\)\. Since this proposal is symmetric, the acceptance ratioαs​t​r​u​c​t\\alpha\_\{struct\}\(Algorithm 1, Line 9\) simplifies to the ratio of unnormalized posteriors:

αs​t​r​u​c​t=min⁡\(1,P​\(𝒟∣𝐒​\(𝜸′\),𝐤​\(𝜸′,𝜽\)\)⋅P​\(𝜸′\)P​\(𝒟∣𝐒​\(𝜸\),𝐤​\(𝜸,𝜽\)\)⋅P​\(𝜸\)\)\\alpha\_\{struct\}=\\min\\left\(1,\\frac\{P\(\\mathcal\{D\}\\mid\\mathbf\{S\}\(\\boldsymbol\{\\gamma\}^\{\\prime\}\),\\mathbf\{k\}\(\\boldsymbol\{\\gamma\}^\{\\prime\},\\boldsymbol\{\\theta\}\)\)\\cdot P\(\\boldsymbol\{\\gamma\}^\{\\prime\}\)\}\{P\(\\mathcal\{D\}\\mid\\mathbf\{S\}\(\\boldsymbol\{\\gamma\}\),\\mathbf\{k\}\(\\boldsymbol\{\\gamma\},\\boldsymbol\{\\theta\}\)\)\\cdot P\(\\boldsymbol\{\\gamma\}\)\}\\right\)\(13\)where𝐒​\(𝜸\)\\mathbf\{S\}\(\\boldsymbol\{\\gamma\}\)is the stoichiometry defined by the active subset\. Crucially, during a structure move, the auxiliary parameters𝜽\\boldsymbol\{\\theta\}are held fixed\.

Parameter Move\.We employ a Random Walk Metropolis kernel for the continuous parameters\. The proposal distribution is a Gaussian centered at the current value:θj′∼𝒩​\(θj\(t−1\),σp​r​o​p2\)\\theta^\{\\prime\}\_\{j\}\\sim\\mathcal\{N\}\(\\theta\_\{j\}^\{\(t\-1\)\},\\sigma\_\{prop\}^\{2\}\)\. The acceptance ratio is:

αp​a​r​a​m=min⁡\(1,P​\(𝒟∣𝜸,𝜽′\)​P​\(𝜽′∣𝜸\)P​\(𝒟∣𝜸,𝜽\)​P​\(𝜽∣𝜸\)\)\\alpha\_\{param\}=\\min\\left\(1,\\frac\{P\(\\mathcal\{D\}\\mid\\boldsymbol\{\\gamma\},\\boldsymbol\{\\theta\}^\{\\prime\}\)P\(\\boldsymbol\{\\theta\}^\{\\prime\}\\mid\\boldsymbol\{\\gamma\}\)\}\{P\(\\mathcal\{D\}\\mid\\boldsymbol\{\\gamma\},\\boldsymbol\{\\theta\}\)P\(\\boldsymbol\{\\theta\}\\mid\\boldsymbol\{\\gamma\}\)\}\\right\)\(14\)If thermodynamic constraints are violated by𝜽′\\boldsymbol\{\\theta\}^\{\\prime\}, the prior termP​\(𝜽′∣𝜸\)P\(\\boldsymbol\{\\theta\}^\{\\prime\}\\mid\\boldsymbol\{\\gamma\}\)evaluates to zero, leading to automatic rejection\.

Algorithm 1Bayesian Inference of Reaction Network Topology and Kinetic Parameters0:Data

𝒟\\mathcal\{D\}, Atomic Matrix

𝐀\\mathbf\{A\}, Candidate Reactions

ℛ\\mathcal\{R\}
0:Posterior samples of Topology

𝐒\\mathbf\{S\}and Parameters

𝐤\\mathbf\{k\}
1:Initialize

γ\(0\),𝐤\(0\)\\gamma^\{\(0\)\},\\mathbf\{k\}^\{\(0\)\}satisfying constraints

2:for

t=1t=1to

TTdo

3:// Structure Move \(Topology Sampling\)

4:Propose new topology

γ′\\gamma^\{\\prime\}by flipping bit

jjof

γ\(t−1\)\\gamma^\{\(t\-1\)\}
5:Construct candidate stoichiometry

𝐒′\\mathbf\{S\}^\{\\prime\}based on

γ′\\gamma^\{\\prime\}
6:if

𝐀⋅𝐒⋅j′≠𝟎\\mathbf\{A\}\\cdot\\mathbf\{S\}^\{\\prime\}\_\{\\cdot j\}\\neq\\mathbf\{0\}then

7:Reject

γ′\\gamma^\{\\prime\}\(Violation of Mass/Electron Conservation\)

8:else

9:Compute Acceptance Ratio

αs​t​r​u​c​t\\alpha\_\{struct\}with Sparsity Prior

10:Accept

γ\(t\)=γ′\\gamma^\{\(t\)\}=\\gamma^\{\\prime\}with probability

min⁡\(1,αs​t​r​u​c​t\)\\min\(1,\\alpha\_\{struct\}\)
11:endif

12:

13:// Parameter Move \(Kinetic Constant Sampling\)

14:Propose

𝐤′\\mathbf\{k\}^\{\\prime\}from proposal distribution

q​\(𝐤′∣𝐤\(t−1\)\)q\(\\mathbf\{k\}^\{\\prime\}\\mid\\mathbf\{k\}^\{\(t\-1\)\}\)
15:ifThermodynamic Loop Constraints \(

Δ​G\\Delta G\) violatedthen

16:Reject

𝐤′\\mathbf\{k\}^\{\\prime\}\(Enforce Thermodynamic Consistency\)

17:else

18:Calculate Likelihood

P​\(𝒟∣𝐒\(t\),𝐤′\)P\(\\mathcal\{D\}\\mid\\mathbf\{S\}^\{\(t\)\},\\mathbf\{k\}^\{\\prime\}\)via stiff ODE integration

19:Accept

𝐤\(t\)=𝐤′\\mathbf\{k\}^\{\(t\)\}=\\mathbf\{k\}^\{\\prime\}based on Metropolis criterion

20:endif

21:endfor

22:returnMCMC chains for

𝐒\\mathbf\{S\}and

𝐤\\mathbf\{k\}

Algorithm 2CIGP\-driven Active Learning0:Initial data

𝒟i​n​i​t\\mathcal\{D\}\_\{init\}, Physical ODE model

fp​h​y​sf\_\{phys\}, Acquisition function

α​\(⋅\)\\alpha\(\\cdot\), Budget

NB​ON\_\{BO\}
0:Optimized process conditions

𝐮∗\\mathbf\{u\}^\{\*\}, Calibrated parameters

𝐖∗\\mathbf\{W\}^\{\*\}
1:

𝒟←𝒟i​n​i​t\\mathcal\{D\}\\leftarrow\\mathcal\{D\}\_\{init\}
2:for

i=1i=1to

NB​ON\_\{BO\}do

3:// Hybrid Model Construction

4:Define Mean Function:

μ​\(𝐮\)=fp​h​y​s​\(𝐮;𝐖\)\\mu\(\\mathbf\{u\}\)=f\_\{phys\}\(\\mathbf\{u\};\\mathbf\{W\}\)
5:Define Kernel:

k​\(𝐮,𝐮′\)=RBF​\(𝐮,𝐮′;Θg​p\)k\(\\mathbf\{u\},\\mathbf\{u\}^\{\\prime\}\)=\\text\{RBF\}\(\\mathbf\{u\},\\mathbf\{u\}^\{\\prime\};\\Theta\_\{gp\}\)
6:

7:// Joint Training \(Constrained Optimization\)

8:Define Objective:

J​\(𝐖,Θg​p\)=ℒG​P−λ​σf2J\(\\mathbf\{W\},\\Theta\_\{gp\}\)=\\mathcal\{L\}\_\{GP\}\-\\lambda\\sigma\_\{f\}^\{2\}
9:Solve:

𝐖∗,Θg​p∗←argmax​J\\mathbf\{W\}^\{\*\},\\Theta\_\{gp\}^\{\*\}\\leftarrow\\text\{argmax \}J
10:subject to

σm​i​n2≤σf2≤ϵt​o​l\\sigma\_\{min\}^\{2\}\\leq\\sigma\_\{f\}^\{2\}\\leq\\epsilon\_\{tol\}
11:

12:// Physics\-Aware Acquisition

13:Select next experiment point

𝐮n​e​w\\mathbf\{u\}\_\{new\}:

14:

𝐮n​e​w=argmax𝐮∈𝒰​α​\(𝐮∣𝒟,𝐖∗\)\\mathbf\{u\}\_\{new\}=\\text\{argmax\}\_\{\\mathbf\{u\}\\in\\mathcal\{U\}\}\\alpha\(\\mathbf\{u\}\\mid\\mathcal\{D\},\\mathbf\{W\}^\{\*\}\)\(e\.g\., GWU, PC\-EI\)

15:

16:// Experimental Feedback

17:Conduct experiment:

yn​e​w←LabOracle​\(𝐮n​e​w\)y\_\{new\}\\leftarrow\\text\{LabOracle\}\(\\mathbf\{u\}\_\{new\}\)
18:Update dataset:

𝒟←𝒟∪\{𝐮n​e​w,yn​e​w\}\\mathcal\{D\}\\leftarrow\\mathcal\{D\}\\cup\\\{\\mathbf\{u\}\_\{new\},y\_\{new\}\\\}
19:endfor

20:returnBest found condition

𝐮∗\\mathbf\{u\}^\{\*\}and Parameters

𝐖∗\\mathbf\{W\}^\{\*\}

Integrated Workflow\.The proposed framework operates in two coupled phases to achieve robust autonomous discovery\. First,PC\-MCMC \(Algorithm 1\)is employed to identify the most probable reaction topology and estimate initial kinetic parameters from preliminary data\. Second,CIGP\-driven Active Learning \(Algorithm 2\)utilizes the identified topology \(embedded infp​h​y​sf\_\{phys\}\) to construct a hybrid surrogate model\. This model guides the experimental design to simultaneously refine the kinetic parameters𝐖\\mathbf\{W\}and optimize process conditions𝐮\\mathbf\{u\}, as detailed below\.

![Refer to caption](https://arxiv.org/html/2606.23757v1/figA1.jpg)Figure 5:Visualization of structural failure modes under ablation settings\.Panels a and b depict the absence of the sparsity prior where the sampler converges to a dense topology in panel a and assigns high posterior probability to kinetically prohibited pathways to minimize residual error in panel b\. Panels c and d illustrate the absence of thermodynamic constraints where the principle of microscopic reversibility is violated in panel c as forward chain initiation is selected with a PIP near 1\.0 while the reverse termination is suppressed with a PIP near 0\.0 resulting in a physically inconsistent mechanism despite accurate trajectory fitting in panel d\.Visual Analysis of Failure Modes\.The structural consequences of relaxing inductive biases are further elucidated through posterior landscape visualization as shown in Figure[5](https://arxiv.org/html/2606.23757#A1.F5)\. In the absence of the sparsity\-inducing prior as illustrated in panels a and b of Figure[5](https://arxiv.org/html/2606.23757#A1.F5)a distinctively dense topology is retrieved\. High posterior inclusion probabilities exceeding 0\.6 are assigned indiscriminately to candidate reactions including kinetically prohibitive pathways such as the direct dissociation of hydrogen\. While reasonable alignment with observational data is maintained in the reconstructed trajectory this is attributed to overfitting wherein spurious structural complexity is recruited to absorb measurement noise rather than isolating the causal signal\. Conversely the violation of microscopic reversibility upon the removal of thermodynamic constraints is starkly illustrated in panels c and d of Figure[5](https://arxiv.org/html/2606.23757#A1.F5)\. A structural asymmetry is observed where forward reactions such as chain initiation are identified with unit probability whereas the corresponding reverse steps are completely suppressed with a posterior inclusion probability approaching zero effectively decoupling the equilibrium relationship\. Although a tight numerical fit is achieved via parameter drift the loss of detailed balance renders the discovered mechanism physically inconsistent and unreliable for extrapolation beyond the training regime\.

#### Implementation Details\.

For the mechanism discovery baselines,SINDywas configured with a polynomial feature library \(degree≤3\\leq 3\) and the STLSQ optimizer \(thresholdλ=0\.1\\lambda=0\.1\)\.PySRwas trained for 100 iterations using default binary operators\(\+,−,∗,/\)\(\+,\-,\*,/\)and a complexity penalty of 0\.001\. For theH2\+B​r2H\_\{2\}\+Br\_\{2\}benchmark, initial concentrations were sampled via LHS from\[0\.1,5\.0\]\[0\.1,5\.0\]M and trajectories were generated by integrating the mass\-action ODE\. TheStyrene Epoxidation Oraclesimulates a competitive reaction network: main epoxidation \(k1=106​e−55000/R​Tk\_\{1\}=10^\{6\}e^\{\-55000/RT\}\) and side hydrolysis \(k2=1010​e−85000/R​Tk\_\{2\}=10^\{10\}e^\{\-85000/RT\}\), subject to standard Arrhenius kinetics\. The released code includes scripts for Fig\. 3, Fig\. 4, the acquisition comparison, and smoke\-test verification; generated CSV summaries are written underexamples/outputs/\.

## Appendix BComputational Complexity and Runtime Analysis

To assess the practical feasibility of the proposed framework, we evaluated the wall\-clock execution time of both the PC\-MCMC \(Structure Discovery\) and CIGP \(Active Learning\) modules\. All experiments were conducted on a standard personal laptop equipped with an Intel Core i7\-12700H CPU \(2\.30 GHz\) and 32 GB RAM, without utilizing any GPU acceleration\.

### B\.1Runtime Breakdown

Table[4](https://arxiv.org/html/2606.23757#A2.T4)summarizes the computational cost compared to standard baselines\.

Table 4:Average wall\-clock time for model inference and optimization tasks\.Task / MethodComputational LoadRuntimeHardwarePhase 1: Mechanism Discovery \(H2\+B​r2H\_\{2\}\+Br\_\{2\}Benchmark\)PC\-MCMC \(Ours\)10,000 MCMC steps \(w/ Stiff ODE\)28\.5±2\.128\.5\\pm 2\.1minCPUSINDy \(Baseline\)Sparse Regression \(STLSQ\)<1\.0<1\.0secCPUPySR \(Baseline\)Symbolic Regression \(100 gens\)5\.2±0\.55\.2\\pm 0\.5minCPUPhase 2: Closed\-Loop Optimization \(Styrene Epoxidation\)CIGP \(Ours\)Hyperparam Opt \+ Acquisition \(per iter\)1\.2±0\.31\.2\\pm 0\.3secCPUStandard GP \(Baseline\)Analytical Inference \(per iter\)<0\.1<0\.1secCPU
### B\.2Analysis of Computational Overhead

Structure Discovery \(PC\-MCMC\)\.The PC\-MCMC sampler requires solving a stiff ODE system at each Metropolis\-Hastings step to evaluate the likelihood\. While this incurs a higher computational cost \(≈30\\approx 30minutes\) compared to regression\-based methods like SINDy, it is negligible in the context of mechanistic discovery, which is typically a one\-time offline process\. Furthermore, the fixed\-dimensional parameter space of our Spike\-and\-Slab formulation avoids the complexities of trans\-dimensional sampling \(e\.g\., RJMCMC\), ensuring convergence within a reasonable timeframe on personal computing devices\.

Active Learning \(CIGP\)\.For the online optimization phase, the CIGP inference entails numerical integration of the mean functionfp​h​y​sf\_\{phys\}\. Although this increases the cost per iteration to≈1\.2\\approx 1\.2seconds \(compared to<0\.1<0\.1seconds for a standard GP\), this overhead is small relative to laboratory reaction and analysis times\. The full PC\-MCMC\-CIGP loop should nevertheless be viewed as a two\-stage workflow: PC\-MCMC is an offline structure\-discovery step, while CIGP acquisition is the online step used between experiments\.