Identifying Neural Dynamics Using Interventional State Space Models
Authors: Amin Nejatbakhsh
Yixin Wang
Venue: ICML 2025 poster (2025)
Keywords: Causal dynamical systems
interventions
state space models
photo-stimulation
micro-stimulation
Primary area: applications->neuroscience_cognitive_science
TL;DR: We develop interventional state space models, a class of causal models that can predict neural responses under causal perturbations and identify neural dynamics.

Abstract: Neural circuits produce signals that are complex and nonlinear. To facilitate the understanding of neural dynamics, a popular approach is to fit state space models (SSM) to the data and analyze the dynamics of the low-dimensional latent variables. Despite the power of SSM to explain the dynamics of neural circuits, these models have been shown to merely capture statistical associations in the data and cannot be causally interpreted. Therefore, an important research problem is to build models that can predict neural dynamics under causal manipulations. Here, we propose interventional state-space models (iSSM), a class of causal models that can predict neural responses to novel perturbations. We draw on recent advances in causal dynamical systems and present theoretical results for the identifiability of iSSM. In simulations of the motor cortex, we show that iSSM can recover the true latents and the underlying dynamics. In addition, we illustrate two applications of iSSM in biological datasets. First, we applied iSSM to a dataset of calcium recordings from ALM neurons in mice during photostimulation. Second, we applied iSSM to a dataset of electrophysiological recordings from macaque dlPFC during micro-stimulation. In both cases, we show that iSSM outperforms SSM and results in identifiable parameters. The code is available at https://github.com/amin-nejat/issm.

Identifying Neural Dynamics Using Interventional State Space Models
Amin Nejatbakhsh Yixin Wang? 2
Abstract
dynamics when certain edges are functionally removed from the circuit. This powerful strategy enables testing sophisti- cated neural hypotheses. For example, O'Shea et al. (2022) used perturbations to understand whether dynamics in the motor cortex are path-following (driven by an upstream brain region), low-dimensional, or high-dimensional. An- other example by Feulner al. (2022) uses a similar strategy to investigate whether feedback drives plasticity for rapid learning in the motor cortex. Another study by Sanzeni et al. (2023) uses optogenetic perturbations uncover the degree of coupling in the visual cortex of mice and monkeys. They demonstrate through modeling that under strong network coupling, perturbations lead to a reshuffling of responses within the circuit.
Neural circuits produce signals that are complex and nonlinear. To facilitate the understanding of neural dynamics, a popular approach is to fit state space models (SSM) to the data and analyze the dynamics of the low-dimensional latent vari- ables. Despite the power of SSM to explain the dynamies of neural circuits, these models have been shown to merely capture statistical associ- ations in the data and cannot be causally inter- preted. Therefore, an important research problem is to build models that can predict neural dynam- ics under causal manipulations. Here, we propose interventional state-space models (iSSM), a class of causal models that can predict neural responses to novel perturbations. We draw on recent ad- vances in causal dynamical systems and present theoretical results for the identifiability of iSSM. In simulations of the motor cortex, we show that iSSM can recover the true latents and the underly- ing dynamics. In addition, we illustrate two appli- cations of iSSM in biological datasets. First, we applied iSSM to dataset of calcium recordings from ALM neurons in mice during photostimula- tion. Second, we applied iSSM to a dataset of elec- trophysiological recordings from macaque dIPFC during micro-stimulltion. In both cases, we show that iSSM outperforms SSM and results in identifi- able parameters. The code is available at https: / //github cwm/ ami ne jat/ issm.
The main insight of these works is that in the absence of perturbations (i.e. observational regime), neural dynam- ies are confined to low-dimensional spaces, and models that are built upon observational data are not able to cap- ture neural dynamics outside of the low-dimensional space. However, during perturbations (i.e. interventional regime), the neural state is driven outside of the task space, providing more information about dynamics in the global neural state space (Jazayeri & Afraz, 2017). This insight allows to build sophisticated hypotheses that can only be tested using perturbations (Fig. 1). Interventional studies are critical for determining the causal contribution of neural dynamics to behavior and perception. For example, a study by Shahbazi et (2022) used targeted electrical stimulation to manipu- late a monkey's perception.
Furthermore, many of the popular models used in neuro- science suffer from identifiability issues (Maheswaranathan et al., 2019). In one line of work, researchers have developed similarity metrics that are agnostic to non-identifiabiliy transformations (see Sucholutsky et al. (2023) for a review). However, in many cases, the model parameters or latent variables are biologically meaningful, and recovering them is desired. Therefore, the need for developing identifiable models for neuroscience data analysis is an overarching goal. For example, Zhou & Wei (2020) used an identifiable VAE as opposed to a vanilla VAE and inferred latent variables that encode the geometry of the task in an unsupervised manner.
1. Introduction
Understanding neural data requires identifying the dynam- ies that underlie it. The principled way to achieve this is through causal perturbations. When a perturbation is deliv- ered, the activity of perturbed neurons is dissociated from their upstream neurons, facilitating inspection of the circuit
"Center for Computation Neuroscience, Flatiron Institute, New York, US University Michigan, Michigan, US. Correspondence to: Amin Nejatbakhsh <anejatbakhsh @ ffatironinsiiute.org..
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Identifying Neural Dynamics Using Interventional State Space Models
Neural Recordings
Oiservaiina Craphical Model
Oboevvaiocal Dijectory
(A.B)
Normal State
x
Latret Activty Space
340
[A.8]
Inaervetional Gruphical Model
Neasi State Space
Time
Intevvationd Trajenndee
Stimulation
(0)
X1
Stindattinn Diliver
y
you
Neunl Stote
Space
Figure 1. Overview. (A) Schematic of neural recordings in observational (top) and interventional (bottom) regimes. (B) Neural dynamics in the state space. The observational data (top) is confined to a tow-dimensional task space, whereas the interventional data (bottom) explores the state space, enabling the testing of causal neural hypotheses. (C) Graphical model of SSM (top) compared with iSSM (bottom).
2. Related Work
Model identification in static data The emerging field of causal representation learning provides statistical treat- ments for recovering the true parameters of statistical mod- els (Schölkopf et al., 2021). Most of the developments correspond to static models and can be broadly catego- interven-
State space models To understand neural circuits, a popu- lar strategy is to build low-dimensional state space models (SSM). Driven by the neural manifold hypothesis, we often assume that neural data lies on a low-dimensional mani- fold. The challenge then becomes discovering the latent manifold and characterizing how the dynamics evolve in the low-dimensional space. Subsequently, suite of SSMs has been developed covering wide range of assumptions and applications. A typical SSM follows a dynamic model and an emission model described by the following equations:
rized into identification using observational tional data. (1) Observational: While early theoretical guarantees have been limited to linear mixing and asym- metric noise (Comon, 1994), these results have been ex- tended to nonlinear mixing (Locatello et al., 2019; Xi & Reddy, 2023), and nonlinear mixing with observa- tion noise (e.g.. VAEs) (Khemakhem et al., 2020a), and multi-environment data (Lachapelle al., 2023). (2) Inter- ventional: With access to interventional data, identifiability results can be extended to broader classes of models. Lippe et al. (2022) shows that with sparse interventions, we can recover latents up to permutation, scaling, and offset. Ahuja et al. (2023) utilize independent support properties (Wang Jordan, 2021) and provide identifiability guarantees. These results have been further extended to nonparametric latents with linear and nonlinear mixing functions (von Küigelgen et al., 2023; Buchholz et al., 2023; Varici et al., 2023).
Et+1 (x())))) Et: == fo(x)) + di, E, p(E:), 8, 1~pple p(ô.).
With this general formulation, models depart based on the specification of ge, fe: p(5,). Linear dynami- cal systems (LDS) assume that both go. fo are linear and p(e:), are multivariate normal distributions. A sep- arate line of work assumes that g is switching linear and develops algorithms that jointly infer switching times as well as latent states (Petreska et al., 2011; Linderman et al., 2017; Fox et al., 2008). These models have been successful, particular when there are abrupt changes in the dynamics.
Model identification in dynamic data More recently, the- oretical results on statistical model identification have been extended Markov models and switching linear dynamical systems Malsells-Rodaa et al., 2023). These results provide the identification of the model parameters up to a class of nuisance transformations (c.g., affine). Most relevant to our work are Yao et al. (2022: 2021). The main shortcoming of these works is that they do not incorporate noise in the observation space, which is crucial for modeling biological
Although SSMs have been primarily used for fitting obser- vational data, there have been a few attempts to apply them to interventional data as well. However, responses to per- turbations are often modeled as additive, which makes the SSM models non-causal. We will elaborate on this further in section 3.1.
2

---

