Modeling Latent Neural Dynamics with Gaussian Process Switching Linear Dynamical Systems

Authors: amber  hu
david  zoltowski
aditya  nair
david  anderson
lea  duncker
scott  linderman
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 08:42:00

Abstract: Understanding how the collective activity of neural populations relates to computation and ultimately behavior is a key goal in neuroscience. To this end, statistical methods which describe high-dimensional neural time series in terms of low-dimensional latent dynamics have played a fundamental role in characterizing neural systems. Yet, what constitutes a successful method involves two opposing criteria: (1) methods should be expressive enough to capture complex nonlinear dynamics, and (2) they should maintain a notion of interpretability often only warranted by simpler linear models. In this paper, we develop an approach that balances these two objectives: the Gaussian Process Switching Linear Dynamical System (gpSLDS). Our method builds on previous work modeling the latent state evolution via a stochastic differential equation whose nonlinear dynamics are described by a Gaussian process (GP-SDEs). We propose a novel kernel function which enforces smoothly interpolated locally linear dynamics, and therefore expresses flexible -- yet interpretable -- dynamics akin to those of recurrent switching linear dynamical systems (rSLDS). Our approach resolves key limitations of the rSLDS such as artifactual oscillations in dynamics near discrete state boundaries, while also providing posterior uncertainty estimates of the dynamics. To fit our models, we leverage a modified learning objective which improves the estimation accuracy of kernel hyperparameters compared to previous GP-SDE fitting approaches. We apply our method to synthetic data and data recorded in two neuroscience experiments and demonstrate favorable performance in comparison to the rSLDS.

