q-bio.NCMar 5, 2026

Behavior-dLDS: A decomposed linear dynamical systems model for neural activity partially constrained by behavior

Authors: Eva Yezerets, En Yang, Misha B. Ahrens, Adam S. Charles

Organizations: Department of Biomedical Engineering Center for Imaging Science Johns Hopkins University · Department of Biology University of North Carolina · Janelia Research Campus Howard Hughes Medical Institute

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Brain-wide recordings of large-scale networks of neurons now provide an unprecedented view into how the brain drives behavior. However, brain activity contains both information directly related to behavior as well as the potential for many internal computations. Moreover, observable behavior is executed not only by the brain, but also by the spinal cord and peripheral nervous system. Behavior is a coarse-grained product of neural activity, and we thus take the view that it can be best represented by lower-dimensional latent neural dynamics. Capturing this indirect relationship while disambiguating behavior-generating networks from internal computations running in parallel requires new modeling approaches that can embody the parallel and distributed nature of large-scale neural populations. We thus present behavior-decomposed linear dynamical systems (b-dLDS) to disentangle simultaneously recorded subsystems and identify how the latent neural subsystems relate to behavior. We demonstrate the ability of b-dLDS to decouple behavioral vs. internal computations on controlled, simulated data, showing improvements over a state-of-the-art model that uses behavior to supervise all dynamics based on behavior. We also demonstrate b-dLDS's interpretability benefits on a task-driven RNN dataset featuring a nonlinear relationship between behavior and activations. We then show that b-dLDS can further scale up to tens of thousands of neurons by applying our model to a large-scale recording of a zebrafish hindbrain during the complex positional homeostasis behavior, wherein b-dLDS highlights asymmetry in behavior-related dynamic connectivity networks.

Figures & tables

Figure 1: behavior-dLDS architecture. A: The latent, lower-dimensional neural manifold represents the processes that generate observable neural activity, mapped to the ambient space via D\bm{D} , and behavior, via Ψ\bm{\Psi} . B: The established dLDS model, inside the dashed lines, is updated here with a mapping from dynamics coefficients c\bm{c} to behavior b\bm{b} via Ψ\bm{\Psi} . This contributes another regularization term to the joint inference of latent states x\bm{x} and dynamics coefficients c\bm{c} .
Figure 2: Simulated system composed of two independent sub-systems with 10 behavior traces. Example trial (of 50 randomly sampled trials); R2R^{2} calculated over all 50 samples. A,B: “Neural” ground truth and reconstruction by b-dLDS. Signal is generated from a set of 6 ground truth dynamics operators and corresponding dynamics coefficients, with D=I\bm{D}=\bm{I} . C,D: Simulated behavior and reconstruction from b-dLDS, where behavior is simulated from the ground truth dynamics coefficients. E: True dynamics coefficients used to generate the data. F: Learned dynamics coefficients (absolute value). G: Learned dynamics operators (1-15). H: Learned Ψ\bm{\Psi} , where each column corresponds to a dynamics operator and dynamics coefficient. Only column 3 is nonzero, and it perfectly reconstructs column 1 of the ground truth Ψ\bm{\Psi} . This is the correct column to reconstruct because the first ground truth dynamics coefficient was used to generate the simulated behavior traces.
Figure 3: b-dLDS vs. CLDS and PSID. CLDS learns only behavior-conditioned dynamics, while b-dLDS learns both behavior-generating and intrinsic ones. PSID learns only a shared latent neural-behavioral space, while b-dLDS learns a hierarchical relationship between latent neural states, latent dynamics, and behavior. (Note: CLDS requires an 80/20 train/test split for optimization, so all models (CLDS, b-dLDS, and PSID) were only trained/learned on 40 out of 50 simulated trials (where CLDS then saw the remaining 10 trials for testing). Results shown here are on the same 40 “train” trials across models.) A-G: b-dLDS, 2 independent systems simulation, 1 trial out of 40 shown, D=I\bm{D}=\bm{I} , with simulated behavior generated as a linear combination of ground truth dynamics coefficients. R2R^{2} calculated across 40 trials. A, B: 1 trace out of 8 shown, data and reconstruction. E: True dynamics coefficients c\bm{c} . Stars match selected time points in H . Note that at the first two time points, one top-left (purple) dynamics operator and one bottom-right (yellow) dynamics operator are active, while at the third time point, only a top-left dynamics operator is active. D: b-dLDS-inferred c\bm{c} . The b-dLDS model is initialized with more than enough dynamics operators (15 learned vs. 6 ground truth operators). Noisy dynamics operators in F correspond to unused (black) dynamics coefficients in D . E: Behavior and reconstruction (1 trace). G: Ψ\bm{\Psi} is a row vector with one dynamics coefficient mapping to one behavior (true value: 1, learned value: 0.8765, other true values 0, learned near 0). H: CLDS learns an 8x8 dynamics matrix A^t\widehat{\bm{A}}_{t} , while b-dLDS learns f^\widehat{\bm{f}} and c^\widehat{\bm{c}} , which can be combined via ∑m=1Mf^mc^mt\sum_{m=1}^{M}\widehat{\bm{f}}_{m}\widehat{\bm{c}}_{mt} to create an equivalent matrix. We show 3 example time points here as described above. b-dLDS learns the block-diagonal dynamics matrices, while CLDS does not. The full comparison across time is shown in Figure 0.4 . PSID learns a dynamics matrix A^\widehat{\bm{A}} that describes the relationship between latent states. I: CLDS neural data reconstruction (1st trace). J: PSID neural data reconstruction (1st trace). K: b-dLDS achieves lower relative mean squared error than CLDS on dynamics reconstruction, regardless of the number of dynamics or the number of behaviors. CLDS used more than 1 TB of RAM when we tried to run it on 10 behavior traces.
ModelIterationsλΨ\lambda_{\Psi}L1L_{1} / L2L_{2} normRb2R^{2}_{b}
b-dLDS50052.60.99
  500.660.99
  1000.390.99
  2000.070.96
