Deep Generative Modeling for Identification of Noisy, Non-Stationary Dynamical Systems

## Abstract

Abstract A significant challenge in many fields of science and engineering is making sense of time-dependent measurement data by recovering governing equations in the form of differential equations. We focus on finding parsimonious ordinary differential equation (ODE) models for nonlinear, noisy, and non-autonomous dynamical systems and propose a machine learning method for data-driven system identification. While many methods tackle noisy and limited data, non-stationarity – where differential equation parameters change over time – has received less attention. Our method, dynamic SINDy, combines variational inference with SINDy (sparse identification of nonlinear dynamics) to model time-varying coefficients of sparse ODEs. This framework allows for uncertainty quantification of ODE coefficients, expanding on previous methods for autonomous systems. These coefficients are then interpreted as latent variables and added to the system to obtain an autonomous dynamical model. We validate our approach using synthetic data, including nonlinear oscillators and the Lorenz system, and apply it to neuronal activity data from C. elegans. Dynamic SINDy uncovers a global nonlinear model, showing it can handle real, noisy, and chaotic datasets. We aim to apply our method to a variety of problems, specifically dynamic systems with complex time-dependent parameters.

## 1 Introduction

Many fields of science and engineering now benefit from unprecedented amounts of data due to increased efforts and technological breakthroughs in data collection. The challenge is to use these measurements to expand our understanding of dynamical systems in areas like climate science, neuroscience, ecology, finance, and epidemiology. Machine learning methods, such as neural networks, are widely used for data-driven modeling, offering high prediction accuracy but limited interpretability. In contrast, traditional techniques that identify ordinary and partial differential equations (ODEs and PDEs) provide interpretable and generalizable insights into the system’s underlying physics. While neural networks may lose accuracy as conditions change, in many systems the governing differential equations remain reliable. The key question is whether we can combine the strengths of deep learning with the clarity and simplicity of data-driven differential equation models.
A key challenge in data-driven system identification is that many systems exhibit nonlinear behavior, such as switching between dynamical regimes [, , , ]. These "hybrid systems" [], where continuous dynamics shift at discrete events, are more challenging to define and simulate than classical systems with smooth vector fields [, ]. Standard methods often assume that the data comes from a system governed by a fixed set of equations and terms, but time-varying hidden variables can further hinder identification of the system’s underlying dynamics. This motivates our focus on non-autonomous (or non-stationary) systems, where sudden shifts or hidden continuous dynamics complicate accurate modeling and prediction.
We introduce dynamic SINDy, a data-driven method for finding non-autonomous dynamic systems with switching or continuously-varying latent variables. These systems are described by:

$$ $\displaystyle\dot{\textbf{x}}=f(\textbf{x}(t),t)$ (1) $$