Identifying Neural Dynamics Using Interventional State Space Models
datasets. Previous work can be broadly categorized into two groups. Some studies consider the transient interven- tional effects, while others investigate the persistent effects in the stationary regime (Schölkopf & von Kügelgen, 2022:
1931). Therefore, the model we consider here theoretically has the capacity to fit complex datasets.
3.1. Interventional State Space Models
Malinsky& & Spirtes, 2018; Besserve & Schölkopf, 2022; Benkō et al., 2018; Malinsky& & Spirtes, 2018; Peters et al., 2022). Hansen & Sokol (2014) uses differential equations as structural equations in dynamical systems. Abuja al. (2021) considers (deterministic) linear dynamics (referred to as mechanism) and nonlinear emissions (referred to as rendering) and proves that the latent space of such model is identifiable from observational data up to mechanism in- variances. Lippe et al. (2023) show that for linear dynamics, we have access to binary interventional data, then the la- tents are identifiable up permutation. Yao et al. (2022); Song et al. (2024b) focus on the identification of latent non- stationary dynamics using observational data. Hyvarinen & Morioka (2016; 2017): Hyvarinen et al. (2019); Hälvä et al. (2021) focus on extending nonlinear ICA and its identifia- bility results to temporal settings. They impose constraints on the mixing function and the latent dynamics to achieve identification using only observational data.
Notation Consider an experiment with N recorded neu- rons over T time steps repeated for K trials. We denote neural responses at time t by yt, where y/ is an N-vector that concatenates the spike counts or calcium activities of all neurons. We assume the existence of a time-dependent latent variable ze E PP where D is the dimension of the latent space. We present the interventional model and elabo- rate on its differences from the observational model.
The first modeling assumption that distinguishes iSSM from SSM is that we assume perturbing neurons affect the latent dynamics directly, which will consequently affect neural responses in the next time point according to the emissions model. The second, more critical assumption is that when- ever a latent node is perturbed, its activity is dissociated from all its upstream nodes. This assumption is easy to incorporate in a linear model, which is achieved by ignoring the columns in the dynamics matrix corresponding to the perturbed node. Denoting the interventional input to indi- vidual channels at time t by u, E RM, we model æ, y, as:
In addition to the statistical literature on model identifica- tion, recent work in dynamical systems theory has utilized the Koopman theory to find conditions, such as sampling frequency, for the exact identification of the continuous-time dynamical systems from sampled data (Zeng et al., 2022).
Lt+1 1(Bu, = 0} @ Ax, + Bu, + Et, y1 P(y))fo(z))).
(3)
(4)
3. Methods
where et *N(00, and @ denotes element-wise mul- tiplication. A E RDx captures latent dynamics, while B E DXM captures the effect of neural perturbations on latent dynamics. We place a centered Laplace prior on B with the scale parameter s encourage its sparsity, which is critical for the identifiability of the results as we discuss in section 3.3. Q E R is the covariance of e and fe is a generic nonlinear function mapping latents to observations. If the intervention u, is zero, the model follows observa- tional dynamics, but in the presence of an intervention, the model decouples the intervened node from its parents. We term this model interventional SSM iSSM. In section 3.3 we theoretically characterize the conditions under which iSSM is identifiable.
Here, we develop new class of latent variable models that aim to capture neural dynamics in both observational and interventional regimes. We base our model on the frame- work of Causal Inference (CT) (Pearl et al., 2016). Instead of directly modeling the joint distribution of the data, CI uses structural equations to describe the generative process of the data. In this framework, interventions are modeled as changing the structural equations. Equivalently, when intervention is performed on a node, it is disconnected from all its parents in the generative model, and its distribution is set to a new distribution. A major benefit of modeling the interventions in this way is that having access to interven- tional data allows us to identify the model parameters, as we will see in section 3.3.
Justification of intervening on latents: iments in neuroscience (both optical and electrophysiologi- cal) still do not perform interventions at the single-neuron level. In some cases, the light- gated proteins (c.g., Chan- nelrhodopsin) are expressed in all the neurons of a certain subtype, and broad-fiell illumination is used to activate the population of neurons of that subtype. In other experiments, light-gated proteins are expressed in a sparse subset of neu- rons, and two-photon lasers are used to activate only those neurons. In these cases, due to the spillover effect of lasers, it is nearly impossible to activate individual neurons with
Many causal exper-
Although our framework is general and can be applied to various types of SSMs, in this paper, we focus on adding in- terventional components to the SSMs with linear dynamics and nonlinear observations. These models are natural ex- tensions of linear models with a higher capacity to express complex datasets (Gao et al., 2016). It has been argued theoretically that linear dynamics in a latent space with sufficiently large dimension, followed by a nonlinear emis- sion, is powerful to model any dynamical system (Koopman,
3

---

Identifying Neural Dynamics Using Interventional State Space Models
enough spatial precision. Similarly, for electrophysiological techniques such micro-stimulation,, each electrode targets several neurons located in its vicinity opposed to a single neuron. While iSSM allows for performing interventions on individual neurons by setting B to the identity matrix (see Supp. Fig. 6 for an example of this), it is still useful to think of latents as providing a natural grouping of the neurons into behaviorally relevant subsets. We demonstrate this in our results in section 4.
tions; (3) further identify So and P((a))r)) up to permu- tation, coordinate-wiss shifting and scaling.
Begin with the first step of identifying ((ffo(z)))ee) from Pef We make the following assumptions on the observation model.
Assumption 3.1 (Bounded completeness of P(Welz;)). The function P(ye|zs) (rnn)-wheree 2; = ( is bounded com- plete in y:- Specifically, a function s(X,Y) ) is bounded complete in Y if ʃ g(X) y) = 0 implies g(X) 0 almost surely for any measurable function g(X) bounded in Li-metric (Yang et al., 2017).
3.2. Inference
Since our model involves a nonlinear emission as well as a non-conjugate noise model, we resort variational infer- ence techniques. Our goal is to infer the posterior distribu- tion Peo ( (1.T. M..T. (I:T) while optimizing the parameters O. We follow the methodology of reparameterization and amortized inference, but adapt some parts to our specific interventional scheme. For a review of variational methods for state space models see Archer et al. (2015). Denoting the approximate posterior distribution by qx (2)r), the ELBO loss function is presented below:
When the observational model satisfies the bounded com- pleteness assumption, we can identify P( (((fo(z:))ter) from P((t)heT) (We detail the proof in Appendix C.) Many common functions P(y x1) satisfy the bounded com- pleteness condition, including exponential families (Newey & Powell, 2003), location-seal scale families (Hu & Shiu, 2018), and nonparametric regression models (Darolles al., 2011). It is common assumption to guarantee the existence and the uniqueness of solutions to integral equations, most com- monly used in nonparametric causal identification in proxy variables and instrumental variables (Miao et al., 2018; Yang et al., 2017; D' D'Haultfocuille, 2011). I). We refer the readers to Chen al. (2014) for a detailed discussion of completeness.
C [( θ) = En- log pe (Yhrr, unrr, x1rr)
- log q (x1: Wirr BBar, up:r)
We next proceed to identifying So and p(( up to affine transformations. We require the following assumption on the mixing function fo-
where is reparameterized as =(n) = μæ + .. The functions Jr, a are typically parameterized by neural net- works (called the recognition network) with an architecture that matches the dataset domain. Here, we choose an LSTM for the recognition network.
ssumption 3.2 (Mixing function). The mixing function fe is piecewise linear, continuous, and injective.
While the piecewise linear assumption may appear restric- tive, we note that it entails flexible choices of sa(-), in- cluding (deep) ReLU networks that can approximate com- plicated functions. We finally leverage the interventional data to achieve coordinate-wise identification of so and ((IxheeT) We make the following assumptions on the latent dynamics.
Another important addition that makes the inference in our model possible to apply the interventional structure di- reetly to the approximate posterior during training. To do this, we replace be with BBBu, = 0} & 44 + Bu, dur- ing each iteration of optimization. This ensures that the interventional data indeed manipulates the causal graph con- sistently in the approximate posterior. We refer the reader to the supplementary C.1 for more details on inference.
Assumption 3.3 (Faithfulness). When the system is unper- turbed, there does not exist non-zero vectors V, W such that V TT 1 W xt, t.
3.3. Theoretical results: On the identifiability of iSSM
Loosely, this assumption guarantees that the latent variables x, are all non-trivially connected. In other words, there does not exist linear transformations such that the two neighbor- ing timestep V x++1 and WT x, will become independent. This assumption guarantees that any independence relation- ships observe in the dynamical system is due to the in- tervention; it will facilitate the alignment of latent variables. We further describe the requirements of the interventions that need to be performed for identifying iSSM. Under these assumptions, we can achieve the identification of iSSM as follows.
We provide sufficient conditions for the identifiability of iSSMs. We show that, given a sufficient set of do- interventions, one can identify both the latent dynamics matrix A and the mixing function f-(- of the iSSM. This identifiability of the latent dynamics enables us to extrapo- late to novel, unseen interventions.
To identify the latent dynamics of iSSM, we proceed in three steps: (1) identify P(.fo(z))her) from the ob- served data distribution P(fycher) (2) identify fe and P((=))er) from (((f(z;))ee up to affine transforma-
4

---

Identifying Neural Dynamics Using Interventional State Space Models
A1
Tree 01
B1
Interne Obs
C1
Tatent Flow Field
m1
Signols
Reconsmrction Acouney
True 100mm
Tiee Infered
Lo-
likmed Sirr
Latens
Dynumt Alroon
G.5
1.0
0.5
-10
SSM
158M
-15
1.0
Recatteral Dysumiis
0.5
0.0
2,00
1.0
0.5
0.0
-15
-13 40
-15
C3
Court
Figure 2. Models of Motor Dynamics. (A1,A2) Flow field underlying dynamic attractor (A1) and rotational dynamics (A2) models of the motor cortex and sampled interventional trajectories. Blue dots represent stimulation times. (BI.B2) Reconstruction of the data using iSSM. (CI,C2) True (black) and inferred (red) linear latent flow fields. The true fields are generated from the evolution of w; in equations 5 and 6. The inferred flow fields are generated according to the fitted latent linear dynamies A. (D1,D2) True (black) and inferred (red) dynamics of the latents æt (top panels) and observations ye (bottom panels) using iSSM. Blue shades in each plot correspond to stimulation times. (E1,E2) Comparison between SSM (observational model) and iSSM (ours). Reconstruction correlation between true and inferred latents (top pancls) and observations (bottom pancls) with an increasing number of interventions is shown. With more interventions, iSSM can better identify the latents.
Theorem 3.4 (Block identifiability of iSSM and generaliza- tion to unseen interventions). Under Assumptions 3. to 3.3, the latent dynamics A and the mixing function of fo(:) can be block identified up to permutation, and coordinate-wiss shifting and scaling; the block corresponds to a separa- tion between intervened and unintervened latents. That is, Å = AAlI + c, where A is an invertible block diago- nal matrix (with two blocks of size 1 1i = 0} and [ * 0}), II is a permutation matrix, and c. is a constant vector. The observations' distribution P(yy:hez) under a novel unseen une" intervention is also identifiable when the intervention only acts on latents in separate blocks, i.e. {i for any i,j E {ko Bu b + 0}.
Below we extend Theorem 3.4 to multiple interventions, for which we can achieve full identifiability if the interventions are sufficiently diverse.
Assumption 3.5 (Unordered pairs condition (Hyttinen et al., 2013)). A set of interventions satisfies the unordered pair condition on the latents, if for any unordered pair 2t,5, Zt.j- there exists an intervention such that Tt,; is intervened on but Ft, is not, or æt.j is intervened on but æt,y-
Corollary 3.6 (Identifiability of iSSM under sufficiently diverse interventions). If the interventions satisfy Assump- tion 3.5, then the iSSM is identifiable up to permutation, along with coordinate-wise scaling and shifting, Á AX* + c for some diagonal A", The observations' distri- bution under any novel interventions is also identifiabl.
The proof of Theorem 3.4 is in Appendix C. This result establishes the block identifiability of iSSM and its predic- tive power for unseen interventions. It says, given single intervention trial, one can separate out the intervened la- tents from the unintervened ones through iSSM. It will also enable us to extrapolate to novel, unseen interventions as long as the novel interventions only touch upon latents that are already separated. The key idea of the proof relies on independence testing among the latents, built on the obser- vation that intervention would make the intervened latent independent of the values of all other latents in previous time steps. This result illustrates how interventions can help identify latent variables by inducing statistical independence among the latents, revealing latent dynamics in non-linest state-space models.
Corollary 3.6 crucially relies on Assumption 3.5 to achieve separation among all pairs of latents, thanks to Theorem 3.4, hence achieving full identifiability and extrapolation to all unseen interventions.
Remark 3.7 (unknown sparse B). We note that the identifi- ability results are primarily driven by the intervention status of Bu, where the matrix B unknown. However, a sparse B would ensure that each dimension of the interventional input u, only affects few latent dimensions, which has two important consequences. First, an interventional input u, that is nonzero only in a few dimensions would typically alter only a small subset of latents, enabling identification of intervened latents from unintervened ones based on the
5

---

Identifying Neural Dynamics Using Interventional State Space Models
Theorem 3.4. Second, designing diverse interventions is usually easier under sparse or known B. Therefore, As- sumption 3.5 implicitly presupposes the presence of some priori structure in B.
A1
Dynami Antreoor
Dynemics
Recnsst uction Accuncc
Roturiomal Dyramizs
Laumecc
Latrib
0
80
166
100
4. Results
0.5
4.1. Identifying motor cortical dynamics in simulations
16-
Inferred
00:
T'me
Pynamic Attoctor
To illustrate how iSSM leads to identification, we take in- spiration from models of the motor cortex. A key obser- vation in the motor cortex is the presence of rotational tra- jectories (Churchland et al., 2012). From computational perspective, it has been argued that rotational trajectories provide a basis for motor neuron activations and muscle movements. has also been argued that the rotational ba- sis provides robustness to noise and interventions (Logiaco et al., 2021). Inspired by these results, multiple dynamical models have been proposed for rotational activities in the motor cortex (Laje& & Buonomano, 2013; Sussillo et al., 2015). The first model, called Rotational Dynamics (RD) proposes that the motor cortex has underlying rotational dy- namics. As a result, in this model, the rotational dynamics are generated within the motor cortex independent of input or feedback activity (Fig. 21; Sussillo al. (2015)), Eq. 5 describes the dynamics and emissions of RD.
Latzth
SSM
1.0
0.5-
00--
0
80
160
0
0
80 160
Time
Obs Dimenaakss
Figure 3. High-dimensional Models Motor Dynamics. To mimic high-dimensional neural responses and low-dimensional latent dynamics, generated latents from equations 5, 6, and noisy observations from equation 7. (A1,A22 Latents (æt)) and observa- tions (g:) for example trial in 5 dimensions for the DA model. (B (B2) Same as panel A for the RD model. accuracy of the latents (lop) and observations (bottom) for the RD model using SSM and iSSM models as a function of observation dimension. Results suggest that iSSM not only reconstructs the observations but also recovers the true latents. (C2) Same as C1 for the DA model.
(C1) Reconstruction
Rotational Dynamics (RD):
dx 0 aæ1
x cos (os())) 22 sin (22)
+ Es, y/ =
+ 8;
(5)
di
where as before æ, y denote the true latents and observations, and ô represent true latent noise and observation noise. We set the generative noise values to zero in all experiments.
an interventional strategy to discover whether the dynamics in the motor cortex follow either of these regimes. Similarly, here we ask if interventional data can distinguish between these models. We generate data from RD and DA to address this. The latent states xxt) in both models follow linear dy- namics, while the observation model in both cases is highly nonlinear. Therefore, recovering the true latents is not a trivial task.
The second model, called Dynamic Attractor (DA), as- sumes that the underlying dynamics of the motor cortex a rounded attractor. In this model, the rotational dynamics in motor neurons are generated by some upstream region moving the state along the attractor (Laje & Buonomano, 2013). Eq. 6 describes the dynamics and emissions of DA.
During data generation, we apply repeated interventions interleaved by resting periods for the network to return its stationary state. The dynamics of latents and observations are shown in Fig. 2A,D. Although in the absence of inter- ventions both models produce the same trajectories, one can observe that the interventional trajectories exhibit distinct characteristics (Fig. 2A).
Dynamic Attractor (DA): dx:
dt
02 x22)
01
+ Ce. y:
(1 - a1) cos( (22)
(1 x2) sin x22)
+ 5
(6)
While these models have distinct characteristics and pro- pose different underlying circuit mechanisms, Galgali al. (2023) show that the trial averages of these models can be precisely the same, limiting our ability to identify the true dynamics of the motor cortex solely from observational data.
Consistent with O'Shea al. (2022), our results suggest that in the presence of interventional data using the iSSM model, one can identify the underlying dynamics and emissions (Fig. 2B,C) and recover the true latent variables (Fig. 2D). This recovery continues to improve as we collect more inter- ventional data, emphasizing the importance of perturbation experiments in causal hypothesis testing (Fig. 2E). The re- producibility details are included in the Appendix B.
O' Shea et al. (2022) refer to these models as low- dimensional vs. fath-following dynamical systems and use
6

---

Identifying Neural Dynamics Using Interventional State Space Models
A1
Observations
A2
Smoothed
B
Latent Activities
Trai Rec. Acc.
Train Rec. Avc.
0.56
Curst Latert aocncee
D2
Test Rec. Acc.
120
0.35
0.48
Latent Index
0.30
0.40
3SM
ISSM
Neurons
Lotest Dixe.
0.40
Consistency
60
30-
0:35
Laton
10
24
0:30
1:25
0.00
1.81
Time
3.61
5:42
0.00
1.81
3.61 Time (6)
5422
Calcium Activity
0.0
Time (5)
2.5
5.0
0.1 Spority
Figure 4. Results on Mice Dataset. (A1) Calcium responses of ALM neurons in an example trial during stimulation (specified by the shaded blue band and the white arrow). (A2) Smoothed responses given by the mean posterior of the model. (B) Mean latents discovered by iSSM for correct (green) vs. incorrect (red) trials. The dynamics of the latents distinguish between correct and incorrect trials without any prior knowledge of the behavior. (C1.C2) Training (CI) and test (C2) reconstruction accuracy of iSSM for range of hyperparameters (varying number of latents and the sparsity of B). Both higher sparsity and a larger number of latents improve the accuracy. (D1.D2) Comparison between SSM (blue) and iSSM (red) for the optimal selection of the hyperparameters for training (DI) and test (D2) data. (E) Consistency score of matrix across random initializations of the model parameters for SSM and iSSM, showing that iSSM leads to better identifiability (see section 4.5).
4.2. Simulations high-dimensional motor activity
4.3. Identifying dynamics underlying short-term memory in mice
To mimic high-dimensional neural responses in the motor cortex, we projected data from the DA and RD models to higher dimensions. For each run, generated random orthogonal matrix C I DxN and used it in the generation of y: in the following way:
Persistent activity is hallmark of short-term memory across species (Romo et al., 1999; Fuster & Alexander, 1971). How can a network of neurons produce activities in response input stimulus that is maintained after the stimulus is re- moved? Multiple network mechanisms are proposed to un- derlie persistent activity. Among those, one popular model is known as Functionally Feedforward (FF) model (Gold- man, 2009). FF assumes that the network consists few smaller subnetworks that are connected in feedforward manner. Since these subnetworks do not necessarily need to form spatial cluster in the brain, experimentally find- ing footprints of this type of connectivity is not feasible. However, the theoretical properties of the model have been well-studied. For example, results robustness to structural noise (Qian et al., 2024). An alternative model for the persistent activity known as Line Attractor (LA) model (Seung, 1996). Under LA circuit model, the activity of an upstream region pushes the state of the circuit along the line attractor, and the dynamics preserve the state until a new input has arrived (see Appendix A. for results on models working memory).
High-dimensional Simulations:
y: C
cos(ar))
+ 8,, y: =C= [1 11) cos(az)) mz) (in()
5,
sin(#2)
(7)
for RD and DA respectively. In this case, 5, E R is nigh-dimensiooal noise added the responses. We set the
standard deviation of noise to 0.1 for all high-dimensional experiments.
commonly argued that FF
An example of the latents 2; and 5-dimensional responses y: shown in Fig. 3A.B. We generated data for varying number of observation dimensions and applied SSM and iSSM models (Fig. 3C). While both SSM and iSSM achieve high reconstruction accuracy, only iSSM is able to identify the true latents (Fig. 3C1, (2). The error bars represent the standard deviation of the reconstruction accuracy across 3 seeds of randomly generated trials for each observation dimension. Other generative and inference parameters are the same as the previous section (see Appendix B).
We applied iSSM public dataset of targeted photostim- ulation in the anterior lateral motor cortex (ALM) of mice

---

Identifying Neural Dynamics Using Interventional State Space Models
during short-term memory task (Daie al., 2021). The task included a sample epoch where auditory cue guided the mice for a left vs. right cue to get water reward. The sample epoch was followed by a delay epoch of 3 seconds, where the mice needed to engage working memory to keep track of the guided cue. Finally, during the response period, the mice received the reward if the lick direction was correct. The photostimulation was delivered during the delay period for short amount of time, started simultaneously with the delay period, or after 1 or 2 seconds.
livering micro-cireuit electrical stimulations (Nejatbakhsh et al., 2023). We analyzed 6 datasets, 3 with only observa- tional data and 3 with a combination of observational and interventional data.
In Fig. 5A, we show firing rates recorded from each of the 96 electrodes for an interventional session. In each interven- tional session, two electrodes were repeatedly stimulated while recordings were performed from all other electrodes. We apply iSSM to this dataset and infer the latents and pa- rameters. Fig. 5A2 shows a sample from the inferred model closely matching the data. The reconstruction accuracy on the training and testing sessions is larger for iSSM com- pared to baseline SSM across a range of latent dimensions (Fig. 5E), suggesting that the model can better generalize to test trials. To find the optimal hyperparameters, we per- formed vross-validation The results are shown in Fig. 5D.
Calcium recordings were done in 179 identified neurons for 77 repeated trials (Fig. 4A). There were 8 photostimulation channels targeted to stimulate neurons according to their response selectivity.
We set the dimension of interventional inputs u, the number of photostimulation channels and fitted the stim- ulation matrix B with a sparsity penalty. The smoothed and denoised neural activities are shown in Fig. 4A2. The reconstruction accuracy of the data for both training and testing trials for iSSM was larger than the baseline SSM model across a range of latent dimensions (Fig. 4D). Further- more, the latents learned by the model show distinct mean trajectories for correct incorrect trials, suggesting that they capture behaviorally meaningful dynamics (Fig. 4B). The optimal latent dimension and sparsity penalty for this dataset were found based on cross-validation (Fig. 4C).
4.5. Performance Metries
Here, we briefly describe the metries used throughout the experiments in the paper.
Reconstruction accuracy: The reconstruction accuracy is defined as the correlation coefficient between the true and reconstructed signals (latents or observations). For latents, it can only be computed for synthetic experiments where ground truth exists. This metric is commonly used in the causal representation learning literature for assessing the identifiability of the latents (Khemakhem et al., 2020a;b; Song al., 2024a).
4.4. Generalizing to test interventions in primates
Understanding network dynamics to control behavior has been a longstanding challenge in neuroscience. The over- arching goal is to deliver targeted stimulation to a network of neurons to steer the dynamics or the behavior towards pre-determined outcome (Haimerl et al., 2023; Jou al., 2023). A first step towards understanding the circuit effects or behavioral influences of network manipulations is to build models that can predict the response to interventions. The space of possible interventions is combinatorial and intractable to cover. Therefore, an alternative approach is to build models that can generalize to unseen interventions.
Consistency score: A hallmark of identification is robust- ness to initialization. To test whether the iSSM results in identifiability, we ran the model several times with different random initializations. We compute the consistency score across random initialization by first aligning the columns of B across runs to account for permutation invariance of the latents, followed by computing the Euclidean distance between aligned matrices. The aligned distances are con- siderably smaller for iSSM compared to SSM. providing evidence for the identifiability (e.g.. Fig.
We showed theoretically in section 3.3 that iSSM has this property. Concretely, if we fit the iSSM model to inter- ventional data, where the dataset consists of small set of canonical interventions, the model can generalize to un- observed interventions. To validate this empirically, we showed results on a synthetic dataset (Fig, 2). Here we want to test whether these results hold in a real biological dataset.
5. Discussion
Summary Here we proposed iSSM, a framework for joint modeling of observational and interventional data. We pro- vided theoretical results showing that the iSSM model, when fitted on interventional data, leads to the identifiability of latents as well dynamics and emissions. To illustrate iSSM's 's applicability, we showed results on 3 different ex- amples covering range of assumptions. The first example was a synthetic dataset with linear dynamics and nonlinear emissions. The second example was calcium recordings from the mouse ALM region with targeted photostimulation
The dataset consisted of electrophysiological recordings us- ing electrode arrays implanted on the prefrontal cortex of macaque monkeys during quiet wakefulness (resting) while the animals were sitting awake in the dark. The electrode array included 96 electrodes that were also used for de-
8

---

Identifying Neural Dynamics Using Interventional State Space Models
Observations
Smoothed
Activities
D1
E1
Train Rec. Acc.
30
Anhhn Stim
9 705-
ISSM
Stim
Latent Index
Laten Din
80
0.700
0.090
W-- w
0.695
E2
0,690
20
6.455
pla
6.650
Neurons
0.8
Time
1.5
G01
Lune Dixs
Rec. Acc. Dist.
D2
Test Rec. Acc.
Consistency of
40
0.8
12
0.545
Consistency Score
Tarat Dan
10
0,540
10
0.535
Time
1,20
1.65
0.30
Time
1.30
Firing Rate
Tole
0.001 6601 0.1 Syarsry
SSM
ÍSSM
Figure 5. Results Primate Dataser. (A1) Unit firing rates for training interventional session (stimulation times are shown by the blue band and white arrow). (A2) Sample responses for the same trial inferred by iSSM. (B) iSSM's latent dynamics show differential responses during stimulation resting periods. (C) The distribution of training and test reconstruction accuracy across trials. (D1,D2) Training (DI) and test (D2) mean reconstruction accuracy is shown with varying hyperparameters (number of latents and sparsity parameter of B). (EI,E2) Comparison between SSM (blue) and iSSM (red) for increasing number of latents. Reconstruction accuracies shown for training (EI) and lest (E2) trials. Both SSM and iSSM benefit from a larger number of latents, with iSSM consistently outperforming SSM. Consistency of the inferred across random initializations showing better identifiability for iSSM (see section 4.5).
delivered by channels that targeted groups of neurons. The third example was electrophysioloical recordings from the macaque monkey's prefrontal cortex with micro-stimulation delivered by the same recording electrodes. In all cases, our results show impressive generalization capabilities and parameter recovery, confirming our theoretical results.
Third, critical modeling assumption that we made that interventions directly affect the latents. Furthermore, our theoretical results require performing interventions on in- dividual latents through neurons that are causally linked to them. Further experimental validation is required to test whether this assumption is valid in neural data.
Fourth, here we rely on variational inference (VI). While VI is commonly used in state space modeling, is an approx- imate method. Furthermore, the choice of the recognition network (in our case, L.TMs) can amplify the approxima- tion gap. Better inference algorithms and more powerful architectures can further improve the identifiability results.
Limitations We identify several limitations of our work. First, in this work, we focused on a generative model that has linear dynamics. While the inference model can still capture nonlinearities through its recognition network, explicitly modeling nonlinearities and providing theoretical results is an important limitation of this work.
Finally, this work establishes a framework for modeling neural data under interventions and is intended to motivate future investigations into its limitations and extensions.
Second, as depicted in Fig. 1, we think of interventions forcing the neural state out of its low-dii www-dimeniionl manifold. Therefore, it is valid question how can we still assume the low dimensionality of dynamics. Our argument is the fol- lowing: (1) While the observational data is low-dimensional, we think that interventional data is higher-dimensional but still lives in much lower dimension than the number of neurons. This is simply because certain configurations of neural states are not biologically possible. (2) The optimal latent dimension depends on the sparsity of the B as well as the amount of data available. With larger datasets and more interventions, we conjecture that increasing the number of latents always improves the results.
Acknowledgmentss
YW was supported by the Office of Naval Research under grant N00014-23-1- 1-25900, the National Science Foundation under grant No. 2231 174, No. 2310831, No. 2428059, No. 2435696, No. 2440954, and a Michigan Institute for Data Science Propelling Original Data Science (PODS) grant. The authors also acknowledge helpful discussions with Kayvon Daie, Roozbeh Kiani, Alex H Williams, Timo- thy Kim, and constructive feedback from the reviewers.
9

---

Identifying Neural Dynamics Using Interventional State Space Models
Impact Statement
Daie, K., Svoboda, K., and Druckmann, S. Targeted photo- stimulation uncovers circuit motifs supporting short-term memory. Nature neuroscience, 24(2):25--265, 2021.
The modeling presented in this work aims to offer a method- ological framework for advancing our understanding of neu- ral computation under interventions. While the analysis of electrophysiological and calcium imaging data holds poten- tial for long-term impact on the understanding of neurologi- cal disorders across species, such implications lie beyond the scope of this theoretical modeling. We therefore antici- pate no immediate societal consequences from this study.
Darolles, S., Fan, Y., Florens, J.-P., and Renault, E. Non- parametric instrumental regression. Econometrica, 79(5): 1541-1565, 2011.
D'Haultfoeuill, X. On the completeness condition in non- parametric instrumental problems. Econometric Theory, 27(3):460-471, 2011.
References
Feulner, B., Perich, M. G., Miller, L. E., Clopath, C., and
Ahuja, K., Hartford, J., and Bengio, Y. Properties from mechanisms: an equivariance perspective on identifiable representation learning. arXiv preprint arXiv:2110.15990 2021.
Gallego, J. A. Feedback- based motor control can guide
plasticity and drive rapid learning. bioRxiv, pp. 2022-10,, 2022.
Fox, E., Sudderth, E., Jordan, M., and Willsky, A. Nonpara- metric bayesian learning of switching linear dynamical systems. Advances neural information processing sys- tems, 21, 2008
Ahuja, K., Mahajan, D., Wang. Y., and Bengio, Y. Inter- ventional causal representation learning. In International Conference Machine Learning, pp. 37- 407. PMLR, 2023.
Fuster, J. M. and Alexander, G. Neuron activity related short-term memory. Science, 173(3997):622-65,, 1971.
Archer, E., Park, I. M., Buesing, L., Cunningham, J., and Paninski, L. Black box variational inference for state space models. arXiv preprint arXiv:1511.07367, 2015.
Galgali, A. R., Sahani, M., and Mante, V. Residual dynamics resolves recurrent contributions to neural computation. Nature Neuroscience, 26(2):326-3) 2023.
Balsells-Rodas, C., Wang, Y., and Li, Y. On the identifiabil- ity of switching dynamical systems, 2023.
Beiran, M. and Kumar, A. Prediction of neural activity in commectome-constraincee recurrent networks. bioRxiv, pp. 02, 2024.
Gao, Y., Archer, E. W., Paninski, L., and Cunningham, J. P. Linear dynamical neural population models through nonlinear embeddings. NeurIPS, 29, 2016.
Benkō, Z., Zlatniczki, Á. Stippinger, M., Fabó, D., Sólyom, A., Eröss, L., Telcs, A., and Somogyvári, Z. Complete inference of causal relations between dynamical systems. arXiv preprint arXiv:1808.. 10806, 2018.
Goldman, M. S. Memory without feedback in a neural network. Neuron, 61(4):621-634 2009.
Haimerl, C., Ruff, D. A., Cohen, M. R., Savin, C., and Simoncelli, E. P. Targeted v1 comodulation supports task- adaptive sensory decisions. Nature communications, 14 (1):7879, 2023.
Besserve, M. and Schölkopf, B. Learning soft interven- tions in complex equilibrium systems. In Uncertainty in Artificial Intelligence, pp. 170-180. PMLR, 2022.
Hälvä, H., Le Corff, S., Lehéricy, L., So, J., Zhu, Y., Gas- siat, E., and Hyvarinen, A. Disentangling identifiable features from noisy data with structured nonlinear ica. Advances in Neural Information Processing Systems, 34: 1624-1633, 2021.
Buchholz, S., Rajendran, G., Rosenfeld, E., Aragam, B., Schölkopf, B., and Ravikumar, P. Learning linear causal representations from interventions under general nonlin- ear mixing. arXiv preprint arxiv:2306.02235. 2023.
Chen, X., Chernozhukov, V., Lee, S., and Newey, W. K. Local identification of nonparametric and semiparametric models. Econometrica, 82(2):785-809, 2014.
Hansen, N. and Sokol, A. Causal interpretation of stochastic differential equations. 2014.
Churchland, M. M., Cunningham, J. P., Kaufman, M. T.,
Hu, Y. and Shiu, J.-L. Nonparametric identification using instrumental variables: sufficient conditions for complete- ness. Econometric Theory, 34(3):659-693, 2018.
Foster, J. D., Nuyujukian, P., Ryu, S. I., and Shenoy, K. V.
Neural population dynamics during reaching. Nature, 487 (7405):51-56, 2012.
Hyttinen, A., Eberhardt, F., and Hoyer, P. O. Experiment selection for causal discovery. The Journal of Machine Learning Research, 14(1330041-071, 2013.
Comon, P. Independent component analysis, a new concept?
Signal processing, 36(3):287-314, 1994.
10

---

Identifying Neural Dynamics Using Interventional State Space Models
Hyvarinen, A. and Morioka, H. Unsupervised feature ex- traction by time-contrastive learning and nonlinear ica. Advances in neural information processing systems, 29, 2016.
Lippe, P., Magliacane, S., Löwe, S., Asano, Y. M., Co- hen, T., and Gavves, E. Biscuit: Causal representa- tion learning from binary interactions. arXiv preprint 1XXi::3306.0644333 2023.
Hyvarinen, A. and Morioka, H. Nonlinear ica of temporally dependent stationary sources. In Artificial Intelligence and Statistics, pp. 460 469. PMLR, 2017.
Lipshutz, D., Nejatbakhsh, A., and Williams, A. H. Disen- tangling recurrent neural dynamics with stochastic repre- sentational geometry. In ICLR 2024 Workshop on Repre- sentational Alignment.
Hyvarinen, A., Sasaki, H., and Turnrr, R. Nonlinear ica ing auxiliary variables and generalized contrastive learn- ing. In The 22nd International Conference on A Artificial Intelligence and Statistics, pp. 859-868. PMLR, 2019.
Locatello, F., Bauer, S., Lucic, M., Ractsch, G., Gelly, S., Schölkopf, B., and Bachem, O. Challenging common assumptions in the unsupervised learning of disentangled representations. In international conference on machine learning, pp. 41 4114 4114-4124, 4124. PMLR, 2019.
Jazayeri, M. and Afraz, A. Navigating the neural space in search of the neural code. Neuron, 9305):1003-1014, 2017.
Logiaco, L., Abbott, L., and Escola, S. Thalamic control of cortical dynamics in model of flexible motor sequencing. Cell reports, 35(9), 2021.
Jou, C., Hurtado, J. R., Carrillo-Segura, S., Park, E. H., and Fenton, A. A. On the results of causal optogenetic engram manipulations. bioRxiv, pp. 2023-02222022 05, 2023.
Maheswaranathan, N., Williams, A., Golub, M., Ganguli, S., and Sussillo, D. Universality and individuality in neural dynamics across large populations of recurrent networks. Advances in neural information processing systems, 32, 2019.
Khemakhem, L., Kingma, D., Monti, R., and Hyvarinen, A. Variational autoencoders and nonlinear ica: A unifying framework. In International Conference on Artificial Intelligence and Statistics, pp. 2207-2217. PMLR, 2020a.
Malinsky, D. and Spirtes, P. Causal structure learning from multivariate time series in settings with unmeasured con- founding. In Proceedings of 2018 ACM SIGKDD work- shop on causal discovery, pp. 23-47. PMLR, 2018.
Khemakhem, I., Monti, R., Kingma, D., and Hyvarinen, A. Ice-beem: Identifiable conditional energy- deep models based on nonlinear ica. Advances in Neural Infor- mation Processing Systems, 13:12768-12777 2020b.
Miao, W., Geng, Z., and Tchetgen Tchetgen, E. Identify- ing causal effects with proxy variables of an unmeasured confounder. Biometrika, 105(4):987-993, 2018.
Koopman, B. O. Hamiltonian systems and transformation in hilbert space. Proceedings of the National Academy of Sciences, 17(5):315-311 -31811 1931.
Nejatbakhsh, A., Fumarola, F., Esteki, S., Toyoizumi, T., Kiani, R., and Mazzucato, L. Predicting the effect of micro-stimulation on macaque prefrontal activity based on spontaneous circuit dynamics. Physical Review Re- search, 544:0113 2023.
Lachapelle, S., Deleu, T., Mahajan, D., Mitliagkas, I., Ben- gio, Y.. Lacoste-Julien. S., and Bertrand, Q. Synergies between disentanglement and sparsity: generalization and identifiability multi-task learning. In International Con- ference on Machine Learning, pp. 18171-18206. PMLR, 2023.
Newey, W. K. and Powell, J. L. Instrumental variable esti- mation of nonparametric models. Econometrica, 71(5): 1565-1578. 2003.
Laje, R. and Buonomano, D. V. Robust timing and motor patterns by taming chaos in recurrent neural networks. Nature neuroscience, 16(7):925-933, 2013.
O'Shea, D. J., Duncker, L., Goo, W., Sun, X., Vyas, S., Trautmann, E. M., Diester, I., Ramakrishnan, C., Deis- seroth, K., Sahani, M., et al. Direct neural perturbations reveal a dynamical mechanism for robust computation. bioRxiv, pp. 12, 2022.
Linderman, S., Johnson, M., Miller, A., Adams, R., Blei, D., and Paninski, L. Bayesian learning and inference in re- current switching linear dynamical systems. In Artificial intelligence and statistics, pp. 914 922. PMLR. 2017.
Pearl, J., Glymour, M.. and Jewell, N. Causal inference in statistics: A primer. John Wiley & Sons, 2016.
Lippe, P., Magliacane, S., Löwe, S., Asano, Y. M., Cohen, T., and Gavves, E. Causal representation learning for instantaneous and temporal effects in interactive systems. In The Eleventh International Conference on Learning Representations, 2022.
Peters, J., Bauer, S., and Pfister, N. Causal models for dy- namical systems. In Probabilistic and Causal Inference: The Works of Judea Pearl, pp. 671-690. 2022.

---

Identifying Neural Dynamics Using Interventional State Space Models
Petreska, B., Yu, B. M., Cunningham, P., Santhanam, G., Ryu, S., Shenoy, K. V., and Sahani, M. Dynamical segmentation of single trials from population neural data. Advances in neural information processing systems. 24, 2011.
Varici, B., Acarturk, E., Shanmugam, K., Kumar, A., and Tajer, A. Score-based causal representation learning with interventions. arXiv preprint arXiv:2301.082230 2023.
von Kügelgen, J., Besserve, M., Liang, W., Gresele, L., Kekić, A., Bareinboim, E., Blei, D. M., and Schölkopf, B. Nonparametric identifiability of causal represen- tations from unknown interventions. arXiv preprint 23X1::230600052,,, 2023.
Qian, W., Zavatone-Vedt Veth, J. A., Ruben, B. S., and Pehlevan, C. Partial observation can induce mechanistic mismatches in data-cont models neural dynamics. bioRxiv, pp. 2024 05, 2024.
Wang. Y. and Jordan, M. I. Desiderata for representa- tion learning: A causal perspective. arXiv preprint arXiv:2109.03995, 2021.
Romo, R., Brody, C. D., Hernández, A., and Lemus, L. Neuronal correlates of parametric working memory in the prefrontal cortex. Nature, (996673) 1999.
Xi, Q. and Bloem-Reddy, B. Indeterminacy in generative models:
Sanzeni, A., Palmigiano, A., Nguyen, T. H., Luo, J., Nassi, J. J., Reynolds, H., Histed, M. H., Miller, K. D., and Brunel, N. Mechanisms underlying reshuffling of visual responses by optogenetic stimulation in mice and mon- keys. Neuron, 11 111(24):102-41155, 2023.
Characterization and strong identifiability. In International Conference on Artificial Intelligence and Statistics, pp. 0912-6939, PMLR, 2023.
Yang. S., Wang, L., and Ding, Identification and estima- tion of causal effects with confounders subject instru- mental missingness. arXiv preprint arXiv:1702.03955: 2017.
Schölkopf, B. and von Kügelgen, J. From statistical to causal learning. arXiv preprint arXiv:2204.0000 2022.
Schölkopf, B., Locatello, F., Bauer, S., Ke, N. R., Kalch- brenner, N., Goyal, A., and Bengio, Y. Towards causal representation learning, 2021. URL hatpps::/// //axxiv.///xv... org/ org/abs/ abs/
Yao, W., Sun, Y., Ho, A., Sun, C., and Zhang, K. Learning temporally causal latent processes from general temporal data. arXiv preprint arXiv:2110.054
2021.
Yao, W., Chen, G., and Zhang, K. Temporally disentangled representation learning. Advances in Neural Information Processing Systems, 35:26492-2663 2022.
Seung, H. S. How the brain keeps the eyes still. Pro- ceedings of the National Academy of Sciences, 93(23): 13339-13344, 1996.
Zeng, Z., Yue, Z., Mauroy, A., Gonçalves, J., and Yuan, Y. A sampling theorem for exact identification of continuous- time nonlinear dynamical systems. In 2022 IEEE 6lst Conference on Decision and Control (CDC), pp. 6686 6692. IEEE, 2022.
Shahbazi, E., Ma, T., Pernus, M., Scheirer, W. J., and Afraz, A. Perceptography: unveiling visual perceptual halluci- nations induced by optogenetic stimulation of the inferior temporal cortex. bioRxiv, 2022.
Song. X., Li, Z., Chen, G., Zheng, Y., Fan, Y., Dong, X., and Zhang, K. Causal temporal representation learning with nonstationary sparse transition. Advances in Neu- ral Information Processing Systems, 37:77098-77131, 2024a.
Zhou, D. and Wei, X.-X. Learning identifiable and inter- pretable latent models of high- high-dimensional neural activity using pi-vae. Advances Neural Information Processing Systems, 33:7234-7247, 2020.
Song. X., Yao. W., Fan, Y., Dong, X., Chen, G., Niebles, C., Xing, E., and Zhang, K. Temporally disentangled representation learning under unknown nonstationarity. Advances in Neural Information Processing Systems, 36, 2024b.
Sucholutsky, I,, Muttenthaler, L., Weller, A., Peng, A., Bobu, A., Kim, B., Love, B. C., Grant, E.,, Achterberg, J., Tenen- baum, J. B., et al. Getting aligned on representational alignment. arXiv preprint arXiv:2310.130/1 2023.
Sussillo, D., Churchland, M. M., Kaufman, M. T., and Shenoy, K. V. A neural network that finds a naturalistic solution for the production of muscle activity. Nature neuroscience, 1877):0025 2015.
12

---

Identifying Neural Dynamics Using Interventional State Space Models
A. Experimental Results
A.I. Models of Working Memory
Various sources of non-identifiabilit can make challenging to recover the truc latents and dynamics in neural data. We elaborate on two of these sources here. First, neural recordings are undersampled, meaning that from a large pool of neurons involved in the computation only a small fraction are recorded. Undersampling or partial observation is indeed a significant origin of non-identifiabiiity as discussed in the literature (Beiran & Litwin-Kumar, 2024). Second, non-identifiability can be caused by the mixing of input-driven and recurrent activity in the network (Galgali et al., 2023). Lipshutz et al. show that noise correlations can be used to disentangle input-driven and recurrent activity. Here, we take complementary approach and investigate whether interventions can help improve identifiability. A recent theoretical study investigates the effect of partial observations in the context of persistent activity (Qian et al., 2024).
Qian et al. (2024) presented two alternative models of working memory, functionally feedforward (FF) and low-rank connectivity (LR). FF characterized by the existence of sub-circuite in the network that sequentially process their inputs and pass the outputs to the next sub-circuits in feedforward manner (Fig. 6A1). On the other hand, LR models perform their computations recurrently by utilizing line-attractor type dynamics (Fig. 6A2). Importantly, this study shows that when the system is partially observed, observational models (c.g. SSM) cannot distinguish between the two models. Specifically, they show that linear dynamical system (LDS) models have a built-in bias for characterizing the dynamics as FF, regardless of the underlying dynamics (Qian et al., 2024).
Here we ask whether interventional models (e.g. iSSM) can help identify the dynamics. To address this, we generated interventional data from FF and LR models with 5 neurons (shown in Fig. 6B,C). The interventions are performed on individual neurons interleaved by resting duration (blue shades in Fig. 6B). We then used iSSM with N = D 5 and set B to identify to check if we can recover the underlying dynamics. Setting the matrix B to identity amounts to assigning one latent per neuron, an important variation of iSSM that enables modeling interventions on individual neurons. In Fig. 6D-F we show that while both SSM and iSSM successfully reconstruct the observational data (Fig. 6D), only iSSM is able to identify the dynamics matrices (Fig. 6E,F). These results suggest that the ability to perform interventions on individual neurons accompanied by models that leverage the interventional data enable uncovering more details of the causal flow of information in circuit of neurons.
A1
191
Sumple Obs.
Sumple States
D1
Roc. Accuracy
EI
F1
SSM
0.98
Sate
C-
0.90
Interre A
0.94
0.92
A2
I22
D22
E2
F2
0.98
0.0
Low Book
Interred
-0.4
0.91
-0.
Time
60c
Thu
0.00
-0.50
-0.25
Thon A
0.00
0.55
Figure 6. Models of Working Memory: Following (Qian et al., 2024) here we generate data from feedforward (FF) and low-rank (LR) models of working memory to test whether iSSM can recover the true underlying flow field parameterized by matrix A. (A1,A2)
Schematies of circuit mechanisms for FF, LR. (B1, True (black) and inferred (red) interventional signals generated from FF (B1)
and I.R (B2) in 5 dimensions. The inference is done using iSSM with linear dynamics and linear observations. Blue shades represent interventions. (C1,C2) Same data shown in the top 2 PC space with blue dots representing interventions. (D1,D2) Comparison of observational reconstruction accuracy using SSM systems, iSSM estimator has lower variance and is more accurate on average. (EI,E2,F1.F2) Scatter plot of the true values of Aej (x-axis) vs. their inferred values (y-axis) for FF (EI,Fl) and LR (E2,F2) models using SSM (E1,E2) and iSSM (FI,F2) across 3 runs. This result suggests that iSSM is able to leverage interventional data to identify the matrix A.
iSSM. While both SSM and iSSM achieve high accuracy for both FF and LR
13

---

Identifying Neural Dynamics Using Interventional State Space Models
B. Reproducibility Details
B.1. Generative Parameters for Synthetic Experiments
In table B.1 we include the parameters used for generating synthetic datasets for Fig. 2. The stimulation signals u, are parameterized by the following parameters: Stim Amp the overall amplitude of the interventional inputs; Stim Noise: the standard deviation of the zero-mean gaussian noise added the interventional inputs satisfy the assumption3.1 Stim to Rest Ratio the ratio of the amount ot time spend delivering stimulation signals vs. resting; Stim Rep: the number of stimulation repetitions.
We generate data for a duration of T seconds, for K trials, using the time increment dt. The number of observation dimensions and latent dimensions are denoted by N, D respectively. The equations governing the generative process the data given interventional inputs u, are included Eq. with parameters a for the Rotational Dynamics and a1 02 for the Dynamic Attractor.
Table 1.
Generative parameters for synthtic experiments.
Model Noise
Stim Amp
Stim Noise
Stim to Rest Ratio
Stim Rep
Params
Initialization
Rotational
20
10
2
2
0.05
1.2
0
1)
0.5
10
Dynamics
Dynamic
a1
20
20
10
2
2
0.05
0
1)
1
0.5
1
10
Attractor
a2
1.2
B.2. Initialization and Inference Parameters
Table B.2 contains the parameters for initializing our generative process well inference and optimization parameters. The Emission model a fully connected (FC) neural network with = 100 hidden units. For Poisson observation model, we include an additional so softplus us transformation to map the emission outputs to positive values. The Obs. column the table determines what observation noise model was used in each experiment. For the synthetic experiments, we set the B matrix to identity without fitting it, therefore we did not use any Sparsity regularization. For real data experiments we used cross-validation to find the optimal value. For the primate dataset, the training data is the first half of each session and the test data is the second half. For the mice dataset the training data is the lick right trials while the test data is lick left trials. The optimization parameters include the learning rate denoted by Optim LR, number of iterations denoted by Optim (Iter), and the number of hidden units for the LSTM that parameterizes the mean and variance of the variational posterior denoted by LSTM (H). Furthermore, we include the initialization values in table B.2. The initialization parameters include the noise covarinace the dynamies' initial step denoted by Inil (wo noise); the covariance of LDS denoted by Init (LDS a); the matrices A, B; and the covariance scaling of the likelihood model referred to as Init (I.I. a)).
Table 2.
Initialization and inference pavameters for synthetic real data experiments.
Sparsity
Optim
Optim
LSTM
Init
Init (L.DS o)
Init
Init
Emission
Obs.
(s)
(LR)
(Iter)
(H)
(ata noise)
(A, B)
(LII o)
Rotational
(O 1)
FC (H=100)
0.05
Normal
N/A
0.01
1000
10
0.05
0.05
0.05
Dynamics
I
Dynamic
FC (H=100)
0.05
Normal
N/A
0.01
1000
10
0.05
0.05
0.05
Attractor
I
Mice
N10, 1)
Dataset
FC (H=100)
0.05
Normal
X-Val
0.001
1000
10
0.05
0.05
N(0.1 1)
0.05
Primate
FC (H=100)
N00,1
Dataset
softplus
Poisson
X- X-Val
0.01
1000
10
0.05
0.05
A/
N/A
14