b-dLDS-x5000.96443.880.95
  54.440.96
Table 1: b-dLDS vs. b-dLDS-x Ψ\bm{\Psi} sparsity and behavior reconstruction across levels of regularization on the Frobenius norm of Ψ\bm{\Psi} . For b-dLDS, we start with 5 as this is the setting used in the model shown in Figure 4 B. For b-dLDS-x, we start with 0.9644 as this is the BayesOpt-identified setting used in the model shown, and then compare to other stronger regularization settings.
Figure 4: b-dLDS vs. b-dLDS-x on delayanti task-driven RNN data . One seed of each model shown. D=I\bm{D}=\bm{I} . For behavior data reconstruction, b-dLDS R2=0.99R^{2}=0.99 , while b-dLDS-x R2=0.95R^{2}=0.95 . The b-dLDS-x behavior reconstruction is constrained to be a linear combination the latent states, while b-dLDS learns a sparser, more interpretable set of dynamics coefficients based on the behavior. A: TaskRNN model activations, “behavior” traces (23 total inputs and outputs), and latent states (equivalent to the activations because there is no noise in the observations; the activations are known directly). B: b-dLDS c^\widehat{\bm{c}} and b^\widehat{\bm{b}} example. C: Corresponding b-dLDS Ψ^\widehat{\bm{\Psi}} . D: b-dLDS-x c^\widehat{\bm{c}} and b^\widehat{\bm{b}} example. E: Corresponding b-dLDS-x Ψ^\widehat{\bm{\Psi}} . F: Comparison of Ψ^\widehat{\bm{\Psi}} across regularization strengths on the Frobenius norm for b-dLDS vs. b-dLDS-x.
Figure 5: Four strongest Ψ\bm{\Psi} connections between dynamics coefficients and motor activity in a zebrafish model with one behavior trace. b-dLDS R2=0.91R^{2}=0.91 for neural data reconstruction (over whole time series concatenated; per neuron varies); R2=0.97R^{2}=0.97 for behavior data reconstruction. Dynamic connectivity maps for dynamics operators and corresponding dynamics coefficients shown. Motor-aligned dynamics are utilized throughout the recording or a mix of earlier and/or later in the recording.
Criteria“Left” side“Right” side% asymmetry
# neuron centroids255024006%
Variance of dF/F over time per tracemean 0.0099, median 0.0052, std 0.0153mean 0.0109, median 0.0057, std 0.01749% (means)
R2R^{2} per neuronmean 0.38, median 0.35, std 0.23mean 0.46, median 0.43, std 0.2717% (means)
# start/end points for largest connections  6 of 25 DOs: 18%-118%
Table 2: b-dLDS connectivity asymmetry in zebrafish hindbrain dynamics. To form a baseline of asymmetry that might stem from other sources (e.g., imaging artifacts) we compute the asymmetry in the geometry and variance of the fluorescence data itself. These baselines indicate that up to 17% asymmetry might be attributed to other sources. Of the 25 DOs, 6 displayed stronger asymmetry, with more start and/or end points of their connectivity on the left side than on the right.
Appendix figures & tables17 assets

