Source-linked AI summary
FluxDisco: Symbolic Regression for Stoichiometric Dynamical Systems via Monte Carlo Graph Search
Cassandra Durr, Alvaro Köhn-Luque, Chris Jewell, Lloyd A. C. Chapman
TL;DR
Existing dynamical symbolic regression methods must balance predictive accuracy with interpretable equations while avoiding physically invalid expressions. FluxDisco jointly discovers fluxes under known stoichiometry using an adapted Monte Carlo Graph Search procedure. Across physical and biological systems, it recovers governing dynamics accurately, though joint flux estimation scales poorly to high-dimensional, complex systems.
Problem
Existing approaches can sacrifice interpretability or produce expressions that violate known physical laws, while assuming equation forms may impose restrictive and incorrect dynamics.
Method
FluxDisco jointly estimates flux expressions under known stoichiometry and adapts Monte Carlo Graph Search for flux discovery.
Results
FluxDisco accurately recovers governing dynamics across diverse physical systems and varying observation noise, outperforming benchmarked methods in trajectory reconstruction and coupled-dynamics identification.
Takeaways & Limitations
Known stoichiometry enables physically adherent, interpretable equation discovery for flux-based systems where existing approaches fall short.
Takeaways & Limitations
Joint flux estimation scales poorly to high-dimensional, complex ODE systems because the probability of sampling a complete correct flux set diminishes rapidly.
Abstract
from arXiv · showhide
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.
1 Introduction
Dynamical symbolic regression learns interpretable differential equations from trajectories, but unconstrained approaches can misrepresent dynamics or violate physical laws. FluxDisco addresses this by jointly discovering fluxes under known stoichiometry using a graph-based search formulation.
- Assuming governing-equation forms in advance can impose restrictive, incorrect assumptions and produce inaccurate forecasts when the form is misspecified.
- Black-box methods can predict dynamical evolution accurately but do not provide explicit governing equations, limiting interpretability.
- Dynamical symbolic regression infers differential-equation forms and parameters directly from observed trajectories.
- FluxDisco assumes known stoichiometry, jointly estimates unknown flux expressions, reduces the expression search space, and promotes physical adherence.
- Flux discovery formulation: The framework formulates flux discovery as a search over complete or partial flux expressions, expanded through context-free grammar actions and evaluated against observed data.
- Search strategy: MCGS improves on tree search by merging states reachable through distinct action sequences, allowing information to propagate across shared graph nodes.
2.1 Search Algorithm
FluxDisco adapts Monte Carlo Graph Search to jointly discover flux expressions for stoichiometric ODE systems. Its workflow combines graph-based expression search, ODE-based reward evaluation, parsimony control, and interval rules designed to focus exploration.
- Search workflow: FluxDisco treats joint flux discovery as an episodic graph search that builds complete expressions from blank states through grammar rules.Each episode traverses partial or complete flux-expression states and ends when expressions are terminal or the computational budget is exhausted.
- Search workflow: The workflow expands and selects graph nodes, rolls out partial fluxes, solves their ODE systems, and evaluates rewards using parsimony and goodness-of-fit.Reward information is then propagated to ancestor nodes to guide subsequent episodes.
- Selection: FluxDisco selects nodes using local and future reward while retaining the maximum sampled reward to avoid prematurely soft-pruning high-performing states.This modifies average-reward selection for equation discovery, where one optimal functional form is sought rather than average performance.
- Selection: Beta credible intervals replace conservative BKL confidence intervals to differentiate states more effectively in expansive equation-discovery search spaces.The authors report fewer Beta bounds remaining at their initialized value, while acknowledging increased risk that the true mean lies outside the intervals.
- Conditional expansion: Completing expressions with terminal rules after κ grammar rules preserves a dense reward signal and avoids leaving most nodes unevaluated when many non-terminal actions are available.This differs from assigning zero reward to paths that have not reached a terminal state within κ applied rules.
- Reward evaluation: ODE-based reward calculation improves robustness to noise and preserves the original data as the search target, but requires ODE solving and iterative multi-dimensional constant optimisation.The computational burden is presented as the cost of avoiding noise-amplifying finite differencing and smoothing-related distortion.
- Reward evaluation: The reward penalizes non-parsimonious equations through ηC, where η controls the penalty and C is the maximum complexity of the individual fluxes.A flux’s complexity counts the operations and substitutions required to construct its mathematical representation.
2.2 Simulation Study
The simulation study benchmarks FluxDisco on diverse stoichiometric systems and noise conditions, evaluating numerical accuracy, coupled-expression recovery, runtime, and efficiency against existing equation-discovery methods.
- Benchmark systems: FluxDisco is evaluated on simulated flux-based systems spanning epidemic, ecological, chemical, and biological dynamics, including oscillations, bifurcations, and fast–slow behavior.Data are generated under noiseless, low-noise, and high-noise conditions.
- Search constraints: Known stoichiometries and system-specific grammar exclusions constrain candidate flux expressions using prior physical or mechanistic knowledge.Examples include excluding chemically or epidemiologically irrelevant reactants from selected fluxes.
- Evaluation metrics: Numerical evaluation compares integrated candidate trajectories with noiseless ground-truth trajectories using NMSE relative to a constant-mean baseline.A naïve mean prediction corresponds to NMSE = 100%.
- Method comparison: The study compares FluxDisco with sparse, symbolic, and deep-learning equation-discovery methods using numerical, symbolic, runtime, and composite efficiency measures.FluxDisco is the only compared method described as imposing known stoichiometries while estimating system fluxes.
- Evaluation metrics: Symbolic evaluation measures whether top-ranked equations contain the relevant coupled terms with sufficiently accurate coefficients, without requiring full-system identification.The coupled-expression criterion is intentionally more lenient because it ignores uncoupled expressions and does not penalize extra terms.
3 Results
FluxDisco closely reconstructs system trajectories and achieves strong numerical and coupled-dynamics identification performance across the benchmark, while its accuracy requires substantially greater computational time.
- Trajectory reconstruction: Estimated trajectories closely align with noiseless and observed trajectories across all systems and regimes, including noisy conditions.The comparison is presented in Figure 6.
- Numerical performance: The maximum observed trajectory reconstruction error remains below 2% across all systems, although errors generally increase with noise.Fairen-Velarde has the highest errors, followed by the unstable Brusselator.
- Ablations: State merging and grammar-rule exclusions improve numerical accuracy and efficiency, but their gains diminish as the search budget increases.The ablations also show benefits from imposing known stoichiometry for numerical and symbolic performance.
- Numerical performance: FluxDisco produces the lowest errors across all systems (< 2%) and remains well below the baseline, whereas the non-stoichiometric variant performs considerably worse.Other methods exceed the baseline in some experiments, sometimes by orders of magnitude.
- Symbolic performance: FluxDisco achieves high coupled-expression identification accuracy except on Fairen-Velarde, which remains entirely unidentified by every method.Removing imposed stoichiometry reduces identification performance to the level of the remaining benchmarked systems.
- Efficiency: FluxDisco has the lowest error and highest identification rate but runtime exceeding alternative methods by orders of magnitude.The paper frames this as a speed–accuracy tradeoff for an offline equation-discovery task.
4 Discussion
FluxDisco achieves strong identification and reconstruction performance, but its computational cost and scaling limitations define important practical boundaries. Known stoichiometry improves efficiency, while complex systems remain difficult for exact equation recovery.
- Performance trade-offs: FluxDisco produces the lowest trajectory reconstruction errors and highest identification rate of coupled dynamics, but has the longest runtimes.The authors frame this as a trade-off because equation discovery is typically an offline task.
- Performance trade-offs: The efficiency score balances NMSE error and runtime equally across ODE systems and noise levels.Figure 10 compares runtime and efficiency across benchmarked models.
- Scope and limitations: On the complex Fairen-Velarde system, FluxDisco fits trajectories but does not precisely recover the ground-truth expressions.None of the benchmarked approaches successfully recovered its governing equations, reflecting the difficulty of exact discovery in complex systems.
- Scope and limitations: Scaling to high-dimensional, complex ODE systems is poor because joint flux estimation requires correct expressions for all fluxes within one graph node.As flux number or complexity increases, the probability of sampling a complete correct set diminishes rapidly.
- Scope and limitations: The hard stoichiometric constraint makes flux estimation inaccurate when the supplied stoichiometry is misspecified.The framework is designed for settings where stoichiometry is known a priori.
5 Conclusion
The paper addresses physically meaningful equation discovery for noisy, multidimensional stoichiometric systems by adapting MCGS to jointly estimate fluxes under known stoichiometry. FluxDisco recovers governing dynamics across physical systems and noise levels, while future work targets partial observability and scalability.
- Motivation: Existing dynamical symbolic regression methods are generally unsuitable for the constraints of flux-based stoichiometric systems.The conclusion frames the problem as extracting physically meaningful governing equations from noisy, multidimensional data.
- Approach: FluxDisco jointly estimates system fluxes while imposing known stoichiometry through a specialized adaptation of the MCGS algorithm.This design targets physically adherent governing equations.
- Findings: FluxDisco accurately recovers governing dynamics across diverse physical systems and varying observation-noise levels.The reported evaluation emphasizes robustness across systems and noise conditions.
- Future work: Future work should address full-system observability and the scalability challenges of joint flux estimation.The paper identifies latent states and high-dimensional complex ODE systems as practical targets for extension.
S1 Methodology Details
The methodology adapts stochastic MCGS for joint flux discovery by replacing sparse terminal rewards with rollout-based estimates and reward-sensitive uncertainty bounds. Thompson sampling is used to make node selection more responsive in large search spaces.
- Method design: The methodology section specifies the stochastic MCGS formulation, credible intervals for mean rewards, and the complete reward calculation procedure.These components support the framework’s search and evaluation process.
- Method design: FluxDisco uses stochastic rather than deterministic MCGS because partial flux expressions cannot receive meaningful data-error evaluations.Deterministic MCGS therefore creates sparse rewards concentrated at terminal complete expressions.
- Method design: In deterministic MCGS, optimistic upper bounds dominate selection when non-terminal states receive zero reward, weakening guidance from discounted terminal rewards.This produces an undirected search in sparse-reward environments.
- Method design: The stochastic variant completes partial expressions with a rollout policy to estimate potential reward, providing dense sampled rewards for navigating the search space.These estimates help distinguish neighboring nodes before expressions are complete.
- Reward uncertainty: BKL confidence bounds can remain near their initial values in large state-action spaces, causing little differentiation between nodes.The issue arises because states are rarely revisited and the relevant sampling ratio remains high.
- Reward uncertainty: Thompson sampling generalizes to scalar bounded rewards by updating a Beta-based pseudo-posterior using each reward as fractional success and failure evidence.Credible intervals are then obtained from posterior quantiles.
- Reward uncertainty: Beta credible intervals are immediately responsive to sampled rewards and do not rely on the total-to-local graph-visit ratio.This motivates their use in large search spaces instead of BKL bounds.
- Implementation: Algorithm S1 gives the complete procedure for evaluating rewards after rollout for sampled candidate flux expressions.The procedure evaluates a sampled set of candidate flux expressions post-rollout.
S2 Hyperparameters
The supplementary material records FluxDisco’s hyperparameter choices in Table S1.
- Hyperparameters: FluxDisco’s hyperparameter choices are listed in Table S1.
S3 Data Simulation
The study simulates noisy trajectories for benchmark ODE systems and evaluates candidate flux expressions by combining them with known stoichiometry, optimizing constants, and measuring error. The simulations enforce physical constraints and focus on a challenging single-realisation setting.
- Data simulation: Simulated data are generated for the benchmark ODE systems, with parameters documented in the supplementary simulation materials.The simulations include initial conditions, constants, time spans, sampling frequencies, and noise levels.
- Flux evaluation: Candidate flux expressions are folded, assigned a complexity measure, joined through the stoichiometry matrix, and optimized with fitted constants.The resulting reward and error calculations use the optimized system equations and sum of squared errors against observations.
- Noise generation: Noise is sampled from a standard normal distribution, scaled by each variable’s range, and added to noiseless trajectories.When noise is absent, the procedure returns clean trajectories directly.
- Constraint processing: Negative noisy values are set to zero because all benchmark systems require non-negative observations, while SIR data are additionally normalized to sum to one at each time step.The non-negativity correction introduces a slight positive bias but is used to preserve physical realism.
- Experimental design: The experiments use one initial condition per system, with SIR simulations using tailored parameters and 80 time points rather than the 150 used for other systems.The SIR setup was chosen to produce comparable infection peaks and capture a full epidemic without an excessively long tail.
S4 Ablation Studies
The ablation studies examine state merging and grammar-rule exclusions in FluxDisco. Graph search improves efficiency and error at lower budgets, while grammar exclusions trade faster runtime for initially poorer solutions; both comparisons converge with sufficient search.
- State merging: Graph search achieves higher median efficiency and lower median NMSE than tree search without state merging, especially at lower episode counts.The advantage is attributed to improved sample efficiency and wider information propagation, while the methods converge as the search budget increases.
- State merging: The tree-search ablation shows a steep early increase in missed state merges, whereas graph search propagates information across shared states.Fairen-Velarde is excluded because its distinct hyperparameters prevent a uniform comparison.
- State merging: State merging adds runtime overhead, but the median runtime difference is marginal and tree search can take longer on some systems at long horizons.Lotka-Volterra and the unstable Brusselator are examples where long-run tree-search runtimes exceed graph-search runtimes.
- Grammar exclusions: Grammar-rule exclusions create more terminal nodes and shorter selection paths, causing episodes to terminate earlier and reducing runtime.The excluded rules are terminal actions, so uniform sampling reaches completed equation sets more quickly.
- Grammar exclusions: Without grammar exclusions, efficiency is lower because poorer candidates are identified, but the two approaches converge after the ablation finds sufficiently low-error solutions.The faster ablation runtime therefore trades search speed for early solution quality.