where x is vector-valued. A simple such example is $\dot{\textbf{x}}=A(t)\textbf{x}$. Importantly, we focus on systems where the time-varying component and the main variables of interest x are separable (e.g., $\dot{\textbf{x}}=f(\textbf{x},t)=\sin(t)\textbf{x}$, but not $\dot{\textbf{x}}=f(\textbf{x},t)=\sin(t\textbf{x})$. Another assumption is that if multiple trajectories of the system are available, these all display the same underlying switching or hidden variable dynamics.
Dynamic SINDy combines the interpretability of differential equations with the power of deep learning. It uses a deep generative model to uncover sparse governing equations directly from data, employing a variational autoencoder (VAE) to generate time series for differential equation coefficients. This enables data-driven discovery of equations for noisy and non-autonomous systems. The paper is organized as follows: Section 2 introduces key concepts, including SINDy, variational autoencoders, and dynamic VAEs. Section 3 describes the methodology, covering the datasets and the dynamic SINDy framework. Section 4 demonstrates dynamic SINDy’s performance on various systems, including non-autonomous oscillators, Lorenz, Lotka-Volterra, and neural activity data from C. elegans. It also compares dynamic SINDy to switching linear dynamical systems [] and group sparse regression methods []. Section 5 concludes the paper.

## 2 Background and previous work

### 2.1 System identification of non-linear dynamical systems (SINDy)

SINDy (Sparse Identification of Nonlinear Dynamics) [] is a data-driven method that uses sparse regression on a library of nonlinear candidate functions to match data snapshots with their derivatives, revealing the governing equations. The method assumes that only a few key terms explain the system’s dynamics. More specifically, consider $\textbf{x}(t)\in\mathbb{R}^{d}$ governed by the ODE: $\dot{\textbf{x}}(t)=f(\textbf{x}(t))$. Given $m$ snapshots of the system $\textbf{X}=[\textbf{x}(t_{1}),\textbf{x}(t_{2}),...,\textbf{x}(t_{m})]^{T}$ and the estimated time derivatives $\dot{\textbf{X}}=[\dot{\textbf{x}}(t_{1}),\dot{\textbf{x}}(t_{2}),...,\dot{
\textbf{x}}(t_{m})]^{T}$ , we construct a library of candidate functions $\Theta(\textbf{X})=[1,\textbf{X},\textbf{X}^{2},...,\textbf{X}^{p},\sin(
\textbf{X}),\cos(\textbf{X}),...]$.
We then solve a sparse regression problem, $\dot{\textbf{X}}=\Theta(\textbf{X})\Xi$, to identify the optimal coefficients $\Xi$ and to reduce the number of terms, enforcing parsimony. A sparsity-promoting regularization function $R$ is added to the final loss to yield:

$$ $\displaystyle\hat{\Xi}=\textrm{argmin}_{\Xi}(\dot{\textbf{X}}-\Theta(\textbf{X })\Xi)^{2}+R(\Xi)$ (2) $$

Several innovations have followed the original formulation of SINDy [, , , ]. For instance, integral and weak formulations [, ] have enhanced the algorithm’s robustness to noise. Of relevance to our study, SINDy’s generalization to non-autonomous dynamical systems has been previously explored using group sparsity norms [] or clustering algorithms [].

### 2.2 (Dynamic) Variational autoencoders for system identification

The Variational Autoencoder (VAE) [, ] combines neural network-based autoencoders with variational inference for probabilistic modeling and data generation. Unlike standard autoencoders, VAEs stand out due to two key features: (i) VAEs encode input data $X$ as a distribution in the latent space, allowing the decoder to generate new data by sampling from this distribution; and (ii) a regularization term ensures the latent space resembles a standard (e.g., normal) distribution, making it continuous (nearby points generate similar outputs) and complete (all points produce meaningful data). Further mathematical details can be found in Supplementary Material (SM) Section 1.1. A related method of interest is HyperSINDy []. It combines VAEs with SINDy to discover differential equations from data. The VAE approximates the probability distribution of equation coefficients, so that once trained, HyperSINDy generates accurate stochastic dynamics and quantifies uncertainty, making it a powerful tool for model discovery.
In order to adapt the VAE/SINDy framework to non-autonomous systems, we would like to implement generative architectures that capture the temporal dependencies in sequential data. Dynamic VAEs (DVAEs) is an approach that extends VAEs to handle time series data []. A number of DVAE architectures are described that use recurrent neural networks or state-space models to address both latent and temporal relationships [, , , , , , , , ]. We specifically use timeVAE [], which has shown strong performance in generating time series data by processing entire sequences with dense and convolutional layers to capture correlations. Our approach is flexible, allowing the VAE architecture to be swapped for other models better suited to the data or system under study (SM Sec. 1.2).

### 2.3 Other machine learning methods for non-autonomous dynamical systems

Traditionally, methods for handling hybrid or switching systems often involve dividing time or space into segments []. For instance, reduced-order models for nonlinear systems segment time intervals into smaller windows, then build a local, reduced approximation space for each segment [, , ]. Clustering methods are also employed for modeling, particularly in complex fluid flows, where clusters represent states that can transition via a Markov model [, ] or via dynamic mode decomposition with control [, ]. Through data-informed geometry learning, authors in [] reconstruct the relevant “normal forms", which are prototypical realizations of the dynamics, providing bifurcation diagram and insights about the parameters even for non-autonomous systems. Yet another method [] applies Koopman operator theory using DMD algorithms to find time-dependent eigenvalues, eigenfunctions, and modes in linear non-autonomous systems.
We compare dynamic SINDy with two existing methods (Section 4.6). First, we look at a method (recurrent SLDS) [] that extends switching linear dynamical systems (SLDS) [, ] by generating transitions through changes in a continuous latent state and external inputs, rather than relying on a discrete Markov model for switching states. This model breaks the data into simpler segments and is interpretable, generative, and efficiently fitted using modular Bayesian inference. Second, we examine a method from [, ] that uses group-sparse penalization for model selection and parameter estimation. This method assumes shared sparsity across parameters by applying group-sparsity regularization to smaller time windows in the data, identifying the system for each segment, and then combining the results.

## 3 Methods

### 3.1 Datasets

We use a synthetic dataset capturing dynamics of a non-autonomous harmonic oscillator:

$$ $\displaystyle\dot{x}$ $\displaystyle=A(t)y$ $\displaystyle\dot{y}$ $\displaystyle=B(t)x$ (3) $$

where $A(t)$ and $B(t)$ are the time-varying coefficients of the ODE. The time dependence of these coefficients renders the system non-autonomous and difficult to discover using classical methods. We test our approach to see if it can handle switching coefficients, as well as explore continuously varying coefficients, such as sinusoidal functions at different frequencies or finite Fourier series (Figure 1A, Suppl. Fig. 2). To ensure robustness against randomness, we add Gaussian noise with varying levels of variance to the time series.
We replace a set of constant coefficients with a set of time series (sigmoidal, switching, sinusoidal, finite Fourier series) for more complex systems, such as the chaotic Lorenz system:

$$ $\displaystyle\dot{x}$ $\displaystyle=\sigma(t)(y-x)$ $\displaystyle\dot{y}$ $\displaystyle=x(\rho(t)-z)-y$ (4) $\displaystyle\dot{z}$ $\displaystyle=xy-\beta(t)z$ $$

We use large-scale neural recordings from whole-brain imaging to model neuronal population dynamics. C. elegans, with its 302 precisely mapped neurons, offers an ideal balance of simple behavior and complex neuronal activity. We analyze calcium imaging data from Kato et al., which includes neural recordings from the head ganglia and manual annotations of seven behaviors: forward movement, reversal, two types of reversal-to-forward turns, and two forward-to-reversal turns []. Previous studies show that high-dimensional neuronal activity simplifies into low-dimensional patterns, with clear clusters in principal component space representing forward and backward movements. This provides a valuable opportunity to study the link between neural activity and behavior.

### 3.2 System identification for non-autonomous dynamical systems

We explore various VAE architectures designed for inference and generation of time series data. The input is the original time series $X$, and the output are time series of ODE coefficients:

$$ $\Xi_{1:t}=V(\textbf{X}_{1:t})$ (5) $$

where $V$ is the (VAE) architecture, and $\Xi_{1:t}$ is the output time series. ‘Autoencoder" is a misnomer because the input is not designed to match the output in this VAE architecture. The ODE coefficients are linearly combined with a pre-determined SINDy library of basis functions to yield $\hat{\dot{\textbf{X}}}$:

$$ $\hat{\dot{\textbf{X}}}(t)=\Theta(\textbf{X}(t),t)\cdot\Xi(t)$ (6) $$

where $\Theta(\textbf{X}(t),t)$ is a row vector comprising of a polynomial basis up to cubic monomials: $\Theta(\textbf{X}(t),t)=[1\hskip 7.22743ptX_{1}(t)\hskip 7.22743pt...\hskip 7.
22743ptX_{n}(t)\hskip 7.22743ptX_{1}^{2}(t)\hskip 7.22743pt...\hskip 7.22743
ptX_{n}^{3}(t)]$, where $X_{i}$ are features of X. Although we choose a polynomial basis for all of our experiments, the basis can change depending on the problem at hand or any prior information [].
Our goal is to match $\tilde{\dot{\textbf{X}}}$, the derivative we estimate from data using numerical methods, to the output $\hat{\dot{\textbf{X}}}$ of our model (Eq. (5-6)). The loss function takes the following form:

$$ $\textrm{loss}=\sum_{t}||\tilde{\dot{\textbf{X}}}(t)-\hat{\dot{\textbf{X}}}(t)| |+\lambda_{1}R_{kld}+\lambda_{2}R(\Xi)$ (7) $$

where $\lambda_{1,2}$ are hyperparameters of the optimization and $R,R_{kld}$ are regularization terms. $R_{kld}$ is the Kullback-Leibler divergence (KLD) loss, part of the ELBO (evidence lower bound) loss in VAEs (see SM Section 1.1). Regularization terms impose that $\Xi(t)$ is sparse (in coefficients) to encourage parsimony and that $\Xi(t)$ has minimal total variation. More details about the loss function and training, specifically the inference and generation models, can be found in the SM, Sec. 1.3.
We focus on two neural network architectures in our experiments. First, timeVAE (SM Sec. 1.2.1, Suppl. Fig. 1A) is simple for proof-of-concept testing []; however, its major drawback is that it requires the entire time series as input, which can be impractical for long sequences, especially in high-dimensional systems due to memory constraints. To address this, we introduce a new architecture called dynamic HyperSINDy (SM Sec. 1.2.2, Suppl. Fig. 1B). Alternatively, we can use DVAE architectures, which allow sequential data input, overcoming timeVAE’s limitations [].

Figure: Figure 1: (A). Synthetic dataset to test dynamic SINDy with non-autonomous harmonic oscillators (Eq. (3)). Top: Example (SINDy) coefficient time series $A(t)$; Bottom: corresponding trajectories in phase space (B). Dynamic SINDy general architecture schematic; two DVAEs shown as example.
Refer to caption: x1.png

## 4 Results

### 4.1 System identification of non-autonomous harmonic oscillators

We begin by identifying noisy, non-autonomous dynamical systems using a simple toy model – a non-autonomous harmonic oscillator with time-varying ODE coefficients (Eq. 3, Figure 1). First, we vary the coefficient $A(t)$ in a switch-like fashion (Figure 2(a)-(c)). The system behaves like a classic harmonic oscillator, but with a frequency switch. The inferred coefficients (Figure 2, Suppl. Fig. 4) and the reconstructed trajectories (Suppl. Fig. 6) align well with the true values. These trajectories are generated during testing, with $z$ sampled from a standard normal distribution.
When varying both $A(t)$ and $B(t)$ as sinusoids with different frequencies, the resulting trajectories generally capture their oscillations, though some higher error and a large outlier appear toward the end (Figure 2(d), Suppl. Fig. 3). We also successfully reproduce coefficients composed of multiple frequencies (a finite Fourier series) in Figures 2(e)-(f). In (f), some error occurs in the first half because the system approaches a fixed point where the derivative is nearly zero. In such cases, system identification becomes difficult, as multiple solutions can produce the same dynamics.

Figure: Figure 2: Dynamic SINDy generates coefficient time series that match ground truth for non-autonomous harmonic oscillators (Eq. (3)). (a)-(f) different examples of time-varying $A(t)$, $B(t)$.
Refer to caption: x2.png

### 4.2 Uncertainty quantification in non-autonomous harmonic oscillators

We use VAEs to quantify uncertainty by estimating the standard deviation of the coefficients over time. Therefore we generate multiple trajectories by sampling $z$ from a standard normal distribution during testing. Figures 3(a)-(c) show examples of trajectories from networks trained on noisy data with two noise levels: low ($0.01$) and high ($0.5$) standard deviations. As expected, trajectories vary more under high noise than low noise. Our results show that the estimated standard deviation generally follows the true coefficient variations. First, we compute the standard deviation across generated samples at each time point and average these deviations (Figure 3 B(a)). Second, we subtract a smooth mean from the trajectory samples and calculate the standard deviation over time (Figure 3 B(b)). Both methods demonstrate that standard deviation aligns with the ground truth, particularly for switch signal coefficients, but is less clear for Fourier series coefficients. Further work is needed to improve standard deviation estimation, considering the VAE architecture and hyperparameters.

Figure: Figure 3: (A) Dynamic SINDy generates coefficient time series for different levels of Gaussian noise in the coefficient. (B) Inferred noise (standard deviation, or std) scales with ground truth Gaussian noise for different time-varying coefficients. (a) std computed over many generated samples, then averaged (b) std computed over time, then averaged over samples (see Sec. 4.2)
Refer to caption: x3.png

### 4.3 System identification in a non-autonomous, chaotic toy dataset

We next modified the Lorenz system by allowing one of its key parameters ($\sigma,\rho,\beta$) to vary over time, similar to the non-autonomous harmonic oscillator examples. The modified Lorenz equations are:

$$ $\displaystyle\dot{x}$ $\displaystyle=\sigma(t)(y-x)$ $\displaystyle\dot{y}$ $\displaystyle=x(\rho-z)-y$ $\displaystyle\dot{z}$ $\displaystyle=xy-\beta z$ (8) $$

Here, $\sigma(t)$ varies over time as a sigmoid, switch function, sinusoid, or as a Fourier series with 7 overlapping frequencies. Despite these changes, the system still converges to a global attractor.
For system identification, we used two dynamic SINDy architectures: the timeVAE, effective for shorter time series (1000–2000 points), and dynamic HyperSINDy (SM Sec. 1.2.2), suitable for longer time series. Training occurs in two stages: first, we apply a sparsity penalty to set small coefficients to zero; second, we fine-tune the remaining coefficients. After training, we remove the encoder and generate time series from the decoder, closely matching the ground truth across different parameters and functions (Figure 4, Suppl. Fig. 5). Hyperparameters are listed in SM, Sec. 1.3.

Figure: Figure 4: Dynamic SINDy generates coefficient time series that match ground truth for Lorenz dynamics (Eq. (8)). (a)-(h) different examples of time-varying $\sigma(t)$, $\rho(t)$, $\beta(t)$.
Refer to caption: x4.png

### 4.4 Dynamic SINDy used for identifying latent variables and their dynamics

Dynamic SINDy is particularly useful for discovering hidden (latent) variables from incomplete datasets. We demonstrate this using a toy model from ecology: the Lotka-Volterra equations, which describe predator-prey dynamics between two species:

$$ $\displaystyle\dot{x}$ $\displaystyle=\alpha x-\beta xy$ $\displaystyle\dot{y}$ $\displaystyle=-\gamma y+\delta xy$ (9) $$

In our example, we only observe the prey population, $x$, and aim to use dynamic SINDy to uncover the hidden predator population, $y$, and reconstruct a full 2D autonomous system in $x$ and $y$.
We apply dynamic SINDy to $x$, using a library with just three terms: $x,x^{2},x^{3}$. As expected, $\dot{x}$ is expressed solely in terms of $x$, with the $x^{2}$ and $x^{3}$ terms vanishing. We derive a time series for the coefficient $\tilde{y}$, where $\dot{x}=x\tilde{y}(t)$. This inferred $\tilde{y}$ correlates with the hidden $y$, where $\tilde{y}=q\cdot(\alpha-\beta y)$, with $q$ being a scaling factor applied to $x$ before using dynamic SINDy. From $\tilde{y}(t)$, we can infer $y$ and compare it to the true population. In noiseless data, we accurately reconstruct the predator dynamics 6A, but with more noise, recovery becomes harder 6B. Using $\tilde{y}$, we form a new 2D system of equations:

$$ $\displaystyle\dot{x}$ $\displaystyle=a\cdot x\tilde{y}$ $\displaystyle\dot{\tilde{y}}$ $\displaystyle=b+c\cdot x+d\cdot\tilde{y}+e\cdot x\tilde{y}$ (10) $$

where $a,b,c,d,e$ are new model parameters. Comparing the inferred coefficients to the original Lotka-Volterra system by changing variables from $y$ to $\tilde{y}$ and using standard SINDy and the pysindy package, we find a close match (Figure 6C).
We applied this same approach to the non-autonomous harmonic oscillator (Eq. 3, SM Sec. 2.1), further confirming that dynamic SINDy can successfully identify hidden variables and form complete autonomous systems.

Figure: Figure 5: Dynamic SINDy can be used for latent variable discovery. (A). Inferred (blue) versus true (orange) $y$ time series, from noiseless Lotka-Volterra. (B). Same as in (A), when noise of different standard deviations $\sigma$ is added to the Lotka-Volterra trajectory. (C). Inferred versus true coefficients in the Lotka-Volterra 2D ODE system, using SINDy for system identification.
Refer to caption: x5.png

### 4.5 Dynamics in the nematode C. elegans during locomotion

Modern neuroscientific data is noisy, nonlinear, and incomplete, with recordings from hundreds or thousands of neurons, yet many network features and neurons remain unmeasured. This makes it a challenging test for dynamic SINDy’s ability to model such complex systems. To demonstrate our method’s potential, we use a dataset of C. elegans neural activity (Sec. 3.1, []). Unlike previous approaches that rely on probabilistic state space models [] or hidden Markov models [, , ], our method uncovers a global nonlinear switching model [, ]. This model captures key features of the neural data: two stable fixed points representing forward and reversal behaviors, transitions between them, and variability in those transitions, reflecting real neural dynamics.
We first apply PCA to the data from one animal to obtain low-dimensional dynamics that cluster according to behavioral states (Figure 6A,B). Notably, only two dimensions are necessary to differentiate between forward, backward, and turning behaviors, although differentiating between various types of turns requires more dimensions. For a minimally complex model, we focus on the neural trajectory described by the dominant PC mode and its derivative. Our goal is to identify a nonlinear, parsimonious, and global model of the form:

$$ $\displaystyle\dot{x}$ $\displaystyle=y$ (11) $\displaystyle\dot{y}$ $\displaystyle=f(x,\beta)+u(t)$ (12) $$

where $x$ is the data projected onto the first principal component, $f$ is an unknown function, $u(t)$ is a potential switching or control signal, and $\beta$ is a vector of parameters we would like to fit to our data.
We apply dynamic SINDy to minimize the error between the model’s derivatives $(\dot{x}(t),\dot{y}(t))$ and the dominant PC derivatives from the data. Following sections 4.1 and 4.3, we identify the sparsity pattern of the SINDy coefficients, enforcing $\dot{x}=y$. The method highlights the terms $1,x,y,x^{2},x^{3}$ for describing $\dot{y}$ and calculates their time-varying coefficients (SM Sec. 3.1, 3.2). To further simplify the model, we set coefficients for all variables to be constant, except the flexible term $u$ which we can also reduce to a time series of switches (Figure 6D, see SM Sec. 3.2.1) without meaningfully affecting global dynamics. Converting $u$ into a switching signal simplifies this term, helping to regularize the model and improve interpretability. This approach aligns with previous studies showing bistability and sudden transitions in behavior.
Our approach identifies a cubic function for the differential equation model: $\ddot{x}=f(x,\beta)+u(t)=0.002\cdot x^{3}+0.0087x^{2}-0.22\cdot y+0.05\cdot x+
u_{i}$, where $u_{i}$ alternates between $u_{0}=-0.266$ and $u_{1}=0.044$ (see SM Sec. 3.1 for details). Each time $u$ switches, the cubic function shifts, altering the fixed point that the trajectory converges to. This switching signal $u$ enables the transitions between the two fixed points, which correspond to forward and reversal behaviors. Overall, the reconstructed data captures key features like fixed points and transitions (Figures 6E-F, 6G). By labeling the trajectory based on behaviors, we align the inferred dynamics with the training data (Figure 6H).
The reconstructions are accurate regardless of whether we use the processed switching term $u$ (Figure 6 D) or the original time series $u$ (Figure 6 C). However, $u$ alone does not adequately explain the data; removing other terms leads to poor fits or instability. By systematically eliminating different terms, we find that all are essential for capturing the dynamics. When we initialize the inferred ODE system from different starting points and use $u$ from training, the resulting dynamics qualitatively match the data. This suggests that our method effectively avoids overfitting. Unlike Morrison et. al., which relies on selecting a model based on human-labeled behavioral states, dynamic SINDy is fully data-driven and does not require labeled data to partition the phase space []. Additionally, unlike Fieseler et. al., our model accommodates nonlinear dynamics with two stable fixed points []. A key advantage of our ODE model is the potential for biologically interpretable parameters (see []. SM Sec. 3.3 offers a more comprehensive discussion of the benefits of our framework, comparing our global nonlinear switching model to previous studies. In summary, we have demonstrated that dynamic SINDy can do automatic data-driven model discovery, generating a nonlinear model with minimal input from the data scientist.

Figure: Figure 6: (A). C. elegans neuronal activity is low-dimensional and clusters according to behavior; (B). Neuronal activity in phase space given by the first principal coordinate and its derivative; (C). Dynamic SINDy inferred constant term; (D). Processing coefficient in (C) as switch; (E) and (F) ODE model’s match to ground truth trajectories; (G) (and (H)) 2D model trajectory (with labeled behavior).
Refer to caption: x6.png

### 4.6 Dynamic SINDy and other methods for system identification

We begin by comparing dynamic SINDy with switching linear dynamic systems (SLDS) [, ] and its extension, rSLDS []. SLDS uses a discrete latent variable, $z_{t}$, to partition the state space between switches (see SM Sec. 4.1), simplifying complex nonlinear dynamics into more manageable linear segments. The rLDS extension allows switches to depend on continuous latent states and external inputs using logistic regression []. We evaluate how well SLDS/rSLDS identifies switching signals in the dynamical systems studied so far, specifically inferring where the latent variable $z$ changes for switching to occur. We use coefficients based on sigmoids and two types of switching signals (refer to the "ground truth" in Figure 7A). Running SLDS or rLDS generates samples of the latent variable $z$ that segment the training trajectory, enabling us to compare this segmentation with the actual ground truth switches.
We find that for a sigmoidally varying coefficient, SLDS identifies the switch fairly well (Figure 7A (a), (b)), as shown by the colored trajectories and the insets comparing the $z_{t}$ time series to the ground truth; although for Lorenz dynamics, the predicted change in the latent state $z_{t}$ is slightly delayed relative to the actual switch (Figure 7A (b)). However, SLDS struggles when there are multiple state switches in the time series (Figures 7A (c) and (d)). For the harmonic oscillator, rSLDS produces a model with too many switches and is more complex than the ground truth. For Lorenz dynamics, both SLDS and rSLDS switch periodically whenever the dynamics change between attractors, but this periodicity does not match the true switches defined by the coefficients. To address these challenges, we added time as a new dimension to the dataset, represented as a simple feature vector $[1,2,…,T]$, where $T$ is the total number of time steps. The goal was for SLDS/rSLDS to recognize that the switches are time-dependent rather than state-dependent. However, this addition did not improve the performance of SLDS or rSLDS.
Another method for identifying non-autonomous systems, discussed in references [, ] and SM Sec. 4.2, involves dividing the trajectory into smaller time windows and applying SINDy to each segment while enforcing a consistent sparsity pattern across all windows. We tested this approach on two toy datasets, using SINDy coefficients modeled as sigmoids, sinusoidal functions, and a Fourier series with seven frequencies. Without the group sparsity regularization, the sparsity patterns varied across the windows, highlighting the importance of group sparsity for achieving a coherent solution. The group sparsity approach worked well when the coefficients were sigmoid functions with varying smoothness (Figure 7 B (a) and (c)). However, it struggled with sinusoidal and Fourier series coefficients, particularly in the Lorenz system (7 B(b) and (d)). In cases of misidentified coefficients, the algorithm also generated incorrect sparsity patterns. We conclude that neither SLDS or rSLDS, nor the group sparsity method are as effective as our method in identifying non-autonomous dynamical systems from data.

Figure: Figure 7: (A) SLDS/rSLDS infers switching behavior for non-autonomous harmonic oscillators and Lorenz dynamics as coefficients vary. (a)-(d) left: ground truth dynamics, labeled switch values colored blue and orange. Inset shows true coefficient. (a)-(d) right: dynamics labeled by inferred switch. Inset: ground truth $z$ and samples of discrete latent values $z$ labeled by switch. (B) Inferred SINDy coefficients versus ground truth using group sparsity method.
Refer to caption: x7.png

## 5 Conclusion

We have developed dynamic SINDy, an extension of SINDy designed for data-driven identification of noisy, nonlinear, and non-autonomous dynamical systems, as well as for discovering latent variables.
We demonstrated the effectiveness of dynamic SINDy on both benchmark synthetic datasets and a real, noisy, chaotic dataset of neuronal activity from C. elegans.
However, our method has some limitations. First, the DVAE architecture has many hyperparameters to tune, and results may not be robust to these settings, especially in noisy datasets. A systematic approach for hyperparameter tuning and addressing multiple solutions is necessary. To prevent overfitting, we should encourage simpler time series through regularization. Although we included a simple approximation of total variation term in our loss function, realistic datasets might require more sophisticated regularization. Future research should explore different DVAE architectures to evaluate their accuracy in reproducing dynamics and their ability to quantify uncertainty. It would also be valuable to apply our method to experiments with different trial dynamics and types of noise than those studied here (e.g., varying switching points across trials). Lastly, we aim to apply dynamic SINDy to realistic data from fields like biology, physics, and engineering to uncover hidden dynamics and fully utilize its potential for discovering latent variables.

## 6 Supplemental information

Supplemental Materials are included with the main manuscript. A repository with the code to generate the figures in the main paper can be found at [https://github.com/dvoina13/Dynamic_SINDy](https://github.com/dvoina13/Dynamic_SINDy).

## 7 Acknowledgements

The authors acknowledge support from the National Science Foundation AI Institute in Dynamic
Systems (grant number 2112085). We benefited from conversations with Preston Jiang, Ryan Raut,
and Sam Otto.
Funding sources have no direct involvement in the study design, in the analysis,
interpretation of results, or in the writing of the manuscript.

## 1 Methods

### 1.1 Variational Autoencoders

In this section, we elaborate on the mathematical foundation of the Variational Autoencoder (VAE) architecture [, ].
Like standard autoencoders, VAEs have an encoder and decoder network to process input data and generate output. However, instead of mapping inputs to fixed points in the latent space, the encoder maps them to a probability distribution. The decoder then samples from this distribution to reconstruct the input. This probabilistic framework reduces overfitting by introducing variability into the latent space. After computing the reconstruction error, the network is trained via backpropagation, with the VAE relying on the reparameterization trick to ensure gradients can be propagated through the network.
Mathematically, we aim to approximate the data distribution $p^{*}(X)$ of some given observations $X$. When direct computation is intractable, we introduce a latent variable $z$ such that $p^{*}(x)$ is decomposed as:

$$ $\displaystyle p^{*}(x)=\int_{z}p^{*}(x|z)p^{*}(z)dz$ (13) $$

where $p^{*}(x|z)$ is the likelihood and $p^{*}(z)$ is a prior, often set to a standard normal distribution. Since this integral is difficult to compute, we approximate $p_{\theta}(x|z)$ with a neural network parameterized by $\theta$. To estimate the posterior distribution $p^{*}(z|x)$, we approximate it with another neural network $q_{\phi}(z|x)$, parameterized by $\phi$. This is the core idea of variational inference: complex distributions are approximated by simpler, parametrized ones through optimization. We arrive at the following objective:

$$ $\displaystyle\log p_{\theta}(x)$ $\displaystyle=\log\int_{z}p_{\theta}(x,z)dz$ $\displaystyle=\log\int_{z}p_{\theta}(x,z)\frac{q_{\phi}(z|x)}{q_{\phi}(z|x)}dz$ $\displaystyle=\log E_{z\sim q_{\phi}(z|x)}[\frac{p_{\theta}(x,z)}{q_{\phi}(z|x )}]\geq E_{z}[\frac{\log p_{\theta}(x,z)}{q_{\phi}(z|x)}]$ (14) $$

by Jensen’s inequality. This leads to the evidence lower bound (ELBO):

$$ $\displaystyle\mathcal{L}=E_{z}\frac{\log p_{\theta}(x,z)}{q_{\phi}(z|x)}$ (15) $$

The ELBO sets a lower bound for the evidence of observations and maximizing $\mathcal{L}$ will increase the log-likelihood of $X$. To find the parameters $\theta,\phi$ so as to maximize the ELBO, it is convenient to re-write $\mathcal{L}$ in the following way:

$$ $\displaystyle\mathcal{L}$ $\displaystyle=\int_{z}q_{\phi}(z|x)\log(\frac{p_{\theta}(x,z)}{q_{\phi}(z|x)})$ $\displaystyle=\int_{z}q_{\phi}(z|x)\log(\frac{p_{\theta}(x|z)p(z)}{q_{\phi}(z| x)})=E_{z\sim q_{\phi}(z|x)}\log(p_{\theta}(x|z))-D_{KL}(q_{\phi}(z|x)||p(z))$ (16) $$

where $D_{KL}$ is the Kullback-Leibler divergence between the approximate posterior $q_{\phi}(z|x)$ and the prior $p(z)$. The ELBO comprises two terms: the expected log-likelihood of the data, and a regularization term that enforces similarity between the posterior and the prior.
In the autoencoder perspective, the encoder network maps inputs to the latent space via $q_{\phi}$, and the decoder maps the latent variables back to the input space via $p_{\theta}$. Both networks are trained jointly using stochastic gradient descent to optimize the ELBO. More details on the VAE framework can be found in Kingma et al.’s excellent review [].

### 1.2 Dynamic VAEs

Generating time series data presents unique challenges due to the intricate temporal relationships and the distribution of features at each time point. One common approach is using generative adversarial networks (GANs), which often incorporate recurrent neural networks (RNNs) for both generation and discrimination. However, despite numerous proposed architectures, GANs have struggled to capture the complex temporal dependencies inherent in time series data. Yoon et al. [] introduced a novel approach that blends the supervised training used in autoregressive models with the unsupervised training of GANs. While we experimented with this method for generating time series, the training proved to be time-consuming and impractical for our datasets (see Section 3.1). Further limitations are discussed by Desai et al. [].
As a result, we shifted our focus to methods based on Variational Autoencoders (VAEs) for time series, leveraging deep learning techniques to model complex temporal patterns more effectively.
An extensive review [] offers a unified framework for several VAE models extended to handle temporal and sequential data. These models, collectively referred to as dynamic VAEs (DVAEs), share common notation, methodology, and a standardized mathematical formalism. The review covers various approaches, including Deep Kalman Filters [, ], Kalman Variational Autoencoders [], Stochastic Recurrent Networks [], Variational Recurrent Neural Networks [, ], Stochastic Recurrent Neural Networks [], Recurrent Variational Autoencoders [], and Disentangled Sequential Autoencoders []. In the following section, we will expand on the mathematical framework common to these methods, as outlined in [].
Briefly, given a time-series $X_{1:T}$, and assuming latent variables $Z_{1:T}$, the goal is to specify the joint distribution of the observed and latent sequential data $p_{\theta}(X_{1:T},Z_{1:T})$, where $\theta$ denotes the parameters of the true distribution’s probabilistic model. DVAEs are hierarchical models in which both observed and latent variables are treated as time-ordered vectors. These models are often causal, meaning the distribution of variables at time $t$ depends only on previous time steps. This causality imposes the following factorization:

$$ $\displaystyle p(X_{1:t},Z_{1:T})=\prod_{t=1}^{T}p(x_{t},z_{t}|x_{1:t-1},z_{1:t -1})=\prod_{t=1}^{T}p(x_{t}|x_{1:t-1},z_{1:t})p(z_{t}|x_{1:t-1},z_{1:t-1})$ (17) $$

The joint distribution of observed and latent variable sequences can be factorized using the chain rule. Crucially, different models proposed in the literature make different conditional assumptions to simplify the dependencies in the conditional distribution. For example, a simple model may make the following simplifications: $p(X_{t}|X_{1:t-1},Z_{1:t})=p(X_{t}|Z_{t})$ and $p(Z_{t}|X_{1:t-1},Z_{1:t-1})=p(Z_{t}|Z_{t-1})$. In addition, different models may implement different network architectures to approximate $p_{\theta}$ and $q_{\phi}$. A detailed account of the kind of assumptions that each model implements to simplify
(17) can be found in [].

#### 1.2.1 timeVAE

TimeVAE is a variational autoencoder designed to generate multivariate time-series data []. It extends the standard VAE framework to model both the latent space and the temporal dependencies of a sequence of data vectors. Supplementary Figure 1A illustrates the basic TimeVAE architecture, which uses dense and convolutional layers without requiring specific time-series knowledge. The decoder allows for customizable distributions by adding layers to capture time-series components like level, trends, and seasonality, though we used the base version that excludes these custom structures in our experiments.
The input to the encoder is a 3D array of size $N\times D\times T$, where $N$ is the batch size, $D$ is the number of feature dimensions, and $T$ is the number of time steps. The encoder processes the data through convolutional layers with ReLU activations, flattens the output, and then applies a fully connected layer. The final encoder layer has $2d$ units, representing the mean and variance of a multivariate Gaussian distribution, where $d$ is the dimensionality of the latent space, a key hyperparameter. The reparametrization trick is used to sample from the Gaussian distribution, parameterized by the encoder’s output.
The decoder reconstructs the data by first passing the sampled latent vector $z$ through a fully connected layer, reshaping it into a 3D array, and processing it through a series of transposed convolutional layers with ReLU activation. The last time-distributed fully connected layer produces the final output that matches the original input dimensions.
Training TimeVAE involves optimizing the ELBO loss function (discussed in Section 1.1) with different weights on the reconstruction error and KL divergence between the approximate posterior $q_{\phi}(z|x)$ and the prior $p_{\theta}(z)$. Hyperparameters are tuned to determine the appropriate balance between reconstruction loss, KL divergence, and additional regularization terms (e.g., sparsity and total variation, for our problem set-up).
TimeVAE has been tested on four multivariate datasets []: (1) a 5-dimensional sinusoidal dataset with varying frequencies, amplitudes, and phases; (2) a 6-dimensional stock market dataset from Yahoo Finance; (3) a 28-dimensional appliances energy prediction dataset from the UCI Machine Learning Repository; and (4) a dataset with 15 features of hourly air quality sensor readings from the UCI Machine Learning Repository. The results show that TimeVAE performs comparably to top generative models across various metrics, is computationally efficient, and outperforms existing methods in next-step prediction tasks, particularly when training data is limited [].

#### 1.2.2 dynamic HyperSINDy

Inspired by previous work [, ], we developed a hierarchical architecture to address non-autonomous systems, illustrated in Supplementary Figure 1B. The main text focuses on the dynamic SINDy framework with a timeVAE architecture, while Supplementary Figures 3, 4, 5 show that incorporating dynamic HyperSINDy results in coefficients and trajectories that closely match the ground truth.
The first level consists of a standard VAE with encoder and decoder modules (SM Sec. 1.1). The decoder generates $X$ with a probability distribution $p_{t}(X)$ at each time step $t$. The next level introduces a hypernetwork, implemented as either a long short-term memory network (LSTM) or multi-layer perceptron (MLP), which updates the decoder’s weights to adjust the probability distribution for the following time step, allowing the system to capture temporal drift in the output:

$$ $\displaystyle W_{\textrm{Decoder}}(t+1)$ $\displaystyle=W_{\textrm{Decoder}}(t)+\sum_{i=1}^{M}\alpha_{i}(t)\cdot D_{i}, \textrm{where }$ $\displaystyle\textrm{LSTM}(t)$ $\displaystyle=\alpha(t)=[\alpha_{1}(t),...,\alpha_{n}(t)]$ (18) $$

Here, ${D_{i}}$ are fixed basis tensors to be learned, and $\alpha_{i}$ are hypernetwork outputs. This architecture adapts to changing dynamics, adjusting the decoder based on reconstruction error and updating the output probability distribution.

We modified this architecture for our problem. Instead of a VAE generating data, our decoder produces SINDy coefficients, which, when combined with the SINDy library, replicate system dynamics. The decoder approximates the true pdf of the SINDy coefficients. This setup builds on [] by adding a hypernetwork that updates decoder weights, forming what we call dynamic HyperSINDy. This extension allows for time-varying SINDy coefficients, processed sequentially rather than requiring the entire time series as input (as in timeVAE).
Two primary training methods are used for dynamic HyperSINDy:

- •
Online Learning: Ideal for switching systems where the network adapts as dynamics change. However, network parameters evolve, requiring tracking of parameter changes and identifying switch points after training. The hypernetwork is not needed in this setup.
- •
Alternate learning:
The hypernetwork is trained first with fixed main module parameters and basis tensors $D_{i}$, followed by adjustment of the main module parameters/$D_{i}$, while fixing the hypernetwork. This method is best for continuously varying SINDy coefficients, with LSTM as the preferred hypernetwork.

Training used the hyperparameters listed in Table 1. We processed one trajectory at a time (trial batch size of 1) and used batch sizes of 1-10 time steps. The latent dimension of the VAE was set to 25 while the starting threshold was $0.1$. Every 50 epochs, we evaluated and set to zero any SINDy coefficients with a mean absolute value below a threshold. A relaxed L0 norm in the loss function encouraged sparsity in the SINDy coefficients, following [, ].
Several hyperparameters were increased progressively during training. The threshold rose by $0.005$ every $50$ epochs until it plateaued, alongside the weight $\lambda_{kl}$ for the KL divergence term, which increased until it reached a maximum value $\lambda_{max}$. The threshold plateaus as well once $\lambda_{max}$ is reached. We fixed the number of basis tensors $D_{i}$ to 10, which combined linearly with hypernetwork outputs to form decoder weights via Eq. (18).
The encoder consisted of four fully connected layers with hidden dimensions of $64$, using ELU activation and an input dimension twice that of $X$, as it takes $X$ and $\dot{X}$ as input. The decoder also had four hidden layers, with a hidden dimension of $256$ and ELU activation. The hypernetwork, either an LSTM or MLP, contained two layers with an input dimension of $25$. We trained using the AdamW optimizer with an initial learning rate of $0.001$, weight decay of $1e-5$, gradient clipping at $1$, and Amsgrad. Additionally, an exponential learning rate scheduler with $\gamma=0.999$ was used. Many hyperparameters match those in [].

**Table 1: Hyperparameters for Dynamic HyperSINDy**
| hyperparameter | value |
| --- | --- |
| batch size (trials) | 1 |
| batch size (time steps) | 1-10 |
| latent variable dimension | 25 |
| threshold | 0.1 |
| threshold interval | 50 |
| threshold increment | 0.005 |
| $\lambda_{kl}$ | 0.01 |
| $\lambda_{kl}$ increment | $\lambda_{kl}/5$ |
| $\lambda_{kl}$ max | 1 |
| M (number of basis tensors) | 10 |
| hidden dim (decoder) | 256 |
| hidden dim (encoder) | 64 |
| input dim (LSTM) | 25 |
| gradient clip | 1.0 |
| cell dimension (LSTM). | 30 |
| optimizer | AdamW |
| weight decay | 1e-5 |
| amsgrad | True |
| learning rate | 0.001 |
| learning rate scheduler | ExponentialLR |
| gamma | 0.999 |

### 1.3 Training dynamic SINDy with timeVAE: methodology and hyperparameters

Training timeVAE requires normalizing the data beforehand. While [] normalizes by subtracting the minimum and dividing by the maximum to scale the data between 0 and 1, we normalize by dividing only by the maximum value. This normalization method affects the SINDy coefficients produced by our method, so we re-scale the resulting time-series before comparing them to the ground truth in synthetic datasets.
The loss function used to train our timeVAE architecture is:

$$ loss $\displaystyle=\lambda_{MSE}\cdot||\hat{\dot{X}}-\dot{X}||_{2}^{2}+\lambda_{KL} \cdot\textrm{KL div}+\lambda_{sp}\cdot<||\xi_{i,j}(t)||_{1}>_{i,j}$ (19) $\displaystyle+\lambda_{tv}\cdot\frac{<||\xi_{i,j}(t+1)-\xi_{i,j}(t)||_{1}>_{i, j,t}}{<||\xi_{i,j}(t)||>_{i,j,t}+\epsilon}$ (20) $$

where $\epsilon$ is the machine precision limit. The hyperparameters $\lambda_{MSE},\lambda_{KL},\lambda_{sp},$ and $\lambda_{tv}$ balance accuracy and complexity by adjusting the weights on the different loss terms: $\lambda_{MSE}$ controls the mean squared error, while the others handle regularization.

- •
The first term represents the mean squared error between the inferred derivative $\hat{\dot{X}}$ using dynamic SINDy and the derivative from the data $\dot{X}$. For all synthetic datasets, the ground truth derivative is the one used to obtain the trajectories $X$.
- •
The second term is the Kullback-Leibler divergence (KL div), a standard term in variational autoencoders (discussed in SM Sec. 1.1). It measures how closely the posterior distribution of $z$, as computed by the encoder given $X$, resembles the standard normal distribution. The KL divergence has an analytic form:
$\displaystyle\textrm{KL div}=-\frac{1}{2}\cdot<(1+2\log(\sigma_{z_{i,j}})-
\sqrt{\mu_{z_{i,j}}}-\exp{(2\log{\sigma_{z_{i,j}}})})>_{i,j}$
(21)
where $\langle\cdot\rangle$ indicates averaging over latent dimensions $i$ and data points $j$, and $\mu_{z_{j}}$ and $\sigma_{z_{j}}$ represent the mean and standard deviation of $z_{j}$, with $\mu_{z_{j}}$ and $\log(\sigma_{z_{j}})$ as the encoder outputs for each input $X_{j}$.
- •
The third term in Eq. (20) is a sparsity penalty that encourages some SINDy coefficients to be zero.
- •
The fourth term is a normalized total variation penalty that prevents drastic changes in the solution over time.

#### 1.3.1 Non-autonomous harmonic oscillators

For the non-autonomous harmonic oscillators, we use the hyperparameters in Table 2 to train the timeVAE architecture. These remain constant across datasets, despite differences in the time-varying coefficients $A(t)$ and $B(t)$. However, key hyperparameters like $\lambda_{sp}$ and $\lambda_{tv}$ vary depending on the dataset, as shown in Table 3. Training is performed using the ADAM optimizer with a weight decay of $1e-5$ and gradient clipping at $1$.

**Table 2: Hyperparameters for timeVAE (non-autonomous harmonic oscillator)**
| hyperparameter | value |
| --- | --- |
| batch size | 1 |
| latent dimension | 2 |
| threshold | 0.01 |
| library size | 3 |
| $\lambda_{MSE}$ | 3 |
| $\lambda_{KL}$ | 1000 |

**Table 3: Hyperparameters for timeVAE at different phases of training (non-autonomous harmonic oscillator)**
| hyperparameters | dataset | at first training phase | at second training phase |
| --- | --- | --- | --- |
| $\lambda_{sp}$ | A(t) sigmoid | 50 | 0 |
|  | A(t) switch signal 1 | 500 | 0 |
|  | A(t) switch signal 2 | 200 | 0 |
|  | A(t) finite Fourier series | 1 | 0 |
| $\lambda_{tv}$ | A(t) sigmoid | 100 | 1000 |
|  | A(t) switch signal 1 | 100 | 1000 |
|  | A(t) switch signal 2 | 100 | 1000 |
|  | A(t) finite Fourier series 1 | 0 | 2 |

#### 1.3.2 Lorenz dynamics

For the results in Sec. 4.3, involving the chaotic system with a time-varying parameter in the Lorenz dynamics, we follow the same steps as before (data normalization/post-processing, loss function, and two-stage training: first for sparsity pattern, then for coefficient recovery), but with different hyperparameters listed in Table 4. These hyperparameters remain constant, regardless of how the Lorenz parameters vary over time.
During the first training stage, when the sparsity penalty is non-zero, the batch size is set to $10$ to ensure the correct sparsity pattern is learned. We use RMSProp with a weight decay of $10^{-5}$ and gradient clipping at $10$. The threshold gradually increases from $0.05$ to $0.1$ in increments of $0.025$ per epoch, while $\lambda_{sp}$ rises from 0 0 to $20$ in steps of $1$. $\lambda_{tv}$ is fixed at $1000$. This gradual increase in hyperparameters follows a successful approach from a related study [].

**Table 4: Hyperparameters for timeVAE (Lorenz dynamics)**
| hyperparameter | value |
| --- | --- |
| latent dimension | 5 |
| library size | 3 |
| $\lambda_{MSE}$ | 3 |
| $\lambda_{KL}$ | 1000 |

#### 1.3.3 Lotka Volterra

The incomplete Lotka Volterra system has only one variable $x$, therefore the library has three terms: $x,x^{2},x^{3}$. For training, we use the hyperparameters listed in Tables 5 and 6.

**Table 5: Hyperparameters for timeVAE (Lotka Volterra)**
| hyperparameter | value |
| --- | --- |
| batch size | 1 |
| latent dimension | 2 |
| library size | 3 |
| $\lambda_{MSE}$ | 3 |
| $\lambda_{KL}$ | 1000 |

**Table 6: Hyperparameters for timeVAE at different phases of training (Lotka Volterra)**
| hyperparameters | at first training phase | at second training phase |
| --- | --- | --- |
| $\lambda_{sp}$ | 0.1 | 0 |
| $\lambda_{tv}$ | 0 | 0 |

## 2 Latent variable discovery

### 2.1 Non-autonomous harmonic oscillator

We can use the same approach with the non-autonomous harmonic oscillator as with the Lotka-Volterra system. We set $A(t)=-4$ and vary $B(t)$ sigmoidally such that $B(t)=2+\frac{1}{1+exp(5+t)}$. After dynamic SINDy identifies a trajectory for $B$, we add it to $(x,y)$ to form a 3D dynamical system. Using SINDy on $(x,y,B)$, we discover the following ODE which is almost exactly identical to the true dynamics, given that $B$ is a sigmoid that can be described by $\dot{B}=-6+5B-B^{2}$:

$$ $\displaystyle\dot{x}$ $\displaystyle=-3.997y$ $\displaystyle\dot{y}$ $\displaystyle=1\cdot Bx$ (22) $\displaystyle\dot{B}$ $\displaystyle=-5.875+4.903B-0.981B^{2}$ $$

## 3 Dynamic SINDy for system identification of neuronal dynamics in the nematode C. elegans

### 3.1 Results

Like in Morrison et. al., we have discovered a dynamical system model that switches between two stable fixed points. The differential equation model is expressed through a cubic function: $\dot{x}=a\cdot x^{3}+b\cdot x^{2}+c\cdot x+d\cdot y+u$ with distinct fixed points corresponding to the different switching states of $u$.
More precisely, the differential model inferred has the form:

$$ $\displaystyle\dot{x}=y$ $\displaystyle\dot{y}=-0.002\cdot x^{3}+0.0087x^{2}-0.22\cdot y+0.05\cdot x+u_{ i},\hskip 7.22743pti=1,2$ (23) $\displaystyle u_{0}\approx-0.266;u_{1}\approx 0.044.$ $$

When $u=u_{0}<0$, the dynamical system has one stable fixed point at $-5.25$ (the other roots of the cubic equation are complex). This fixed point corresponds to the reversal behavior. Then, when $u=u_{1}>0$, the dynamical system has two stable and one unstable fixed point: $-2.32$ and $7.88$ stable fixed points and $-1.19$ unstable. Therefore, varying $u$ can generate a bifurcation. In practice, the trajectory shifts between $-5.25$ and $-2.32$ in tandem with behavioral switches between reversal and forward states.
The two-dimensional model is a simple model that fits the first principal component and captures stable state clusters and turning trajectory variability. Once the low-dimensional coordinates are identified, dynamic SINDy effectively enables data-driven model discovery. Future work will extend this approach to multiple animals to test its generality across individuals.

### 3.2 Training dynamic SINDy on C. elegans data

We apply timeVAE to infer dynamics for a single worm as a proof-of-concept, demonstrating dynamic SINDy’s capability for data-driven discovery. Since there is only one trajectory per worm, uncertainty quantification isn’t possible. As with synthetic datasets, we normalize the trajectories and train timeVAE to infer the differential equation’s sparsity pattern using the hyperparameters from Tables 7 and 8.
With a threshold of $0.01$, only the terms $1$, $x$, and $x^{2}$ are considered important. We add $y$ and $x^{3}$ for comparison with the Morrison et al. model and retrain with a fixed sparsity pattern, omitting sparsity regularization from the loss. Our SINDy coefficients vary over time, matching behavioral transitions between forward and reversal locomotion. Increasing total variance regularization was not effective, so we averaged non-constant SINDy coefficients over time to simplify the model. This is part of the training process, where we take the average of all non-constant coefficients at the last layer of the network to yield the model output and backpropagate. The constant SINDy coefficient is not constrained, but all other coefficients do not change in time.
The resulting differential equation model, detailed in Eq. 3.1, includes the time series $u(t)$, which is shown in Figure 6C (main text). We interpret $u(t)$ as a switching variable and hypothesize that even a simple switching time series can qualitatively capture the neural activity data. To test this hypothesis, we post-process the $u(t)$ time series to generate a switch-like signal.

#### 3.2.1 Post-processing the switching signal

Starting with the $u(t)$ time series inferred using dynamic SINDy, we perform the following steps:

- •
Subtract the mean of $u(t)$ over time. We also note the approximate minimum and maximum values, which will be used later.
- •
Scale the data by a large factor (1000) and apply a pointwise sigmoid function across time, producing a time series of switches between 0 and 1.
- •
Finally, re-scale the time series to vary between the previously determined minimum and maximum values, and then add back the mean $<u(t)>_{t}$ to obtain the final post-processed switching time series.

To evaluate the accuracy of our model, we integrate the differential equation from Eq. 3.1 using the post-processed switching signal $u(t)$ and compare the resulting trajectory to the real trajectory (Figure 6E-G, main text). Since the C. elegans data has low time resolution ($\Delta t=0.35749752$), we interpolate the data using the CubicSpline function from the scipy.interpolate library. We reduce the time step to $\Delta t/100$ and perform numerical integration using the Euler method.

**Table 7: Hyperparameters for timeVAE (C. elegans data)**
| hyperparameter | value |
| --- | --- |
| batch size | 1 |
| latent dimension | 2 |
| threshold | 0.01 |
| library size | 3 |
| $\lambda_{MSE}$ | 3 |
| $\lambda_{KL}$ | 1000 |

**Table 8: Hyperparameters for timeVAE at different phases of training (C. elegans data)**
| hyperparameters | at first training phase | at second training phase |
| --- | --- | --- |
| $\lambda_{sp}$ | 10 | 0 |
| $\lambda_{tv}$ | 100 | 1000 |

### 3.3 Background: related studies and comparisons to our model

#### 3.3.1 Comparison with state space models

Our findings with dynamic SINDy reveal a key similarity with the probabilistic state space model proposed by Linderman et al. []: both models switch between different dynamical regimes. However, despite being nonlinear, our model is more parsimonious in several ways.
Linderman et al. propose a hierarchical recurrent state space model that switches between simple linear models, using Bayesian inference to fit the model at scale []. This model decomposes complex nonlinear neural activity into discrete states with simple linear dynamics, which correspond to behaviorally relevant aspects of the worm’s behavior. The transition probabilities depend on both the preceding state and the position in continuous state space, with each discrete state largely tied to the activation of specific neuron clusters.
While this model provides insights into C. elegans neural dynamics, its linear state space models are local. The model switches between eight discrete states, each representing a smaller linear system fitted to the data to explain local dynamics []. In contrast, dynamic SINDy discovers a global nonlinear ODE model that switches between only two states. Thus, we simplify the model by replacing eight local linear regimes with a more compact nonlinear system switching between two states.
A future direction is to develop a generative model for the switching behavior of $u(t)$, possibly using a probabilistic model or differential equation linked to the variables $x$ and $y$. This would allow us to eliminate the dynamic SINDy network post-training, retaining a global nonlinear switching differential equation with just four parameters, compared to the many more parameters required by the hierarchical recurrent SLDS for its eight linear systems (even when considering only the continuous variable dynamics for a fair comparison with our approach in PC space).
Moreover, the hierarchical recurrent SLDS is a statistical model that doesn’t directly map onto network dynamics or account for biologically realistic state transitions. While further research is needed to validate our model’s connection to biological measurements of neural activity, nonlinear differential equations like ours are potentially more interpretable. For example, a single parameter change in a global nonlinear model similar to ours can reproduce different long-timescale behaviors observed in C. elegans [] (see Sec. 3.1.4 below). This modulation mirrors distinct changes in state distribution and switching frequencies seen in experiments, which are linked to specific neuromodulators and neurons [, ].
The challenges discussed here also apply to simpler models based on Markov dynamics, such as hidden Markov models (HMMs) [, , ].

#### 3.3.2 Comparison with a nonlinear global model with control

Our C. elegans neural activity modeling is inspired by Fieseler et al. and Morrison et al. [, ]. Unlike state space models, Morrison et al. discovered a minimally parameterized global nonlinear model with control that mimics Hidden Markov model state transitions within a single dynamical system. This model captures key features of the C. elegans calcium imaging data, including two stable fixed points for forward and reversal behaviors, state transitions triggered by control signals, and variability in transition trajectories that match neural activity data [].
The model is represented as:

$$ $\displaystyle\dot{\textbf{x}}=F(\textbf{x},\beta)+\textbf{u}(t)$ (24) $$

where $\beta$ governs longer timescale dynamics, and $u(t)$ is a control signal operating on faster timescales. $u(t)$ is a one-dimensional signal that may integrate multiple local and non-local processes. This separation of intrinsic dynamics and control inputs increases the model’s interpretability.
Nonlinear control has been used in other biological networks to describe switching between multiple stable states [, , ] A significant advantage is that a nonlinear model can have multiple fixed points corresponding to different behavioral states – in the case of C. elegans, forward and reversal motion. A heuristic model capturing C. elegans behavior is:

$$ $\displaystyle\dot{\textbf{x}}$ $\displaystyle=\textbf{y}$ $\displaystyle\dot{\textbf{y}}$ $\displaystyle=-(\textbf{x}-1)(\textbf{x}-\beta)(\textbf{x}+1)+\lambda\textbf{y }+\textbf{u}(t)$ (25) $$

When $u=0$ this cubic system has two stable fixed points at $x=\pm 1$ corresponding to forward and reversal motions, as well as an unstable point at $\beta$. The fixed points correspond to locations in the state space where $F(x,\beta)=0$ and $u(t)=0$. Transitions to the other stable fixed point occur when $u(t)\neq 0$, corresponding to a shift in the behavior of the animal. To capture the stochasticity of the data, additional noise terms are added to $\dot{x},\dot{y}$.
The model is fit to reproduce the dominant PCA mode of the neural activity data. Importantly, the optimization is done by using the manually annotated behavioral labels to determine when the control switches values. For forward and reversal motion, $u=0$, while each type of turn (reversal to forward and vice versa) corresponds to a different $u$ value. These distinct models are fitted to the corresponding time series segments based on the annotations. The resulting nonlinear control model and the parameters found through optimization is fully described in [].
Both our model and the global nonlinear model with control described above employ nonlinear terms in the dynamics. These are global models with few parameters that capture the most important qualitative features in the C. elegans data.
A key difference is that the control variable in Morrison et al. takes $3$ values, including 0 0 during stable states, while dynamic SINDy’s switching variable takes $2$ values that influence the fixed points and are longer-lasting than the transient controls.
A key advantage of our method is that it is entirely data-driven, requiring no behavioral annotations or manual fitting. We directly input the low-dimensional neural activity time series, allowing dynamic SINDy to automatically discover the governing equations. This reduces the effort required from the data scientist while still capturing the system’s essential dynamics.

#### 3.3.3 Comparison with a linear model with control

A related study proposed a global linear model with control
whereby a linear dynamical system is actuated by temporally sparse control signals []. Denoting $\textbf{x}_{j}=\textbf{x}(t_{j})$, neural activity across neurons at time $t_{j}$, and X a matrix of neuronal data at different snapshots in time, $\textbf{X}=[\textbf{x}_{1},...,\textbf{x}_{m}]$, dynamic mode decomposition (DMD) provides a linear model for the dynamics of the state space:

$$ $\displaystyle\textbf{X}^{\prime}=\textbf{A}\textbf{X}$ (26) $$

where $\textbf{X}^{\prime}=[\textbf{x}_{2},\textbf{x}_{3},...\textbf{x}_{m+1}]$ is offset by one time step compared to X. Since a linear model alone cannot capture the neural activity data, DMD with control (DMDc, []) is employed to distinguish between the underlying dynamics and control signals $\textbf{U}=[\textbf{u}_{1},\textbf{u}_{2},…\textbf{u}_{m}]$, where $\textbf{u}_{j}=\textbf{u}(t_{j})$ are actuation signals at a snapshot in time. DMDc regresses to the linear control system:

$$ $\displaystyle\textbf{X}\textquoteright=\textbf{A}\textbf{X}+\textbf{B}\textbf{U}$ (27) $$

The control signal can either be fixed using manually annotated behavioral onsets in a supervised setting or learned jointly with $A$ and $B$ (Algorithm 1, []). To avoid trivial solutions, the control signals are constrained to be sparse, meaning transitions between states should be infrequent. The following loss function, incorporating an $l_{0}$ regularization, is minimized using the sequential least squares thresholding algorithm:

$$ $\displaystyle\textrm{loss}=\min_{\textbf{A},\textbf{B},\textbf{U}}=||\textbf{A }\textbf{X}+\textbf{B}\textbf{U}-\textbf{X}^{\prime}||_{2}+\lambda||\textbf{U} ||_{0}$ (28) $$

If control signals are internally generated, they are either random or encoded within the network. Sparse variable selection and time-delay embeddings test the influence of present and past data to determine which neurons predict the controls:

$$ $\displaystyle\textbf{u}_{k}=\textbf{K}_{1}\textbf{x}_{k}+\textbf{K}_{2}\textbf {x}_{k-1}+…$ (29) $$

Key findings from this study include that the unsupervised algorithm produces control signals somewhat correlating with manually annotated behavioral onsets, though it fails to capture forward motion onsets. This suggests that neurons involved in forward motion (^1^11This comparison between model and data at the single neuron level is only possible due to identification of neurons with stereotyped identities, as described in []) require non-trivial nonlinearities throughout the time series for full reconstruction, not just control signals at the onset. Reversal neurons are well-modeled by the supervised control signals, implying fewer required nonlinearities other than the transition signal itself. Turns are also largely captured by the neuronal activity in specific cells, although there is much more variability.
The global linear framework, with internally generated control signals, partially explains neuronal activity but produces a weak qualitative and quantitative match with the data. For instance, the correlation between real neuronal activity and the model is at or only slightly above $0.5$, even for neurons not implicated in forward motion, that the model supposedly is successful in capturing. Incorporating manual annotations significantly improves the model’s accuracy, but it shifts the approach from data-driven to supervised.
The inability to capture forward motion and the ineffectiveness of forward control signals demonstrate that nonlinearities are necessary for many interactions in the system. A linear model with control can produce only one fixed point, with all other states simply being longer-lived, conflicting with empirical evidence that both forward and reversal behaviors in C. elegans are stable states []. Empirical studies have instead shown that multiple behavioral states appear to be stable [].
These results support our approach using dynamic SINDy, which automatically identifies nonlinear systems with multiple fixed points, without requiring manual annotations. Dynamic SINDy discovers a dynamical system with switching signals as opposed to controls, while incorporating nonlinearities that could be critical in capturing complex neural dynamics. Furthermore, our framework, as described in Sec. 4.4, could discover how transition signals depend on the data to create a fully closed-loop feedback system. Future work should confirm that our low-dimensional model can explain neural activity according to cell class.

## 4 Other methods for system identification of non-autonomous dynamic systems

### 4.1 Switching Linear Dynamical System (SLDS)

The generative model for switching linear dynamical systems (SLDS) is as follows: for each time $t=1,2,...,T$, there is a discrete latent state $z_{t}\in\{1,2,...,K\}$ that follows Markov dynamics:

$$ $z_{t+1}|z_{t},\{\pi_{k}\}_{k=1}^{K}\sim\pi_{z_{t}},$ $$

where $\{\pi_{k}\}_{k=1}^{K}$ is the Markov transition matrix and $\pi_{k}\in[0,1]^{K}$ is the k-th matrix row. In addition to $z_{t}$, there is a continuous latent state $x_{t}\in\mathbb{R}^{M}$ following linear dynamics that depend on $z_{t}$:

$$ $\displaystyle x_{t+1}=A_{z_{t+1}}x_{t}+b_{z_{t+1}}+v_{t},v_{t}\sim N(0,Q_{z_{t }})$ (30) $$

where $A_{z_{t}},Q_{z_{t}}$ are matrices and $b_{z_{t}}$ is a vector depending on the latent state $z_{k}\in{1,...,K}$. In addition, we have access to observables $y_{t}$, generated from the continuous latent state $x_{t}$:

$$ $\displaystyle y_{t}=Cx_{t}+d+w_{t},w_{t}\sim N(0,S)$ (31) $$

where $C,S,d$ are shared matrices and a vector across different discrete states $z_{t}$. We denote the complete set of parameters as $\theta=\{\pi_{k},A_{k},Q_{k},b_{k},C,S,d|k=1,...,K\}$ and learn SLDS using Bayesian inference and a set of convenient priors as detailed in [].
An extension of SLDS – rsLDS – allows the discrete switches to depend on the continuous latent state and exogeneous inputs through a logistic regression []. Specifically, when a discrete switch occurs whenever a continuous state enters a particular region of state space, SLDS is unable to learn this dependence, while rSLDS designed to address this state dependence. An important contribution of [] is an inference algorithm leveraging Polya-gamma auxiliary variable methods to make inference fast, scalable, and, easy. We make use of this algorithm through the open-source rSLDS libraries [].

#### 4.1.1 Training

To train the SLDS/rSLDS models, we first select the number of latent states in advance. For switch-like underlying coefficients, we choose two states to maximize SLDS’s ability to identify the switching dynamics. We then perform principal component analysis (PCA) on the data, followed by fitting a simpler autoregressive hidden Markov model (AR-HMM), which lacks continuous latent states. The discrete latent variables $z$ inferred from the AR-HMM, along with the data projected onto a small number of principal components, are used to initialize the SLDS/rSLDS algorithms. The SLDS/rSLDS training algorithms are implemented from the following github repository: [https://github.com/slinderman/recurrent-slds](https://github.com/slinderman/recurrent-slds).

### 4.2 A method based on group sparsity

We focus on two studies [, ] that address system identification in non-autonomous systems. Both studies propose a method for identifying ODEs with time-varying SINDy coefficients by dividing the trajectory into smaller time windows and applying SINDy to each segment while maintaining the same sparsity pattern across all segments. In [], this sparsity pattern is enforced using group sparsity regularization, while [] introduces a novel algorithm based on sequential thresholding least squares (STLQ) from the original SINDy paper []. This STLQ adaptation averages the SINDy coefficients across time windows and compares the average to a threshold, setting coefficients below this threshold to zero.
The group sparse penalized method for model selection and parameter estimation is used with datasets of multiple trajectories that share the same physical laws, but differ in bifurcation parameters []. This framework is subsequently adapted to switching systems, whereby in a Lorenz system, the parameter $\alpha$ changes from $-1$ to $6.66$ at some unknown time. The framework matches our problem, therefore we adapt the algorithm proposed in this analysis to the non-autonomous dynamical systems we study (Sec. 3.3).
Adapting the notation in [] to our own, we have a total of $M$ time windows that partition the trajectory, and we denote time windows by $i$. Data points from specific time windows are indexed by superscript $i$, while different variables of the system are denoted by subscript $j$, $j=1,...,n$. For instance, in the Lorenz system variables $x,y,z$ correspond to $x_{j}$, $j=1,2,3$, where $x_{j}^{i}$ corresponds to variable $x_{j}$ within time window $i$. We can then define variable $\Xi_{j}$ in terms of $\xi_{j}^{i}$ which are SINDy coefficients corresponding to variable $j$ within time window $i$:

$\Xi_{j}=\begin{bmatrix}|&|&&|\\
\xi_{j}^{1}&\xi_{j}^{2}&...&\xi_{j}^{M}\\
|&|&&|\end{bmatrix}$
Next we can define the data matrix $X^{(i)}$, the velocity matrix $\dot{X}^{(i)}$, and the dictionary matrix $\Theta^{(i)}$ as:
$X^{i}=\begin{bmatrix}|&|&&|\\
x_{1}^{i}&x_{2}^{i}&...&x_{n}^{i}\\
|&|&&|\par\end{bmatrix}=\begin{bmatrix}x_{1}(t_{1};\lambda^{(i)})&x_{2}(t_{1};
\lambda^{(i)})&...&x_{n}(t_{1};\lambda^{(i)})\\
x_{1}(t_{2};\lambda^{(i)})&x_{2}(t_{2};\lambda^{(i)})&...&x_{n}(t_{2};\lambda^
{(i)})\\
...&...&...&...\\
x_{1}(t_{l_{i}};\lambda^{(i)})&x_{2}(t_{l_{i}};\lambda^{(i)})&...&x_{n}(t_{l_{
i}};\lambda^{(i)})\end{bmatrix}$
$\dot{X}^{i}=\begin{bmatrix}|&|&&|\\
\dot{x}_{1}^{i}&\dot{x}_{2}^{i}&...&\dot{x}_{n}^{i}\\
|&|&&|\par\end{bmatrix}=\begin{bmatrix}\dot{x}_{1}(t_{1};\lambda^{(i)})&\dot{x
}_{2}(t_{1};\lambda^{(i)})&...&\dot{x}_{n}(t_{1};\lambda^{(i)})\\
\dot{x}_{1}(t_{2};\lambda^{(i)})&\dot{x}_{2}(t_{2};\lambda^{(i)})&...&\dot{x}_
{n}(t_{2};\lambda^{(i)})\\
...&...&...&...\\
\dot{x}_{1}(t_{l_{i}};\lambda^{(i)})&\dot{x}_{2}(t_{l_{i}};\lambda^{(i)})&...&
\dot{x}_{n}(t_{l_{i}};\lambda^{(i)})\end{bmatrix}$
and
$\Theta^{i}=\begin{bmatrix}\textbf{1}_{l_{i},1}&X^{i}&(X^{i})^{2}&(X^{i})^{3}&.
..\par\end{bmatrix}$

Using this notation, the optimization problem can be rewritten as a least-square fitting :

$$ $\displaystyle\min_{\Xi_{j}}\sum_{i=1}^{n}||\Theta^{i}\xi_{j}^{i}-\dot{X}_{j}^{ i}||_{2}^{2}$ (32) $$

for each $j=1,...,n$.

To help prevent overfitting, we add regularization to this cost function, by including a penalty on the number of active candidate functions. The main assumption of this method is that coefficients $\xi^{i}$ have the same support set (in $j$) for each $i$, but can differ in value. Thus we can group each row $j$ together to be either zero or nonzero, therefore the number of active (nonzero) rows in $\Xi_{j}$ is sparse. The cost function now implements the following group-sparse optimization problem:

$$ $\displaystyle\min_{\Xi_{j}}\sum_{i=1}^{m}||\Theta^{i}\xi_{j}^{i}-\dot{X}_{j}^{ i}||_{2}^{2}+\lambda||\Xi_{j}||_{2,0}$ (33) $$

where the $l_{2,0}$ penalty is defined as:

$$ $\displaystyle||A||_{2,0}=\#\{k:(\sum_{l=1}^{m}|a_{kl}|^{2})^{1/2}\neq 0\}$ (34) $$

for any matrix $A=[a_{k,l}]$. Although the problem is non-convex, the authors in [] propose to solve it numerically using an iterative hard thresholding algorithm, the group hard-iterative thresholding algorithm for dynamical systems:
Group Hard-Iterative Thresholding Algorithm for Dynamical Systems:

where $supp(x)$ is the support set of $x$, i.e. the indices of $x$ that correspond to the nonzero values.
While we implement this algorithm and test it on the data, we have found, surprisingly, that a simpler algorithm is more effective in many cases:
Simple sequential thresholding algorithm:

Comments starting with $``\%"$ are provided throughout the code to clarify its meaning. We use this algorithm throughout to showcase our results using the group sparsity method.
For training, we have varied the total time for the trajectories, the number of batches used, the time window, as well as the precise algorithm used. Throughout these experiments we have found that the algorithm was highly sensitive to whether the data was normalized or not, specifically we have found that not normalizing the data yielded superior results.

Supplementary Figures

Figure: Figure 1: (A) Schematic of timeVAE architecture from []; (B) Schematic of dynamic HyperSINDy architecture described in SM Sec. 1.2.2
Refer to caption: x8.png

Figure: Figure 2: Time-varying coefficients (above) and corresponding dynamics (below) for the non-autonomous harmonic oscillator of Eq. (3) (main text). Complementary to Figure 1A (main text)
Refer to caption: x9.png

Figure: Figure 3: Data-driven discovery of sinusoid SINDy coefficients (a and b)and trajectory reconstruction (c) of a non-autonomous harmonic oscillator with dynamic HyperSINDy.
Refer to caption: x10.png

Figure: Figure 4: Data-driven discovery of sigmoid and switch signal coefficients (a and c) and the corresponding trajectory reconstruction (b and d) of a non-autonomous harmonic oscillator with dynamic HyperSINDy.
Refer to caption: x11.png

Figure: Figure 5: Data-driven discovery of sigmoid (a), constant (b, c and e) and sinusoid (e) SINDy coefficients and trajectory reconstruction (d) of the Lorenz dynamics with dynamic HyperSINDy. Time series of coefficients that correspond to constants in the real dynamics sometimes inherit frequency content from the dynamics, as in (e).
Refer to caption: x12.png

Figure: Figure 6: Trajectory reconstructions with timeVAE for a non-autonomous harmonic oscillator with different coefficients (sigmoid, switch signal, finite Fourier series) and different levels of noise in the coefficients ($0.01$, $0.1$, $1$).
Refer to caption: x13.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 2
questions: - For the switch signals (Fig. 2 a-c, also Fig. 3A low noise setting), the inferred ODE parameter time series seem to exhibit high frequency oscillations on top of the correct switch-like dynamics. Is there an intuitive explanation why the encoder-decoder architecture struggles in inferring the correct switching dynamics and how this could be addressed?
- Results of Fig. 3B look rather weak to me, can the authors report Pearson’s $r$ of noise lvl vs. std?
- I’m confused by section 4.6 & Fig. 7; How exactly does the dynamic SINDy approach compare now to the proposed baseline methods based on SLDS and (vanilla?) SINDy with a group sparsity norm? I think Fig. 7 would be much clearer if the authors would find a design to compare all comparison methods side-by-side.
- ll. 409-411: Can the authors provide references for the mentioned studies?
- How do other methods like reservoir computing compare to the dynamic SINDy approach qualitatively and quantitatively in the settings discussed in the manuscript (see e.g. [1])?
- How does the approach perform on e.g. benchmarks used in [2], which exhibit different bifurcations than the ones discussed in this paper?

I am very happy to increase my score if the authors adequately address my concerns and questions.

References:

[1] Köglmayr, Daniel, and Christoph Räth. "Extrapolating tipping points and simulating non-stationary dynamics of complex systems using efficient machine learning." Scientific Reports 14.1 (2024): 507.

[2] Patel, Dhruvit, and Edward Ott. "Using machine learning to anticipate tipping points and extrapolate to post-tipping dynamics of non-stationary dynamical systems." Chaos: An Interdisciplinary Journal of Nonlinear Science 33.2 (2023).
soundness: 3
strengths: - Learning a parsimonious representation of non-autonomous DS is extremely important and relevant in many scientific disciplines, which makes the approach very promising.
- I think it is highly interesting that the encoder-decoder architecture is able to predict the ODE parameters with this level of fidelity in an unsupervised fashion (as there is no direct reconstruction loss for the ODE parameter time series involved in the loss function (7)).
- The method is tested against other baselines and also on a real-world dataset (C. elegans).
summary: The paper proposes ‘dynamic sparse identification of nonlinear dynamics’ (dynamic SINDy), a deep learning framework for identifying governing equations in noisy, non-stationary and nonlinear dynamical systems (DS). By combining variational autoencoders (VAEs) and previous work on SINDy, it enables unsupervised inference of the underlying ODE systems’ parameters while extracting a global and parsimonious nonlinear dynamical model. The approach is validated on both synthetic and real-world data and is compared to other methods in the field, demonstrating great potential for scientific machine learning community.
weaknesses: - The authors should stick to the ICLR style guide and use the 'author, year' reference style instead of mere numerical numbers (i.e. APA style instead of IEEE). This increases readability and helps the reader to understand the train of thought of the authors, as one directly sees on which work the authors base certain statements.
- Center box in Fig. 1B is in parts hard to read as (font) sizes vary a lot. I think it would be better to shrink down Fig. 1A a touch and to increase size of Fig. 1B, especially as it describes the main framework of the manuscript.
- I also think the figure group titles ('suptitles') above Fig. 1, 2, and 4 are superfluous and their message should be put into the figure caption. This would create additional space (e.g. to compensate for the change in referencing style).
- I think ‘dynamic HyperSINDy’ deserves a bit more attention in the main text, which lacks explanation on how this approach really works. Explaining this method in the supplement makes the corresponding results a bit hard to read and almost forces the reader to read the supplement section 1.2.2.
- All of the employed (benchmark) datasets are fairly low dimensional (2-3D). The authors do not address the scalability of the method to high dimensional systems (which can not be sufficiently described by the first few PCA components). I think this is a major drawback, as this setting is highly relevant to many real-world systems.

Minor details:
- typo: Fig. 3B y-axis label say “approxiate std”
- l. 353 It just says 6A and 6B, while the authors probably reference Fig. 5A and 5B? Also l. 360 it says 6C instead of 5C.
- typo: supplement l. 262 it says weight decay of 1e5 (I assume 1e-5?)

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: have the authors applied their method to any datasets requiring a higher dimensional latent space to see how quality of the learned dynamical system scales with dimensionality of the latent space? 

have the authors considered comparisons to a method more adept at handling more smooth like transitions of dynamics such as [1] which considers a smoothly switching latent system.  

[1] Kurle, Richard, et al. "Deep rao-blackwellised particle filters for time series forecasting." Advances in Neural Information Processing Systems 33 (2020): 15371-15382.
soundness: 3
strengths: dynamical system reconstruction or system identification is an important topics, and learning more interpretable models of system dynamics, such as time varying ode representations like the authors, has broad application to many fields.  additionally, the authors consider a comprehensive amount of toy experiments — i appreciate that the authors considered experiments with several time varying motifs (i.e. switching, sigmoid, switch, fourier) and show that their method can recover the time varying coefficients with a calibrated measure of uncertainty.
summary: the authors introduce dynamic sindy — an extension of the sindy for uncovering nonstationary dynamical systems from data.  the authors use a time-series VAE architecture to map their data to time-varying coefficients that are linearly combined with a fixed library of basis functions to produce an estimate of the data derivative.   they conduct experiments using several toy datasets showing that their dynamic sindy can recover the coefficients of time-varying dynamical systems, even in the case that the entire system is not observed.  finally, they show on low-dimensional representations of c. elegans neural recordings, that their method  recovers a representative dynamical system of the first principal component.
weaknesses: i found the separation between what has previously been done in the literature and what are the authors exact main novel contributions to be unclear; for example, more precise statements about the differences with hypersindy (around line 096) would have been very helpful.  in its current form, it is hard to parse from the manuscript what their exact technical advances are.

a lot of real-estate in the paper is taken up by the experimental plots.  often i found the amount of information conveyed by the plots disproportionate to the amount of space they take up — making more compact figures seems like it would work to the authors advantage.  additionally, information could be conveyed better i.e. thick lines and their ordering (i.e. green/blue lines in Fig 2d are not clear).  Fig 1B has a lot of small labels and the zoomed out view of timeVAE does not feel like it helps much.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: - Variational inference is known to provide overconfident uncertainty estimates because of the KL term encouraging mode-covering behavior. Can the authors discuss this further and the impact of this on the method?
- How robust are the identified parameters? Is there a quantitative relationship between the robustness and (1) the size of the library, (2) the level of noise in the system?
- It is stated that the method can only deal with non-stationarity arising from separable time-varying variables. Can the authors elaborate more on why this is the case?
soundness: 2
strengths: - The authors address an important and often overlooked problem in modeling time series data: learning interpretable models from non-stationary dynamical systems.
- The problem and the proposed method are presented clearly, and the paper is generally easy to follow.
summary: The authors propose DynamicSINDY, an approach that aims to learn sparse ODEs with time-varying coefficients from noisy, non-stationary time series data. The authors achieve this by combining SINDy with sequential VAEs to probabilistically infer the coefficients and their time-varying values. The method is evaluated on three synthetic datasets and a calcium imaging dataset of C. elegans.
weaknesses: - While the motivation of the problem is clear, I find that the conducted experiments fail to convince the reader of the method's impact in real-world settings. The experiments mainly focus on synthetic datasets that are artificially created to fit this problem and avoid most challenges often encountered in real-world datasets (high dimensionality, non-Gaussian noise, large space of possible coefficients). I think the paper would be greatly benefit if the authors can demonstrate the method in such contexts. 

- The C. elegans dataset used to demonstrate the method is quite simple, and the method is only applied to low-dimensional representations obtained via PCA (which in this case is enough to explain the data). The impact of using DynamicSINDY in this case is not well motivated, and the obtained results don't add any much scientific insights, especially since identifiability is not discussed.

- While this is not necessarily always a weakness, the proposed method is a straightforward combination of two existing approaches. Taken together with the limited experiments section and the lack of significant technical innovations, I find the overall contribution of the paper in its current form rather limited.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: - What can one do with the uncertainty? Is it necessary for accurate estimate? Would the uncertainty provide insights for scientific questions? Suggestion: elaborate more or showcase scenarios where the uncertainty is useful vs. method.
- How does the method perform for noise driven dynamics? e.g. system with multiple meta-stable points or line attractors.
- What scientific insights the proposed method could offer for the systems that don't have particular a priori form of ODEs? If we don't know u(t) switches, would it discover that?
- Line 461 typo: rLDS
- Suggestion: put Dynamics SINDy reconstruction together with SLDS result in Fig. 7.
soundness: 3
strengths: - Modeling Time-Varying Coefficients: Dynamic SINDy employs deep generative models, specifically variational autoencoders (VAEs), to learn the time-varying nature of ODE coefficients. This enables the identification of non-autonomous systems that exhibit complex dynamics.
- Uncertainty Quantification: The use of VAEs allows dynamic SINDy to quantify uncertainty in the estimated ODE coefficients. 
- Latent Variable Discovery: Dynamic SINDy can effectively uncover hidden (latent) variables that influence the observed dynamics. This is demonstrated through an example with the Lotka-Volterra equations, where only prey population data is available.
- Application to Real-World Data: The paper validates dynamic SINDy's capabilities on both synthetic and real-world data, i.e. the C. elegans data.
summary: This paper presents dynamic SINDy, a machine learning method designed to identify the governing equations of noisy, nonlinear, and non-autonomous dynamical systems from data. Dynamic SINDy combines variational inference with sparse identification of nonlinear dynamics to identify time-varying coefficients in sparse ordinary differential equations. The method is particularly valuable for non-stationary systems with changing parameters over time. The contributions include
- Modeling Time-Varying Coefficients: Dynamic SINDy employs deep generative models, specifically variational autoencoders (VAEs), to learn the time-varying nature of ODE coefficients. This enables the identification of non-autonomous systems that exhibit complex dynamics.
- Uncertainty Quantification: The use of VAEs allows dynamic SINDy to quantify uncertainty in the estimated ODE coefficients. This is crucial for understanding the reliability of the identified model.
weaknesses: - Uncertainty Quantification: though sec 4.2 shows the estimated standard deviation follows the truth, it only confirms the accuracy of estimation. Nonetheless, it is unclear what use the uncertainty could provide.
- The examples are mostly using systems that driven by mean/input processes. It's unclear how the proposed method would perform for noise driven dynamics.
- In the C. elegans example, the assumed form depends heavily on prior knowledge (dimensionality, input) on the dynamics. Though the proposed method has shown very good performance, it's unclear what scientific insights the proposed method could offer especially for the systems that people have limited knowledge about.