FluxDisco: Symbolic Regression for Stoichiometric Dynamical Systems via Monte Carlo Graph Search

Cassandra Durr; Alvaro Köhn-Luque; Chris Jewell; Lloyd A. C. Chapman

stat.MLSep 4, 2026

FluxDisco: Symbolic Regression for Stoichiometric Dynamical Systems via Monte Carlo Graph Search

Authors: Cassandra Durr, Alvaro Köhn-Luque, Chris Jewell, Lloyd A. C. Chapman

Organizations: Lancaster University · University of Oslo

Abstract

Dynamical symbolic regression methods identify governing differential equations from noisy data, balancing interpretability and predictive accuracy. However, standard methods often produce expressions that violate known physical laws. To address this, we propose FluxDisco, a physics-informed framework tailored for flux-based, stoichiometric ODE systems. By leveraging a known stoichiometry, we reduce the expression search space and ensure physical adherence. Our framework adapts the Monte Carlo Graph Search algorithm for the unique challenges associated with joint flux discovery of stoichiometric systems. We evaluate our method across a range of physical and biological systems, demonstrating its ability to accurately recover governing dynamics through interpretable equations.

Figures & tables

Panel of Figure 1 — Figure 1: Flow diagram of the susceptible-infectious-recovered (SIR) epidemic model with state variables $\mathbf{x}(t)$ and fluxes $\mathbf{v}(\mathbf{x}(t))$ .

Panel of Figure 3 — Figure 3: The four phases of the MCTS and MCGS algorithm applied to symbolic regression. Unlike MCTS, MCGS allows nodes to be reached via multiple paths, enabling more efficient information propagation across the search structure. For simplicity, a single flux expression per state is illustrated. Highlighted nodes indicate the path traversed during a single iteration of the algorithm. Within-node text shows the flux expression, between-node text shows the applied grammar rule, and arrow directions indicate information flow.

Panel of Figure 4 — Figure 4: Workflow of FluxDisco, our modified MCGS algorithm.

Figure 5: BKL confidence intervals compared with Beta credible intervals (Equation 7 ). Both panels illustrate that the more frequently a node is visited, the more its intervals contract towards the reference line, $l(s)=u(s)$ .

Panel of Figure 6 — Figure 6: Illustration of the complexity calculation of a flux, where the flux expression is visualised as an expression tree. Blue nodes represent non-terminal tokens, orange nodes represent applied terminal actions, and grey nodes represent applied non-terminal actions.

ODE System	$\dot{\mathbf{x}}$	Number of coupled expressions
SIR (Standard)	$\begin{bmatrix}{\color[rgb]{0.043,0.2305,0.5547}-c_{0}x_{0}x_{1}}\\ {\color[rgb]{0.043,0.2305,0.5547}c_{0}x_{0}x_{1}}{\color[rgb]{0.8945,0.4141,0.0977}-c_{1}x_{1}}\\ {\color[rgb]{0.8945,0.4141,0.0977}c_{1}x_{1}}\end{bmatrix}$	2
SIR (Squared)	$\begin{bmatrix}{\color[rgb]{0.043,0.2305,0.5547}-c_{0}x_{0}x_{1}^{2}}\\ {\color[rgb]{0.043,0.2305,0.5547}c_{0}x_{0}x_{1}^{2}}{\color[rgb]{0.8945,0.4141,0.0977}-c_{1}x_{1}}\\ {\color[rgb]{0.8945,0.4141,0.0977}c_{1}x_{1}}\end{bmatrix}$	2
SIR (Square root)	$\begin{bmatrix}{\color[rgb]{0.043,0.2305,0.5547}-c_{0}x_{0}\sqrt{x_{1}}}\\ {\color[rgb]{0.043,0.2305,0.5547}c_{0}x_{0}\sqrt{x_{1}}}{\color[rgb]{0.8945,0.4141,0.0977}-c_{1}x_{1}}\\ {\color[rgb]{0.8945,0.4141,0.0977}c_{1}x_{1}}\end{bmatrix}$	2
Lotka-Volterra	$\begin{bmatrix}c_{0}x_{0}{\color[rgb]{0.043,0.2305,0.5547}-c_{1}x_{0}x_{1}}\\ {\color[rgb]{0.043,0.2305,0.5547}c_{1}x_{0}x_{1}}-c_{2}x_{2}\end{bmatrix}$	1
Brusselator	$\begin{bmatrix}{\color[rgb]{0.043,0.2305,0.5547}x_{0}^{2}x_{1}}{\color[rgb]{0.8945,0.4141,0.0977}-c_{0}x_{0}}+c_{1}-x_{0}\\ {\color[rgb]{0.043,0.2305,0.5547}-x_{0}^{2}x_{1}}{\color[rgb]{0.8945,0.4141,0.0977}+c_{0}x_{0}}\end{bmatrix}$	2
Fairen-Velarde	$\begin{bmatrix}{\color[rgb]{0.043,0.2305,0.5547}-\frac{x_{0}x_{1}}{1+c_{0}x_{0}^{2}}}+c_{1}-x_{0}\\ {\color[rgb]{0.043,0.2305,0.5547}-\frac{x_{0}x_{1}}{1+c_{0}x_{0}^{2}}}+c_{2}\end{bmatrix}$	1

