Organizations: School of Mathematics, Southeast University, Nanjing, 210096, P. R. China. · School of Electronic Engineering and Computer Science, Queen Mary University of London, London, E1 4NS, UK. · School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK. · Dipartimento di Fisica ed Astronomia, Universit`a di Catania and INFN, 95123, Catania, Italy. · Complexity Science Hub Vienna (CSHV), Vienna, Austria. · Research Institute of Intelligent Complex Systems and MOE Frontiers Center for Brain Science, Fudan University, Shanghai, 200433, P. R. China. · School of Mathematical Sciences, LMNS, and SCMS, Fudan University, Shanghai, 200433, P. R. China.
Predicting the behavior of ultra-large complex systems, from climate to biological and technological networks, is a central unsolved challenge. Existing approaches face a fundamental trade-off: equation discovery methods provide interpretability but fail to scale, while neural networks scale but operate as black boxes and often lose reliability over long times. Here, we introduce the Sparse Identification Graph Neural Network, a framework that overcome this divide by allowing to infer the governing equations of large networked systems from data. By defining symbolic discovery as edge-level information, SIGN decouples the scalability of sparse identification from network size, enabling efficient equation discovery even in large systems. SIGN allows to study networks with over 100,000 nodes while remaining robust to noise, sparse sampling, and missing data. Across diverse benchmark systems, including coupled chaotic oscillators, neural dynamics, and epidemic spreading, it recovers governing equations with high precision and sustains accurate long-term predictions. Applied to a data set of time series of temperature measurements in 71,987 sea surface positions, SIGN identifies a compact predictive network model and captures large-scale sea surface temperature conditions up to two years in advance. By enabling equation discovery at previously inaccessible scales, SIGN opens a path toward interpretable and reliable prediction of real-world complex systems.
Figures & tables
Figure 1: SIGN pipeline. SIGN facilitates accurate and scalable equation discovery through a two-phase process. a , An example of the network of coupled Rössler oscillators that we aim to reconstruct from its time-series observations. b , Candidate libraries for intrinsic terms F and coupling terms C , each generated from a set of nonlinear basis functions. c , Phase I: Sparse regression on a small subset of nodes, followed by DBSCAN clustering to determine a global support mask. d , Phase II: Message passing on the full graph using the learned mask to estimate the shared coefficients of F and C . e , Masked intrinsic and coupling messages, used to form the recovered equations. f , Coefficient recovery error (sMAPE) across network topologies and sizes. g , Reconstructed versus true trajectories from integrating the inferred equations, showing the model’s predictive accuracy.
Figure 2: Equation discovery across diverse networked dynamical systems using the SIGN framework. a , Kuramoto phase-oscillator network. Here xi(t) denotes the phase of oscillator i , and the interaction between connected oscillators follows the sinusoidal coupling sin(xj−xi) , where j ranges over neighbors of i defined by the network adjacency Aij (i.e., edges indicate which pairs interact). The left panel reports the inference error of the coupling-term coefficient across synthetic scale-free networks with sizes of 103 and 105 nodes, respectively, and two large empirical networks (GitHub and Catster), while the right panel shows an example of true versus reconstructed phase trajectories xi(t) for three representative nodes. b , SIS model on small-world and empirical networks: coefficient inference errors for the epidemic-spreading dynamics across network types and sizes. c , Michaelis–Menten (MM) regulatory model: coefficient inference errors for gene-regulatory dynamics across multiple network topologies and scales. d , FitzHugh–Nagumo model: coefficient inference errors for synthetic simulations and a large-scale empirical brain network. e , Hindmarsh–Rose model: coefficient inference errors across model terms, highlighting the distribution of errors across coefficients.
Figure 3: SIGN supports equation inference under noise, limited observations, complex (out-of-basis) dynamics, and imperfect network information. All experiments are repeated 10 times. Panels report coefficient errors (sMAPE) and representative trajectory reconstructions where applicable. a , Noise robustness: sMAPE under additive Gaussian noise as the signal-to-noise ratio (SNR) decreases (down to 60 dB) across the benchmark systems. b,c , Sampling resolution on ultra-large networks ( 105 nodes): b shows sMAPE as the number of observed time points is reduced (down to 200 ); c shows sMAPE as the sampling interval Δt is varied (with Δt=0.01∼0.05 ). d,e , Non-canonical dynamics (basis mismatch): representative trajectory reconstructions for systems whose generating functions are not contained in the candidate library, including a mutualistic model with fractional nonlinearities and the chaotic Chua circuit; reconstruction error is quantified by mean squared error (MSE) across nodes. f , Phase heterogeneity: Kuramoto oscillators with natural frequencies drawn from N(0.3,σ) ; sMAPE as the dispersion σ increases. g , Strong coupling: FitzHugh–Nagumo networks as coupling strength increases toward synchronized dynamics, quantified by the order parameter ⟨R⟩→1 . h , Structured topologies: sMAPE on low-rank random graph models (RGM) and stochastic block models (SBM) with 1,000 nodes. i , Structural incompleteness: sMAPE under missingness, with up to 20% node removal and 40% edge deletion prior to inference. j , Baselines and scalability: coefficient-error comparison with Two-stage and LaGNA under structural incompleteness, and runtime scaling with network size from 10 to 105 nodes.
Figure 4: Scalable prediction of large-scale FitzHugh–Nagumo neural dynamics under observational noise using SIGN. a , Representative phase-space trajectories obtained by integrating the inferred equations, shown for multiple noise conditions; solid and dashed curves indicate the training and testing intervals, respectively, and nodes are randomly sampled with noise intensities centered around the median. b , Distribution of node-wise prediction error (MSE) at SNR =50 dB (left) and =30 dB (right). c , Two-dimensional density plots comparing true versus predicted values at SNR =50 dB (left) and =30 dB (right); density indicates how frequently pairs fall into the same region. d , Prediction error (MSE) evaluated across temporal segments, with the training and testing intervals indicated. e , Comparison on a 1,000 -node simulated dataset under SNR =50 dB and =30 dB, reporting predictive errors for SIGN and neural-network-based predictors.
Figure 5: SIGN equation inference and integration-based prediction of large-scale sea surface temperature (SST) dynamics. a , Spatial map of node-wise MAPE over the 8-year training period (June 2002–June 2010), summarizing regional prediction error (mean MAPE =1.93% ). b , Spatial map of node-wise MAPE over the 2-year testing period (July 2010–June 2012) (mean MAPE =3.55% ). c , Representative SST trajectories at four randomly selected grid points: observed SST (blue) versus trajectories obtained by integrating the SIGN-inferred equations (purple). Solid and dashed segments indicate training and testing periods, respectively. d , Distribution of MAPE aggregated over all nodes and time steps. e , Two-dimensional density plot comparing predicted versus observed SST values (shown over 24 – 29∘C ); density indicates the frequency of value pairs. f , Relationship between local SST variability (standard deviation) and node-wise prediction error (MAPE), evaluated on 10,000 randomly sampled nodes. g , MAPE as a function of prediction horizon over the full 10-year window, with the training/testing boundary indicated.
Our ability to predict, control, and ultimately understand complex systems rests on discovering the equations that govern their dynamics. Identifying these equations directly from noisy, limited observations has therefore become a central challenge in data-driven science, yet existing library-based sparse regression methods force a compromise between automation, statistical rigor, and computational efficiency. Here we develop Bayesian-ARGOS, a hybrid framework that reconciles these demands by combining rapid frequentist screening with focused Bayesian inference, enabling automated equation discovery with principled uncertainty quantification at a fraction of the computational cost of existing methods. Tested on seven chaotic systems under varying data scarcity and noise levels, Bayesian-ARGOS outperforms two state-of-the-art methods in most scenarios. It surpasses SINDy in data efficiency for all systems and noise tolerance for six out of the seven, with a two-order-of-magnitude reduction in computational cost compared to bootstrap-based ARGOS. The probabilistic formulation additionally enables a suite of standard statistical diagnostics, including influence analysis and multicollinearity detection that expose failure modes otherwise opaque. When integrated with representation learning (SINDy-SHRED) for high dimensional sea surface temperature reconstruction, Bayesian-ARGOS increases the yield of valid latent equations with significantly improved long horizon stability. Bayesian-ARGOS thus provides a principled, automated, and computationally efficient route from scarce and noisy observations to interpretable governing equations, offering a practical framework for equation discovery across scales, from benchmark chaotic systems to the latent dynamics underlying global climate patterns.
Yuzheng Zhang, Weizhen Li, Rui Carvalho
Department of Engineering, Durham University, Lower Mountjoy, South Road, Durham, DH1 3LE, United Kingdom. · College of Control Science and Engineering, Zhejiang University, 38 Zheda Rd.(Yuquan Campus), Hangzhou, 310027, P.R. China. · Institute for Data Science, Durham University, South Road, Durham, DH1 3LE, United Kingdom.
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
Discovering the governing equations of a dynamical system from observed trajectories provides deeper insight into its structure than mere prediction of future states. We present a data-driven approach to model discovery based on complex-valued product-unit networks, in which each unit represents a complex monomial and the network output is a sparse linear combination of such monomials. In contrast to established library-based methods such as SINDy, our approach does not require a predefined set of candidate functions: the relevant monomials, including those with fractional or negative exponents, are learned directly from data. Across four chaotic benchmark systems (Lorenz63, Lorenz84, the Four-Wing attractor, and a fractional variant of Lorenz63), we recover the exact governing equations in 90% of trials for the first three systems, and in 70-90% of trials for the fractional case, using at least 3000 training points. Applied to real-world human-gait accelerometer signals, the model produced stable trajectories with bounded prediction errors, corresponding to an RMSE of approximately 12-14% of the signal amplitude range over a test horizon three times longer than the training interval, demonstrating its potential for high-dimensional systems in which analytic equations are unavailable.
Martin Brückmann, Babette Dellen, Uwe Jaekel
Department of Mathematics, Informatics, and Technology, RheinAhrCampus Remagen University of Applied Sciences Koblenz Joseph-Rovan-Allee 2, 53424 Remagen, Germany