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).
Department of Machine Learning and Neural Computing, Donders Institute for Brain, Cognition and Behaviour, Radboud University, Nijmegen, the Netherlands
Data Science and Analytics Institute, University of Oklahoma, Norman, OK, USA · School of Industrial and Systems Engineering, University of Oklahoma, Norman, OK, USA