Table 2: Coupled expressions within each ODE system considered

Method	F.D.	Joint Disc.	Fluxes	Category	Ref
Ours	✗	✓	✓	MCGS
Ours (No Stoic.)	✗	✓	✗	MCGS
SPL	✓	✗	✗	MCTS	[ 23 ]
SINDy	✓	✓	✗	Sparse regression	[ 7 , 18 ]
PySR	✓	✗	✗	Genetic programming	[ 10 ]
ProGED	✗	✓	✗	Monte Carlo, grammar-based	[ 22 ]

Table 3: Overview of methods included in our comparative study. Traits include whether finite differencing of state trajectories are required (F.D.) , whether the method estimates all ODEs in a system jointly or independently (Joint Disc.) , whether the method accounts for the flux-based nature of the dynamical systems (Fluxes) , and what general category the method falls under.

Figure 7: Comparison of estimated, noiseless, and observed trajectories.

Figure 8: Normalised MSE (NMSE) across all benchmarked systems and levels of noise. Errors are averaged across fluxes per system and noise level.

Figure 9: Comparison of Normalised Mean Squared Error (NMSE) across benchmarked models. Error lines indicate the range of NMSE across varying noise levels. The vertical line denotes the baseline error associated with predicting the data mean across all time points (NMSE = 100%).

Figure 10: Coupled term identification accuracy in top-predicted equations across benchmarked methods for each ODE system. Results are averaged across all noise levels.

Panel of Figure 11 — Figure 11: Comparison of runtime and efficiency across benchmarked models across ODE systems and noise levels. The efficiency score equally balances error (NMSE) and runtime (seconds).

\hat{\dot{\mathbf{x}}}(t)=\mathbf{S}\hat{\mathbf{v}}_{\text{fold}}(\mathbf{x}(t))

Algorithm S1 Reward calculation

Symbol	Description	Value	Fairen-Velarde Modifications
$\gamma$	Discount factor for future rewards	0.9
$\epsilon$	Accuracy measure for propagation (Algorithm 4, [ 19 ] )	0.01
$\tau$	Temperature in reward function	0.005	Increased to 0.05 to increase reward scale.
$\eta$	Parsimony coefficient in reward function	0.99	Increased to 0.995 to allow for more complex expressions.
$\kappa$	Maximum flux expression depth	6	Increased to 12 to allow for more complex expressions.
	Episodes	100	Decreased to 40 to compensate for increased expression depth.

Table S1 : Selected hyperparameter values

\mathbf{x}(t)=\text{ODESolver}(\dot{\mathbf{x}}(t),\mathbf{x}(0))

Algorithm S2 Data Simulation

