Scientific Machine Learning of Chaotic Systems Discovers Governing Equations for Neural Populations

TL;DR

This work presents PEM-UDE, a framework that synergistically combines the prediction-error method with universal differential equations to learn interpretable, symbolic governing equations from chaotic data, even when observations are noisy or partial. By first fitting a UDE under PEM and then extracting symbolic forms, the method yields analytically tractable models that generalize across brain regions and connectivity patterns. The authors demonstrate success on chaotic physical systems (Rössler, Petrzela-Polak circuit) and on neural-population dynamics under sparse connectivity, predicting emergent relationships between connectivity, oscillation frequency, and synchrony that align with intracranial data. The approach advances mechanistic, multi-scale brain modeling by linking microscale neuronal parameters to macroscale population behavior while offering robust handling of noise and partial observability. This has broad implications for constructing interpretable, generalizable models of neural dynamics across diverse architectures.

Abstract

Discovering governing equations that describe complex chaotic systems remains a fundamental challenge in physics and neuroscience. Here, we introduce the PEM-UDE method, which combines the prediction-error method with universal differential equations to extract interpretable mathematical expressions from chaotic dynamical systems, even with limited or noisy observations. This approach succeeds where traditional techniques fail by smoothing optimization landscapes and removing the chaotic properties during the fitting process without distorting optimal parameters. We demonstrate its efficacy by recovering hidden states in the Rossler system and reconstructing dynamics from noise-corrupted electrical-circuit data, in which the correct functional form of the dynamics is recovered even when one of the observed time series is corrupted by noise 5x the magnitude of the true signal. We demonstrate that this method can recover the correct dynamics, whereas direct symbolic regression methods, such as STLSQ, fail to do so with the available data and noise. Importantly, when applied to neural populations, our method derives novel governing equations that respect biological constraints such as network sparsity - a constraint necessary for cortical information processing yet not captured in next-generation neural mass models - while preserving microscale neuronal parameters. These equations predict an emergent relationship between connection density and both oscillation frequency and synchrony in neural circuits. We validate these predictions using three intracranial electrode recording datasets from the medial entorhinal cortex, prefrontal cortex, and orbitofrontal cortex. Our work provides a pathway to develop mechanistic, multi-scale brain models that generalize across diverse neural architectures, bridging the gap between single-neuron dynamics and macroscale brain activity.

Figures (11)

• Figure 1: Schematic of the approach presented here. We begin with observations of a physical or biological system that have some time-varying activity (data). From these, we adapt prior work to form an initial description of the system under observation (known constraints). We then fit a UDE to describe the unknown portion of the dynamics that are not described by this initial estimate of the physical system, adapting the prediction-error method to assist in training UDEs on chaotic systems (learn dynamics). Finally, we learn a symbolic form of the dynamics described by the UDE (generate equations) and show that this symbolic form makes correct predictions of future activity beyond the initial observations in the system (validate).
• Figure 2: The combined PEM-UDE approach is able to learn the dynamics of chaotic systems that traditional UDEs are unable to capture. a. Comparison of training error for the UDE system (fails to accurately capture dynamics) to the PEM-trained UDE (succeeds). The difference in accuracy is already stark 10% into the full training session (inset), where unassisted UDEs are not progressing. b. The application of PEM smooths out an otherwise intractable parameter loss landscape, with the hyperparameter $K$ tuning the steepness of the landscape. c. Comparison of the mean-square-error of the fit dynamics of the UDE to the PEM-trained UDE. d. Using STLSQ to replace the PEM-trained UDE results in an expression that not only captures the original dynamics (first 200 time points) but continues to accurately describe them well into the future (last 100 time points). e. The surface of the Rössler attractor demonstrates chaotic filling of phase space, with a basin near the x-y plane and occasional bursts along the z-axis. f. Relative error of the PEM-trained UDE at each point of the attractor surface.
• Figure 3: The PEM-UDE approach is capable of recovering the correct form of a chaotic circuit even in the presence of significant observational noise. a. The analog Petrzela-Polak circuit that produces the dynamics examined in this section. The observed states are the dimensionless forms of $V_1$, $V_2$, and $I$, where $V_2$ and $I$ (green boxes) are observed with only a small amount of noise, and $V_1$ (blue box) is observed with large observational noise. b. In this example system, we assume a faulty observation of one state ($V_1$), where the true dynamics (top) are masked by observational noise of 5x the true signal's magnitude (bottom). c. Original attractor surface, with true form of the observed dimension $x$ in color with magnitude of error from observational noise shown. d. Attractor surface that is fit by the PEM-UDE method with noise in the third dimension (notice the scattered colors and the noise has significantly corrupted the dynamics of the state $x$). e. The true form of the state to be recovered is shown at the top, with the form recovered by STLSQ applied to the PEM-UDE dynamics is shown at the bottom. The functional form - although not the exact parameters - is recoverd from the PEM-UDE fit (table compares parameters fit by STLSQ to original parameters).
• Figure 4: Learning novel neural mass models to account for regional heterogeneity of neural connection sparsity allows for broader insights into the origin of different frequency bands within the brain. Single neurons produce spike trains (a), which are encoded as a raster plot of individual dots for each neuron/spike. In a population, this allows for easy viewing of the group activity, which in synchronized regimes has clearly visible bursts (b). In larger behavioral experiments, electrode arrays can be placed to record population-level activity from multiple brain regions (c). The relative degree of local connectivity within brain regions varies across the brain (d), with deeper regions typically exhibiting greater connectivity and more cortical regions having sparse connections within regions. To account for this sparsity, we train a novel next-generation neural mass model to learn terms that can the emergent population dynamics as connections become more sparse (e). The PEM-UDE learned dynamics exhibit phenomena that are typical of real brain regions, with a shift from theta to alpha as the dominant frequency as network connections become increasingly sparse (f). Increased sparsity also constrains the synchronization that is possible within the network (g). Figure created with BioRender.com.
• Figure 5: Intracranial EEG shows synchrony variation across regions as predicted by symbolic PEM-UDE. As shown in Fig. \ref{['fig:4']}g, the typical synchrony in regions with greater degrees of connectivity will be higher than that in regions with lower degrees of connectivity. We also show that this holds in biological brains by examining three datasets. In the one with greatest connectivity (rat medial entorhinal cortex; left column) there is a greater degree of synchrony (measured by phase-locking value; PLV) over the duration of the experiment compared to the prefrontal cortex (middle panel). This difference is even more noticeable compared to the human orbitofrontal cortex at rest (right panel), where overall synchronization tends to be very low during spontaneous activity. As all cumulative distributions (upper row) are based on synchrony values from many individuals, synchrony distributions from each subject are also shown (bottom row) to illustrate the variability across individuals. Figure created with BioRender.com.