---

Identifying Neural Dynamics Using Interventional State Space Models
C. Proof of Theorem 3.4
We consider the interventional state space model (iSSM),
Yt ~ (())))
(8)
2++1 = 1(Bu: : 0) 8 Ax, + Bu, + 6x-
(9)
Step I: Identifying the distribution of _ fo( We begin with identifying the marginal distribution of P(zz) from P(y.). The core assumption we rely on in this step is bounded completeness, which we define in Assumption 3.1
The bounded completeness of (x)) implies that P(z:) is identifiable from P(y)). It is because P(z)) must be the unique solution to the integral equation f' P(g))x))(())dzz = P(y:). Specifically, there are two solutions to this equation (z,z), Pz))) then they must be equal. It is due to the bounded completeness of P(yr|z)): the two solutions must satisfy int P(y ( - P2(z))]zzz= = O, which implies Pr = P2(2)).
Step 2: Affine identification of (a()) and P(( In this step, we establish the affine identification of the mixing function fo(-) by invoking Theorem 3.5 of Balsells-Rodas et (2023): identifying fo(-) from (x))( is a special case of identifying the mixing function in switching dynamical system.
To enable identification, we require Assumption 3.2. In particular, the mixing function should be a piece-wise linear function. Lemma C.1 (Theorem 3.5 of Balsells-Rodas et al. (2023)). Under Assumption 3.2, the mixing function sa(-) and the latent distribution P( (((2t) heT) can be identified from P (()( up to affine transformation.
This lemma is an instantiation of Theorem 3.5 in Balsells-Rodas et al. (2023) in the special case of linear transition dynamics.
Step 3: Identification of s; via interventions. The previous step shows that we can identify ø, up to affine transformation. In this step, we show that, if two solutions of are affine transformations of each other, they must coincide up to block affine transformation, if they agree on the interventional distributions, under Assumption 3.3. This argument implies that the interventional distributions can identify æ; (up to block-wise permutation, and shifting and scaling.)
Concretely, consider two sets of latent variables (at tEeT and (#t where they are affine transformations of each other
i = Ma;- + c, VI.
(10)
Suppose both sets satisfy Equation (23) across all intervention environments, namely,
xt+1 [ = 0} 8 Ax; + Bu: + c.+1 =i 0} Aæ; Bu: + Ec,
(11)
(12)
where both e,, é; are i.i.d over time. Then we will prove that M = AII, where A is an invertible block-diagonal matrix (with al 11(B, = sized block and 1 TT ((Bu, * O)-sized block), and IT is a permutation matrix.
Single intervention case. Below we first consider the simplest case where only one latent variable is intervened, i.e. 1 Bu, = 0}- = 0. We will entend to general interventions later.
We achieve identification using the following observation. Suppose the jth latent mty was intervened in an environment, namely 1((Bu,, = 0} 0. Then we have
Tetij = (Bu,), + a.j for E (t:Bu, + 0}.
(13)
and thus 24+15 D, for all I. The reason is that the intervention set X4+1, to be (But), plus a random noise component, hence independent of all components of xx.
Below argue that, if we also find component j' of 2141 such that f.t1, 1 i1, then My.-s 0. i.e. Retiy must be an affine transformation of 2441.j-
To make this argument, we write
E1415' = MJ.- MM + My MyJI1 +op,
(14)
Lts MJ.-.
MH.-, It.- j Cgy
(15)
15

---

Identifying Neural Dynamics Using Interventional State Space Models
Then since I 1 Bt we have that
Co0($2415;;31) O.
(16)
This implies
0 = Cou( 2x)
(17)
y- (,+1, + My J. My. -,tt + My + -,- ;0t,)))
(18)
- 1141,-5
(M} , ,0 54
+1.) +Cou(My 22 M - + Cov( My 2 1. M1 ((rs)
(19)
-Cou(M, M - () + Cou(M) / - My gF
(20)
The last equation is due to Equation (13). implies that My.-j 0. The reason is that Assumption 3.3 guarantees one of
V, W must be zero, if we write Equation (20) as Cov(V Ft+1, W xx.) = 0, where W = My and is matrix where VJ and V-j- = My-.j. which implies My = Thus we have that the j'th dimension of æ; that achieves the independence property is mapped the jth dimension of æ; up to scaling and shifting; one can separate out the intervened latent from the unintervened ones up to permutation, and coordinate-wise shifting and scaling.
We note that the use of covariance is not essential to the proof. Rather, any statistic that can characterize independence relationships (e.g. moment generating function) apply here and can be leveraged to make this identification argument.
Multiple intervention case. We next extend to the case when multiple latents j E I were intervened, which gives Tt-1, 1 2, in the true latent space.
We find a latent space (i.e. 2.) for the iSSM model that is compatible with the interventional data while maximizing the cardinality of I'' that satisfies a similar independence relationship as the true intervened latents 1t+1, J. ât. In other words, among all candidate ât that can fit the interventional data, we find the one that has the most number of dimensions independent of it, The faithfulness assumption (Assumption 3.3) implies that one cannot use linear transformations to induce additional independence relatinoship. hence [I'] = [Z].
We then leverage the linear identifiability from step 2 that gives
( = Mz, Mz + MI + Cz.
(21)
This step mimics Equation (14) except we use an additional permutation matrix Ii to align the I' entries in é with the I entries in x. This alignment is for notational convenience.
Then following the same argument as in Equation (17), Cov( (iL04+1 .t) = 0 implies that Mr = 0. That is, we can successfully separate out the intervened latents from the un-intervened latents: the intervened latents were mapped to I' in the estimated latent space up to permutation, and block-wise shifting and scaling.
As consequence of identifying the parameters of the iSSM up to blockwise affine transformation, we can predict the observation distributions for novel unseen interventions as long as the interventions touch upon only those latent nodes that are in different blocks; these latents are already separated out: {AII)s = 0 for any i,j E ( Bu,ew) }& * 0}. In other words, the submatrix restricted to the intervened latents in the new intervention {ATT) I x must be diagonal, ensuring that interventions performed on these latents are properly modeled and predicted.
Sufficiently diverse interventions. We finally extend Theorem 3.4 to set of sufficiently diverse interventions. The unordered pairs condition (Assumption 3.5) ensures that Myj 0 for any pair of i,j latents, following the same argument as in Equation (17), since we can separate out any pair of latents as long they are not always either being intervened at the same time or being un-intervened at the same time. Hence, one can identify the latents up to permutation, and coordinate-wist shifting and scaling, i.c. Á = II c for some diagonal A'. It further implies that one can identify the latent dynamics matrix A also up to permutation, and coordinate-wisg shifting and scaling. It also enables us to predict for all unseen interventions given the full identifiability of the dynamics,
Yt P(y. fo(::)),
(22)
x1+1 = 1(B = 0} @ Ax. + Bu: + Et.
(23)
16

---

Identifying Neural Dynamics Using Interventional State Space Models
C.1. Inference Details
Our statistical model is described by the following probabilistic decomposition.
(a MI., u1r)) = p(B) pe II pe (x (p)p(()))))))
where Q [, R. Muu Qo) contains all generative parameters of the model and Mox Qo are the mean and covariance of the initial condition distribution pe(coo) We may include an independent Laplace prior over the elements of B to encourage its sparsity.
Notice that some of these parameters depend on the specific instantiation of the model. For example, for normal observations R E N is the observation covariance whereas the observation mean given by the nonlinear transformation (o For the model with Poisson observations there are no observational parameters since the rate of the Poisson distribution generated by applying a positive transform (e.g. softplus) the output of the nonlinear transform fo()).. In addition, for some instantiations of the model we set B to predefined matrix and do not optimize it. For example, when want to explicitly model the causal interactions between neurons, we set D = N and set B to the identity matrix.
Given a dataset of K trials of length denoted by { (A) (k) K
yiirr " 1:T k=1
our goal is to approximate the posterior distribution
pe(#) pe M1:7) while finding the optimal generative parameters that fit the data the best denoted by Ô Since marginal likelihood optimization and analytical posterior inference in this model is not tractable, we follow variational inference (VI), where the goal is to find a distribution qa (di.Tly1: uI:T) that best approximates the posterior distribution while optimizing the parameters @ via maximizing the evidence lower bound (ELBO).
[( θ = Egg [log (a M1.T, u1r)) log qu (x1.7 W1:T, u1:T)
Using the reparameterizatiion trick, achieve empirical estimate of the gradient of ELBO using samples from an independent noise distribution.
(k) Fir ha (n; Vi.r uir)
(k)
nut (0,,1)
(k)
V Vo.e(( = E. V Po pe (h UI.T, u.r)), VIrr, MIr) log (ho (n: Firr, u1r)]]
In practice, we use Monte Carlo estimates of the expectation above by sampling from an independent noise distribution and evaluating the gradient inside the expectation for those samples. The number of samples used for estimating the gradient introduces trade-off between the accuracy of the estimation computation time. We empirically find that single sample is often enough to achieve accurate estimates. We choose Adam optimizer to perform our stochastic gradient- based optimization.
Crucially, the function ho which takes in the data at each trial and combines with independent noise to generate the latent trajectories an LSTM. We choose LSTM to match the data domain (i.e. time series) and ensure that the variational distribution respects the time causality while being expressive enough to fit our datasets. Specifically, we have the following decomposition of the variational distribution.
T N (uo (Yi:t, ura), o (Vi:t. u1:))
g (XIT)
(k) h4 (m): y VI.T uITT) μ* (3 14 + 17 x T (v) u (.))
(K)
(k)
(k)
where both μlg ah are LSTM functions and the inference is done in an amortized fashion. An important distinction in our variational inference scheme the application of interventional inputs directly in the variational mean f1. == a (y1.r, u. This is achieved via replacing the variational mean at time point t denoted by μr with its interventional counterpart 1 iBBu 0} 3 pt + Bu
17

---

Identifying Neural Dynamics Using Interventional State Space Models
C.2. Computational Complexity
In order to examine the computational complexity of our variational framework, we first expand the gradient calculation.
V C (
By V e log (he (n; VI.T, U1rr), VI.r, u1:r) log qo (ho (1): yI.T, uur)))
V v pog h1 (î: (k) (.r) Vi.ro u (.7)) V log qo (h (7); (k) (1:T))
(1)
where fj is sample from i..d gaussian noise. The estimator above involves (1) generating h4 by running two LSTM forward calls and running backprop to compute the gradients wrt Φ; (2) passing the generated latent trajectory to the generative model with linear dynamics and nonlinear emissions and computing gradients of the log probability w.rt the parameters O Both these scale linearly with time T, and quadratically with the latent dimension D and N. addition, running backprop on the LSTM functions as well as the emission function depends on the size of the architecture used to realize those functions.
18


Reviews:
[{"review_id": "SOvZ8chK8f", "reviewer": "Reviewer_3PQB", "summary": "This paper provides an extension of the state space model (SSM) to an interventional SSM (iSSM). iSSM is able to causally identify and infer external inputs as an intervention to the neural dynamics, while also inferring the latent and reconstruct the observation accurately. Methods, assumptions, and derivations are clearly presented. These good properties are supported by one synthetic and three real-world experiments.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate you and the other reviewers for your time and thoughtful evaluation of our work. We found the feedback to be highly constructive, as well as both fair and encouraging. Below we address your specific questions:\n\n**Definition of reconstruction accuracy:** We noticed that the definition of the reconstruction accuracy is missing in the paper. The reconstruction accuracy is defined as the correlation coefficient between the true and reconstructed latents, which can only be computed for synthetic experiments with ground truth latents. It’s a common metric used in the causal representation learning literature for assessing the identifiability of the latents [1,2,3].\n\n[1] Khemakhem, I., Kingma, D., Monti, R., & Hyvarinen, A. (2020, June). Variational autoencoders and nonlinear ica: A unifying framework. In International conference on artificial intelligence and statistics (pp. 2207-2217). PMLR.\n\n[2] Khemakhem, I., Monti, R., Kingma, D., & Hyvarinen, A. (2020). Ice-beem: Identifiable conditional energy-based deep models based on nonlinear ica. Advances in Neural Information Processing Systems, 33, 12768-12778.\n\n[3] Song, X., Li, Z., Chen, G., Zheng, Y., Fan, Y., Dong, X., & Zhang, K. (2024). Causal temporal representation learning with nonstationary sparse transition. Advances in Neural Information Processing Systems, 37, 77098-77131.\n\n\n**Other variational distributions:** This is a great suggestion, we thank the reviewer for pointing this out. LSTMs have been used in previous SSM literature [1]. We tried vanilla RNNs and LSTMs and LSTMs were significantly better than RNNs. The performance enabled by LSTMs is already sufficient for identifiability experiments we show in the paper. Unfortunately due to the limited time we will not be able to provide new results using new variational families but if we benefit from new variational families we will incorporate the new results in the final manuscript upon acceptance.\n\n[1] Krishnan, Rahul G., Uri Shalit, and David Sontag. \"Deep Kalman Filters.\" arXiv preprint arXiv:1511.05121 (2015).\n\n**Compare against other SSMs:** In all experiments where we compare SSM with iSSM, we use the exact same model, i.e. linear dynamics and nonlinear emissions. The only difference between SSM and iSSM is the interventional formulation. In addition, the focus of our experiments are identifiability as opposed to mere reconstruction of the data. That said, one can develop new iSSMs with different types of dynamics and emissions using the same framework, such as nonlinear dynamics and nonlinear emissions. Notice that the our inference method will still be effective in other scenarios, and only the generative model should account for the new components. This work attempts to lay the foundation of applying interventional models to neural datasets. While our theoretical results currently do not support nonlinear dynamics we’ll consider developing theoretical results for more complex models in the future work. \n\nThanks again for your time; we look forward to your final assessment.\n\n&nbsp;\n\n### **Further responses to reviewer xGnX**\n\n**Question 1**: Our main intuition is that the perturbations kick the state of the system outside of the \"observational manifold\" and force the system to explore regions in the state space that are not visited in the observational regime (we depict this in the schematic in Fig. 1). In addition, our theoretical results suggest that if the true generative model of the data follows our modeling assumptions, using interventional models allows for identifying all model components which theoretically leads to out of distribution generalization. However, in reality the true generative model of the data is far more complex and our model components only provide an approximation to their true counterparts.\n\n**Question 2:**  iSSM is equivalent to SSM in the absence of interventional data. If $\\boldsymbol{u}_t = 0$ then the term that involves $\\boldsymbol{u}_t$ is ignored and we recover the observational SSM model. It’s only in the presence of interventional inputs where the latent nodes dissociate from their parents in the graphical model that enable all the benefits that we exploited (enabling the model to visit new states, changing the graphical model of the data providing more information about the model, etc.).\n\n**Question 3:** Yes. Using LSTMs as the variational family has become standard in the literature of SSMs. Some benefits of LSTMs include (1) it’s causal, therefore it respects the arrow of time (2) it’s context dependent, therefore it integrates information from far away time points to inform the model predictions (3) it’s data efficient and scales well with dimensionality and time. Due to these benefits and specifically since we observed sufficiently good results we did not experiment with alternatives. The only alternative we tried was vanilla RNNs which drastically underperformed LSTMs.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "Have checked in detail before Theorem 3.4.", "experimental_designs_or_analyses": "Yes.", "relation_to_broader_scientific_literature": "Related to a series of SSMs, which is a class of commonly used neural data and neural dynamics analysis tools in computational neuroscience.", "essential_references_not_discussed": "/", "other_strengths_and_weaknesses": "* The presentation is great.\n* The experiment is comprehensive and scientifically interesting.", "other_comments_or_suggestions": "/", "questions_for_authors": "* What is the definition of reconstruction **accuracy**. People usually use test likelihood, rather than accuracy, to evaluate observation data, at least in a lot of machine learning applications and especially the neuroscience field.\n* Is the reconstruction accuracy based on the optimal latent $\\boldsymbol x_t$ inferred by the variational distribution?\n* Have the authors think more about the choice of variational distribution and its validity? LSTM is only unidirectional information flow, which might not be powerful enough to serve as the variational distribution. See, [this](https://proceedings.neurips.cc/paper/2007/hash/2bcab9d935d219641434683dd9d18a03-Abstract.html) and [this](https://openreview.net/pdf?id=2FKzbEE24s).\n* Since the SSM model used in this paper is essentially a nonlinear LDS, there are a lot of LDS as mentioned at the beginning of this paper, such as SLDS. Is it possible to compare iSSM with SLDS? Comparing with the simple LDS seems a bit weak.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Xc7JGWeqTT", "reviewer": "Reviewer_xGnX", "summary": "The authors propose iSSM, which is a linear dynamical systems model that accounts for causal perturbations to the neural population activity. The authors apply this model to two synthetic datasets inspired by literature on motor cortex dynamics, and apply their model to a variety of real neural population datasets that include perturbations.", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate you and the other reviewers for your time and thoughtful evaluation of our work. We found the feedback to be highly constructive, as well as both fair and encouraging. Below we address your specific questions:\n\n**Generalizing to other SSMs:** Thank you for bringing this up. Indeed we present our theoretical results for linear dynamics and nonlinear observations. The algorithmic part of the paper is readily applicable to more complex generative models. Since we’re using variational inference which is agnostic to the underlying model, changing the components of the generative model–i.e. dynamics, emissions, and noise model–will not impact the inference. For identifiability results, the key idea of the proof lies in independence testing: The intervened latents at time t+1 are independent of the others at previous time step t, and this independence leads to the identifiability of the latents. The linearity assumption in the latent space does make this argument easy to make, since we only need to consider the covariance. That said, this argument is readily generalizable to other settings where $x_{t+1}$ and $x_{t}$ are nonlinearly related (e.g. polynomial function or ReLU function up to sparsity constraints), as long as the statistical independence of the independent latents is still sufficient to identify the latents. That said, our goal for this paper is laying the foundation of applying interventional models to neural datasets and we leave these developments for future work. \n\n**Hallmark of short-term memory:** We thank the reviewer for noticing this. While we agree that multiple mechanisms are proposed for working memory, here we’re using \"persistent activity\" loosely without any implications of the underlying mechanism. Notice that even in Goldman 2009 the proposed mechanism still explains the observed persistent activity. We will clarify in the revised manuscript that the observed persistent activity should not be confused with the underlying mechanism generating it.\n\n**Reproducibility:** Thank you for pointing this out. While we have not detailed the specifics of our inference scheme, we’re using standard variational inference tools. That said, we added a new section to the supplementary about the details of inference and computational complexity for completeness. Additionally, we realize that we’ve omitted some important reproducibility details, therefore we added a new section to the supplementary on this. Table [1](https://f.uguu.se/qDwlATSB.png), [2](https://d.uguu.se/jCNcSZzj.png) summarize these results.\n\n**Recovering true latents:** Thank you for this great observation. The arrows presented are all from the ground truth model, we have not shown the inferred dynamics by the model. We will include this in the revised manuscript. That said, notice that it’s not quite straightforward to visualize the inferred dynamics. The dynamics inferred by the model are encoded in the variational LSTM and the fitted emission. There are two ways for visualizing the inferred dynamics. (1) Using the inferred A: This is not ideal because of the variational gap. (2) Using the emission & variational posterior. This approach requires us to generate a grid in the observation space. However, due to the use of LSTMs the variational posterior is context dependent. We’re open to the reviewer’s suggestions for visualizing the dynamics, but so far we haven’t found an appropriate visualization.\n\n\n**Definition of iSSM:** The SSM and iSSM in all experiments share the exact same characteristics–i.e. same latent dimension, emission & inference architecture, step size, etc.–except for the interventional components where SSM models interventions as an additive term. We did this to focus on the utility of interventional models while accounting for all other confounding factors and bringing the two models to the same footing. In the codes, there is an `interventional: bool` flag which if set to `False` will model the interventions in an additive way.  For the Poisson experiment (Fig. 4) the model is equivalent to Gao et. al. 2016, except that they did not include any inputs to the model. For the Normal experiment (Fig. 2, 3) the model departs from Gao et. al. 2016 in terms of the observational noise model.\n\n**Low-d synthetic experiments:** The synthetic experiment in Fig. 2 is mainly presented to provide intuition. It’s an example where latents are clearly defined from a behavioral standpoint and the nonlinear emission is required to produce rotational trajectories. We think it’s an intuitive example for the computational neuroscience audience because it’s a toy version of two prevalent hypotheses about rotational dynamics in the motor cortex. The synthetic example in Supp. Fig. 7 uses higher dimensional latents.\n\nWe used the remaining space in the responses to reviewer **3PQB** and **kBw8** to address your remaining concerns. Thanks again for your time; we look forward to your final assessment.", "claims_and_evidence": "The authors present theoretical results for the identifiability of iSSM. The authors claim that the “framework is general and can be applied to different types of SSMs, in this paper we focus on adding interventional components to the SSMs with linear dynamics and nonlinear observations.” However, I don’t think there was enough information in the paper to conclude to what extent the authors’ results apply. For example, can we replace the dynamics matrix A with a neural network and would the results still hold? How about if both matrix A and B are replaced with a neural network such that we have a completely nonlinear dynamics x_{t+1} = g(x_t, u_t), where g is an MLP?\n\nI think the authors shouldn’t say persistent activity is a hallmark of short-term memory (Line 325), especially when they reference Goldman, 2009. The idea in Goldman, 2009 is precisely that we can memorize without persistent neural activity.", "methods_and_evaluation_criteria": "I couldn’t find sufficient details of how inference was done in the model, either in the main text or the Appendix. For example, having in the Appendix the details of the LSTM architecture, how the number of units in the LSTM was chosen, whether it is bidirectional or unidirectional, and whether it is an LSTM that runs forward in time or backward in time would be helpful. I also couldn’t find sufficient details on how the synthetic datasets were generated. For example, it was not clear what a_1 and a_2 were set to be. How many trajectories were generated for this synthetic dataset? What was the initial condition of the dynamics? How was the dynamical system perturbed in this simulation? Without these details (sufficient enough for the readers to reproduce the authors’ results), it is difficult to evaluate the authors’ claim.\n\nThe authors claim that in the synthetic datasets, iSSM can recover the true latents and the underlying dynamics. However, the analyses didn’t have important details. For example, the authors show gray arrows (flow fields) in Figure 2A and B. Are the arrows both from the ground truth data or are they model-inferred? The flow fields for true and inferred seem identical. A clarification on this would be helpful. I would like to see how well iSSM recovers the flow fields, not just the trajectories.\n\nIt was also not mentioned in the paper (both main text and Appendix) how SSM is defined. I think the definition of what SSM is should be in the main text. Is it a linear dynamical system (LDS) model with nonlinear observations similar to Gao et al., 2016? What inference method was used for SSM? Is it similar to the inference procedure for iSSM?\n\nFurthermore, I am uncertain whether this is the right type of synthetic datasets because it seems that the number of observations here is only two, whereas in the real datasets that the authors look at, the number of neurons tends to be higher.\n\nIs “Dynamic Attractor” a ring attractor? The flow field shown in Figure 2A1 and 2B1 seem to suggest this. However, if this is a ring attractor, I am confused because in Equation (6), dx/dt = 0 if and only if x_1 = 0 and x_2 = 1. This means that y_t = [cos(1); 0] for all t. However, the trajectory y_t is not a circle, contrary to what is shown in Figure 1A1 and 2B1. I think it is not clear to me how the flow field plots are generated. For “Dynamic Attractor”, is there a way for the neural state to move without perturbations? Rotational dynamics are seen in neural data even when there are no optogenetic or electric stimulations. \n\nIn Figure 2B, the observation states are sometimes negative, but in Figure 2C the observations and latents are always positive. Is the trajectory in Figure 2A-B not the same as C? A clarification on how many trajectories are generated in total in this dataset would be helpful here as well.\n\nIt was not mentioned in the paper how B with sparsity penalty was implemented. How was the scale parameter s chosen? What was s in the experiments?", "theoretical_claims": "I did not check the proofs of the theoretical claims. However, I would like to note that one of the assumptions is that the function f is an injective function. It was not clear how this was enforced during training of the authors’ model.", "experimental_designs_or_analyses": "It was not clear how reconstruction accuracy was calculated in Figure 2. I am also surprised that using this metric, reconstruction accuracy is close to 0 for both iSSM and SSM in Figure 2D. It was not clear how the error bar is computed. For this analysis, are the authors assuming that the latent dimensionality is 2?\n\nI think it is concerning that there is not enough details in the experiments to reproduce the results. For example, for the real datasets, how was the latent dimensionality chosen? What was the hyperparameter space? How was the hyperparameter optimized? How was the training and test split made? When the authors mention they do cross-validation, how many folds? (e.g., “To find the optimal hyperparameters, we performed cross-validation.” (Line 379) How many folds were used? What fold was used for showing the result in Fig 4D?)\n\nConsistency score is defined only with respect to B. Have the authors evaluated consistency for A?\n\nIn Figure 3, the mean latents for correct and incorrect trials seem very similar. Given that the reconstruction accuracies for both SSM and iSSM are less than ~0.5, I’m not sure if the latent trajectories here should be trusted. To have some confidence, I think it would make sense to fit this dataset to an established model such as LDS with linear observations to see how SSM and iSSM perform with respect to this LDS baseline.\n\nIs data relevant to Figure 4 spike trains? It is odd that the plots in Figure 4A are in terms of firing rates, and not binned spike counts. Do the authors assume Poisson observations or Gaussian?", "relation_to_broader_scientific_literature": "A major strength of the paper is that the idea of incorporating interventional data into SSMs is interesting, and, if it works, would be an important contribution to the field. As the authors have pointed out, inferring latent dynamics from neural population activity is difficult, and if we have interventional data that can further constrain the model to learn the correct dynamics, this would help us understand neural computation. However, the experimental evidence supporting this idea is weak.", "essential_references_not_discussed": "I can’t think of essential references not cited.", "other_strengths_and_weaknesses": "Please see above.", "other_comments_or_suggestions": "I have no other comments or suggestions. The paper was written clearly.\n\n## update after rebuttal\n\nI sincerely thank the authors for the detailed response. I could access the link and look at the details on estimated flow fields, 6d synthetic experiments, and details on initialisation, hyperparameters and inference (though I believe full details on the inference should be available in the Appendix for this to be published).\n\n**Dynamic Attractor:** if the equation should be $\\frac{dx}{dt} = [a_1x_1, a_2(1-x_1)]^\\top$ instead of $\\frac{dx}{dt} = [a_1x_1, a_2(1-x_2)]^\\top$, then I don't think there is a fixed point in this system if $a_1 = -20$ and $a_2 = 1.2$ as stated in the shared link. This is inconsistent with the authors' claim that the system is a dynamic attractor.\n\n**Sparsity on B:** I re-read the Appendix of the original submission and could not find information about the loss function. The loss shows up in 3.2, but it doesn't have any details on regularization. The shared link says the regularization coefficient was found using cross-validation, but no more details could be found.\n\n**Rotational Dynamics:** Should we be worried that the red flow field looks more like a line attractor, whereas for the black arrows, there doesn't seem to be a line of fixed points? Have the authors compared the eigenvalues?\n\n**Reconstruction accuracy:** It makes sense that if the latents are not identifiable, the correlation could be close to 0 for the latents (D1), but it is still puzzling to me why the neural reconstruction (D2) is also close to 0. This is also partly why I think having benchmarks to models like LDS (with linear emission) could be helpful. Does the LDS show that the neural reconstruction accuracy is close to 0?\n\nThese remaining unresolved questions raise concerns about the validity of the analysis, and I will keep my score as is.", "questions_for_authors": "1. The authors mention that “models that are built upon observational data are not able to capture neural dynamics outside of the low-dimensional space.” However, it is not clear to me whether this method can model dynamics that are outside the latent x. How do the authors model dynamics that are outside the intrinsic manifold?\n\n2. Can iSSM model inputs that are not interventional as well? If so, how does that affect its performance in comparison to SSM?\n\n3. Could we have used something other than LSTM for the posterior? I wonder if there are any alternative ways of doing inference.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "uEEFME8KEe", "reviewer": "Reviewer_kBw8", "summary": "The paper introduces Interventional State Space Models (iSSM), a class of causal state-space models designed to identify neural circuit dynamics and predict responses to causal manipulations, addressing the limitations of traditional state-space models (SSMs) which lack causal interpretability. By explicitly modeling perturbations as interventions on latent neural dynamics, the authors demonstrate theoretically that iSSMs achieve identifiability of latent dynamics under suitable conditions. Empirically, they validate their model using simulated motor cortex dynamics and biological datasets involving calcium imaging from mouse ALM neurons during targeted photostimulation, and electrophysiological recordings from macaque dlPFC during micro-stimulation. Across these examples, iSSM consistently outperforms conventional observational SSMs in reconstructing neural activity, demonstrates robustness to initialization, and effectively generalizes to predict neural responses under novel interventions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate you and the other reviewers for your time and thoughtful evaluation of our work. We found the feedback to be highly constructive, as well as both fair and encouraging. Below we address your specific questions:\n\n**Correct vs. incorrect trials:** Please notice the latents closer to the end of the trial where the correct and incorrect latents start to diverge. We noticed that that’s the time when the animal is performing its action selection. This is relevant for us because the area where the recordings are performed ALM is thought to involve action selection. Therefore our results suggest that the causal latents discovered by the model distinguish between the action selection mechanism for the correct vs. incorrect trials.\n\n**Generalization to unseen interventions:** As described in the last paragraph of section 4.3, we divide the data into train and test trials with different interventions performed in each set. Our results show that (Fig. 4 C,D,E) that the test reconstruction accuracy is on par with training reconstruction accuracy, supporting the argument that the model is able to generalize to test interventions.\n\n\n**Compare against other baselines:** Thank you for pointing at the weakness of the paper. We want to clarify that our contribution is not developing a new latent variable model, rather it’s adding the interventional components to the existing SSMs. To our best knowledge, this is the first causal dynamical systems model that incorporates latents and observational noise. Therefore we are not aware of other alternatives to use for our comparisons, we are happy to include new comparisons if the reviewer has specific methods in mind.\n\nThanks again for your time; we look forward to your final assessment.\n\n&nbsp;\n\n### **Further responses to reviewer xGnX**\n**Explanation of dynamic attractor:**  We borrowed the term \"dynamic attractor\" from the previous literature on this toy example [1]. In our understanding, a ring attractor in the neuroscience literature implies a specific latent mechanism where the spatial information is encoded in the connectivity pattern of the network. Here the dynamic attractor has non-rotational latents that are nonlinearly transformed to generate rotational dynamics in the observational space. Indeed the dynamical structure generated in the observational space is a ring attractor but we remained consistent with the terms in the literature. \n\nIndeed you are correct that in the non-noisy setting $\\boldsymbol{x}=(0,1)$ or equivalently $\\boldsymbol{y}=(0,0)$ is an unstable fixed point of the model. So the model will not move away from it. However, except for that pathological case, all other settings of $\\boldsymbol{x}$; $\\frac{d\\boldsymbol{x}}{dt}$ will be nonzero and the model does not require perturbation to move along the circle.\n\n[1] Galgali, A. R., Sahani, M., and Mante, V. Residual dynamics resolves recurrent contributions to neural computation. Nature Neuroscience, 26(2):326–338, 2023.\n\n**Positive observations in Fig. 2C**: Thank you for your observation. This is a plotting issue, and the y-axis in the plots is not meaningful. We will fix that in the final revision (if accepted). We min-subtract the data in each dimension to stack different axes of the latents on top of each other. What’s shown in the state plots (Fig. 2 A, B) shows the exact numbers corresponding to the trajectories. \n\n**Sparsity penalty on $\\boldsymbol{B}$**: We mention in line 132-134 (right column) that the sparsity penalty is implemented by placing a Laplace prior on the elements of $\\boldsymbol{B}$ with a scale parameter $s$. In our real data experiments we do cross-validation over $s$ (the columns of Fig. 3C and Fig. 4D). For synthetic data we don’t use a sparsity prior.\n\n**Injective $f$**: While we don't explicitly enforce that the function f is an injective function, this function maps from a low-dimensional latent space for x to a high-dimensional observational space. So if this function fits the data well, then this function is often injective without explicitly enforcing it, as long as none of the latent dimensions collapses and are not used --- namely many different observations are mapped to the same latent variable value. We can often empirically check that this issue does not occur.\n\n**Reconstruction accuracy:** The reconstruction accuracy is defined as the correlation coefficient between true and reconstructed signals (latents or observations). Notice that except for the case of stimulation count is 0, the observational reconstruction accuracy is always close to 1 for both models. However, the latent reconstruction accuracy which directly measures the identifiability of the latents is only close to one for iSSM as the stimulation count increases. This is the main signature that the added interventional component to the model indeed leads to identifiability. The latent dimension for all the results in this figure is set to the ground truth 2.", "claims_and_evidence": "The main claims are well-supported by both theoretical arguments and empirical results.", "methods_and_evaluation_criteria": "Yes, the paper clearly demonstrates alignment between the proposed methods, evaluation approaches, and the key problems it addresses. The authors provide theoretical guarantees showing that iSSMs achieve identifiability under suitable assumptions (e.g., bounded completeness, piecewise linearity, sufficient interventions), thus enabling accurate predictions under novel perturbations. They empirically validate these claims using synthetic motor cortex simulations and biological datasets (calcium imaging in mice and electrophysiology in macaques), showing that iSSM consistently outperforms standard observational models (SSM) in reconstructing neural activity, recovering latent dynamics, and ensuring robustness of learned parameters across initializations.", "theoretical_claims": "Yes.", "experimental_designs_or_analyses": "Yes. I have the following questions:\n- In Figure 3.B, what's the difference of latent activities between correct and incorrect trials?\n- The goal of Section 4.3 is to evaluate iSSM's ability to generalizetounseeninterventions. Could you explain more about how this is achieved via Figure 4?", "relation_to_broader_scientific_literature": "The paper’s key contributions build upon and directly address limitations highlighted in prior literature: traditional SSMs, effectively capture statistical correlations but fail to offer causal interpretability.  This work specifically incorporates theoretical insights on the identifiability of nonlinear causal dynamical systems and extends them by allowing observational noise, a crucial step for realistic biological modeling", "essential_references_not_discussed": "Yes.", "other_strengths_and_weaknesses": "Strengths:\n- Combining causal inference concepts with traditional state-space modeling techniques, bridging a critical gap between purely statistical and causal modeling approaches.\n- Advances theoretical understanding by providing formal identifiability conditions, clearly extending prior theoretical work to the specific and challenging domain of neural dynamics.\n\nWeakness:\n- Only a standard SSM is used as a baseline. The authors should consider including additional baseline models specifically designed for causal dynamical systems to better demonstrate the superiority of iSSM.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "N/A", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "sw5j10hTbW", "reviewer": "Reviewer_KJn6", "summary": "The paper introduces interventional state space models (iSSM), a novel framework designed to predict neural responses to novel perturbations, addressing the limitations of traditional state space models (SSM) which capture statistical associations without causal interpretation. The authors establish the identifiability of iSSM. The model is validated through simulations of the motor cortex and applied to biological datasets, including calcium recordings from mouse ALM neurons and electrophysiological recordings from macaque dlPFC. The results demonstrate that iSSM outperforms traditional SSMs in terms of parameter identifiability and predictive power.\n\n## update after rebuttal\nThanks for the detailed response and the additional sections on reproducibility, inference, and computational complexity. I also appreciate the short discussion on the optimal latent dimension. These changes address my feedback and strengthen the paper. My evaluation remains positive, and I’ll keep my original score of 4 (Accept) unchanged. I encourage the authors to further explore testing the optimal latent dimension in future work.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate you and the other reviewers for your time and thoughtful evaluation of our work. We found the feedback to be highly constructive, as well as both fair and encouraging. Below we address your specific questions:\n\n**Nonlinearity limitation:** We completely agree with the reviewer that linear dynamics introduce a limitation into the model. However, developing theoretical results for more complex models is quite challenging. To our knowledge this is the first paper that develops such results for a non-LDS noisy SSM. We hope this brings the causal inference and neuroscience communities closer together and inspires the development of more theoretical results for more complex models.\n\n**Complexity of implementation:** Thanks for pointing this out, we have included a new section in the supplementary material detailing the inference and its computational complexity. To summarize the section we find that our algorithm scales as $\\mathcal{O}(TN^2 + TD^2)$ where $T$ is the length of the trajectories, $D$ is the latent dimension, and $N$ is the number of neurons. Therefore our inference is linear in time and quadratic in the number of neurons.\n\n**Limited score of interventions:** The existing evidence for interventions forcing the neural state outside of the observational trajectories is discussed in both prior theoretical and experimental papers [1, 2, 3]. That said, we believe that collectively as a field our understanding of how interventions affect the neural states is quite limited. The model we present here is among the first models that incorporate interventions into the modeling framework and therefore we believe it’s critical to develop interventional models to better understand how interventions influence neural dynamics.\n\n[1] Jazayeri, M. and Afraz, A. Navigating the neural space in search of the neural code. Neuron, 93(5):1003–1014, 2017.\n\n[2] O’Shea, D. J., Duncker, L., Goo, W., Sun, X., Vyas, S., Trautmann, E. M., Diester, I., Ramakrishnan, C., Deisseroth, K., Sahani, M., et al. Direct neural perturbations reveal a dynamical mechanism for robust computation. bioRxiv, pp. 2022–12, 2022.\n\n[3] Daie, K., Svoboda, K., and Druckmann, S. Targeted photostimulation uncovers circuit motifs supporting short-term memory. Nature neuroscience, 24(2):259–265, 2021.\n\n**Expressivity of the emissions:** While we agree that linear dynamics is a limiting factor of the model, we would like to point out that adding nonlinear emissions to the linear dynamics changes its expressivity quite drastically. We reference the Koopman theory in the paper, which posits that a (possibly infinite dimensional) linear dynamical system followed by a nonlinear emission is sufficient to express \"any\" autonomous dynamical system. The catch here is of course \"infinite dimensional\" but intuitively as the latent dimension grows larger the expressivity of the model also increases.\n\n**Increasing latent dimension:** You are absolutely correct that if we assume low-dimensionality of the data then increasing the latent dimension is counter-intuitive. We would like to clarify that while we think the optimal latent dimension crucially depends on the task and dataset, we think of interventions as kicking the state of the system outside of its low-dimensional manifold increasing the dimensions spanned by the trajectories. Therefore the conjecture we've included in the discussion section is only valid for large datasets where diverse and numerous interventions are performed. But we have only presented this argument as a conjecture and we agree that experimental validation is needed for confirming it.\n\n**Discussion of the implementation:** Thank you for the suggestion, we added a new section in the supplementary material on the details of our variational inference algorithm. We have also included our code package in the submission which follows standard code development practices. \n\n**LSTMs, hyperparameters:**\nUsing LSTMs in the variational posterior has become a standard practice in the field (e.g. [1]). We did experiment with vanilla RNNs but they drastically underperformed LSTMs. We will include these results if the paper gets accepted. Regarding the hyperparameter search, our two hyperparameters are the sparsity of $\\boldsymbol{B}$ denoted by $s$ and the latent dimension $D$. The cross-validation and the range of parameters are included in Fig. 3C, 4D. We used a logarithmic axis for s to cover a wide range of possibilities and for $D$ we increased it up to a point where we see a drop of validation performance. Regarding how the hyperparameters change w.r.t. the task at hand, this question is highly problem specific. We think these value highly affected by the details of the experimental setup such as the task, type of recording, animal model, etc.\n\n[1] Krishnan, Rahul G., Uri Shalit, and David Sontag. \"Deep Kalman Filters.\" arXiv preprint arXiv:1511.05121 (2015).\n\nThanks again for your time; we look forward to your final assessment.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "I checked the proof of Theorem 3.4. It seems correct.", "experimental_designs_or_analyses": "Yes", "relation_to_broader_scientific_literature": "The paper could be related to causal discovery and control of dynamical systems.", "essential_references_not_discussed": "-", "other_strengths_and_weaknesses": "## Strengths\n- Theoretical Contribution: The paper provides a theoretical foundation for the identifiability of iSSM, which is crucial for causal inference in general especially when the goal is to design some sort of control strategy to steer the system towards a desired state.\n- Empirical Validation: The model is tested on both synthetic and real-world datasets, showing generalization capabilities and parameter recovery.\n- Practical Relevance: The application of iSSM to real biological datasets demonstrates its potential utility in neuroscience research, particularly in understanding neural dynamics under perturbations. Although the current state of technology may not allow a granular intervention on latent dimensions, the result of this work could be of higher practical use when such interventional capabilities become possible.\n\n## Weaknesses\n- Nonlinearity Limitation: The current model focuses on linear dynamics, and while the inference model can capture nonlinearities, explicitly modeling nonlinear dynamics in the latent space remains a limitation.\n- Complexity of Implementation: The implementation details and computational complexity of iSSM are not thoroughly discussed, which could be a barrier for adoption by other researchers. \n- Limited Scope of Interventions: The paper primarily assumes that interventions force the neural state out of its low-dimensional attractor manifold. I think this assumption and its consequences need to be deeper understood as they seem to be critical for the results of this work.", "other_comments_or_suggestions": "-", "questions_for_authors": "- Although the emission model introduces nonlinearity to the model, the dynamics in the latent space is linear and cannot express qualitative behaviours that only nonlinear dynamical systems may have, e.g. limit cycles. This might be a guide for where not to use this model. Can you provide examples of neural circuits for which the linear dynamics is insufficient to exhibit the underlying dynamics regardless of how expressive the emission model is?\n- It seems the latent dimension is crucial for this model as those dimensions are supposed to correspond to biologically meaningful behaviour and to be physically intervenable. The stated observation in the paper “increasing the number of latents always provides better test accuracy” seems counterintuitive as there should be a number of dimensions that corresponds to the physical properties of the system and beyond that there is a risk of overfitting. Could you elaborate on this observation and any reasoning behind it?\n- I recommend including a pseudocode in the appendix and provide a more detailed discussion of the implementation. An ablation study would also be useful to better understand the author’s design choice, e.g. why choosing LSTM as the recognition network over other options. The authors have mentioned they have done cross-validation to find the optimal set of hyper-parameters, but it’s unclear which hyperparameters were involved in that search and with which range. Also how the range of search for hyper-parameters changes with the problem at hand.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper introduces the interventional state space model (iSSM), a principled extension of SSMs designed to handle external interventions in neural dynamics. The key strengths of the submission lie in its clear alignment between problem formulation, proposed method, and evaluation. Theoretical results establishing identifiability under appropriate assumptions are compelling, and the paper provides experimental evaluation on both synthetic and real-world datasets.

**However, several issues limit the paper’s overall strength**. The evaluation could be improved: the paper does not include benchmarks against other models (e.g., even a simple LDS). The comparison would have been helpful, especially since the models' performance on the simple synthetic datasets is concerning (e.g., when the number of interventions is 0, neural reconstruction is close to 0).

Their main methods figure remains unclear: the dynamics equation -- particularly in the Dynamic Attactor example -- does not appear to be consistent with the authors' figure depicting the flow field, which shows the ring attractors. Inferred flow fields were not presented initially and remain only partially analyzed post-rebuttal, which must be in the paper. There are also concerns about the choice of evaluation metrics for Figure 4 (e.g., Pearson’s r vs. log-likelihood for spike trains), and incomplete clarity around inference details despite an appendix update.

In sum, the contribution is promising with several reviewers recommending accept, yet key issues do remain. In summary, the authors should strengthen the empirical analysis, clarify model-inference details (and code), and improve quantitative evaluation of learned dynamics prior to publication.
Decision: Accept (poster)