System	$\mathbf{x}(0)$	$\mathbf{c}$	Time Span	Time Points	Noise
SIR (Standard)	$[0.98,0.02,0]$	$[0.4,0.1]$	$[0,80]$	$80$	$[0,0.02,0.05]$
SIR (Squared)	$[0.88,0.12,0]$	$[1.3,0.08]$	$[0,80]$	$80$	$[0,0.02,0.05]$
SIR (Square root)	$[0.98,0.02,0]$	$[0.15,0.06]$	$[0,80]$	$80$	$[0,0.02,0.05]$
Lotka-Volterra	$[10,5]$	$[1,0.1,1.5]$	$[0,30]$	$150$	$[0,0.03,0.07]$
Brusselator (Stable)	$[0.5,2]$	$[1.5,1]$	$[0,30]$	$150$	$[0,0.03,0.07]$
Brusselator (Unstable)	$[0.5,2]$	$[0.5,2]$	$[0,30]$	$150$	$[0,0.03,0.07]$

Table S2 : Simulation parameters for the benchmarked ODE systems, including initial conditions

\mathbf{x}(0)

, constant values

\mathbf{c}

, time spans, sampling frequencies, and noise levels

\epsilon

Figure S1 : Graph search vs. tree search ablation study results. The left panels show median scores aggregated across all conditions, and the right panels break down median scores by ODE system.

Figure S2 : Search results with grammar rule exclusions vs. ablation without exclusions. The left panels show median scores aggregated across all conditions, and the right panels break down median scores by ODE system.

Explore similar work

Symbolic Discovery of Stochastic Differential Equations with Genetic Programming
Mar 10, 2026Sigur de Vries, Sander W. Keemink, Marcel A. J. van GervenStochastic Differential Equations Symbolic Regression
Department of Machine Learning and Neural Computing, Donders Institute for Brain, Cognition and Behaviour, Radboud University, Nijmegen, the Netherlands
Data-driven Progressive Discovery of Physical Laws
Mar 14, 2026Mingkun Xia, Weiwei ZhangScientific Discovery Symbolic Regression
School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China · International Joint Institute of Intelligent Fluid Mechanics, Northwestern Polytechnical University, Xi’an 710072, China · National Key Laboratory of Aircraft Configuration Design, Xi’an 710072, China
LLM-ODE: Data-driven Discovery of Dynamical Systems with Large Language Models
Mar 21, 2026Amirmohammad Ziaei Bideh, Jonathan GryakODE Discovery Dynamical Systems
The Graduate Center, CUNY New York, New York, USA · Queens College and The Graduate Center, CUNY New York, New York, USA

Mar 10, 2026J/K move · Enter open · S save

Symbolic Discovery of Stochastic Differential Equations with Genetic Programming

Sigur de Vries, Sander W. Keemink, Marcel A. J. van Gerven

Department of Machine Learning and Neural Computing, Donders Institute for Brain, Cognition and Behaviour, Radboud University, Nijmegen, the Netherlands

Stochastic Differential Equations Symbolic Regression

Open paper AskPDF arXiv

Automated scientific discovery aims to improve scientific understanding through machine learning. A central approach in this field is symbolic regression, which uses genetic programming or sparse regression to learn interpretable mathematical expressions to explain observed data. Conventionally, the focus of symbolic regression is on identifying ordinary differential equations. The general view is that noise only complicates the recovery of deterministic dynamics. However, explicitly learning a symbolic function of the noise component in stochastic differential equations enhances modelling capacity, increases knowledge gain and enables generative sampling. We introduce a method for symbolic discovery of stochastic differential equations based on genetic programming, jointly optimizing drift and diffusion functions via the maximum likelihood estimate. Our results demonstrate accurate recovery of governing equations, efficient scaling to higher-dimensional systems, robustness to sparsely sampled problems and generalization to stochastic partial differential equations. This work extends symbolic regression toward interpretable discovery of stochastic dynamical systems, contributing to the automation of science in a noisy and dynamic world.

Abstract

Explore similar work

Symbolic Discovery of Stochastic Differential Equations with Genetic Programming

Data-driven Progressive Discovery of Physical Laws

LLM-ODE: Data-driven Discovery of Dynamical Systems with Large Language Models