Main paper: 
 introduction: 
 Computations in the brain are thought to be implemented through the dynamical evolution of neural activity. Such computations are typically studied in a controlled experimental setup, where an animal is engaged in a behavioral task with relatively few relevant variables. Consistent with this, empirical neural activity has been reported to exhibit many fewer degrees of freedom than there are neurons in the measured sample during such simple tasks [1] . These observations have driven the use of latent variable models to characterize low-dimensional structure in high-dimensional neural population activity [2, 3] . In this setting, neural activity is often modeled in terms of a low-dimensional latent state that evolves with Markovian dynamics [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] . It is thought that the latent state evolution is related to the computation of the system, and therefore, insights into how this evolution is shaped through a dynamical system can help us understand the mechanisms underlying computation [15] [16] [17] [18] [19] [20] [21] .In practice, choosing an appropriate modeling approach for a given task requires balancing two key criteria. First, statistical models should be expressive enough to capture potentially complex and nonlinear dynamics required to carry out a particular computation. On the other hand, these models should also be interpretable and allow for straightforward post-hoc analyses of dynamics. One model class that strikes this balance is the recurrent switching linear dynamical system (rSLDS) [8] . The rSLDS approximates arbitrary nonlinear dynamics by switching between a finite number of linear dynamical systems. This leads to a powerful and expressive model which maintains the interpretability of linear systems. Because of their flexibility and interpretability, variants of the rSLDS have been used in many neuroscience applications [19, [22] [23] [24] [25] [26] [27] [28] and are examples of a general set of models aiming to understand nonlinear dynamics using compact and interpretable components [29, 30] .However, rSLDS models suffer from several limitations. First, while the rSLDS is a probabilistic model, typical use cases do not capture posterior uncertainty over inferred dynamics. This makes it difficult to judge the extent to which particular features of a fitted model should be relied upon when making inferences about their role in neural computation. Second, the rSLDS often suffers from producing oscillatory dynamics in regions of high uncertainty in the latent space, such as boundaries between linear dynamical regimes. This artifactual behavior can significantly impact the interpretability and predictive performance of the rSLDS. Lastly, the rSLDS does not impose smoothness or continuity assumptions on the dynamics due to its discrete switching formulation. Such assumptions are often natural and useful in the context of modeling realistic neural systems.In this paper, we improve upon the rSLDS by introducing the Gaussian Process Switching Linear Dynamical System (gpSLDS). Our method extends prior work on the Gaussian process stochastic differential equation (GP-SDE) model, a continuous-time method that places a Gaussian process (GP) prior on latent dynamics. By developing a novel GP kernel function, we enforce locally linear, interpretable structure in dynamics akin to that of the rSLDS. Our framework addresses the aforementioned modeling limitations of the rSLDS and contributes a new class of priors in the GP-SDE model class. Our paper is organized as follows. Section 2 provides background on GP-SDE and rSLDS models. Section 3 presents our new gpSLDS model and an inference and learning algorithm for fitting these models. In Section 4 we apply the gpSLDS to a synthetic dataset and two datasets from real neuroscience experiments to demonstrate its practical use and competitive performance. We review related work in Section 5 and conclude our paper with a discussion in Section 6.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: In our abstract and introduction, we propose a new method for interpretable modeling of neural dynamics. In our paper, we apply our method and demonstrate its interpretability and competitive performance on a synthetic dataset and two real datasets.• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
gaussian process stochastic differential equation models: 
 Gaussian processes (GPs) define nonparametric distributions over functions. They are a popular choice in machine learning due to their ability to capture nonlinearities and encode reasonable prior assumptions such as smoothness and continuity [31] . Here, we review the GP-SDE, a Bayesian generative model that leverages the expressivity of GPs for inferring latent dynamics [10] .Generative model In a GP-SDE, the evolution of the latent state x ∈ R K is modeled as a continuous-time SDE which underlies observed neural activity y(t i ) ∈ R D at time-points t i ∈ [0, T ]. Mathematically, this is expressed asEQUATIONThe drift function f : R K → R K describes the system dynamics, Σ is a noise covariance matrix, and dw ∼ N (0, dtI) is a Wiener process increment. Parameters C ∈ R D×K and d ∈ R D define an affine mapping from latent to observed space, which is then passed through a pre-specified inverse link function g(•).A GP prior is used to model each output dimension of the dynamics f (•) independently. More formally, if f (•) = [f 1 (•), . . . , f K (•)] T , thenEQUATIONwhere κ Θ (•, •) is the kernel for the GP with hyperparameters Θ.Interpretability GP-SDEs and their variants can infer complex nonlinear dynamics with posterior uncertainty estimates in physical systems across a variety of applications [32] [33] [34] . However, one limitation of using this method with standard GP kernels, such as the radial basis function (RBF) kernel, is that its expressivity leads to dynamics that are often challenging to interpret. In Duncker et al. [10] , this was addressed by conditioning the GP prior of the dynamics f (•) on fixed points f (x * ) = 0 and their local Jacobians J(x * ) = ∂ ∂x f (x)| x=x * , and subsequently learning the fixed-point locations x * and the locally-linearized dynamics J(x * ) as model parameters. This approach allows for direct estimation of key features of f (•). However, due to its flexibility, it is also prone to finding more fixed points than those included in the prior conditioning, which undermines its overall interpretability.The rSLDS models nonlinear dynamics by switching between different sets of linear dynamics [8] . Accordingly, it retains the simplicity and interpretability of linear dynamical systems while providing much more expressive power. For these reasons, variations of the rSLDS are commonly used to model neural dynamics [19, [22] [23] [24] [25] [26] [27] [28] [29] .The rSLDS is a discrete-time generative model of the following form:EQUATIONwhere dynamics switch between J distinct linear systems with parameters {A j , b j , Q j } J j=1 . This is controlled by a discrete state variable s t ∈ {1, . . . , J}, which evolves via transition probabilities modeled by a multiclass logistic regression,EQUATIONThe "recurrent" nature of this model comes from the dependence of eq. ( 4) on latent space locations x t .As such, the rSLDS can be understood as learning a partition of the latent space into J linear dynamical regimes seprated by linear decision boundaries. This serves as important motivation for the parametrization of the gpSLDS, as we describe later.Interpretability While the rSLDS has been successfully used in many applications to model nonlinear dynamical systems, it suffers from a few practical limitations. First, it often produces unnatural artifacts of modeling nonlinear dynamics with discrete switches between linear systems. For example, it may oscillate between discrete modes with different discontinuous dynamics when a trajectory is near a regime boundary. Next, common fitting techniques for rSLDS models with non-conjugate observations typically treat dynamics as learnable hyperparameters rather than as probabilistic quantities [26] , which prevents the model from being able to capture posterior uncertainty over the learned dynamics. Inferring a posterior distribution over dynamics is especially important in many neuroscience applications, where scientists often draw conclusions from discovering key features in latent dynamics, such as fixed points or line attractors.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: We discuss limitations of our work in the discussion section of the paper.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
gaussian process switching linear dynamical systems: 
 To address these limitations of the rSLDS, we propose a new class of models called the Gaussian Process Switching Linear Dynamical System (gpSLDS). The gpSLDS combines the modeling advantages of the GP-SDE with the structured flexbility of the rSLDS. We achieve this balance by designing a novel GP kernel function that defines a smooth, locally linear prior on dynamics.While our main focus is on providing an alternative to the rSLDS, the gpSLDS also contributes a new prior which allows for more interpretable learning of dynamics and fixed points than standard priors in the GP-SDE framework (e.g., the RBF kernel). Our implementation of the gpSLDS is available at: https://github.com/lindermanlab/gpslds. The core innovation of our method is a novel GP kernel, which we call the Smoothly Switching Linear (SSL) kernel. The SSL kernel specifies a GP prior over dynamics that maintains the switching linear structure of rSLDS models, while allowing dynamics to smoothly interpolate between different linear regimes.For every pair of locations x, x ′ ∈ R K , the SSL kernel with J linear regimes is defined as,EQUATIONwhere c j ∈ R K , M ∈ R K×K is a diagonal positive semi-definite matrix, σ 2 0 ∈ R + , and π j (x) ≥ 0 with J j=1 π j (x) = 1. To gain an intuitive understanding of the SSL kernel, we will separately analyze each of the two terms in the summands.The first term, κ (j) lin (x, x ′ ), is a standard linear kernel which defines a GP distribution over linear functions [31] . The superscript j denotes that it is the linear prior on the dynamics in regime j. M controls the variance of the function's slope in each input dimension, and c j is such that the variance of the function achieves its minimum value of σ 2 0 at x = c j . We expand on the relationship between the linear kernel and linear dynamical systems in Appendix A.The second term is what we define as the partition kernel, κ (j) part (x, x ′ ), which gives the gpSLDS its switching structure. We interpret π(x) = [π 1 (x) . . . π J (x)]T as parametrizing a categorical distribution over J linear regimes akin to the discrete switching variables in the rSLDS. We model π(x) as a multiclass logistic regression with decision boundaries {w T j ϕ(x) = 0 | j = 1, . . . , J -1}, where ϕ(x) is any feature transformation of x. This yields random functions which are locally constant and smoothly interpolate at the decision boundaries. More formally,EQUATIONwhere w J = 0. The hyperparameter τ ∈ R + controls the smoothness of the decision boundaries. As τ → 0, π(x) approaches a one-hot vector which produces piecewise constant functions with sharp boundaries, and as τ → ∞ the boundaries become more uniform. While we focus on the parametrization in eq. ( 6) for the experiments in this paper, we note that in principle any classification method can be used, such as another GP or a neural network.The SSL kernel in eq. ( 5) naturally combines aspects of the linear and partition kernels via sums and products of kernels, which has an intuitive interpretation [35] . The product kernel κ(j) lin (x, x ′ )κ (j)part (x, x ′ ) enforces linearity in regions where π j (x) is close to 1. Summing over J regimes then enforces linearity in each of the J regimes, leading to a prior on locally linear functions with knots determined by π(x). We note that our kernel is reminiscent of the one in Pfingsten et al. [36] , which uses a GP classifier as a prior for π(x) and applies their kernel to a GP regression setting. Here, our work differs in that we explicitly enforce linearity in each regime and draw a novel connection to switching models like the rSLDS.Figure 1A depicts samples from each kernel in 1D. Figure 1B shows how a SSL kernel with J = 2 linear regimes separated by decision boundary x 2 1 + x 2 2 = 4 (top) produces a structured 2D flow field consisting of two linear systems, with x 1 -and x 2 -directions determined by 1D function samples (bottom).The full generative model for the gpSLDS incorporates the SSL kernel in eq. ( 5) into a GP-SDE modeling framework. Instead of placing a GP prior with a standard kernel on the system dynamics as in eq. ( 2), we simply plug in our new SSL kernel so thatEQUATIONwhere the kernel hyperparameters are Θ = {M , σ 2 0 , {c j } J j=1 , {w j } J-1 j=1 , τ }. We then sample latent states and observations according to the GP-SDE via eq. (1) . A schematic of the full generative model is depicted in fig. 1C .Incorporating inputs In many modern neuroscience applications, we are often interested in how external inputs to the system, such as experimental stimuli, influence latent states. To this end, we also consider an extension of the model in eq. ( 1) which incorporates additive inputs of the form,EQUATIONwhere v(t) ∈ R I is a time-varying, known input signal and B ∈ R K×I maps inputs linearly to the latent space. The latent path inference and learning approaches presented in the following section can naturally be extended to this setting, with updates for B available in closed form. Further details are provided in Appendix B.3-B.6.For inference and learning in the gpSLDS, we apply and extend a variational expectationmaximization (vEM) algorithm for GP-SDEs from Duncker et al. [10] . In particular, we propose a modification of this algorithm that dramatically improves the learning accuracy of kernel hyperparameters, which are crucial to the interpretability of the gpSLDS. We outline the main ideas of the algorithm here, though full details can be found in Appendix B.As in Duncker et al. [10] , we consider a factorized variational approximation to the posterior,EQUATIONwhere we have augmented the model with sparse inducing points to make inference of f tractable [37] . The inducing points are located at {z m } M m=1 ⊂ R K and take values (f k (z 1 ), . . . , f k (z M )) ⊤ = u k . Standard vEM maximizes the evidence lower bound (ELBO) to the marginal log-likelihood log p(y | Θ) by alternating between updating the variational posterior q and updating model hyperparameters Θ [38] . Using the factorization in eq. ( 9), we will denote the ELBO as L(q(x), q(u), Θ).For inference of q(x), we follow the approach first proposed by Archambeau et al. [39] and extended by Duncker et al. [10] . Computing the ELBO using this approach requires computing variational expectations of the SSL kernel, which we approximate using Gauss-Hermite quadrature as they are not available in closed form. Full derivations of this step are provided in Appendix B.3. For inference of q(u), we follow Duncker et al. [10] and choose a Gaussian variational posterior for q(u k ) = N (u k | m k * u , S k * u ), which admits closed-form updates for the mean m k * u and covariance S k * u given q(x) and Θ. Duncker et al. [10] perform these updates before updating Θ via gradient ascent on the ELBO in each vEM iteration.In our setting, this did not work well in practice. The gpSLDS often exhibits strong dependencies between q(u) and Θ, which makes standard coordinate-ascent steps in vEM prone to getting stuck in local maxima. These dependencies arise due to the highly structured nature of our GP prior; small changes in the decision boundaries {w j } J-1 j=1 can lead to large (adverse) changes in the prior, which prevents vEM from escaping suboptimal regions of parameter space. To overcome these difficulties, we propose a different approach for learning Θ: instead of fixing q(u) and performing gradient ascent on the ELBO, we optimize Θ by maximizing a partially optimized ELBO,EQUATION)Due to the conjugacy of the model, the inner maximization can be performed analytically. This approach circumvents local optima that plague coordinate ascent on the standard ELBO. While other similar approaches exploit model conjugacy for faster vEM convergence in sparse variational GPs [37, 40] and GP-SDEs [41] , our approach is the first to our knowledge that leverages this structure specifically for learning the latent dynamics of a GP-SDE model. We empirically demonstrate the superior performance of our learning algorithm in Appendix C.It is straightforward to obtain the approximate posterior distribution over f * := f (x * ) evaluated at any new location x * . Under the assumption that f * only depends on the inducing points, we can use the approximationq(f * ) = K k=1 p(f * k | u k , Θ)q(u k) du k which can be computed in closed-form using properties of conditional Gaussian distributions. For a batch of points {x * i } N i=1 , this can be computed in O(N M 2 ) time. The full derivation can be found in Appendix B.4. This property highlights an appealing feature of the gpSLDS over the rSLDS. The gpSLDS infers a posterior distribution over dynamics at every point in latent space, even in regions of high uncertainty. Meanwhile, as we shall see later, the rSLDS expresses uncertainty by randomly oscillating between different sets of most-likely linear dynamics, which is much harder to interpret.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: We have made our codebase on GitHub publicly available. We provide both data and code for our main synthetic data result, and we also provide references to the dataset and instructions for fitting models to reproduce the second real data result. We do not provide data or code for the first real data result since that data has not been released to the public.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
synthetic data: 
 We begin by applying the gpSLDS to a synthetic dataset consisting of two linear rotational systems, one clockwise and one-counterclockwise, which combine smoothly at x 1 = 0 (fig. 2A ). We simulate 30 trials of latent states from an SDE as in eq. ( 1) and then generate Poisson process observations given these latent states for D = 50 output dimensions (i.e. neurons) over T = 2.5 seconds (Fig. 2B ). To initialize C and d, we fit a Poisson LDS [4] to data binned at 20ms with identity dynamics. For the rSLDS, we also bin the data at 20ms. We then fit the gpSLDS and rSLDS models with J = 2 linear states using 5 different random initializations for 100 vEM iterations, and choose the fits with the highest ELBOs in each model class.We find that the gpSLDS accurately recovers the true latent trajectories (fig. 2C ) as well as the true rotation dynamics and the decision boundary between them (fig. 2D ). We determine this decision boundary by thresholding the learned π(x) at 0.5. In addition, we can obtain estimates of fixed point locations by computing the posterior probabilityK k=1 P q(f ) (|f k (x)| < ϵ) for a small ϵ > 0;the locations x with high probability are shaded in purple. This reveals that the gpSLDS finds high-probability fixed points which overlap significantly with the true fixed points, denoted by stars. In comparison, both the rSLDS and the GP-SDE with RBF kernel do not learn the correct decision boundary nor the fixed points as accurately (fig. 2E-F ). Of particular note, the RBF kernel incorrectly extrapolates and finds a superfluous fixed point outside the region traversed by the true latent states. Figure 2G illustrates the differences in how the gpSLDS and the rSLDS express uncertainty over dynamics. To the left of the dashed line, we sample latent states starting from x 0 = (7, 0) and plot the corresponding true dynamics. To the right, we simulate latent states x samp from the fitted model and plot the true dynamics (in gray) and the learned most likely dynamics (in color) at x samp . For a well-fitting model, we would expect the true and learned dynamics to overlap. We see that the gpSLDS produces smooth simulated dynamics that match the true dynamics at x samp (top). By contrast, the rSLDS expresses uncertainty by oscillating between the two linear dynamical systems, hence producing uninterpretable dynamics at x samp (bottom). This region of uncertainty overlaps with the x 1 = 0 boundary, suggesting that the rSLDS fails to capture the smoothly interpolating dynamics present in the true system.Next, we perform quantitative comparisons between the three competing methods (fig. 2H ). We find that both continuous-time methods consistently outperform the rSLDS on both metrics, suggesting that these methods are likely more suitable for modeling Poisson process data. Moreover, the gpSLDS better recovers dynamics compared to the RBF kernel, illustrating that the correct inductive bias can yield performance gains over a more flexible prior, especially in a data-limited setting.Lastly, we note that the gpSLDS can achieve more expressive power than the rSLDS by learning nonlinear decision boundaries between linear regimes, for instance by incorporating nonlinear features into ϕ(x) in eq. ( 6). We demonstrate this feature for a 2D limit cycle in Appendix D.In this section, we revisit the analyses of Nair et al. [27] , which applied dynamical systems models to neural recordings during aggressive behavior in mice. To do this, we reanalyze a dataset which consists of calcium imaging of ventromedial hypothalamus neurons from a mouse interacting with two consecutive intruders. The recording was collected from 104 neurons at 15 Hz over ∼343 seconds (i.e. 5140 time bins). Nair et al. [27] found that an rSLDS fit to this data learns dynamics that form an approximate line attractor corresponding to an aggressive internal state (fig. 3A ). Here, we supplement this analysis by using the gpSLDS to directly assess model confidence about this finding.For our experiments, we z-score and then split the data into two trials, one for each distinct intruder interacting with the mouse. Following Nair et al. [27] , we choose J = 4 linear regimes to compare the gpSLDS and rSLDS. We choose K = 2 latent dimensions to aid the visualization of the resulting model fits; we find that even in such low dimensions, both models still recover line attractor-like dynamics. We model the calcium imaging traces as Gaussian emissions on an evenly spaced time grid and initialize C and d using principal component analysis. We fit models with 5 different initializations for 50 vEM iterations and display the runs with highest forward simulation accuracy (as described in the caption of fig. 3 ).In fig. 3A-B , we find that both methods infer similar latent trajectories and find plausible flow fields that are parsed in terms of simpler linear components. We further demonstrate the ability of the gpSLDS to more precisely identify the line attractor from Nair et al. [27] . To do this, we use the learned q(f ) to compute the posterior probability of slow dynamics on a dense (80×80) grid of points in the latent space using the procedure in Section 3.4. The gpSLDS finds a high-probability region of slow points corresponding to the approximate line attractor found in Nair et al. [27] (fig. 3C ). This demonstrates a key advantage of the gpSLDS over the rSLDS: by modeling dynamics probabilistically using a structured GP prior, we can validate the finding of a line attractor with further statistical rigor. Finally, we compare the gpSLDS, rSLDS, and GP-SDE with RBF kernel using an in-sample forward simulation metric (fig. 3D ). All three methods retain similar predictive power 500 steps into the future. After that, the gpSLDS performs slightly worse than the RBF kernel; however, it gains interpretability by imposing piecewise linear structure while still outperforming the rSLDS.In neuroscience, there is considerable interest in understanding how neural dynamics during decision making tasks support the process of making a choice [26, 34, [42] [43] [44] [45] [46] [47] . In this section, we use the gpSLDS to infer latent dynamics from spiking activity in the lateral intraparietal (LIP) area of monkeys reporting decisions about the direction of motion of a random moving dots stimulus with varying degrees of motion coherence [42, 43] . The animal indicated its choice of net motion direction (either left or right) via a saccade. On some trials, a 100ms pulse of weak motion, randomly oriented to the left or right, was also presented to the animal. Here, we model the dynamics of 58 neurons recorded across 50 trials consisting of net motion coherence strengths in {-.512, -.128, 0.0, .128, .512}, where the sign corresponds to the net movement direction of the stimulus. We only consider data 200ms after motion onset, corresponding to the start of decision-related activity. To capture potential input-driven effects, we fit a version of the gpSLDS described in eq. ( 8) with K = 2 latent dimensions and J = 2 linear regimes over 50 vEM iterations. We encoded the input signal as ±1 with sign corresponding to pulse direction. In Figure 4A , we find that not only does the gpSLDS capture a distinct visual separation between trials by motion coherence, but it also learns a precise separating decision boundary between the two linear regimes. Our finding is consistent with previous work on a related task, which found that average LIP responses can be represented by a 2D curved manifold [48] , though here we take a dynamical systems perspective. Additionally, our model learns an input-driven effect which appears to define a separating axis. To verify this, we project the inferred latent states onto the 1D subspace spanned by the input effect vector. Figure 4B shows that the latent states separate by coherence (top) and by choice (bottom), further suggesting that the pulse input relates to meaningful variation in evidence accumulation for this task. Fig. 4C shows an example latent trajectory aligned with pulse input; during the pulse there is a noticeable change in the latent trajectory. Lastly, in Fig. 4D we find that the gpSLDS expresses high confidence in learned dynamics where latent trajectories are present and low confidence in areas further from this region. 
related work: 
 There are several related approaches to learning nonlinear latent dynamics in discrete or continuous time. Gaussian process state-space models (GP-SSMs) [49] [50] [51] [52] [53] can be considered a discrete-time analogue to GP-SDEs. In a GP-SSM, observations are assumed to be regularly sampled and latent states evolve according to a discretized dynamical system with a GP prior. Wang et al. [49] and Turner et al. [50] learned the dynamics of a GP-SSM using maximum a posteriori estimation. Frigola et al. [51] and Eleftheriadis et al. [52] employed a variational approximation with sparse inducing points to infer the latent states and dynamics in a fully Bayesian fashion. In our work, we use a continuous-time framework that more naturally handles irregularly sampled data, such as point-process observations commonly encountered in neural spiking data. Neural ODEs and SDEs [54, 55] use deep neural networks to parametrize the dynamics of a continuous-time system, and have emerged as prominent tools for analyzing large datasets, including those in neuroscience [33, [56] [57] [58] . While these methods can represent flexible function classes, they are likely to overfit to low-data regimes and may be difficult to interpret. In addition, unlike the gpSLDS, neural ODEs and SDEs do not typically quantify uncertainty of the learned dynamics.In the context of dynamical mixture models, Köhs et al. [59, 60] proposed a continuous-time switching model in a GP-SDE framework. This model assumes a latent Markov jump process over time which controls the system dynamics by switching between different SDEs. The switching process models dependence on time, but not location in latent space. In contrast, the gpSLDS does not explicitly represent a latent switching process and rather models switching probabilities as part of the kernel function. The dependence of the kernel on the location in latent space allows for the gpSLDS to partition the latent space into different linear regimes. While our work has followed the inference approach of Archambeau et al. [61] and Duncker et al. [10] , other methods for latent path inference in nonlinear SDEs have been proposed [32, 41, 62, 63] . Verma et al. [41] parameterized the posterior SDE path using an exponential family-based description. The resulting inference algorithm showed improved convergence of the E-step compared to Archambeau et al. [39] . Course and Nair [32, 62] proposed an amortization strategy that allows the variational update of the latent state to be parallelized over sequence length. In principle, any of these approaches could be applied to inference in the gpSLDS and would be an interesting direction for future work. 
discussion: 
 In this paper, we introduced the gpSLDS to infer low-dimensional latent dynamical systems from high-dimensional, noisy observations. By developing a novel kernel for GP-SDE models that defines distributions over smooth locally-linear functions, we were able to relate GP-SDEs to rSLDS models and address key limitations of the rSLDS. Using both simulated and real neural datasets, we demonstrated that the gpSLDS can accurately infer true generative parameters and performs favorably in comparison to rSLDS models and GP-SDEs with other kernel choices. Moreover, our real data examples illustrate the variety of potential practical uses of this method. On calcium imaging traces recorded during aggressive behavior, the gpSLDS reveals dynamics consistent with the hypothesis of a line attractor put forward based on previous rSLDS analyses [27] . On a decision making dataset, the gpSLDS finds latent trajectories and dynamics that clearly separate by motion coherence and choice, providing a dynamical systems view consistent with prior studies [42, 48] .In our experiments, we demonstrated the ability of the gpSLDS to recover ground-truth dynamical systems and key dynamical features using fixed settings of hyperparameters: the latent dimensionality K and the number of regimes J. For simulated data, we set hyperparameters to their true values; for real data, we chose hyperparameters based on prior studies and did not further optimize these values. However, for most real neural datasets, we do not know the true underlying dimensionality or optimal number of regimes. To tune these hyperparameters, we can resort to standard techniques for model comparison in the neural latent variable modeling literature. Two common evaluation metrics are forward simulation accuracy [23, 27] and co-smoothing performance [6, 64, 65] .While these results are promising, we acknowledge a few limitations of the gpSLDS. First, the memory cost scales exponentially with the size of the latent dimension due to using quadrature methods to approximate expectations of the SSL kernel, which are not available in closed form. This computational limitation renders it difficult to fit the gpSLDS with many (i.e. greater than 3) latent dimensions. One potential direction for future work would be to instead use Monte Carlo methods to approximate kernel expectations for models with larger latent dimensionality. In addition, while both the gpSLDS and rSLDS require choosing a discretization timestep for solving dynamical systems, in practice we find that the gpSLDS requires smaller steps for stable model inference. This allows the gpSLDS to more accurately approximate dynamics with continuous-time likelihoods, at the cost of allocating more time bins during inference. Finally, we acknowledge that traditional variational inference approaches -such as those employed by the gpSLDS-tend to underestimate posterior variance due to the KL-divergence-based objective [38] . Carefully assessing biases introduced by our variational approximation to the posterior would be an important topic for future work.Overall, the gpSLDS provides a general modeling approach for discovering latent dynamics of noisy measurements in an intepretable and fully probabilistic manner. We expect that our extension will provide a useful addition to the rSLDS and related methods on future analyses of neural data.Here, we draw a mathematical connection between functions sampled from a GP with a linear kernel,EQUATION)and functions of the formEQUATION)This connection is useful for understanding the relationship between our model formulation and the typical linear dynamical systems formulation as in the rSLDS. In particular, we will show that the linear kernel in eq. ( 11) equivalently places a prior on β and β 0 in eq. ( 12), in a similar manner to Bayesian linear regression.By definition of a GP, for any finite set of input locations {x i } N i=1 , we haveEQUATIONwhere Φ ij = κ lin (x i , x j ). Equivalently, for every pair i, j = 1, . . . , N ,EQUATION)Under eq. ( 12), treating β and β 0 as random, this would becomeEQUATIONEquation ( 15) and eq. ( 16) are equivalent and eq. ( 13) is satified if and only ifβ β 0 ∼ N 0, M -M c -c T M c T c + σ 2 0 (17)This shows that a GP with a linear kernel is equivalent to a Bayesian linear regression model with a prior on coefficients of the form in eq. ( 17).We present the full inference and learning details for the gpSLDS with additive inputs,EQUATION)Our approach primarily follows that of Duncker et al. [10] . We use the following notation from Duncker et al. [10] throughout this section:• ⟨h(•)⟩ p(•) denotes the expectation of h(•) with respect to the distribution p(•).• K xz ∈ R N1×N2 denotes the covariance matrix defined by our kernel κ ssl (•, •) for two batches of points {x i } N1 i=1 and{z i } N2 i=1 . Specifically, [K xz ] ij = κ ssl (x i , z j ).If one of these batches has only one point, e.g. N 1 = 1, we denote the corresponding covariance vector as k xz ∈ R 1×N2 . Note, these quantities depend on Θ through the kernel; we omit writing this dependence for brevity.Following Duncker et al. [10] , we first augment our generative model with sparse inducing points{u k } K k=1 ⊂ R M at input locations {z m } M m=1 ⊂ R D .Inducing points can be seen as pseudoobservations of f (•) at input locations {z m } M m=1 . They allow for tractable inference of f at any new batch of N input locations by reducing computational complexity from O(N 3 ) to O(N M 2 ), where typically M is chosen to be much smaller than N [37] . After this augmentation, the joint likelihood of our model becomesEQUATION)Treating u k as pseudo-observations, we assume the following augmented prior:EQUATIONThen we can viewEQUATIONwhereµ f k |u k (x) = k xz K -1 zz u k κ f k |u k (x, x ′ ) = κ ssl (x, x ′ ) -k xz K -1 zz k zx .As in Duncker et al. [10] , we consider a variational approximation to the posterior of the formEQUATION)Using this factorization, we derive the ELBO to the marginal log-likelihood of our model. By Jensen's inequality,log p(y | Θ) = log p(y | x)p(x | f )p(f | u, Θ)p(u | Θ)dxdf du ≥ q(x, f , u) log p(y | x)p(x | f )p(f | u, Θ)p(u | Θ) q(x, f , u) dxdf du = q(x, f , u) log p(y | x)p(x | f ) K k=1 p(u k | Θ) q(x) K k=1 q(u k ) dxdf du = ⟨log p(y | x)⟩ q(x) -⟨KL[q(x)||p(x | f )]⟩ q(f ) - K k=1 KL[q(u k )||p(u k | Θ)] := L(q(x), q(u), Θ),whereEQUATION)To perform inference over the posterior of latent paths q(x), we follow a method first proposed in Archambeau et al. [39] and extended by Duncker et al. [10] .As in Archambeau et al. [39] , we choose a posterior distribution q(x) characterized by a Markov Gaussian process,EQUATION)This distribution satisfies posterior marginals q(x t ) = N (x t | m t , S t ) which satisfy the differential equationsEQUATIONArchambeau et al. [39] maximize the ELBO with respect to q(x) subject to the constraints in eq. ( 25) and eq. ( 26) using the method of Lagrange multipliers. They show that after applying integration by parts, the Lagrangian becomesEQUATIONwhereEQUATIONAs in Archambeau et al. [39] , we assume that λ(T ) = 0 and Ψ(T ) = 0.To find the stationary points of the Lagrangian, we first take derivatives of L with respect to m(0), S(0), m(t), S(t), A(t), b(t) and set them to 0. The derivatives with respect to m(0) and S(0) lead to the updatesEQUATIONwhere we assume a prior on initial conditions p(x 0 ) = N (x 0 | µ(0), V (0)).The derivatives with respect to m(t) and S(t) lead to the stationary equationsEQUATION)EQUATIONwith P ij = 1 2 for i ̸ = j and 1 otherwise. The inclusion of P was proposed by Duncker et al. [10] to adjust for taking derivatives with respect to a symmetric matrix.To take derivatives with respect to A(t) and b(t), we first extend a result from Appendix A of Archambeau et al. [39] to our affine inputs model. This allows us to rewrite the KL-divergence term between the posterior and prior latent paths in the ELBO asEQUATION+ T 0 v(t) T B T ⟨∆f ⟩ q(x),q(f ) dt + 1 2 T 0 v(t) T B T Bv(t)dtwhere ∆f := ff q . The integrand of the first term in eq. ( 34) can be expanded asEQUATION+ Tr A(t) T A(t)(S(t) + m(t)m(t) T + 2m(t) T A(t) T ⟨f ⟩ q(x),q(f ) + b(t) T b(t) -2b(t) T ⟨f ⟩ q(x),q(f ) -2b T A(t)m(t)where we have used the identity ⟨f (x)⟩ q(f ) (x -m) T q(x)= ∂f ∂x q(x),q(f ) S, which can be derived from Stein's lemma. Note that computing eq. ( 35) relies on three quantities, ⟨f ⟩ q(x),q(f ) , f T f q(x),q(f ) , ∂f ∂x q(x),q(f ) which can be written as terms which depend on expectations of the kernel with respect to q(x). We derive these as follows:⟨f ⟩ q(x),q(f ) = k xz K -1 zz u q(u),q(x) = ⟨k xz ⟩ q(x) K -1 zz m u f T f q(x),q(f ) = K k=1 f k (x) 2 q(f k ) q(x) = K k=1 Var q(f k ) [f k (x)] + ⟨f k (x)⟩ 2 q(f k ) q(x) = K k=1 Var p(f k |u) [f k (x)] q(x),q(u) Term 1 + Var q(u) [⟨f k (x)⟩ p(f k | u) ] q(x) Term 2 + ⟨f k (x)⟩ 2 p(f k | u),q(u) q(x) Term 3 = K k=1 k xx -k xz K -1 zz k zx q(x) Term 1 + k xz K -1 zz S k u K -1 zz k zx q(x) Term 2 + k xz K -1 zz m k u (m k u ) T K -1 zz k zx q(x) Term 3 ∂f ∂x q(x),q(f ) = m T u K -1 zz ∂kzx dx q(x) = m T u K -1 zz ∂kzx dx q(x)The above three function expectations can thus be expressed in terms of the following four kernel expectations with respect to q(x) = N (x | m, S): ⟨κ(x, x)⟩ q(x) , ⟨κ(x, z)⟩ q(x) , ⟨κ(z 1 , x)κ(x, z 2 )⟩ q(x) , ∂κ(z, x) ∂x q(x).For our SSL kernel these kernel expectations are not available in closed form, so in practice we approximate them using Gauss-Hermite quadrature.Next, differentiating L with respect to A(t) and b(t) yields the updatesEQUATIONEQUATIONNote that these updates have one key difference from the previously derived updates in Archambeau et al. [39] and Duncker et al. [10] : the input-dependent term Bv(t) in eq. (37) . Intuitively, this is because the posterior bias term b(t) is fully time-varying, so it captures changes in the latent states due to input-driven effects in the posterior.In summary, the inference algorithm for updating q(x) is as follows. In each iteration of vEM, we repeat the following forward-backward style algorithm:1. Solve for m(t), S(t) forward in time starting from m(0), S(0) via eq. ( 25) and eq. ( 26). 2. Solve for λ(t), Ψ(t) backward in time starting from λ(T ), Ψ(T ) = 0 via eq. ( 32) and eq. ( 33). 3. Update A(t) and b(t) via eq. ( 36) and eq. ( 37).After solving these stationary equations, we update m(0) and S(0) via eq. ( 31).Computational details Solving for m(t), S(t), λ(t), and Ψ(t) requires integrating continuoustime ODEs. In practice we use Euler integration with a small discretization step ∆t relative to the sampling rate of the data, though in principle any ODE solver can be used. We found that the ELBO usually converges within 20 forward-backward iterations.For the experiments in this paper, we considered two observation models: Gaussian observations and Poisson process observations.We consider the observation modelEQUATIONwhere R ∈ R D is a diagonal covariance matrix. The expected log-likelihood is available in closed form and is given byEQUATIONClosed-form updates for C, d and R are also available:EQUATIONEQUATIONEQUATIONwhere n ti is the number of observed time points, R * d is the d-th entry of R, and c d is the d-th row of C.The second observation model we consider is Poisson process observations of the formEQUATIONwhere either g(a) = exp(a) (exponential inverse link) or g(a) = log(1 + exp(a)) (softplus inverse link). For the exponential inverse link, the expected log-likelihood is available in closed form and is given byEQUATIONFor the softplus inverse link, the expected log-likelihood is not available in closed-form, but can be approximated by Gauss-Hermite quadrature or a second-order Taylor expansion around m(t).For both Poisson process models, we update C and d using gradient ascent on the expected loglikelihood with the Adam optimizer.Here we derive a closed-form update for B, which linearly maps external inputs to the latent space. The only term in the ELBO which depends on B is ⟨KL[q(x)||p(x | f )]⟩ q(f ) . We differentiate this term as written in eq. ( 34) and arrive at the updateEQUATIONNote that the termT 0 v(t)v(t) T dt -1can be pre-computed since v(t) is known.𝑥 ! 𝑥 " 𝑥 "err($ 𝑤, 𝑤!"#$ )Figure 5 : Comparison between the standard vEM approach in Duncker et al. [10] and our modified vEM approach. A. Estimation error between the true and learned decision boundaries, computed as described in Appendix C. For each vEM approach, we fit 5 gpSLDS models with different random initializations. The estimation errors are denoted in light blue/purple dots. The runs that we display in the next two panels are denoted by a solid blue/purple dot. B. The standard vEM approach fails to learn the true decision boundary of x 1 = 0. C. The modified vEM approach precisely learns this decision boundary.In this section, we empirically compare the standard vEM approach from Duncker et al. [10] to our modified vEM approach in which we learn kernel hyperparameters on a partially optimized ELBO.For this experiment, we use the same synthetic dataset from our main result in Section 4.1. We fit the gpSLDS with 5 different random initializations using both standard vEM and modified vEM. For these fits, we fix the values of C and d throughout learning to ensure that the resulting models are anchored to the same latent subspace (in general, they are not guaranteed to end up in the same subspace due to rotational unidentifiability). Each run was fit with 50 total vEM iterations; each iteration consisted of 15 forward-backward solves to update q(x) and 300 Adam gradient steps with a learning rate of 0.01 to update kernel hyperparameters.To compare the quality of the learned hyperparameters, we quantitatively assess the error between the learned and true decision boundaries. In this simple example with J = 2, the decision boundary can be parametrized as w 0 + w 1 x 1 + w 2 x 2 = 0 for some w = (w 0 , w 1 , w 2 ) T . The true decision boundary is characterized by w true = (0, 1, 0) T . We denote the learned decision boundary as ŵ. Next, we compute an error metric between the learned and true decision boundaries as follows. We first normalize the learned decision boundary and define ŵnorm = ŵ ∥ ŵ∥2 . We do not need to do this for w true since it is already normalized. Then, we use the error metric err( ŵ,EQUATION)Including both terms in the minimum is necessary due to unidentifiability of the signs of ŵ.Figure 5A compares this error metric across the 5 model fits for each vEM method. It is clear that the modified vEM approach consistently outperforms the standard vEM approach in terms of more accurately estimating the decision boundary. In addition, the error for standard vEM has much higher variance, since the algorithm is prone to getting stuck in local maxima of the ELBO. Figures 5B-C display the learned versus true decision boundaries in the latent space for two selected runs. We select the run from each vEM approach which achieved the lowest decision boundary error metric, as denoted by solid dots in fig. 5A . We find that the model fit with standard vEM learns a decision boundary which noticeably deviates from the true boundary. On the other hand, the model with with modified vEM recovers the true boundary almost perfectly. This illustrates that our modified approach dramatically improves kernel hyperparameter estimation in practice and enables the gpSLDS to be much more interpretable in the latent space. To further demonstrate the expressivity of the gpSLDS over the rSLDS, we apply the gpSLDS to a synthetic dataset where the true decision boundary between linear regimes is nonlinear. The rSLDS can only model linear decision boundaries in order for its inference algorithm to remain tractable.For this example, we generate a synthetic dataset consisting of an unstable linear system and a stable linear system separated by the decision boundary x 2 1 + x 2 2 = 4. Both of the linear systems have fixed points at (0, 0). The smooth combination of these linear systems results in a 2D limit cycle (fig. 6A ). We simulate 30 trials of Poisson process observations from D = 50 neurons over T = 2 seconds (fig. 6B ). We initialize the observation model parameters C and d using a Poisson LDS with data binned at 20ms. Then, we fit a gpSLDS with J = 2 regimes, and with π(x) modeled using the feature transformationEQUATION)The results of this experiment are shown in fig. 6C-E . In fig. 6C we find that the gpSLDS successfully recovers the the true latent trajectory with accurate posterior credible intervals for an example trial. Furthermore, fig. 6D demonstrates that by using the quadratic feature transformation in eq. ( 54), the gpSLDS accurately learns the true flow field and true decision boundary x 2 1 + x 2 2 = 4. In addition, the values of π(x) smoothly transition between 0 and 1 near this boundary, highlighting the ability of our method to learn smooth dynamics if present. Lastly, in fig. 6E we plot the inferred posterior variance of our method. We find that the gpSLDS is more confident in regions of the latent space with more data (e.g. at the decision boundary), and less confident in regions of latent space with little to no data.We fit all of our models on a NVIDIA A100 GPU on an internal computing cluster. A breakdown of approximate compute times for the main experiments in this paper includes:• Synthetic data results in Section 4.1: 1.5 hours per model fit, ∼40 hours for the entire experiment. • Real data results in Section 4.2: 1.5 hours per model fit, ∼8 hours for the entire experiment.• Real data results in Section 4.3: 1 hour per model fit, ∼5 hours for the entire experiment.We note that these estimates do not include the full set of experiments we performed while carrying out this project (such as preliminary or failed experiments).Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: We provide all details for the experimental setting that are needed to understand the results for all three sets of experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: Yes, we report thorough quantitative comparisons with error bars for the experiments in Section 4.1 and 4.2. Throughout the paper, we fit multiple (i.e. at least 5) models with different initializations for our experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: By taking a product of a linear kernel and a "partition" kernel, gpSLDS elegantly improves Duncker & Sahani's work on estimating locally linear latent dynamical systems. The paper is a key advance in the topic, and it is clearly written with detailed experiments.

Reviews: review: paper_summary: paper_summary: This paper introduces the Gaussian Process Switching Linear Dynamical System (gpSLDS), a novel approach for modeling latent neural dynamics. This model extends the Gaussian process stochastic differential equations framework by incorporating a new kernel function that supports smoothly interpolated locally linear dynamics. This innovation allows the gpSLDS to maintain the expressiveness needed for complex systems while enhancing interpretability, addressing the limitations of the rSLDS. The paper's contributions include the development of the gpSLDS model, introduction of a novel kernel function to balance expressiveness and interpretability, and a new learning algorithm that improves the accuracy of kernel hyperparameter estimation. The model's effectiveness is shown through applications to both synthetic and real neuroscience data, showing superior performance compared to existing methods like the rSLDS. This advancement provides a more robust framework for understanding neural computation and dynamics.
strength_weakness: strength_weakness: 1. The proposed model effectively transforms a linear dynamical system into a Gaussian Process with a well-designed kernel. While this adaptation supports smoothly interpolated locally linear dynamics, it does compromise computational efficiency, shifting from linear time to cubic time cost.
questions: questions: N/A
limitations: limitations: The authors should explore and discuss the limitations of the gpSLDS, including specific scenarios or conditions under which the model may underperform or fail.
paper_summary: paper_summary: The paper explores latent state inference and parameter learning within a switching stochastic dynamical system. In this context, the dynamics are represented by a stochastic differential equation, with the drift function modeled as a Gaussian process. Notably, the paper introduces a novel kernel for this Gaussian process—a mixture of linear functions that captures the system’s switching behavior. The inference process employs a variational inference framework. To validate its effectiveness, the paper conducts evaluations using both synthetic and real neuroscience data.
strength_weakness: strength_weakness: My primary concern with this work lies in how the latent dynamics, denoted as $f$, are modeled as independent draws from the same Gaussian process (GP) prior. In my view, this implies that the posterior over latent components should be independent in each component. Consequently, learning becomes challenging when only single components of the system are observed. This limitation could pose a substantial problem, especially in classical tracking scenarios where only specific state components are observable.

Additionally, I perceive the contribution as somewhat marginal. Essentially, the paper introduces a new GP kernel, but beyond that, the impact seems limited.

Choice of Linear Kernel:

I wonder about the rationale behind using a linear kernel. Why not explore switching between nonlinear kernels? 

Inference Algorithm Analysis: 

Lastly, a detailed examination of the variational inference algorithm would provide valuable insights. Additionally, a synthetic example for comparison could help demonstrate the algorithm’s performance.

For instance, the authors could select a linear stochastic differential equation (SDE) and perform closed-form parameter estimation using methods like the Expectation Maximization (EM) algorithm. In this scenario, latent state inference could be achieved through RTS smoothing.
questions: questions: - Could the authors provide further details on how the presented approach could incorporate a multi-dimensional correlated Gaussian process prior?
limitations: limitations: -
paper_summary: paper_summary: This paper proposes a new model called Gaussian Process Switching Linear Dynamical System (gpSLDS). The model is more interpretable and infers more stable latent compared with the alternative rSLDS. Particularly, the weird oscillations of the latent can be avoided due to the newly proposed Smoothly Switching Linear (SSL) kernel. Extensive experimental results on one synthetic and two real-world datasets shows the benefit of gpSLDS compared with the one with the traditional RBF kernel and the rSLDS.
strength_weakness: strength_weakness: * Some of the presentations can be improved.
* Lacks some complexity analysis or results. This might be an important factor to be considered for long sequences.

See questions.
questions: questions: * Typo: Line 129, function's slope of
* Figure 1A, if different colors are different samples (should be clearly state if so), then it is meaningless to show these samples unless kernel parameters are told.
* Line 177, $\boldsymbol z_m \in \mathbb{R}^D$? And it seems like this $\boldsymbol z$ is different from the $z_j$.
* What does the conjugacy of the model mean in line 193?
* How about the effect of choosing $K$ and $J$. At least some explanations or guidance and the corresponding effects are necessary for them.
* How about the algorithm complexity. Computing GP is time consuming, which might hinder the computing efficiency for long sequences.
limitations: limitations: /
paper_summary: paper_summary: This paper introduces a new approach to model the low-dimensional latent dynamics of a collection of neurons over time. Their approach balances the two desiderate of capturing complex nonlinear dynamics while remaining interpretable. Specifically, they introduce the Gaussian Process Switching Linear Dynamical System, which models the latent state as a Gaussian process with a specific novel kernel function, which interpolates between different linear dynamics akin to a switching linear dynamical system. This mitigates some of the pathologies of rSLDS and also provides uncertainty estimates for the dynamics. The authors validate their method both on synthetic and experimental data from neural recordings.
strength_weakness: strength_weakness: The biggest weakness the paper currently has in my opinion is that it is hard to judge what the computational trade-offs are between the different methods. While the authors provide estimates for the overall compute they use in all experiments (see Appendix E), there is no comparison between different methods. I would have liked to see an analysis that compares GP-RBF, rSLDS and gpSLDS not only as a function of the number of trials or number of steps simulated forward, but also in terms of their predictive performance as a function of training time / budget. Otherwise, it is hard to judge whether these methods were compared fairly.

One aspect that would make the paper stronger, but is not strictly necessary in my opinion, is to add another baseline, which is different from the two approaches that were combined in gpSLDS, e.g. something like LFADS (Sussillo et al, 2016).

- Sussillo, D., Jozefowicz, R., Abbott, L., & Pandarinath, C. (2016). LFADS - Latent Factor Analysis via Dynamical Systems. In Advances in Neural Information Processing Systems (NeurIPS)
questions: questions: - Figure 3: GP with RBF kernel seems to perform better with respect to forward simulation, why do you think that is?
- Could you comment on how expensive the different approaches are relative to each other in terms of training and prediction?
limitations: limitations: While the paper discusses limitations of their method throughout the paper, this is a bit sparse. I think the paper could be improved by a dedicated limitations paragraph in the discussion in Section 6.

Comments: comment: I have read the rebutal and thank the authors for their answers. I do not have any further questions and will keep my score.
Thank you for the response. I have no further questions and will keep my score.
Thanks for the response. I don't have further questions and I would like to keep my score.
Thank you for answering my questions on computational cost and adding an explicit limitations section. I believe the changes in the rebuttal strengthen the paper and have raised my confidence in my earlier assessment that this paper should be accepted. I've increased my confidence score to reflect this.