Supplementary material from the paper’s appendix.

Appendix

List of Supplementary Figures 0.1 Two-independent-systems simulation with 10 behaviors generated from the first two ground truth states, with some overlap in their combined effects. A: Ground truth signal fed into the model. B: Signal reconstruction from b-dLDS. C: True dynamics coefficients used to generate the data. D: Learned dynamics coefficients (absolute value). E: Learned dynamics operators (1-15). F: Learned Ψ\bm{\Psi} . G: Best reconstruction of the true Ψ\bm{\Psi} (2 columns, other columns all 0). H: Highlighted dynamics coefficient used to simulate behavior data. I: Highlighted dynamics coefficient corresponding to best match of learned Ψ\bm{\Psi} . J: Highlighted corresponding dynamics operator. K: Simulated behavior ground truth. L: Behavior reconstruction.
List of Supplementary Figures 0.2 Two-independent-systems simulation with 10 behaviors generated from a combination of all 6 ground truth states. A: Ground truth signal fed into the model. B: Signal reconstruction from b-dLDS. C: True dynamics coefficients used to generate the data - one was randomly initialized as all zeros. D: Learned dynamics coefficients (absolute value). E: Learned dynamics operators (1-15). F: True Ψ\bm{\Psi} . G: Best reconstruction of the true Ψ\bm{\Psi} (6 columns). H: True behavior. I: Behavior reconstruction.
List of Supplementary Figures 0.3 Two-independent-systems simulation with 10 behaviors generated from the first ground truth state, with added Gaussian noise in the “neural” data. A: Ground truth signal fed into the model, with Gaussian noise scaled by 0.1\sqrt{0.1} . B: Signal reconstruction from b-dLDS. C: True dynamics coefficients used to generate the data - two were randomly initialized as all zeros. D: Learned dynamics coefficients (absolute value). E: Learned dynamics operators (1-15). F: Learned Ψ\bm{\Psi} . G: Best reconstruction of the true Ψ\bm{\Psi} (1 column, other columns all 0). H: Highlighted dynamics coefficient used to simulate behavior data. I: Highlighted dynamics coefficient corresponding to best match of learned Ψ\bm{\Psi} . J: Highlighted corresponding dynamics operator. K: Simulated behavior ground truth. L: Behavior reconstruction.
List of Supplementary Figures 0.4 CLDS vs. b-dLDS learned dynamics, continued. CLDS learns an 8x8 dynamics matrix At\bm{A}_{t} , while b-dLDS learns f\bm{f} and c\bm{c} , which can be combined via ∑m=1Mfmcmt\sum_{m=1}^{M}f_{m}c_{mt} to create an equivalent matrix. We compare each block of the matrix (top left 4x4 through bottom right 4x4) to the ground truth dynamics from the simulation. b-dLDS closely matches the ground truth and correctly learns that the two off-diagonal corners are 0, indicating 2 independent systems; CLDS instead incorrectly follows the behavior in all dynamics. Each row of each heatmap is an entry of the 4x4 block of the overall 8x8 matrix (A12 means row 1 column 2, 16 rows per heatmap).
ModelSeedRy2R^{2}_{y}Rb2R^{2}_{b}MaxIdxMaxIdxMaxIdxTime (s)RAM (KB)
    Ψ^\widehat{\bm{\Psi}}maxR2R^{2} ,maxR2R^{2}max  
    valueΨ^\widehat{\bm{\Psi}}c^\widehat{\bm{c}} vs.c^\widehat{\bm{c}}x^\widehat{\bm{x}} vs.x^\widehat{\bm{x}}  
      b\bm{b} b\bm{b}   
