Modeling real-world systems requires accounting for noise - whether it arises from unpredictable fluctuations in financial markets, irregular rhythms in biological systems, or environmental variability in ecosystems. While the behavior of such systems can often be described by stochastic differential equations, a central challenge is understanding how noise influences the inference of system parameters and dynamics from data. Traditional symbolic regression methods can uncover governing equations but typically ignore uncertainty. Conversely, Gaussian processes provide principled uncertainty quantification but offer little insight into the underlying dynamics. In this work, we bridge this gap with a hybrid symbolic regression-probabilistic machine learning framework that recovers the symbolic form of the governing equations while simultaneously inferring uncertainty in the system parameters. The framework combines deep symbolic regression with Gaussian process-based maximum likelihood estimation to separately model the deterministic dynamics and the noise structure, without requiring prior assumptions about their functional forms. We verify the approach on numerical benchmarks, including harmonic, Duffing, and van der Pol oscillators, and validate it on an experimental system of coupled biological oscillators exhibiting synchronization, where the algorithm successfully identifies both the symbolic and stochastic components. The framework is data-efficient, requiring as few as 100-1000 data points, and robust to noise - demonstrating its broad potential in domains where uncertainty is intrinsic and both the structure and variability of dynamical systems must be understood.
Figures & tables
Figure 1: Overview of our methodology for discovering SDEs with structured diffusion. Here the methodology is applied to a Duffing oscillator, from data acquisition to model optimization. A Collect the raw data. B The time series data is fed into a Gaussian process to extract the deterministic dynamics; GPs can robustly differentiate the data, as shown on the acceleration as compared to a finite finite difference scheme. C.1-2 The preprocessed data is fed to DSR to discover the mathematical expression of the drift. C.3-4 We build a custom dataset Θ from the basis functions of the discovered drift and extract the noise data using the Euler–Maruyama approximation. C.5 . We assume the noise is the sum of contributions from different components, each represented by a basis function from Θ with its noise variance. These noise variances are optimized using MLE to best fit the stochastic process data. After optimization, variances below a certain threshold are removed to obtain a more parsimonious model without compromising predictive accuracy. D . The estimated drift and diffusion functions μ^,σ^ fully determined.
Figure 2: Linear oscillator with additive noise. Predicted coefficients for the drift (left column) and diffusion (right column) as a function of noise level. Red lines indicate the model predictions, while black dashed lines denote the ground truth. Error bars (where shown) correspond to the 2.5%–97.5% quantiles across realizations, reflecting the framework’s uncertainty when applied to five independent datasets obtained from integrations of the same system. The framework consistently preserves structural accuracy across noise levels and yields precise estimates of the diffusion amplitude, although spurious diffusion activations appear at higher noise intensities. In the diffusion panel, C denotes the additive noise constant σa (see ( 7 )).
Figure 3: Multiplicative frequency noise. Estimation of the natural frequency distribution for the linear system with 5% frequency noise. Left: Ground truth frequency variations (blue) compared with estimated variations (black) using the Euler–Maruyama approximation with Δt=0.01 . Right: Ground truth probability density function of the frequency parameter (blue), approximated distribution from the stochastic process data (black dashed), and MLE-predicted distribution (red shaded).
Figure 4: Linear and Duffing oscillators with noise in all parameters. Violin plots show distributions of recovered coefficients at 1% and 5% noise. Blue bars denote ground truth, while red bars indicate predictions. Each distribution is modeled as Gaussian, with the mean corresponding to the drift and the standard deviation to the diffusion of the coefficients. In both panels, C denotes the additive noise constant σa (see ( 8 )).
Figure 5: The experiment setup. Microscope image of two E. coli bacteria trapped in separate cavities connected by a narrow channel [ 40 ] .
Figure 6: Comparison between experimental data and discovered models for two synchronization datasets. Left: Violin plots show the distributions of coupling strength ( k ) and angular velocity ( ω ) for the two oscillators. Blue bars represent empirical estimates obtained with the method of Japaridze et al. [ 40 ] , while red bars indicate predictions from our algorithm. Each distribution is modeled as Gaussian, with the mean corresponding to the drift and the standard deviation to the diffusion of the coefficients. Right: Cosine of the phase difference (ϕ1−ϕ2) over time, comparing simulations of the discovered model with experimental data. The gray-shaded region marks the interval during which synchronization occurs.
Appendix figures & tables10 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 7: The extracted noisy process (blue dots) is plotted against the fit of our tailored GP (shaded grey area). Our method is capable of fitting this heteroscedastic noise with remarkable precision.
Figure 8: Schematic overview of Deep Symbolic Regression (DSR). A recurrent neural network (RNN) generates symbolic expressions token by token, building expression trees from a library of primitives. Each expression is evaluated on the input data and assigned a reward; the RNN is updated via a risk-seeking policy gradient to favour the highest-reward expressions.
Soft Length Prior
Entropy
library
epochs
location
scale
w
γ
learning rate
reward
Linear Systems
{+,−,×,C}
50
4
2
0.02
0.85
5×10−4
neg-mse
Nonlinear Systems
{+,−,×,C}
100
14
2
0.02
0.85
5×10−4
neg-mse
Biological Oscillator
Appendix
Table 1: Hyperparameter settings used for each of the problems studied.
Figure 9: Results from applying the GP-DSR framework to four benchmark SDEs under broad DSR regularization settings (soft length-prior location: 20 tokens, scale: 10 tokens; noise level: 5%). The equations at the top of each panel are the DSR-predicted drift μ^ and diffusion σ^ ; the black curve shows the GP-predicted mean trajectory; the blue dots are the extracted stochastic residuals; and the grey band is the GP posterior uncertainty. No reference trajectory is superimposed because the predicted symbolic form differs structurally from the true equation. (A) Linear oscillator with additive noise; (B) Duffing oscillator with noise affecting all parameters; (C) Linear oscillator with noise affecting all parameters; (D) Van der Pol oscillator with noise on both the nonlinear and linear damping terms.
Figure 10: Bar plots showing the occurrence frequency of each diffusion token in the model predictions, for various systems and noise levels. ( A ) Linear oscillator with frequency noise, ( B ) Van der Pol oscillator with noise in linear and nonlinear damping, and ( C ) Linear oscillator with additive noise. Each bar is computed from 5 independent realizations for each noise level. The framework consistently captures the correct tokens across all realizations; minor spurious activations are observed in some realizations at higher noise levels.
Figure 11: Frequency noise estimation comparison across various time steps Δt (illustrated in red) against the ground truth distribution (in blue).
Figure 12: σL (Lorenz coupling parameter) noise estimation comparison across various time steps Δt (illustrated in red) against the ground truth distribution (in blue).
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.
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
Discovering governing equations from observational data remains a fundamental challenge in scientific modeling, particularly when the underlying mathematical structure is unknown. Traditional sparse identification methods like SINDy excel at discovering parsimonious models but require researchers to specify candidate basis functions a priori, a limitation that often leads to model failure when critical terms are omitted or when systems exhibit unconventional dynamics. Purely symbolic regression approaches offer unlimited flexibility but struggle with noise sensitivity and frequently produce overly complex, unstable equations. We present AutoSINDy, a hybrid Discovery-then-Solve framework that combines the exploratory power of symbolic regression with the robust sparsity-promoting capabilities of SINDy. Our method operates in three stages: (1) PySR-based symbolic regression discovers candidate functional forms from bootstrapped data chunks; (2) a curation pipeline decomposes, expands, and filters these expressions using collinearity analysis to construct a minimal yet comprehensive library; and (3) SINDy identifies sparse governing equations from this custom-tailored library. Extensive experiments across canonical nonlinear systems demonstrate that AutoSINDy consistently recovers ground-truth equations even under high observational noise, achieving a ground-truth recovery rate of 92.8% across all trials. Compared with standard SINDy using enriched libraries and standalone symbolic regression, AutoSINDy achieves higher predictive accuracy, superior generalization to unseen trajectories, and substantially lower symbolic complexity.
Mohammad Amin Basiri, Charles Nicholson
Data Science and Analytics Institute, University of Oklahoma, Norman, OK, USA · School of Industrial and Systems Engineering, University of Oklahoma, Norman, OK, USA
Model-based approaches for (bio)process systems often suffer from incomplete knowledge of the underlying physical, chemical, or biological laws. Universal differential equations, which embed neural networks within differential equations, have emerged as powerful tools to learn this missing physics from experimental data. However, neural networks are inherently opaque, motivating their post-processing via symbolic regression to obtain interpretable mathematical expressions. Genetic algorithm-based symbolic regression is a popular approach for this post-processing step, but provides only point estimates and cannot quantify the confidence we should place in a discovered equation. We address this limitation by applying Bayesian symbolic regression, which uses Reversible Jump Markov Chain Monte Carlo to sample from the posterior distribution over symbolic expression trees. This approach naturally quantifies uncertainty in the recovered model structure. We demonstrate the methodology on a Lotka-Volterra predator-prey system and then show how a well-designed experiment leads to lower uncertainty in a fed-batch bioreactor case study.