b-dLDS10.910.960.68250.82250.3636.11E+04-2.17E+08
 20.90.950.29200.84250.36436.09E+04-2.18E+08
 30.910.970.7430.8430.33214.74E+047.91E+07
Appendix
Table 3: b-dLDS, dLDS, and b-dLDS-x model comparison, with runtimes and memory utilization. b-dLDS-x represents an alternative to b-dLDS where the behavior maps to the latent states rather than to the dynamics coefficients. With additional optimization, b-dLDS-x can match b-dLDS in overall behavior reconstruction. However, b-dLDS provides better interpretability: its max Ψ\bm{\Psi} values are larger, and so are its max R2R^{2} values between the dynamics coefficients and behavior, compared to the max R2R^{2} values between the dynamics coefficients and behavior for b-dLDS-x. In the two versions of b-dLDS shown here, two times out of three, these two indices (max Ψ\bm{\Psi} and max c\bm{c} vs. b\bm{b} R2R^{2} also corresponded; in contrast, they only corresponded in one version of b-dLDS-x,“b-dLDS-x with reg on Ψ\bm{\Psi} ”, which also had worse behavior reconstruction (see “Max Ψ\bm{\Psi} value” and “Idx max Ψ\bm{\Psi} ” vs. “Max R2R^{2} , c\bm{c} vs. b\bm{b} ” and “Idx max c\bm{c} ” for b-dDLS models, and “Max Ψ\bm{\Psi} value” and “Idx max Ψ\bm{\Psi} ” vs. “Max R2R^{2} , x\bm{x} vs. b\bm{b} ” and “Idx max x\bm{x} ” for b-dLDS-x models). Please note that runtime and memory utilization are only meant to be a rough estimate, as these were dependent on other processes happening in parallel on the server, and there is no built-in function for calculating memory utilization for MATLAB on Linux machines; instead we rely on subtracting start from end memory used, sometimes resulting in negative values ( Dvorak, 2026 ) .
ParameterMeaningHow to setSettings used for zebrafish models (in order of appearance in Table 3 )
    
nFnumber of dynamics operators (DOs)Based on models used in zebrafish figures. Also tried matching nF to N for fairness of x^\widehat{\bm{x}} and c^\widehat{\bm{c}} vs. b\bm{b} R2R^{2} comparisons.25 (first 7 model types); 50 (last 3 model types, “nF=50”)
Nnumber of latent dimensionsBased on models used in zebrafish figures.50
lambda_valregularization: sparsity of states ( x^\widehat{\bm{x}} )Based on models used in zebrafish figures.0.05
lambda_historyregularization: smoothness of states ( x^\widehat{\bm{x}} ), but more specifically, reconstruction of states from dynamics operators and dynamics coefficientsBased on models used in zebrafish figures (with exceptions for b-dLDS-x BayesOpt).0.23; 0.1154 (b-dLDS-x BayesOpt model types)
lambda_bregularization: sparsity of dynamics coefficients ( c^\widehat{\bm{c}} )Based on models used in zebrafish figures.0.05
Appendix
Table 4: Meaning of key b-dLDS parameters, heuristics for setting them, and settings used for model comparisons in Table 3 : one motor signal input as behavior to b-dLDS and b-dLDS-x.
List of Supplementary Figures 0.5 Ablation of regularization on Ψ^\widehat{\bm{\Psi}} update step. b-dLDS with these settings optimized (same simulation setup as shown in Figure 3 K) achieves lower relative mean squared error on dynamics reconstruction than ablated models, regardless of the number of behavior traces (1 or 10) or number of dynamics tied to the behavior traces (behavior traces are generated as a linear combination of the simulated dynamics traces).
ParameterMeaningHow to setSettings used for 4 simulation types
    
nFnumber of dynamics operators (DOs)First set to 2x the number of PCs required for 95% variance explained. Then decrease until the median number of active dynamics coefficients starts decreasing. (Goal: minimize the number of DOs while maintaining the inherent dimensionality of the data.) Alternative heuristic: use the number of known behavioral states as a benchmark and double it to get the starting point for the number of DOs, and then start decreasing from there.15
Nnumber of latent dimensionsKnown based on data, given that D=I\bm{D}=\bm{I} .8
lambda_valregularization: sparsity of states x^\widehat{\bm{x}}Increase to reduce the number of states active at the same time (moderately sensitive)1.0000e-04
lambda_historyregularization: smoothness of states x^\widehat{\bm{x}} , but more specifically, reconstruction of states from dynamics operators and dynamics coefficientsIncrease to encourage states to change smoothly from time point to time point (less sensitive)1.0000e-04
lambda_bregularization: sparsity of dynamics coefficients (c)Increase to reduce the number of DOs active at the same time (sensitive)0.25
Appendix
Table 5: Meaning of key b-dLDS parameters, heuristics for setting them, and settings used for 4 simulation model types: one dynamics coefficient maps to all 10 behaviors; 2 dynamics coefficients map to all 10 behaviors; all dynamics coefficients map to all 10 behaviors; and one dynamics coefficient maps to all 10 behaviors, with added Gaussian noise in the simulated data.
ParameterMeaningHow to setSettings used (if differ, first b-dLDS setting shown; then b-dLDS-x)
    
nFnumber of dynamics operators (DOs)Set to 2x the number of PCs required for 95% variance explained.16
Nnumber of latent dimensionsKnown based on data, given that D=I\bm{D}=\bm{I} .128
lambda_valregularization: sparsity of states ( x^\widehat{\bm{x}} )Optimize with BayesOpt.0.9772;0.0293
lambda_historyregularization: smoothness of states ( x^\widehat{\bm{x}} ), but more specifically, reconstruction of states from dynamics operators and dynamics coefficientsOptimize with BayesOpt.0.8702;0.0636
lambda_bregularization: sparsity of dynamics coefficients ( c^\widehat{\bm{c}} )Optimize with BayesOpt.0.25908;0.0315
Appendix
Table 6: Meaning of key b-dLDS and b-dLDS-x parameters, heuristics for setting them, and settings used for the TaskRNN experiments.
ParameterMeaningHow to setSettings used for zebrafish models (1 vs. 4 behavior inputs)
    
nFnumber of dynamics operators (DOs)First set to 2x the number of PCs required for 95% variance explained. Then decrease until the median number of active dynamics coefficients starts decreasing. (Goal: minimize the number of DOs while maintaining the inherent dimensionality of the data.) Alternative heuristic: use the number of known behavioral states as a benchmark and double it to get the starting point for the number of DOs, and then start decreasing from there.25
Nnumber of latent dimensionsSet to at least the number of PCs required for 80-90% variance explained.50
lambda_valregularization: sparsity of states x^\widehat{\bm{x}}Increase to reduce the number of states active at the same time (moderately sensitive)0.05
lambda_historyregularization: smoothness of states ( x^\widehat{\bm{x}} ), but more specifically, reconstruction of states from dynamics operators and dynamics coefficientsIncrease to encourage states to change smoothly from time point to time point (less sensitive)0.23
lambda_bregularization: sparsity of dynamics coefficients c^\widehat{\bm{c}}Increase to reduce the number of DOs active at the same time (sensitive)0.05
Appendix
Table 7: Meaning of key b-dLDS parameters, heuristics for setting them, and settings used for two zebrafish model types in Figure 5 and Supplementary Figures 0.6 , 0.7 , 0.8 , 0.9 , and 0.10 : one motor signal input as behavior, or one motor signal plus three trial types as one-hot encodings.
List of Supplementary Figures 0.6 Four weakest Ψ\bm{\Psi} connections between dynamics coefficients and motor activity in a model with one behavior. b-dLDS R2=0.91R^{2}=0.91 for neural data reconstruction (over whole time series concatenated; per neuron varies); R2=0.97R^{2}=0.97 for reconstruction of smoothed behavior data input to the model. Dynamic connectivity maps and dynamics coefficients shown. The weakest relationships to behavior are those dynamics operators that are rarely utilized throughout the recording.
List of Supplementary Figures 0.7 Four strongest Ψ\bm{\Psi} connections between dynamics coefficients and motor activity in a model with four behaviors (motor and 3 trial types based on stimulation patterns, i.e., fish offset by experimenter). Overall model neural data reconstruction R2=0.93R^{2}=0.93 ; reconstruction of smoothed behavior data input to the model R2=0.997R^{2}=0.997 (all 4 behavior traces concatenated). Dynamic connectivity maps and dynamics coefficients shown. Motor-aligned dynamics are utilized throughout the recording or later in the recording.
List of Supplementary Figures 0.8 Four strongest Ψ\bm{\Psi} connections between dynamics coefficients and swim trial type 1 in a model with four behaviors (motor and 3 trial types based on stimulation patterns between swims). Overall model neural data reconstruction R2=0.93R^{2}=0.93 ; reconstruction of smoothed behavior data input to the model R2=0.997R^{2}=0.997 (all 4 behavior traces concatenated). Dynamic connectivity maps and dynamics coefficients shown. Swim trial-aligned dynamics cover different time segments of the recording.
List of Supplementary Figures 0.9 Four strongest Ψ\bm{\Psi} connections between dynamics coefficients and swim trial type 2 in a model with four behaviors (motor and 3 trial types based on stimulation patterns between swims). Overall model neural data reconstruction R2=0.93R^{2}=0.93 ; reconstruction of smoothed behavior data input to the model R2=0.997R^{2}=0.997 (all 4 behavior traces concatenated). Dynamic connectivity maps and dynamics coefficients shown. Swim trial-aligned dynamics cover different time segments of the recording.
List of Supplementary Figures 0.10 Four strongest Ψ\bm{\Psi} connections between dynamics coefficients and swim trial type 3 in a model with four behaviors (motor and 3 trial types based on stimulation patterns between swims). Overall model neural data reconstruction R2=0.93R^{2}=0.93 ; reconstruction of smoothed behavior data input to the model R2=0.997R^{2}=0.997 (all 4 behavior traces concatenated). Dynamic connectivity maps and dynamics coefficients shown. Swim trial-aligned dynamics cover different time segments of the recording.
List of Supplementary Figures 0.11 Hyperparameter sensitivity experiments related to number of dynamics coefficients and model size. Data from Table 8 . 3 seeds shown for each point. A-D: Zebrafish model parameter sweeps: N, nF, regularization for sparsity of c^\widehat{\bm{c}} , regularization for b\bm{b} reconstruction. E,F: Simulation (Fig. 3 ) parameter sweeps: regularization for sparsity of c^\widehat{\bm{c}} , regularization for b\bm{b} reconstruction.
ModelParam.Val.Median c^\widehat{\bm{c}} useMeanRb2R_{b}^{2}MeanRy2R_{y}^{2}Mean* Setting used
ZebrafishN2576.670.980.970.730.72 
hindbrain 256 0.98 0.70  
  257 0.95 0.72  
  3566.000.970.970.780.78 
  356 0.99 0.78  
  356 0.97 0.79  
Appendix
Table 8: Hyperparameter sensitivity experiments. * indicates settings used in paper figures. Values to test were chosen based on values used in figures, plus 0.01, 0.1, and 1, at least. The results shown are not necessarily converged to their maximum potential R2R^{2} , for example, because only 500 iterations were run for each of these tests (rather than the 5000 used in zebrafish figures).
More in this topic

Latent Dynamics Modeling

Latent dynamics modeling learns compact hidden-state representations and transition rules to predict, explain, and control evolving systems.

See all 253 papers

Explore similar work

CardsList
Jun 9, 2026

Machine Learning Methods for Studying Latent Neural Activity Dynamics

Recent developments in brain recording are driving a demand for machine learning tools capable of decoding the latent structure of large populations of neurons. In this paper, we provide a comprehensive survey that outlines the trajectory of Latent Variable Models (LVMs) from early state-space models to more recent deep generative models. We organize the literature into three closely related domains: (1) Single-Region Latent Dynamics, which includes models such as linear dynamical systems to more complex dynamics represented by Recurrent Neural Networks (RNNs) and Neural Ordinary Differential Equations (ODEs); (2) Multi-Region Communication, which employs probabilistic as well as subspace methods to study how information is transferred across different brain areas considering synaptic propagation delays and network connectivity; and (3) Behavior-Aligned Modeling, which seeks to disentangle neural activity related to task performance from other internal states via supervised or contrastive learning. This survey also includes large-scale neural foundation models, such as Transformers and diffusion models, that rely on large-scale pre-training for optimal performance across subjects. Finally, we conclude and discuss benchmarks, evaluation criteria, and open challenges, such as the ability to identify causal links or directionality of communication, to facilitate future research for bridging interpretable brain dynamics with reliable neural decoding.

Shufeng Kong, Fumei Deng, Xinyi Dong +7

School of Software Engineering, Sun Yat-sen University, Zhuhai, China · School of Computer Science and Artificial Intelligence, Foshan University, Foshan, China · Department of Ecology and Evolutionary Biology, Cornell University, Ithaca, NY, USA +2

Dynamical SystemsLatent Variable ModelsLatent Dynamics ModelingNeural ProcessesNeural Decoding
AskPDF
Oct 1, 2026

Inferring Multi-Timescale Neural Dynamics with Switching Linear Dynamical Systems

Neural activity often exhibits multiple timescales that can vary with behavioral states and task conditions. Identifying these timescales from neural recordings is important for better understanding neural computation and function. However, traditional approaches based on autocorrelation fitting are difficult to scale to high-dimensional population recordings and can become unreliable when neural dynamics change with behavior. State-space models have been a powerful framework for modeling high-dimensional neural population activity through latent dynamical systems, but standard formulations and inference methods do not explicitly account for multiple timescales and therefore do not guarantee accurate recovery of the underlying temporal structure. Motivated by these questions, we introduce the Multi-Timescale Switching Linear Dynamical System (MTS-SLDS), a framework for identifying regime-specific latent timescales from continuous or spiking neural observations. MTS-SLDS combines a multi-lag moment initialization, which captures temporal structure across multiple observation lags, with regime-conditioned Laplace-EM inference, which reduces mixing of dynamical statistics across uncertain regimes. Characteristic timescales can then be extracted directly from the eigenvalues of the learned latent transition matrices. In synthetic and neural experiments with Gaussian and Poisson spike observations, MTS-SLDS accurately recovers timescales and switching structure over multiple datasets.

Lulu Gong, Yongxu Zhang, Shreya Saxena

Center for Neurocomputation and Machine Intelligence, Wu Tsai Institute, Department of Biomedical Engineering, Yale University, New Haven, CT 06510, USA

Stochastic Linear Dynamical SystemsLatent Dynamics ModelingRegime-Switching Dynamical SystemsMultiscale Time-Series Modeling
AskPDF
Apr 26, 2026

Integrative neurocybernetic modeling in the era of large-scale neuroscience

Large-scale neuroscience is generating rich datasets across animals, brain areas, and behavioral contexts, yet our modeling efforts remain fragmented across isolated experiments. We argue that understanding behavior requires integrative neurocybernetic models: understandable dynamical models that capture the closed-loop coupling of brain, body, and environment, treat the brain as a controller pursuing latent objectives, represent structured variation across scales, and scale to heterogeneous datasets. Such models shift the goal from predicting neural recordings in isolation to inferring the organizing principles that govern neural and behavioral dynamics. We outline a practical route toward this goal by combining nonlinear state-space models and meta-dynamical extensions with scalable inference, knowledge distillation, and mixed open- and closed-loop training. By pooling complementary constraints from recordings, behavior, perturbations, and anatomy, integrative neurocybernetic models can provide statistical amplification, few-shot generalization, and mechanistic insight into shared dynamical structure, individual variation, and the control objectives that govern behavior. This agenda offers a model-centric path from fragmented data to a mechanistic science of how brains produce behavior.

Il Memming Park, Gonzalo G. de Polavieja, Ayesha Vermani +12

Champalimaud Centre for the Unknown, Champalimaud Foundation, Lisbon, Portugal · Ryv​​ivyR Inc, NY, USA · Harvard University, Cambridge, MA, USA +3

Computational Cognitive ModelingLatent Dynamics Modeling
AskPDF