Dynamical modeling of nonlinear latent factors in multiscale neural activity with real-time inference
Authors: Eray Erturk
Maryam M. Shanechi
Venue: NeurIPS 2025 poster ()
Keywords: Multimodal Deep Learning
Missing Data
Dynamical Models
Neuroscience
Primary area: neuroscience_and_cognitive_science

Abstract: Real-time decoding of target variables from multiple simultaneously recorded neural time-series modalities, such as discrete spiking activity and continuous field potentials, is important across various neuroscience applications. However, a major challenge for doing so is that different neural modalities can have different timescales (i.e., sampling rates) and different probabilistic distributions, or can even be missing at some time-steps. Existing nonlinear models of multimodal neural activity do not address different timescales or missing samples across modalities. Further, some of these models do not allow for real-time decoding. Here, we develop a learning framework that can enable real-time recursive decoding while nonlinearly aggregating information across multiple modalities with different timescales and distributions and with missing samples. This framework consists of 1) a multiscale encoder that nonlinearly aggregates information after learning within-modality dynamics to handle different timescales and missing samples in real time, 2) a multiscale dynamical backbone that extracts multimodal temporal dynamics and enables real-time recursive decoding, and 3) modality-specific decoders to account for different probabilistic distributions across modalities. In both simulations and three distinct multiscale brain datasets, we show that our model can aggregate information across modalities with different timescales and distributions and missing samples to improve real-time target decoding. Further, our method outperforms various linear and nonlinear multimodal benchmarks in doing so.

Dynamical modeling of nonlinear latent factors in multiscale neural activity with real-time inference
Eray Erturk
Ming Hsieh Department of Electrical and Computer Engineering University Southern California Los Angeles, CA eerturkQusc.edu
Maryam M. Shanechi"
Ming Hsieh Department of Electrical and Computer Engineering
Thomas Lord Department of Computer Science
Alfred E. Mann Department of Biomedical Engineering
Neuroscience Graduate Program University of Southern California
Los Angeles, CA shanechi @du
Abstract
Real-time decoding of target variables from multiple simultaneously recorded neural time-series modalities, such as discrete spiking activity and continuous field potentials, important across various neuroscience applications. However, a major challenge for doing so is that different neural modalities can have different timescales (i.e., sampling rates) and different probabilistic distributions, or can even be missing at some time-steps. Existing nonlinear models of multimodal neural activity do not address different timescales or missing samples across modal- ities. Further, some of these models do not allow for real-time decoding. Here, we develop a learning framework that can enable real-time recursive decoding while nonlinearly aggregating information across multiple modalities with different timescales and distributions and with missing samples. This framework consists of 1) multiscale encoder that nonlinearly aggregates information after learning modainyyodllity dynamics to handle different timescales and missing samples real time, 2) multiscale dynamical backbone that extracts multimodal temporal dynamics and enables real-time recursive decoding, and 3) modality de- coders account for different probabilistic distributions across modalities. In both simulations and three distinct multiscale brain datasets, we show that our model can aggregate information across modalities with different timescales and distributions and missing samples to improve real-time target decoding. Further, our method outperforms various linear and nonlinear multimodal benchmarks doing so.
Introduction
Real-time continuous decoding of target time-series, such as various brain or behavioral states from neural tie--seriss data is of interest across many neuroscience applications. A popular approach for doing so to develop dynamical latent factor models that describe neural dynamics in terms
"Corresponding author: thanechi@usc.cer edu
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

of the temporal evolution of latent variables that can be used for downstream decoding. To date, dynamical latent factor models of neural data have mostly focused on a single modality of neural data, for example, either spiking activity or local field potentials (I.PP) [+ However, brain and behavioral target states are encoded across multiple spatial and temporal scales of brain activity that are measured with different neural modalities. Furthermore, some of these dynamical models have a non-causal inference procedure, which hinders real-time decoding. Therefore, inference of target variables could be improved by developing nonlinear dynamical models of multimodal neural time-series that can, at each time-step, aggregate information across neural modalities and do so in real-time.
natural challenge in developing such multimodal models arises when modalities are not aligned due to their different recording timescales that can be caused by various factors such as fundamental biological differences across modalities-witt some modalities evolving slower than others [ differences in recording devices 6 77, or measurement failures or interruptions [8]00. Further, modalities could have different distributions. For example, spiking activity is a binary-valued time- series that indicates the presence of action potential events from neurons at each time. As such, it has a fast millisecond timescale and is often modeled as count processes, such as Poisson. In comparison, LFP activity is a continuous-valuee modality that measures network-level neural processes, has a slower timescale, and is typically modeled with a Gaussian distribution [5] 7 [I]. We refer to multimodal data with different timescales multiscale data. Thus, to fuse information across spiking and LFP modalities and improve downstream target decoding tasks, their dynamics should be modeled by incorporating their cross-modality probabilistic and timescale differences through careful encoder design.
Existing neural dynamical modeling approaches do not address the nonlinear modeling of multimodal data with different timescales and/or with real-time decoding capability. Specifically, most dynamical models do not capture multimodal neural dynamics and instead focus a single modality of neural activity either by using finear/wwitching-Be
approaches [ [ or by utilizing nonlinear deep learning approaches [2 14). There are also some methods to jointly model unimodal neural activity together with behavior [15 23), but again, their latent factor inference is performed by processing unimodal neural time-series and does not aggregate multimodal neural data. While there has been some recent work on dynamical modeling and decoding of multimodal neural data, many of these works have been linear [24 Motivated by this gap, recent studies have developed nonlinear models of multimodal neural data [8] [1] However, these recent works have not addressed modalities with different timescales; further, the latent factor inference in such dynamical models has been done non-causally over time (1
Beyond neural time-series data, many approaches in other domains have been proposed to combine multiple modalities. However, these are again not focused on addressing the challenge of different timescales and missing samples over time, and their applicability to joint modeling of Poisson and Gaussian distributed modalities encountered in neuroscience has not been investigated (Section 2).
Contributions We introduce Multiscale Real-time Inference of Nonlinear Embeddings (MRINE). a nonlinear dynamical modeling approach designed to nonlinearly fuse multimodal neural time-series with different timescales, distinct distributions, and/or missing samples over time, while supporting inference both in real-time and non-causally. To achieve these capabilities, we: 1) design a multiscale encoder that performs nonlinear information fusion via neural networks after learning modality- specific dynamical models that account for timescale differences and missing samples in real-time by learning the temporal dynamics (Section 3.2) and 2) impose smoothness priors on the latent dynamics via smoothness regularization objectives that also prevent learning trivial identity neural network transformations (Section 3.3).
Through stochastic Lorenz attractor simulations and two independent nonhuman primate (NHP) spiking and LFP neural datasets, we show that MRINE infers latent factors that are more predictive of target variables such as the NHP's arm kinematics. Further, we compare MRINE with various recent linear and nonlinear multimodal methods and show that MRINE outperforms them in downstream decoding from multiscale spike-LFP time-series in the NHP datasets, as well as on a high-dimensional visual stimuli dataset containing neuropixel spikes and calcium imaging data.
2

---

2
Related work
Single-Scale Models of Neural Activity Numerous dynamical models of single-scale neural activity have been developed. Some of these models are in linear generalized linear form [1 12 15
32 34
Linear dynamical models (LDMs) are widely used in real-time applications because they provide real- time and recursive inference algorithms. To enable nonlinear modeling, there has been an increased interest in switching linear models [13] and deep learning architectures including recurrent neural architectures autoencoder-basse
networks (RNN) with nonlinear temporal dynamics 3 19 35
that utilize Markovian linear dynamics to learn a smoothing distribution 22 36., transformer
encoder based models optimized with masked training [37 38 and neural ordinary differential
equations [39]. These models have shown great promise in improving behavior decoding compared
to linear models (3]. A recent work [4] has also developed nonlinear neural network framework termed DFINE that supports flexible inference-i.e. enables both real-time filtering and noncausal smoothing, and accounts for missing observations simultaneously by jointly training an autoencoder
with linear state-space models and utilizing Kalman filtering 40 Another work has proposed a
low-rank structured variational encoding framework for Gaussian state-space models to capture dense covariances with predictive capabilities, while supporting real-time parameterization of their inference network [41). While LDMMbased approaches and some other nonlinear approaches [41) can handle missing samples- either via Kalman filtering or in a similar spirit within probabilistic state-space formulations- aal the above linear and nonlinear methods are designed for a single modality of neural activity and do not address multiscale modeling.
Multimodal Information Fusion Outside neuroscience, fusing multiple data modalities has been extensively researched across many areas including natural language processing (NLP) and computer vision. However, these approaches do not address nonlinear modeling of multiple time-series modalities with different timescales and with real-time inference capability, which we consider here. Specifically, in computer vision, many studies focused on variational autoencoders (VAE), and approximated the joint posterior distribution by factorization it into modality-specific posteriors 42.46 to handle missing modalities, by concatenating the modality-specifte representations 48 Instead of learning the common embedding space via factorization or concatenation, some studies have also employed cross- modality generation 49 50]. Even though some of these methods can fuse time-series modalities with different timescales, they need to do so using separate networks for each modality first non-causally encode each modality into a single vector, which means information fusion cannot be done causally over time-steps (i.e., in real time). However, many neuroscience applications require aggregating information at each time step in a causal manner to perform continuous real-time decoding of target variables. Finally, they are not designed to handle time-series signals with missing samples in time and would rely on data augmentations such as zero-padding, which can be suboptimal by changing the value of missing samples (51533)
Multimodal Models in Neuroscience One line of work in neuroscience aims to jointly model single-scale neural activity and behavior as multimodal signals [15H17] 19). However, in these works, latent factor inference is performed using single-scale neural activity without any information fusion, similar to the single-scale models discussed above. Another line of work aims to model multimodal neural signals. However, these methods are either linear/generalized- -areeiiiiidddleearr or are designed for offline reconstruction without addressing distinct timescales or enabling real-time inference. Specifically, some approaches proposed linear multimodal modeling frameworks [26 27 and learned the model parameters via further do not address different timescales. To partly address this gap. a linear dynamical modeling framework was introduced in 24 in which model parameters are also learned with EM but this time inference can aggregate modalities with different timescales. With a similar lincar formulation of multiscale dynamics, recent work in (28 proposed a more computationally efficient learning framework compared to [24] by using subspace identification. However, both of these approaches are still linear and cannot characterize nonlinearities.
(EM). But these methods do not capture nonlinearity and
To account for nonlinearity, recent studies have developed nonlinear latent factor models for multi- modal neural data [18] 29 31], but their formulation assumes the same timescale for the different modalities and also does not consider missing samples. Thus, such situations, they would need to rely on indirect approaches such as augmentations with zero-padding, which can be suboptimal by changing the value of missing samples 4 511 Furthermore, such multimodal dynamical models have non-causal inference networks and thus do not enable real-time inference of latent
3

---

factors [29 31] Instead, here we develop nonlinear dynamical modeling framework for multimodal neural time-series data that supports real-time and efficient recursive inference, and handles both timescale differences and missing samples by directly leveraging the learned dynamical model to predict these missing samples.
3 Methodology
We assume that we observe discrete neural signals (e.g., spikes) 8, E {0, 1) for T where
T= = [1,,,...., .T} and continuous neural signals (c.g., LFPs) yr E R for t' E T' where T TT. Note that the two different sets T and T' allow for the timescale differences of 8; and Br via different time-indices. As shown in Fig. T and expanded on below, we describe the neural processes generating 8 and ye through multiscale latent and embedding factors, which in turn can be extracted by nonlinearly aggregating information across these multiple neural modalities with multiscale encoder.
b.
stpurveion
Parson AESETAATIONS
Cousiion
Muitiscale IDM
ditribuzon parametoos
MIP
Marrcas
Clanddon
Figure I: a. MRINE model architecture. Multiscale encoder nonlinearly extracts multiscale embed- ding factors (a,) by fusing discrete Poisson and continuous Gaussian neural time-series in real-time (as indicated by thick dashed lines) while accounting for timescale differences and missing samples. Temporal dynamics
then explained with multiscale LDM whose states are the multiscale latent factors (I))). As an example, when y/ is missing, ag and 2 are inferred by processing only $2. b. MRINE multiscale encoder design. Modality-specinff
LDMs learn within-modality dynamics and account for timescale differences or missing samples via Kalman filtering. Then, filtered modality- specific embedding factors (af and . are fused and processed by another fusion network to obtain the multiscale embedding factors aç. As an example, when y/2 is missing, 22/2 is predicted only by the dynamics of modality-specific LDM in Eq. 7 Note that MLP networks are time- inve-invariant, and information about temporal dynamics is incorporated through LDMss
3.1 Model formulation
The model architecture is shown in Fig. La with its main components being a multiscale encoder network for nonlinear information fusion, a multiscale LDM backbone facilitate real-time inference, and nonlinear decoder networks for each modality. We write the generative model as follows:
= Az, + w.
= Cx, + rt
(2)
Po. (8) a,) = (ais A(a:)),
(3)
where A(.) == (.)
pe, (yy ar) = N (y (.
(4)
where u(-) fo, (-)
Here, Eq. 1 and 2 form multiscale LDM, and / E T where T (1,2,.... T} and l' T' where T' C T. this I.DM, a, E RR",, termed multiscale embedding factors, are the observations. We obtain a₂ as the nonlinear aggregation of multimodal information through the multiscale encoder

---

designed in Section 3.2 Further, H, E R"x, termed multiscale latent factors, are the LDM states and model the dynamics as linear in the nonlinear embedding space, which facilitates real-time inference (Section 3.2 Note that from this point onward, we will only use time-index & (except observation models) for embedding and latent factors as T T. Correspondingly.
E IR x and C' E x# are the state transition and observation matrices in the multiscale LDM, respectively; we E IR= and r; E RT are zero-mean Gaussian dynamics and observation noises with covariance matrices W E 1 and R E R Xn..
The observations from different modalities 8 and ye are then generated from at (or ar) as in Eqs. 3 and respectively, with the likelihood distributions denoted by Pe, (st at) and pe, (Ve ar).. Assummng that data modalities are conditionally independent given a, (or ay for yrr). we modeled the discrete spiking activity 8, with a Poisson distribution and the continuous neural signals y with a Gaussian distribution, given that these distributions have shown success in modeling each of these
modalities 2 3 24 54) The means of the corresponding distributions X(-) := So, (.) E R and (x(+) == fo, (+) E RM are parametrized by neural networks with parameters 0g and 8y. Practically, we observed that learning the variance of the Gaussian likelihood yielded suboptimal performance,
thus we set to a constant value, i.e., unit variance, as in previous works [55] 56).
3.2 Encoder design and inferring multiscale factors
To infer a and x, from s and ye we first construct the mapping from 8, and yı to a as:
at = =f ($t. yrr
(5)
where fo(:) represents the multiscale encoder network parametrized by a neural network with parameters 6.
One obstacle in the design of the encoder network is accounting for the different timescales without using augmentation techniques such as zero-padding as commonly done [53] (57) since they can yield suboptimal performance and distort the information during latent factor inference [51 52 58). Thus, it is important to account for timescale differences and the possibility of missing samples when designing the multiscale encoder. Additionally, our goal is to perform multimodal fusion at each time-step while also allowing for real-time inference of factors. We address these problems, which remain unresolved in prior methods, with a multiscale encoder network design shown in Fig. Lpp.
In our multiscale encoder (Fig. [b). at each time-step, first, each modality (s; and ye)) is processed by separate time-invariant multilayer perception (MLP) networks with parameters d, and Øy to obtain modality-speciff embedding factors, a E Rh= and af E R"=,, respectively. Then, we construct modality-specifie LDMs for each modality, whose observations each time-step are their corresponding embedding factors:
2 Asæ; + wi ai C,æi + rē #x A,= w
(6)
(7)
af Cyat + re
where w;; x E RR are modality-specific latent factors, As, Ay E R XNa are the state transition matrices, C. Cy E Reaxn, are the emission matrices, wT, w. E Rh the zero-mean dynamics noises with covariances W a. W, y E R XMa, and rū and r are the zero-mean Gaussian observation noises with covariances R,, Ry IR X for modalities 8; and yy, respectively. We denote the modality-specific LDM parameters by sb. = [A,, C., W. o, R,} and t; = ( Cy, W, y, R for 8; and ye.. respectively.
In our design, we use the sppelit LDMs because they allow us to account for missing samples (whether due to timescale differences missed measurements) by using the learned within- modality state dynamics to predict these samples forward in time, while also maintaining the operation fully real-time/causal. Specifically, given the modality-speciffe LDMs in Eq, and the modality-speciff latent factors, a and x with Kalman filtering, which is real-time and constitutes the optimal minimum mean-squared error estimator for these models [40]. We use the subscript tlj to denote the factors inferred at time i given all observations up to time j. As such, subscripts t/t, t/T and t + klft denote causal/real-time filtering, non-cause smoothing and k-step-ahead prediction, respectively. At this stage, if / an intermittent time-step such that ye is missing (i.e., T and t. & T ), a is obtained with forward prediction using the Kalman
can obtain
5

---

predictor as 2x1 = Ay 159, and similarly for Tij: if s; is missing. Having done this for each modality, we perform information fusion in real-time by concatenating the modality embedding factors and passing them through a fusion network with parameters @m to obtain the initial representation for as, which later becomes the noisy observations of the multiscale LDM formed by Eq. T and 2 (also see Fig. Ta). We denote the learnable multiscale encoder network parameters by { Dy: thss Vy: Om
The multiscale LDM now allows us to infer Wejt with real-time (causal) Kalman filtering, or infer FefT with non-causal Kalman smoothing [60]. Similarly, k-step-ahead predicted multiscale latent factors R-t can be obtained by forward propagating 2 k-times into the future with Eq. [ i.e.,
Te+kkte A A- We denote the parameters of the multiscale LDM by tm {A, C, W, R). We can now obtain the filtered, smoothed, and k-step-ahead predicted parameters of the likelihood functions in Eq.
3 and 4nd by first using Eq. to compute the corresponding multiscale embedding factors-i.e. 0gj Caa where ilj is t/t, 77, and + kijt, respectively- and then forward passing these factors through fo, (+) and fo, (.) (Fig. I.).
3.3 Learning the model parameters
k-step ahead prediction To learn the MRINE model parameters and in order to encourage learning the dynamics, part of the loss, we employ the multi-horizon k-step-ahead prediction loss defined as:
L
Σ AEK
Σ
Flg (po, (e]=-))] + Σ
log (ye, dele-=)))
(8)
F'T 122
where T and T' denote the time-steps when 8 and ye are observed, respectively. 7 is the scaling parameter the log-likelihood values of different modalities are of different scales (see Appendix A.2.2). and IC is the set of future prediction horizons. We note that k-step-ahead prediction is performed by computing k-step-ahead predicted multiscale latent factors, k2t rather than modality- specific ones.
Smoothed reconstruction In addition to the k-step-ahead prediction, we also optimize the recon- struction from smoothed multiscale factors:
Symooth
log (s ((r)) +
Σ V"T
log (pe, arrr)))
where 7' is the last time-step that any modality is observed.
Smoothness regularization To impose smoothness prior learned dynamics and to prevent the model from overfitting by learning trivial identity encoder/decoder transformations, in our loss, we also apply smoothness regularization on po. (s: a1:1), Poy (3 a1:r) and p(æ; (1:r) by minimizing the KL-divergence between the distributions in consecutive time-steps as introduced in [61] for Gaussian-distriuutedd modalities. Here, extend this technique also Poisson-distriuuted modalities as below. Let L.m be the smoothness regularization penalty, defined as:
IT By
Cxm
(sl Pa. (aT...I)) lat...a.) 17g
d Pr. (vrr (e7p)) Po, WT lar,..m))
Smootthness on d (plar agy) plal (4417))
Smocthness
(10)
+ 7x
Smonthneas on
where -(.,-) the KL-divergence between given distributions, subscript i denotes the 7th element of the set, superscript j is the 5th component of the vector (e.g., 8T,), and po, (s, a.jr) and pa, (yrr (rrr) are as in Eq. []ddd]] respectively. Here 7s. 7y and 7ar are the scaling hyperparameters. The smoothness penalties on 8, and yı are computed over the time-steps that they are observed. The penalty on D; is obtained over all time-steps 2; is inferred for all time-steps. After extracting a1r with multiscale encoder, p(x, d1.T) can be obtained with the Kalman smoother, which provides the posterior distribution for the multiscale LDM [40 62]:
p(æ; (1.T) N (W Egpr
(1I)

---

where EijT and Enj are the smoothed multiscale latent factors and their error covariances, respec- tively. To allow the model learn both fast and slow dynamics, we put the smoothness regularization on 2e on half of its dimensions.
To assess the impact of incorporating smoothness regularization terms in Eq. 10 and smoothed reconstruction in Eq. | we performed an ablation study (see Appendix[A.6.2]) which demonstrates that each term contributes the improved performance. Further, in another ablation study on the effect of multiscale modeling (see Appendix [A.64). we show that MRINE's multiscale encoder design is an important contributing factor improved performance, compared the case where missing samples are imputed by zeros and removed from the training objectives above.
Finally, we form the loss as the sum of the above elements and regularization terms, and minimize via mini-batch gradient descent using the Adam optimizer [63 to learn the model parameters {db vh 0z 0y):
MMRINE = Lk Lsmooth + Lsm + 7r (0 By. Os Oy. Pm)
(12)
where (2))) is the L2 regularization penalty on the MLP weights and "r is the scaling hyperparameter.
Moreover, we employ a dropout technique termed time-dropout during training to increase the robustness of MRINE to missing samples even further. See Appendix for more information, Appendix[A.6.. for an ablation study on the effect of time-dropout, and Appendix A.2 for training details and hyperparameters.
4 Results
4.1 Stochastic Lorenz attractor simulations
We first validated that MRINE can successfully aggregate information across multiple modalities by performing simulations with the stochastic Lorenz attractor dynamics defined in Eq. T8. To do so, we generated Poisson and Gaussian observations with 5, 10 and 20 dimensions as described in Appendix A.4.2 We generated systems with different random seeds and performed 5-fold cross- validation for each system. Then, we trained MRINE as well as single-scale networks with either only Gaussian observations only Poisson observations (see Appendix [A.2.T). To assess MRINE's ability to aggregate multimodal information, we compared its latent reconstruction accuracies with those of single-scale networks. For each model, these accuracies were obtained by computing the average correlation coefficient (CC) between the true and reconstructed latent factors (see Appendix A.4.3 for details). We refer to each observation dimension channel for simplicity.
We first considered the Poisson modal- ity as the primary modality and fused with it 5, 10, and 20 Gaussian channels as the secondary modality. As shown in Fig. (2h, for all different numbers of Pois- son channels, gradually fusing Gaus- sian channels significantly improved the latent reconstruction accuracies (p < 10-5, 20, one-sided Wilcoxon signed-rank test). As expected, these improvements were higher in the low in- formation regime where fewer primary Poisson channels were available (5 and 10) compared to the high information regime (20 Poisson channels). Likewise, when the Gaussian modality was the primary modality, latent reconstruction accuracies showed consistent improve- ment across all regimes as Poisson channels were fused (Fig. 2.. < 10-5, n = 20, one-sided Wilcoxon signed-rank test).
a.
1.00
aaomrooomtteett accuacy
0.96
200
Poisson channels
5 Gaues channels
Posson channess
10 Sauss channels
0.50
20 Poessan dhannet male
20 Dausseee sooee chemels
Number Gausson channels
Number Pessean charmele
Figure 2: Latent reconstruction accuracies for the stochas- tie Lorenz attractor simulations. a. Accuracies when 5, 10, or 20 Poisson channels were the primary modality to which Gaussian channels were gradually added. Solid lines show the mean and shaded areas show the standard error of the mean (SEM). b. Similar to a, but with the Gaussian channels being the primary modality.
7

---

4.2 MRINE fused multiscale information in behavior decoding for the NHP datasets
To test MRINE's information aggregation capabilities in real datasets, we used two independent publicly available NHP datasets (6466), For both datasets, discrete spiking activity and continuous LFP signals were recorded while subjects performed 2D reaches either on a grid via a cursor (grid- reaching dataset), from center target random outer target via manipulandum (center-out reaching dataset). We considered the 2D cursor and manipulandum velocities in x and y directions as our target behavior variables (see Appendices A.5.1 and A.5.2 for details), For both datasets, we trained single-scale models with 5, 10, and 20 channels of either spike or LFP signals alone, as well as MRINE models of spike-LFP signals for every combination of these multimodal channel sets, and decoded the target behavior from inferred latent factors (see Appendix [A.333] for details). In our analyses, we used and 3 experimental sessions recorded on different days for the NHP grid reaching
and center-out reaching datasets, respectively, and performed 5-fold cross-validation for each session. In our analyses, LFP and spike signals had different timescales: spikes were observed every 10 ms
5 Stoine chaernea Spito shannels Sphvi charees
L117 ceeeee LFP charnels charnels
00000 00 20
Scharic dood
40
5 Spak cha mecbs
5 LFP channsls 10 FP diannes
344
Sp chemnes
[FP channee
0.1
hart box of LPP charmeda
Number of chltechannee
deouSho SEDURETy (5, (1- 31
Inpeaament Sehutee
Besener ducurta
0.
3
23
Number LFP LFP charnsta
N.vre of Spiia statrols
Figure 3: Behavior decoding accuracies for the NHP grid reaching dataset (top row) and NHP center- reaching dataset (bottom row). Shaded areas and error bars represent SEM. a. Accuracies, when 10, or 20 spike channels were the primary modality and an increasing number of LFP channels, were fused. b. Similar when LFP channels were the primary modality. Percentage improvements in decoding accuracy when 20 LFP channels were added to 5, 10, and 20 spike channels. Asterisks indicate significance of comparison (***; p < 0.0005, one-sided Wilcoxon signed-rank test). d. Similar to c, when 20 spike channels were added to 5, 10, and 20 LFP channels. e-h. Same as a-d but for NHP center-out reaching dataset.
and LFPs every 50 ms (abbreviated as 50 ms LFPs) 6 28). We found that when spiking signals were taken as the primary modality, fusing them with increasing numbers of LFP channels improved the behavior decoding accuracy for all different numbers of primary spike channels (Fig. 33 for NHP grid reaching dataset, < 0.009, = 20; Fig. 3 for NHP center-out reaching dataset, p < 10 = 15, one-sided Wilcoxon signed-rank test). Similar to simulation results, improvements in behavior decoding accuracy were higher in the low-information regimes, i.e., for 5 and 10 primary spike channels. For instance, adding 20 LFP channels to 5, 10, and 20 spike channels improved behavior decoding accuracy by 14.7%. 10.4%, and 5.3%, respectively, for the NHP grid reaching dataset (Fig. Be), and 47.3%, 24. 2..4%. and 9.8% for the NHP center-out reaching dataset (Fig. 3
When LFP signals were considered as the primary modality, behavior decoding accuracies again improved significantly when spike channels were fused (Fig. 3p for grid reaching dataset, p < 10 for 20: Fig. 3t for center-out reaching dataset, < 10 n 15, one-sided Wilcoxon signed-rank test). For example, adding 20 spike channels 5, 10 and 20 LFP channels improved behavior decoding accuracy by 49.8%, 34%, and 29.6% in the grid reaching dataset (Figs. Bdd and by 30, 21.6%, and 20.6% in the center-out reaching dataset (Fig. 3h). Overall, the bidirectional improvements- eengrrdl of whether spikes or L.PPs were the primary modality- for both datasets suggest that spike and LFP modalities may encode non-redundant information about behavior variables that can be fused nonlincarly with MRINE. Interestingly, adding spike channels to LFP
8

---

channels helped more than adding LFP channels to spike channels in the grid grid-reaching dataset, but the reverse held in the center-out reaching dataset. This suggests that while these modalities may have different amounts of behaviorally relevant information depending on the experiment (e.g., task or electrode types/locations), MRINE can combine their complementary strengths and be robust to such differences. Finally, we also tested MRINE when both modalities had the same timescale (10 ms) and found that MRINE again enabled consistent bidirectional improvements in behavior decoding performance for both datasets (Fig. 4h.
4.3 MRINE improved behavior decoding compared with prior multimodal modeling methods
We compared MRINE with prior multimodal methods, namely the recent MSID in 28 mmmLRNN in [29]. MMGPVAE in BT, and MVAE in [4 for both NHP datasets. We trained all models with 5, 10, or 20 channels of 10 ms spikes and 50 ms LFPs. Since mmPLRNN and MMGPVAE do not support training and inference with different timescale signals, we imputed the missing LFP timesteps by their global mean, i.e., zeros due to z-scoring, as is common practice. Further, mmPLRNI and MMGPVAE decodings were performed non-causally unlike MRINE and MSID since mmPLRNN and MMGPVAE do not support real-time recursive inference. For all numbers of primary channels (i.c., all information regimes), MRINE significantly outperformed all baselines in behavior decoding (Table 0 p 10-5, n 20 for NHP grid reaching dataset; p 10-5, 15 for NHP center-out reaching dataset, one-sided Wilcoxon signed-rank test).
NIIP grid reaching
NIIP center-out reaching
Spike
Spike
20 Spike
Spike
10 Spike
20 Spike
Method
LLFP
10
2.PF
L
10
LFP
MVAE
0.326
0.011
0,386
0.009
0.425 0.009
0.392:4000
0,462 0.014
0.544 0.018
MSID
0.380
0.021
0,440
1 0.015
0.519 0.012
0,442
0.022
0.521 1 0.023
0.561 0.020
mmPLRNNN
0.450
0.010
0.477
0.013
0.540 0.011
0.468
0.022
0.470 + 0.029
0.538
0.032
MMGIVVEE
0,276
0.021
0,438 0.019
0.479 0.017
0.495
0.020
0.567
0.020
0.601
0.021
MRINE
0.487 0.007
0.555 0.011
0.611 = 0.012
0.547 + 0.020
0.624 = 0,022
0.649 = 0.021
Table 1: Behavior decoding accuracies for the NHP grid reaching and center-out reaching datasets with 5, 10, and 20 channels of 10 ms spikes and 50 ms LFP for MVAE, MSID, mmPLRNN, MMGPVAE, and MRINE. The pest-performing method is in bold, the second best-performing method is underlined, + represents SEM.
Next, using ablation studies on MRINE, we found that our multiscale encoder design was key in achieving higher performance than baselines. particular, accounting for different timescales using the common technique of zero- cerotimmutation instead of using our multiscale encoder significantly degraded performance (Appendix(AAe) (A.6.4). Also, we performed the same baseline comparisons for the same timescale scenario (10 ms) and found that MRINE again outperforms all baseline methods for both datasets (Table 77. The degradation in performance due to having different timescales vs. having the same timescale was lower in MRINE than in MMGPVAE and mmPLRNN, again highlighting the importance of MRINE's design, especially when modalities have different timescales.
Finally, we trained single-scale models, MSID, MMGPVAE, and MRINE models on high-channel information regimes beyond 20 channels of spike and LFP signals, i.c., 30 and 60 channels. As shown in Table 8 MRINE again outperformed all methods. highlighting that MRINE's s aggrega- tion capabilities extend effectively to nigher-dimensional recordings and that MRINE maintains performance advantages even as the information content increases. Please also see Appendix (ppendix(A.3.9) for trial- averaged latent factor visualizations, Appendix [ for more details on MVAE, MSID, mmPLRNN, and MMGPVAE benchmarks, and Table instead of CC, showing that MRINE's improved performance consistent across metrics.
9 for decoding accuracies using the R22 metric
4.4 MRINE'S information aggregation was robust to missing samples
Next, we studied the robustness of all methods' inferences to randomly missing samples. Here, in addition to having different timescales, we used various sample-dropping probabilities to drop spike or LFP samples in both datasets. In both datasets, for models that were trained with 20 channels of 10 ms spikes and 50 ms LFP, we either fixed the dropping probability for LFP samples as 0.2 while varying that of spike samples (Fig. (p.c) or vice versa (Fig. (h.d) during inference. For the NHP grid reaching dataset, behavior decoding accuracies of MRINE remained robust, for example decreasing

---

by only 5.4% and 20.4% when 40% and 80% of spike samples were missing, respectively, in addition to 20% of LFP samples missing (Fig. 5p. Also, MRINE outperformed all baseline methods across all sample dropping regimes (Fig. bb)) Similarly for the NHP center-out reaching dataset, MRINE again outperformed all baseline methods across all sample dropping regimes (Figs. Fp.d).
In addition to behavior decoding, we also evaluated all methods' information aggregation capabilities in reconstructing spike and LFP signals in the same missing sample scenarios as for the previous case. As shown in Table TO MRINE achieved competitive performance in reconstructing neural modalities as it achieved the best total average rank in terms of how well reconstructed the neural modalities compared to baseline methods (see Appendix 8.5.8 for details; averaged over the sample dropping probability regimes, modalities, and datasets in Fig. 6 This shows MRINE's capabilittes beyond behavior decoding for neural reconstruction.
4.5 MRINE's performance improvements generalized to distinct high-dimensional multimodal neural dataset with visual stimuli
To provide evidence for MRINE'S generalization beyond neural datasets during motor tasks, we evaluated MRINE on a multimodal high-dimensional (i.e., 800-D) neural dataset that contains neuropixel spiking activity and calcium imaging data recorded from the visual cortex during visual stimuli presentation [67 [68] (see Appendix AA.5.10 for details). This dataset contained neuropixel and calcium imaging data with different timescales (i.e., sampled at 120 and 30 Hz, respectively), which were subsequently smoothed with causal Gaussian kernels with standard deviations of 8 ms. In this analysis, we also included CEBRA 17 in our evaluations (see Appendix A.2.8 for CEBRA training details), as well as single-scale models trained on either modality. As the downstream task, we used the frame ID decoding task (0-900, 30Hz, 30s movie), similar [17.. Since the downstream task was of slower frequency (30 Hz rather than 120 Hz), we performed average pooling on the latent factors inferred by MRINE for downsampling. As shown in Table 11 MRINE successfully aggregated information across neuropixel spiking activity and calcium imaging modalities, which had different timescales, and improved frame IDD prediction performance over single-scale models. Furthermore, MRINE outperformed the CEBRA baseline, potentially due to its explicit dynamical backbone that can capture temporal dependencies in sequential time-series data, whereas CEBRA's convolutional architecture does not explicitly model state evolution.
5 Discussion
In this study, we presented MRINE that can dynamically, nonlinearly and in real-time aggregate information across multiple time-series modalities with different timescales or even with missing samples. To achieve this, we proposed a novel multiscale encoder design that first extracts modality- specific representations in real-time while accounting for their timescale differences and missing samples, and then performs nonlinear fusion in real-time to aggregate multimodal information. Through stochastic Lorenz attractor simulations and real NHP datasets. we show MRINE's ability to fuse information across modalities with different (or the same) timescales and with missing samples to achieve better downstream target decoding performance. We show that MRINE outperforms recent linear and nonlinear multimodal methods. Further, using these comparisons in addition to several ablation studies, we show the importance of MRINE's multiscale nonlinear encoder design and training objective outlined in Section 3.3 to enable accurate real-time -time fusion for multiscale data. A current limitation of MRINE is the assumption of time-invariant multiscale dynamics, which may not hold in non-stationary cases, and MRINE models may need to be intermittently retrained across days/sessions. Extending MRINE to track temporal variability, such with switching dynamics or adaptive approaches, is an important future direction [13] (39773. Another promising direction is to incorporate automatic hyperparameter tuning methods, since MRINE introduces additional scaling terms for smoothness regularization penalties that we showed are important contributing factors to its improved downstream decoding performance. A further direction is explore alternative training objectives for accurate cross- moosl generation, which could allow MRINE and other baseline multimodal methods to reconstruct an entirely missing modality from the available one. Finally, incorporating external inputs in the MRINE model to disentangle input dynamics from intrinsic neural dynamics, as done in prior single-scale models 20 21, is an interesting direction. Overall, MRINE can be especially important for applications where target variables are encoded across time-series modalities with different timescales and distribution, and/or where decoding needs to be done in real-time, for example, in brain-computer interfaces 764.
10

---

Acknowledgements
This work was partly supported by the National Institutes of Health (NIH) grants RFIDA056402. R01MH123770, and R61MH135407, the NIH Director's New Innovator Award DP2-MHI26377 and the National Science Foundation (NSF) CRCNS program award I1S-2113271. We sincerely thank all members of the Shanechi lab for helpful discussions.
References
[1] Jakob H Macke, Lars Buesing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking neural populations. In Advances in Neural Information Processing Systems, volume 24, pages 1350-1355 Curran Associates, Inc., 2011 1.
[2] Yuanjun Gao, Evan W Archer, Liam Paninski, and John P Cunningham. Linear dynamical neural population models through nonlinear embeddings. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, vol- ume 29, pages 171. Curran Associates, Inc., 2016.
[3] Chethan Pandarinath, Daniel J. O' Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew Kaufman, Stephen Ryu, Leigh R. Hochberg, Jaimie M. Henderson, Krishna V. Shenoy, L. F. Abbott, and David Sussillo. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10))805-815 October 2018.
[4] Hamidreza A Abbaspourazad, Eray Erturk, Bijan Pesaran, and Maryam M. Shanechi. Dynamical flexible inference of nonlinear latent factors and structures neural population activity, Nature Biomedical Engineering, 8(1):85-108, January 2024. Publisher: Nature Publishing Group.
[5] Bijan Pesaran, Martin Vinck, Gaute T. Einevoll, Anton Sirota, Pascal Fries, Markus Siegel, Wilson Truccolo, Charles E. Schroeder, and Ramesh Srinivasan. Investigating large-scale brain dynamics using field potential recordings: analysis and interpretation. Nature Neuroscience, 21(7):903-919, July 2018. Publisher: Nature Publishing Group.
[6] Han-Lin Hsieh, Yan T. Wong, Bijan Pesaran, and Maryam M. Shanechi. Multiscale modeling and decoding algorithms for spike-field field activity. Journal of Neural Engineering, 16(1):00))))))) (6(1):016018, December 2018.
[7] Hung-Yun Yun Lu, Elizabeth S Lorene, Hanlin Zhu, Justin Kilmarx, James Sulzer, Chong Xie, Philippe Tobler, Andrew J Watrous, Amy L Orsborn, Jarrod tha R Santacruz. Multi-scale neural decoding and analysis. Journal of neural engineering, 18(4):10.1088/ 1088/1741 25552acc 16ff August 2021.
and Saman-
[8] Schalk Burger, Searle Silverman, and Gary Vuuren. Deriving correlation matrices for missing financial time-series data. International Journal of Economics and Finance, 10:105, 09 2018.
[9] Steven Cheng-Xias Xian Li and Benjamin Marlin. Learning from irregulatly-sampled time series:
A missing data perspective. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 5937-5946. PMLR, 13-18 Jul 2020.
[10] Michael Berger, Naubahar Shahryar Agha, and Alexander Gail. Wireless recording from unrestrained monkeys reveals motor goal encoding beyond immediate reach frontoparictal cortex. eLiff,, 9:e51322, May 2020.
[11] Gaute T. Einevoll, Christoph Kayser, Nikos K. Logothetis, and Stefano Panzeri. Modelling and analysis of local field potentials for studying the function of cortical circuits. Nature Reviews Neuroscience, 14(11):770-785 November 2013.
[12] Tze Hui Koh, William E. Bishop, Takashi Kawashima, Brian B. Jeon, Ranjani Srinivasan, Yu Mu, Ziqiang Wei, Sandra J. Kuhlman, Misha B. Ahrens, Steven M. Chase, and Byron M. Yu. Dimensionality reduction of calcium-imaged neuronal population activity. Nature computational science, 3(1):31-1))) 1-885, January 2023.

---

[13] Scott Linderman, Matthew Johnson, Andrew Miller. Ryan Adams, David Blei, and Liam Paninski. Bayesian Learning and Inference Recurrent Switching Linear Dynamical Systems. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, pages 912-922. PMLR, April 2017.
[14] Evan Archer, II Memming Park, Lars Buesing, John Cunningham, and Liam Paninski. Black box variational inference for state space models, November 2015.
[15] Omid G. Sani, Hamidreza Abbaspourazad, Yan T. Wong, Bijan Pesaran, and Maryam M. Shanechi. Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification. Nature Neuroscience, 24(1): 1e01144049999 January 2021.
[16] Cole Hurwitz, Akash Srivastava, Kai Xu, Justin Jude, Matthew Perich, Lee Miller, and Matthias Hennig, Targeted Neural Dynamical Modeling. In Advances in Neural Information Processing Systems, volume 34, pages 19379-293322 Curran Associates, Inc., 2021.
[17] Steffen Schneider, Jin Hwa Lee, and Mackenzie Weygandt Mathis. Learnable latent embeddings for joint behavioural and neural analysis. Nature, 617-7960, May 2023. Publisher: Nature Publishing Group.
[18) Ding Zhou and Xue-Xin Wei. Learning identifiable and interpretable latent models of high- dimensional neural activity using pi-vae. In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS '20, Red Hook, NY, USA, 2020. Curran Associates Inc.
[19] Omid G. Sani, Bijan Pesaran, and Maryam M. Shanechi. Dissociative and prioritized modeling of behaviorally relevant neural dynamics using recurrent neural networks. Nature Neuroscience, 27(10):2033-2455 October 2024. Publisher: Nature Publishing Group.
[20] Parsa Vahidi, Omid G. Sani, and Maryam M. Shanechi. Modeling and dissociation of intrinsic and input-driven neural population dynamics underlying behavior. Proceedings of the National Academy of Sciences, 121(7):e22288871 1221 2024.
[21] Parsa Vahidi, Omid G. Sani, and Maryam Shanechi. BRAID: Input-driven nonlinear dynamical modeling of berrvoora data. In The Thirteenth International Conference on Learning Representations, 2025.
[22] Sayed Mohammad Hosseini and Maryam Shanechi. Dynamical modeling of behaviorally rele- vant spatiotemporal patterns in neural imaging data. In Forty-second International Conference on Machine Learning. 2025.
[23] Lucine L Oganesian, Omid G. Sani, and Maryam Shanechi. Spectral learning of shared dynamics between generalized-linee processes. In The Thirty-eighth Annual Conference Neural Information Processing Systems, 2024.
[24] Hamidreza Abbaspourazad, Mahdi Choudhury, Yan T. Wong, Bijan Pesaran, and Maryam M. Shanechi. Multiscale low-dimensional motor cortical state dynamics predict naturalistic reach- and-grasp behavior. Nature Communications, 12(1):607, January 2021.
[25] Hamidreza Abbaspourazad, Han-Lin Hsieh, and Maryam M. Shanechi. A multiscale dynamical modeling and identification framework for spike-field activity. IEEE Transuctions on Neural Systems and Rehabilitation Engineering, 27(6): 128-11288111 138, 2019.
[26] Mohammad R. Rezaei, Haseul Jeoung, Ayda Gharamani, Utpal Saha, Venkat Bhat, Milos R. Popovic, Ali Yousefi, Robert Chen, and Milad Lankarany. Inferring cognitive state underlying conflict choices in verbal Stroop task using heterogeneous input discriminative-geeeraai decoder model. Journal of Neural Engineering, 20(5):050016,, September 2023. Publisher: IOP Publishing.
[27] Todd P. Coleman, Marianna Yanike, Wendy A. Suzuki, and Emery N. Brown. 3a mixed-filte algorithm for dynamically tracking learning from multiple behavioral and neurophysiological measures. In The Dynamic Brain: An Exploration of Neuronal Variability and Its Functional Significance. Oxford University Press, 01 2011.
12

---

[28] Parima Ahmadipour, Omid G Sani, Bijan Pesaran, and Maryam M Shanechi. Multimodal spiking and field potential population activity. Journal of Neural Engineering, 21(2):26001, March 2024. Publisher: IOP Publishing.
subspace identification for modeling
[29] Daniel Kramer, Philine Bommer, Carlo Tombolini, Georgia Koppe, and Daniel Durstewitz. Reconstructing nonlinear dynamical systems from multi-modal time series. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, edi- tors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 11613-11633, PMLR, 17-23 Jul 2022.
[30] Manuel Brenner, Georgia Koppe, and Daniel Durstewitz. Multimodal Teacher Forcing for Reconstructing Nonlinear Dynamical Systems, December 2022.
[31] Rabia Gondur, Usama Bin Sikandar, Evan Schaffer, Mikio Christian Aoi, and Stephen L., Keeley. Multi-modal Gaussian Process Variational Autoencoders for Neural and Behavioral Data. October 2023.
[32] Lars Buesing, Jakob H Macke, and Maneesh Sahani. Spectral learning of linear dynamics from generalised- linear observations with application to neural population data. In Advances in Neural Information Processing Systems, volume 25, pages 1682-1690. Curran Associates, Inc., 2012.
[33] John P. Cunningham and Byron M. Yu. Dimensionality reduction for large-scale neural recordings. Nature Neuroscience, 17(11):1500-1599 November 2014.
[34] Jonathan C. Kao, Paul Nuyujukian, Stephen I. Ryu, Mark M. Churchland, John P. Cunningham, and Krishna V. Shenoy. Single-trial dynamics of motor cortex and their applications to brain- machine interfaces. Nature Communications, 6(1):7759, July 2015.
[35] Qi She and Anqi Wu. Neural Dynamics Discovery via Gaussian Process Recurrent Neural Networks. In Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, pages 454-464. PMLR, August 2020.
[36] Daniel Durstewitz. A state space approach for piecewise-linear recurrent neural networks for identifying computational dynamics from neural measurements, PLOS Computational Biology, 13(6);1-33, 06 2017.
[37] Joel Ye and Chethan Pandarinath. Representation learning for neural population activity with Neural Data Transformers, August 2021.
[38] Trung Le and Eli Shlizerman. Stndt: Modeling neural population activity with spatiotemporal
transformers. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 17926-17939, Curran Associates, Inc., 2022.
[39] Timothy D. Kim, Thomas Z. Luo, Jonathan W. Pillow, and Carlos D. Brody, Inferring Latent Dynamics Underlying Neural Population Activity via Neural Differential Equations. In Pro- ceedings of the 38th International Conference on Machine Learning, pages 5551-5561. PMLR, July 2021.
[40] R. F. Kalman. A New Approach to Linear Filtering and Prediction Problems. Journal of Basic Engineering, 82(1):35-45 45, March 1960.
[41] Matthew Dowling, Yuan Zhao, and II Memming Park. eXponential FAmily Dynamical Systems (XFADS): Large-scale nonlincar Gaussian state-space modeling. November 2024.
[42] Mike Wu and Noah Goodman. Multimodal Generative Models for Scalable Weakly-Supervist Learning. arXiv:1802.05333 [cx, stat], November 2018.
[43] Yuge Shi, Siddharth N, Brooks Paige, and Philip Torr. Variational Mixture-of-Expers Au- wonccoders for Multi-Modal Deep Generative Models. In Advances Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
13

---

[44] Thomas M. Sutter, Imant Daunhawer, and Julia E. Vogt. Generalized Multimodal ELBO. arXiv:2105.02470 [es, stat], June 2021.
[45] Yao-Hung Hubert Tsai, Paul Pu Liang, Amir Zadeh, Louis-Philippe Morency, and Ruslan Salakhutdinov. Learning Factorized Multimodal Representations. 17Xfv:1886.06117 [cs, stat],
May 2019.
[46] Mihee Lee and Vladimir Pavlovic. Private-Shared Disentangled Multimodal VAE for Learning of Latent Representations. pages 1692-1700, 2021.
[47] Masahiro Suzuki, Kotaro Nakayama, and Yutaka Matsuo. Joint Multimodal Learning with Deep Generative Models, November 2016.
[48] Ramakrishna Vedantam, Ian Fischer, Jonathan Huang, and Kevin Murphy. Generative Models of Visually Grounded Imagination. 17Xiv:1705.10700 [cs, stat], November 2018.
[49] Tom Joy, Yuge Shi, Philip H. S. Torr, Tom Rainforth, Sebastian M. Schmon, and N. Siddharth. Learning Multimodal VAEs through Mutual Supervision. arXiv:2106.12570 [cs]. July 2021.
[50] Gaurav Pandey and Ambedkar Dukkipati. Variational methods for conditional multimodal deep learning. In 2017 International Joint Conference on Neural Networks (IJCNN), pages 308-315, May 2017.
[51] Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recur- rent Neural Networks for Multivariate Time Series with Missing Values. Scientific Reports, 8(1):6085, April 2018.
[52] Brian Wells, Kevin M Chagin, Amy S Nowacki, and Michael W Kattan. Strategies for handling missing data in electronic health record derived data. EGEMS (Washington, DC), 103): 2013.
[53] Feng Zhu, Andrew Sedler, Harrison A Grier, Nauman Ahad, Mark Davenport, Matthew Kauf- man, Andrea Giovannucci, and Chethan Pandarinath. Deep inference of latent dynamics with spatio-temporal super-resolution using selective backpropagation through time. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 2331 1244. Curran Associates, Inc., 2021.
[54] Sergey D. Stavisky, Jonathan C. Kao, Paul Nuyujukian, Stephen I. Ryu, and Krishna V. Shenoy, high performing machine interface driven by lowffrquency local field potentials alone and together with spikes, Journal of Neural Engineering, 12(3):036009, May 2015.
[55] Marco Fraccaro, Simon Kamronn, Ulrich Paquet, and Ole Winther. A disentangled recognition and nonlinear dynamics model for unsupervised learning. In Guyon, U. Von Luxburg, S. Bengio, H. Wallach. R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30, pages 3601-3610. Curran Associates, Inc., 2017.
[56) Mikael Henaff, Alfredo Canziani, and Yann LeCun. Model-Predictive Policy Learning with Uncertainty Regularization for Driving in Dense Traffic. In International Conference on Learning Representations (ICLR). September 2018.
[57] Zachary Lipton, David Kale, and Randall Wetzel. Directly modeling missing data in sequences with rnns: Improved classification of clinical time series. In Finale Doshi- Velez, Jim Fackler, David Kale, Byron Wallace, and Jenna Wiens, editors, Proceedings of the 1st Machine Learning for Healthcare Conference, volume 56 of Proceedings of Machine Learning Research, pages 253-270, Northeastem University, Boston, MA. USA, 18-19 Aug 2016. PMLR.
[58] Yonghong Luo, Xiangrui Cai, Ying Zhang, Jun Xu, and Yuan Xiaojie. Multivariate time series imputation with generative adversarial networks. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Biande Bianchi, and Garnett, editors, Advances Neural Information Pro- cessing Systems, volume 31, pages 1596-1607. Curran Associates, Inc., 2018.
[59] Karl Johan Åström. Introduction to stochastic control theory, volume 70 of Mathematics in science and engineering. Academic Press, United States, 1970.
14

---

[60] H. E. Rauch, F. Tung, and C. T. Striebel. Maximum likelihood estimates of linear dynamic systems. AIAA Journal, 3(8):1445-1450, 1965.
[61] Longyuan Li, Junchi Yan, Haiyang Wang, and Yaohui Jin. Anomaly Detection of Time Series with Smoothness- amnothhnss--nduucing Sequential Variational Auto-Encoder, February 2021.
[62] Roozbeh Dehghannasiri, Xiaoning Qian, and Edward R. Dougherty. A Bayesian robust Kalman smoothing framework for state-space models with uncertain noise statistics. EURASIP Journal on Advances in Signal Processing, 2018 (008((1))) (()555 September 2018.
[63] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), 2015.
[64] Joseph G. Makin, Joseph E. O' Doherty, Mariana M. B. Cardoso, and Philip N. Sabes. Superior arm-movement decoding from cortex with a new, unsupervised-learing algorithm. Journal of Neural Engineering, 15(2):226010, January 2018.
[65] Joseph E. O'Doherty, Mariana M. B. Cardoso, Joseph G. Makin, and Philip N. Sabes. Nonhuman Primate Reaching with Multichannel Sensorimotor Cortex Electrophysiology. May 2020.
[66] Robert D. Flint, Eric W. Lindberg, Luke R. Jordan, Lee E. Miller, and Mare W. Slutzky.Accurate Accurate decoding of reaching movements from field potentials in the absence of spikes. Journal of Neural Engineering, (4):046000 June 2012. Publisher: IOP Publishing.
[67] Saskia E. J. de Vries, Jerome A. Lecoq, Michael A. Buice, Peter A. Groblewski, Gabriel K. Ocker, Michael Oliver, David Feng, Nicholas Cain, Peter Ledochowitsch, Daniel Millman, Kate Roll, Marina Garrett, Tom Keenan, Leonard Kuan, Stefan Mihalas, Shawn Olsen, Carol Thomp- son, Wayne Wakeman, Jack Waters, Derric Williams, Chris Barber, Nathan Berbesque, Brandon Blanchard, Nicholas Bowles, Shiella D. Caldejon, Linzy Casal, Andrew Cho, Sissy Cross, Chinh Dang, Tim Dolbeare, Melise Edwards, John Galbraith, Nathalie Gaudreault, Terri Gilbert, Fiona Griffin, Perry Hargrave, Robert Howard, Lawrence Huang, Sean Jewell, Nika Keller, Ulf Knoblich, Josh D. Larkin, Rachael Larsen, Chris Lau, Eric Lee, Felix Lee, Arielle Leon, Lu Li. Fuhui Long, Jennifer Luviano, Kyla Mace, Thuyanh Nguyen, Jed Perkins, Miranda Robertson, Sam Seid, Eric Shea-Brown, Jianghong Shi, Nathan Sjoquist, Cliff Slaughterbeck, David Sullivan, Ryan Valenza, Casey White, Ali Williford, Daniela M. Witten, Jun Zhuang, Hongkui Zeng, Colin Farrell, Lydia Ng, Amy Bernard, John W. Phillips, R. Clay Reid, and Christof Koch. A large-scale standardized physiological survey reveals functional organization of the mouse visual cortex. Nature Neuroscience, Nature Publishing Group
January 2020. Publisher:
[68] Joshua H. Siegle, Xiaoxuan Jia, Séverine Durand, Sam Gale, Corbett Bennett, Nile Graddis, Greggory Heller, Tamina K. Ramirez, Hannah Choi, Jennifer Luviano, Peter A. Groblewski, Ruweida Ahmed, Anton Arkhipov, Amy Bernard, Yazan N. Billeh, Dillan Brown, Michael A. Buice, Nicolas Cain, Shiella Caldejon, Linzy Casal, Andrew Cho, Maggie Chvilicek, Timothy C. Cox, Kael Dai, Daniel J. Denman, Saskia F. J. de Vries, Roald Dietzman, Luke Esposito, Colin Farrell, David Feng, John Galbraith, Marina Garrett, Emily C. Gelfand, Nicole Hancock, Julie A. Harris, Robert Howard, Brian Hu, Ross Hytnen, Ramakrishnan Iyer, Erika Jessett, Katelyn Johnson, India Kato, Justin Kiggins, Sophic Lambert, Jerome Lecoq, Peter Ledochowitsch. Jung Hoon Lee, Arielle Leon, Yang Li, Elizabeth Liang, Fuhui Long. Kyla Mace, Jose Melchior, Daniel Millman, Tyler Mollenkopf, Chelsea Nayan, Lydia Ng, Kiet Ngo, Thuyahn Nguyen, Philip R. Nicovich, Kat North, Gabriel Koch Ocker, Doug Ollerenshaw, Michael Oliver, Marius Pachitariu, Jed Perkins, Melissa Reding, David Reid, Miranda Robertson, Kara Ronellenfitch, Sam Seid, Cliff Slaughterbeck, Michelle Stoecklin. David Sullivan, Ben Sutton, Jackie Swapp, Carol Thompson, Kristen Turner, Wayne Wakeman, Jennifer D. Whitesell, Derric Williams, Ali Williford, Rob Young, Hongkui Zeng, Sarah Naylor, John W. Phillips, R. Clay Reid, Stefan Mihalas, Shawn R. Olsen, and Christof Koch. Survey of spiking in the mouse visual system reveals functional hicrarchy. Nature, 592(7852):86-92, April 2021. Publisher: Nature Publishing Group.
[69] Christian Song, Han-Lin Hsieh, Bijan Pesaran, and Maryam M Shanechi. Modeling and inference methods for switching regime-dependent dynamical systems with multiscale neural observations. Journal of Neural Engineering, 19660666,,, November 2022. Publisher: IOP Publishing.
15

---

[70] Christian Y Song and Maryam M Shanechi. Unsupervised learning of stationary and switch- ing dynamical system models from Poisson observations. Journal of Neural Engineering. 20(6):066029, December 2023. Publisher: IOP Publishing.
[71] Han-Lin Hsieh and Maryam M. Shanechi. Optimizing the learning rate for adaptive estimation of neural encoding models. PLOS Computational Biology, 14(5):e 10061600066 16,, May 2018. Publisher: Public Library of Science
[72] Parima Ahmadipour, Yuxiao Yang, Edward Chang, and Maryam M Shanechi. Adaptive tracking of human ECoG network dynamics. Journal of Neural Engineering, 18(1):01611, February 2021. Publisher: IOP Publishing.
[73] Yuxiao Yang, Parima Ahmadipour, and Maryam M Shanechi. Adaptive latent state model- ing of brain network dynamics with real-time learning rate optimization. Journal of Neural Engineering. 18(3):036013, 3, March 2021. Publisher: IOP Publishing.
[74] Krishna V Shenoy and Jose M Carmena. Combining decoder design and neural adaptation in mrchne interfaces. Neuron, 84(4):665-688 2014.
[75] Maryam M Shanechi. Brain-machin interfaces from motor to mood. Nature neuroscience. 22(10): 1554-1164, 2019.
[76] Lucine L Oganesian and Maryam M Shanechi. Brain computr interfaces for neuropsychiatric disorders. Nature Reviews Bioengineering, 2:653 2.653-670, 2024.
[77] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, llya Sutskever, and Ruslan Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15(56):1929-195) 2014.
[78]
Leslie N. Smith. Cyclical Learning Rates for Training Neural Networks, June 2015.
[79] Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Yee Whye Teh and Mike Titterington, editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pages 249-256, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR.
[80] Alexey Dosovitskiy and Thomas Brox. Generating images with perceptual similarity metrics based on deep networks. In D. Lee. M. Sugiyama, U. Luxburg. I. Guyon, and R. Garnett, editors, Advances Neural Information Processing Systems, volume 29, pages 658-666. Curran Associates, Inc., 2016.
[81] Samuel R. Bowman, Luke Vilnis, Oriol Vinyals, Andrew Dai. Rafal Jozefowicz, and Samy Bengio. Generating sentences from a continuous space. In Stefan Riezler and Yoav Goldberg, editors, Proceedings of the 20th SIGNIL Conference on Computational Natural Language Learning, pages 10-21, Berlin, Germany, August 2016. Association for Computational Linguis- tics.
[82] Ali Razavi, Aäron van den Oord, Ben Poole, and Oriol Vinyals. Preventing posterior collapse with delta-vaes, 2019.
[83] Huajie Shao, Shuochao Yao, Dachun Sun, Aston Zhang, Shengzhong Liu, Dongxin Liu, Jun Wang, and Tarek Abdelzaher. ControlVAE: Controllable variational autoencoder. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 8655-8664. PMLR, 13-18 Jul 2020.
[84] James P. Harris and Dustin J. Tyler. Biological, Mechanical, and Technological Considerations Affecting the Longevity of Intracortical Electrode Recordings. Critical Reviews& in Biomedical Engineering, (1), 2013. Publisher: Begel House Inc.
[85] James C Barrese, Naveen Rao, Kaivon Paroo, Corey Triebwasser, Carlos Vargas-Irwin, Lachlan Franquemont, and John P Donoghue. Failure mode analysis of silicon-based intracortical microelectrode arrays in non-human primates. Journal of Neural Engineering, 1066)446 November 2013, Publisher: IOP Publishing.
16

---

[86] Robert D. Flint, Michael R. Scheid, Zachary A. Wright, Sara A. Solla, and Mare W. Slutzky. Long-Term Stability of Motor Cortical Activity: Implications for Brain Machine Interfaces and Optimal Feedback Control. Journal of Neuroscience, 36(12):3623-366 22322,3,,,,, March 2016. Publisher: Society for Neuroscience Section: Articles.
[87] Dong Wang, Qiaosheng Zhang, Yue Li, Yiwen Wang. Junming Zhu, Shaomin Zhang, and Xiaoxiang Zheng. Long-term decoding stability of local field potentials from silicon arrays in primate motor cortex during a 2D center out task. Journal of Neural Engineering, 1 11(3):03600,, May 2014. Publisher: IOP Publishing.
[88] Dustin A. Heldman and Daniel W. Moran. Chapter 20 Local field potentials for BCI control. In Nick F. Ramsey and José del R. Millán, editors, Handbook of Clinical Neurology, volume 168 of Brain-Corme Comip-Commputr Interfaces, pages 279-288. Elsevier, January 2020.
[89] Gaurav Sharma, Nicholas Annetta, David Friedenberg, Tony Blanco, Daphne Vasconcelos, Ammar Shaikhoumi, Ali R. Rezai, and Chad Bouton. Time Stability and Coherence Analysis of Multiunit, Single-Unit and Local Field Potential Neuronal Signals Chronically Implanted Brain Electrodes. Bioelectronic Medicine, 2(1))63-71, June 2015. Publisher: BioMed Central.
[90] Arjun K. Bansal, Wilson Truccolo, Carlos E. Vargas-Irwin, and John P. Donoghue. Decoding 3D reach and grasp from hybrid signals in motor and premotor cortices: spikes, multiunit activity, and local field potentials. Journal of Neurophystology, 10765):1377-1355 March 2012.
[91] Arda Uran, Kerim Ture, Cosimo Aprile, Alix Trouillet, Florian Fallegger, Emilie C. M. Revol, Azita Emami, Stéphanie P. Lacour, Catherine Dehollain, Yusuf Leblebici, and Volkan Cevher. 16-Channel Neural Recording System-on-Chie With CHT Feature Extraction Processor in 65-nm CMOS. IEEE Journal of Solid-Statee Circuits, 579)22752-2766, September 2022.
[92] Wilson Truccolo, Leigh R. Hochberg, and John P. Donoghue. Collective dynamics in human and monkey sensorimotor cortex: predicting single neuron spikes. Nature Neuroscience, 1 13(1):105-11,, 1110011-1,,100 January 2010.
[93] Chethan Pandarinath, Vikash Gilja, Christine H Blabe, Paul Nuyujukian, Anish A Sarma, Brittany L Sorice, Emad N Eskandar, Leigh R Hochberg, Jaimic M Henderson, and Krishna Shenoy. Neural population dynamics in human motor cortex during movements in people with ALS. eLife, 4:c07436, June 2015. Publisher: eLife Sciences Publications, Ltd.
[94] Hari Teja Kalidindi, Kevin P Cross, Timothy P Lillicrap, Mohsen Omrani, Egidio Falotico, Philip N Sabes, and Stephen H Scott. Rotational dynamics in motor cortex consistent with feedback controller. eLife, 10:e67256, November 2021. Publisher: eLife Sciences Publications, Ltd.
17

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: All claims made the abstract and introduction are supported by our figures and tables in the main body and the Appendix. Our contributions clearly detailed throughout the Methods section and the Appendix.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We discuss several limitations of our work in the Discussion section. Our results are obtained over multiple session recordings and statistics are gathered across 5-fold cross-validation.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged
the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
18

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification:
Our work does not include any theoretical result.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide the data and the link to code of our implementation in the Appendix. We also provide the hyperparameters and training details in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
19

---

reproducible, they 5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We provide the data and the link to the code of our implementation in the Appendix, in which we also share the configuration files we used to obtain our results. The hyperparameters and training details can be found in the Appendix.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co cc/
haagli/c/guudes//deeebbb.m...aa wttpsuubbiiisssooooooooo.o....ee for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
cuipss.c//ublic//uuie guwdes/ CodeeSumiiiseeeeee onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide the hyperparameters and other training-relatn related details our model and benchmark methods in the Appendix.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We provide error bars and confidence bounds for one standard error of the mean (SEM) in our figures, as well statistical significance levels. SEMs are also provided in our benchmark comparison table for all methods.
Guidelines:
The answer NA means that the paper does not include experiments.
20

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We provide the details of the cluster and device specifications on which the model trainings are performed.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer:
[Yes]
Justification:
All aspects of our study conform with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
21

---

We discuss the impact of our work that can improve the accuracy and robustness of neurotechnologies, e.g.,
Justification:
interfaces, for the treatment of brain disorders in millions of patients. We do not anticipate any negative societal impacts of our work.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation. generating fake profiles, surveillance), faimess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work poses no such risks for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We gave the proper references the creators of the data we used in our study. The developers of the benchmark methods and their implementations are also cited throughout the paper.
Guidelines:
The answer NA means that the paper does not existing assets.
22

---

The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license ge.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, pappessitthheeeocc thoode com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-pase re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We explained all components of our model training and architecture in detail throughout the paper, in addition to releasing our codebase with comments that further explains the implementation details.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Our study does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
23

---

Justification:
Our study does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality the research, declaration is not required.
Answer: [NA]
Justification:
The core development of our method does not involve LLM usage.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy http://///oopppsee//oof/ereree/e for what should or should not be described.
24

---

A Appendix
A.1 Time-dropout
To improve the robustness of the MRINE against missing samples, we developed a regularization technique denoted as sime-dropout" Before training MRINE, based on the availability of the observations denoted by T and Tt' (see Section 3 we define mask vectors m m E IR" for and y, respectively:
[1;
if y, is observed at t (i.e., if ET'
m
(13)
0,
otherwise
1,
if s is observed t (i.e., if E T
m =
(14)
0, otherwise
where subscript is the thh component the vector. The general assumption that T' C T, such that is available for all time-steps where y, is observed, motivated by the faster timescale of spikes compared with field potentials. However, this scenario may not always hold as recording devices can have independent failures leading to dropped samples at any time. To mimic the partially missing scenario where either modality can be missing, as well as the fully missing scenarios where both can be missing, randomly replaced (dropped) elements of m and m by 0 at every training step, with the same dropout probabilities p; for both modalities. Masked time-steps (time-steps with mask value) were not used either during the latent factor inference deseribed in Section 3.2 in the computation of loss terms in Eqs. g and so that the inference and model learning were not distorted by missing samples. Instead, note that inference procedure uses the learned model of dynamics to account for missing samples during both inference and learning. We note that we used the original masks (before applying time-dropout) while computing Cam in Eq. [0] as we wanted to obtain smooth representations for all available observations. Note that time-dropowl differs from masked training that is commonly used for training waanformer-basee networks. Masked training aims to predict masked samples from existing samples unlike our training objective (see Section 3.3 for details). Here, the goal of time-dropout is increase the robustness to missing samples by artificially introducing partially and fully missing samples during model training, rather than being the training objective itself. See the ablation study on the effect of time-dropout in Appendix 6. which shows that MRINE models trained with time-dropout had more robust behavior decoding performance and the effect of time-dropout was more prominent with more missing samples.
In addition to the dropout we also applied regular dropout 77 in the encoder's input and output layers with probability pa-
A.2
Training details
A.2.1
Training single-scale networks
For the case of single-scale networks, we use special case of the architecture in Fig. Ta by replacing the multiscale encoder with an MLP encoder, and by just using a single-scale LDM with one modality's decoder network. In particular, for field potentials, we use Gaussian decoder, and for the spikes we use Poisson decoder to obtain the corresponding likelihood distribution parameters. Also, instead of the multiscale LDM of MRINE, this MLP is then followed by a single-scale LDM to perform inference both in real-time and non-causally. During learning of the single-scale networks, we dropped the terms related to the discarded modality from loss functions in Eqs. and 10
A.2.2 Setting the likelihood scaling parameter
The log- likelihood values of modalities with different likelihood distributions may be of different scales, e.g., Poisson ligeliied has a smaller scale than Gaussian log-likelihood in our case (and this can change arbitrarily by simply multiplying the Gaussian modality with any constant). To prevent the model from putting more weight on one modality vs. the other due to higher log-likelihood scale while learning the dynamics, we scaled the log- likelihood of the modality with a smaller scale by parameter r. To set this parameter, we first computed the time averages of each
25

---

modality over T and T' for 8 and yy, respectively, as:
A
1
Σ 1ET
$1
(15)
T
Dr
(16)
PET
where À E R" and μ E 1 Then, we computed the corresponding log-likelihoods of s and yo by using these time averages the means of their likelihood distributions (assuming unit variance for Gaussian distribution). Finally, we set T as the ratio between the higher and smaller seale log-likelihoods:
T (og
(17)
|T| log(Possonl X))
where 1 E R is the I's vector denoting the unit variances for the Gaussian likelihood. The above allows us to balance the contribution of both modalities in our loss function during learning.
A.2.3 Hyperparameters and codebase
Tables 2 and provide the hyperparameters and network architectures used for training single- scale networks and MRINE on stochastic Lorenz simulations and analyses of the real NHP datasets. Our codebase is available at htttps://githut com/ Shwnechiaab/mmnnnnnn... Labl/minmm//
For all models, we used cyclical learning rate scheduler (78) starting with a minimum learning rate of 0,001, and reaching the maximum learning rate of 0.01 in 10 epochs. The maximum learning rate is exponentially decreased by a scale of 0.99. Across all experiments, batch size was set 32, MLP weights were initialized by Xavier-normal initialization [79], and tanh function was used as the activation function hidden layers. All models were trained on CPU servers (AMD Epyc 7513 and 7542, 2.90 GHz with 32 cores) with parallelization.
For scaling hyperparameters of smoothness regularization penalty, the hyperparameter search each dataset is performed using random inner-training and inner-validation split from the training set of the first fold on the first available session over a small grid. For scaling hyperparameters of the Poisson and Gaussian modalities' smoothness regularization penalties (7 and 79), we used a grid of [o, 30, 50, 100, 250] and [0, 5, 10, 50], respectively. For initial experiments, we recommend using smaller scaling hyperparameters, such as 30 and for Poisson and Gaussian modalities, respectively. For the scaling hyperparameter of the smoothness penalty on I; (7+), we used a grid of [0, 30, 50] after finding % and Ty (in order to control the depth of hyperparameter search). Similarly, we recommend using Yar 30 for initial experiments.
Hyperparameters
Models
d,
8,
B,
GC
TE
SS-Poisson
3.128
3,128
32
32
12,34
0.3
0.4
0.1
100
30
0.0001
200
SS.Caussian
3,128
12,344
200
MRINE
3,188
3,128 1.28
3.128
3,128
32
32
12,34
0.3
0.4
0.1
250
0.001
200
Table 2: Hyperparameters used for the stochastic Lorenz attractor simulations. SS denotes single- scale
network. We represent the architecture's ' various MLP encoders and decoders with their parameter notations, and for each, provide the number of hidden layers and hidden units in order, separated by commas. Specifically. d. and - MIP blocks through which and yv are passed in Fig. [ represent the modality-specific encoder networks for Poisson and Gaussian modalities, respectively. om is the fusion network (the last MLP block in Fig. ()). Modality-speciffe decoder networks are A and By for Poisson and Gaussian modalities, respectively.nn, and n. represent dimensions of a, and respectively. K. is the set of future prediction horizons in Eq. g Pt is the time dropout probability on the mask vectors, and Pa denotes the dropout probability applied in the input and output layers of the encoder network (see Appendix l.. GC represents the global gradient clipping norm on learnable parameters. Ye. "/y and 7/2 are scaling parameters for smoothness regularization penalty in Eq. 10 7 is the L2 penalty MLP weights of the encoder and decoder networks. TE denotes the number of training epochs.
26

---

Hyperparaneters
Models
9,
9v
Pr
94
%
7
SS Poisson
3,128
3,128
64
64
1,2.3,
0.3
0.1
0.1
100
0.0001
500
SS-Gaussian
3.128
3,128
64
64
12,3,4
0.3
0.1
0.1
10
30
0.0001
500
MRINE
3.128
3.128
1.28
3,128
3,128
64
1223,4
0.3
0.1
0.1
0.001
500
Table Hyperparameters used for the NHP grid reaching dataset analysis with same timescales for
both modalities. Hyperparameter definitions the same in Table 2
Hyperparameters
Models
d,
6v
Om
8,
By
Px
94
GC
7/4
7
7
SS-Poin
3.128
3,128
12.3,4
0.1
100
0.0001
500
88 Gaussian
3,128
3,128
64
64
1,2,3,4
0.3
0.1
0.1
30
0.0001
500
MRINE
3,128
3,188
1,188
3,128
3,128
64
64
1,2,34
0.3
0.1
0.1
250
30
0.001
500
Table 4: Hyperparameters used for the NHP grid reaching dataset analysis with different timescales
for the different modalities. Hyperparameter definitions are the same as in Table 2
H Hyperparameters
Models
0,
Qy
Om
0,
0y
Pg
R₃
K
GC
7
%
TE
SS-Poisson
3,128
3,128
64
1,23344
0.0001
SS-Gaussian
3,128
3,128
64
12.3,
03
0.1
0.1
30
0.0001
200
MRINE
3,128
3,128
1,128
3.128
3,128
64
64
12.3,4
0.3
0.1
0.1
50
30
0.001
200
Table 5: Hyperparameters used for the NHP center-out reaching dataset analysis with the and
different timescales for the different modalities. Hyperparameter definitions are the same as in Table
Hyperparamettrs
Models
0,
Uy
A
GC
TE
SS-Gaussian
3,128
3,128
64
64
1,23444
033
0.
0.1
0 0.001 500
MRINE
3.128
3,128
1,128
3.128
3.128
64
12,34
30
0.001
Table 6: Hyperparameters used for the visual stimuli dataset analysis. Hyperparameter definitions are
the same in Table Note that both modalities were treated with Gaussian observation model for
MRINE, and single-scale models for both modalities were trained using the hyperparameter denoted for SS-Gaussian.
A.2.4 MSID training
Multiscale subspace identification (MSID) recently proposed linear multiscale dynamical model of neural activity that assumes linear dynamics. MSID can handle different timescales and allow for real-time inference of latent factors 28 We compared MRINE with MSID and showed that compared with the linear approach of MSID, the nonlinear information aggregation enabled by MRINE can improve downstream behavior decoding while still allowing for real-time recursive inference and for handling different timescales. To train MSID, we used the implementation provided by the authors and set horizon hyperparameters with the values provided in their manuscript, i.e. by f2 = 10. For the MSID) results reported in this work, recommended by the developers in their manuscript (28), we fitted MSID models with various latent dimensionalities consisting of [8, 16, 32, 64] and picked the with the best behavior decoding accuracy found with inner-cross validation done on the training data.
A.2.5 mmPLRNN training
Multimodal piecewise-linee recurrent neural networks (mmPLRNN)) is recent multimodal frame- work that assumes piccewise-lincar latent dynamics coupled with modality-specinf observation models [2)]. As discussed in Section 2 mmPLRNN has shown great promise in reconstructing the underlying dynamical system but its inference network operates non non-causally in time and assumes the same timescale between modalities, unlike MRINE. Therefore, train mmPLRNN models with different timescale signals in Table we performed global-meat imputation for the missing signals as done common practice. Also, for mmPLRNN results reported this work, the behavior decoding with mmPLRN PRNNN was performed non-causally. To train mmPLRNN, we used the variational inference training code provided by the authors in their manuscript. However, the default implementation of mmPLRNN only supports Gaussian and categorical distributed modalities. Thus, implemented
27

---

the Poisson observation model by following Appendix C of [29]. Further, we replaced the default linear decoder networks with nonlinear MLP networks for a fair comparison and better performance.
Finally, we trained all mmPLRNN models for 100 epochs. We also performed hyperparameter searches for latent state dimension, learning rate, and number of neurons in encoder/decoder layers over grids of [16, 32, 64), [lc-4, 5e-4, le-3, 1c-2] and [32, 64, 128], respectively.
As shown in Tables I and MRINE outperformed mmPLRNN in behavior decoding for both NHP datasets across all information regimes. Further, we observed that mmPLRNN experienced higher decline in behavior decoding performance when trained with different timescale signals compared to MRINE and MSID (which support training and inference with multiscale signals and missing samples), indicating the importance of multiscale modeling. We believe that the future-step ahead prediction training objective along with new smoothness regularization terms and smoothed reconstruction are also important elements contributing such performance gap (see Appendix (.662) whereas mmPLRNN is trained on optimizing evidence lower bound (HLBO), whose optimization can be challenging due to KL-divergence term [3).
A.2.6 MVAE Training
Multimodal variational autoencoder (MVAE) is variational nutooencoder-basee architecture proposed in 42 that can account for partially paired multimodal datasets by a mixture of experts posterior distribution factorization. However, the notion of partial observations in our work and MVAE are different. In MVAE, partial observations refer to having partially missing data tuples in each element of the batch, which would translate having completely missing LFP spike signals for given trial/segment of multimodal neural activity. However, as we detailed in Section 3 we are interested in modeling multimodal neural activities with different sampling rates, i.e., partially missing time-steps rather than missing either of the signals completely. Further, as discussed in Section 2 the latent factor inference in MVAE is designed to encode each modality to a single factor, whereas MRINE is designed to infer latent factors for each time-step so that behavior decoding can be performed at each time- To account for all these differences, we trained MVAE models without dynamic backbone by treating each time-step as a different data point that allowed us to train MVAE models with partially missing time-steps as done for MRINE. As shown in Table MVAE showed the lowest overall performance among all methods which could be caused by lacking dynamical backbone, unlike other methods.
A.2.7 MMGPVAE training
Multimodal Gaussian process variational autoencoder (MMGPVAE) is another recent multimodal framework that utilizes Gaussian process to model latent distribution underlying multimodal observa- tions [3]. Distinct from other approaches discussed this work, MMGPVAE inference network extracts the frequency content of the latent factors followed by conversion to time domain representa- tions, rather than direct estimation on the time domain. This approach allows MMGPVAE to prune high-frequency content in the latent factors that help in obtaining smooth representations. Similar to mmPLRNN, MMGPVAE does not allow training on modalities with different timescales, thus, for the comparisons provided in Table we performed global-mean imputation for missing samples as common practice. Further, behavior decoding for MMGPVAE is performed non-causally for all MMGPVAE results reported in this work, as it does not support real-time latent factor inference. To train MMGPVAE, we used the authors' official implementation provided in their manuscript. To provide a fair comparison, we trained MMGPVAE models with 64-dimensional latent factors for each modality, where 32 dimensions of 64-dimensional latent factors were shared across modalities (i.e., MMGPVAE models in this work had 96 dimensional latent factors unlike other methods). All MMGPVAE models are trained for 100 epochs as behavior decoding performances reached their peak performance in around 100 epochs, then started degrading due overfitting. The default encoder/decoder architecture for Gaussian modality in the MMGPVAE codebase was implemented for an image dataset, which resulted in poor performance in our dataset. Therefore, for the en- coder/decoder architectures, we used the same architectures as MRINE. Further, we scaled the likelihood of Poisson modality with the same 7 value as done for MRINE (see Section Section|A.222) We also performed a hyperparameter search for the learning rate in a grid of [le-3, 5e-4, Te-4] since default values resulted in poor convergence. Despite MMGPVAE being the best competitor of MRINE as shown in Table 7 it also experiences a significant decline in its performance when trained
28

---

with different timescale signals (Table D which indicates the importance of multiscale modeling. Further, we believe that MRINE's training objective another important factor contributing to its improved performance over MMGPVAF.
A.2.8 CEBRA training
CEBRA is a contrastive learning-based model for neural activity with a convolutional encoder architecture [17). To train CEBRA models for the NHP grid reaching task and the visual stimuli To train multimodal models of CEBRA, we followed the multisession CEBRA training protocol as done in the tutorial. For the visual stimuli dataset, we used the same hyperparameters as in the tutorial. For the NHP grid-reaching task, we used offser-36 encoder models as they achieved better performance, and set the output dimensionality to 64 for fair comparison with other baseline methods.
dataset, we followed the tutorial provided by the authors 2
A.2.9 LFADS training
LFADS is a sequential uutoencoder-based model of unimodal neural activity proposed by B We used the authors' codebase to train LFADS models. Also, we used the hyperparameters in the second row of Supplementary Table (3 with a factor dimension of 64 and used the default learning rate scheduler and early stopping criterion used in the codebase. We followed two different approaches to train LFADS models. First, we trained LFADS models on concatenated multimodal data, which corresponds to treating multimodal data single modality. Second, implemented an by incorporating coodality-speciffe
extension of multimodal LFADS (denoted as LEADS-multimon decoders with appropriate Poisson and Gaussian observation models. Similar MRINE, we applied scaling on the Poisson likelihood loss term (i.e., scaled by 3, which was roughly the automatically selected scale for MRINE) for balancing.
A.3 Latent factor inference times
In addition to enabling real-time inference of latent factors through model design, it is also important to report per-timestep inference times, as these provide a direct measure of the practical efficiency of the approach. Thus, we computed per-timestep inference times of MRINE and bascline methods with an Intel 17-10700K processor, averaged over more than 25,000 timesteps. MRINE's s latent inference time per timestep is 1.82 + 0.13 ms, which is smaller than the timestep considered in this study (10 ms), and smaller than common train-computer interface (BCI) timesteps (e.g.. 50 ms in (54). thus suggesting real-time BCI applicability. Also, the latent inference times per timestep are 1.83 + 0.12 ms for MMGPVAE, 0.03 + 0.008 ms for MSID, and 25.91 + 3.44 ms for mmPLRNN, The multi-modal GP approach had similar inference time compared to MRINE, and MSID achieved the fastest per-timestep latent factor inference due to its simple linear form.
A.4 Stochastic Lorenz attractor simulations
A.4.1 Dynamical system
The following set of dynamical equations defines the stochastic Lorenz attractor system:
dri (x2 - x1 xz)ddl + 91
dr (px) x1 x2) + Q2
(18)
dr3 (x1 B3a3)dt + q3
where 11, .2, and Vy are the latent factors of the Lorenz attractor dynamics, denotes the discretiza- tion time-step of the continuous system and d is the change of variables in di time. We used 10, and dt = 0.006 as in B] q1, q2. and q3 are zero-mean Gaussian dynamic noises with
p == 28, 3 variances of 0.01. We generated 750 trials each containing 200 time-steps, and the initial condition of each trial was obtained by running the system for 500 burn-in steps starting from a random point. Then, normalized trajectories to have zero mean and maximum value of I across time for each latent dimension.
httmpss//cchra.a//oess//// wtt.d//wssdeemooooooooo//oooo notebooks/Demo_ Allen,hhmllnmmml https://laads. wwtp:/aadd.iiiubiiooooooooo io/ mun-mengger/
29

---

A.4.2 Generating high dimensional observations
To generate the Gaussian-distriuteed modalities, we multiplied the normalized latent factors by random matrix C₁ IR x3 and added zero mean Gaussian noise with variance of 5 to generate noisy observations.
To generate the normalized trajectories by another random matrix C, E R and added log baseline firing rate of 5 spikes/ sec with bin-size of ms, followed by exponentiation. We then generated spiking activity by sampling from the Poisson process whose mean is the simulated firing rates.
modalities, we first generated firing rates by multiplying the
A.4.3 Computing the latent reconstruction accuracy
For MRINE and each single-scale network trained only on Poisson or Gaussian modalities, we obtained the smoothed single-scale or multiscale latent factors Fyrr. E {1, 21, T) for each trial in the training and test sets. To quantify how well the inferred latent factors can reconstruct the true latent factors, we fitted a linear regression model from the inferred latent factors of the training set to the corresponding true latent factors. Using the same linear regression model, we reconstructed the true latent factors from the inferred latent factors of the test set. Then, we computed the Pearson correlation coefficient (CC) between the true and reconstructed latent factors for each trial and latent dimension. The reported values are averaged over trials and latent dimensions.
A.5
Real dataset analyses
A.5.1 Nonhuman primate (NHP) grid reaching dataset
In this publicly available dataset 66 65 a macaque monkey performed 2D target-reaching task by controlling a cursor in 2D virtual environment. All experiments were performed in accordance with the US National Research Council's Guide for the Care and Use of Laboratory Animals and were approved by the UCSF Institutional Animal Care and Use Committee. Monkey I was trained to perform continuous reaches to cireular targets with 5 mm visual radius randomly appearing on an 8-by-8 square or an 8-by-17 rectangular grid. The cursor was controlled by the monkey's ss fingertips, and the targets were acquired if the cursor stayed within a 7.5 mm-by-7.5 mm target acceptance zone for 450 ms. Even though there was no inter-tria interval between sequential reaches, there existed a 200 ms lockout interval after target acquisition during which no target could be acquired. After the lockout interval, a new target was randomly drawn from the set of possible targets with replacement. Fingertip position was recorded with six-axis electromagnetic position sensor (Polhemus Liberty, Colchester, VT) at 250 Hz and non-causally low-pass filtered to reject the sensor noise (4th order Butterworth, with 10 Hz cut- frequency). The cursor position was computed by a linear transformation of the fingertip position, and we computed 2D cursor velocity using discrete differentiation of the 2D cursor position in the x and directions. In our analysis, we used the 2D cursor velocity as the behavior variable to decode.
One 96-channel silicon microelectrode array (Blackrock Microsystems) was chronically implanted into the subject's 's right hemisphere primary motor cortex. Each array consisted of 96 electrodes, spaced at 400 μm and covering 4mm-by-44m area. We used multi- spiking activity obtained at a 10 ms timescale, and LFP signals were extracted from the raw neural signals by low-pass filtering with 300 Hz cut-off frequency, and downsampling to either 100 Hz (10 ms timescale) or 20 Hz (50 ms timescale). In our study, we picked the top spiking and I.PP channels based on their individual behavior prediction accuracies and considered a maximum of 20 channels for each modality. As this dataset consists of continuous recordings without a clear trial structure, we created 1-second non- overlapping segments from continuous recordings to form trials so that we could utilize mini-batch gradient descent during model learning.
A.5.2
NHP center-out reaching dataset
In this publicly available dataset [66], a macaque monkey performed a 2D center-out reaching task while grasping a two-link manipulandum. All experiments were performed with approval from the Institutional Animal Care and Use Committee of Northwestern University. Monkey C was trained to perform reaches from a center position to 2-cm square outer targets in an 8-target environment, where outer targets were spaced at 45 45-degree intervals around 10-cm radius circle. Each trial of the task
30

---

started with the illumination of the center target where the monkey had to hold the manipulandum for a random hold time of 0.5-0.6 seconds. After, the center target disuppeared and an outer target was randomly selected from the pool of possible 8 targets, which signaled the monkey to start the reach. To obtain the reward, the monkey had to reach the outer target within 1.5 seconds and hold the manipulandum at the outer target for random time of 0.2-0.4 seconds. Then, the monkey returned back to the center target position and the next trial started. In our analysis, we used 2D manipulandum velocity as the behavior variable to decode.
One 96-channel silicon microelectrode array (Blackrock Microsystems) was chronically implanted into the subject's 's proximal arm area of primary motor (M1) and premotor (PMd) cortices contralateral to the arm used to perform the task. We used multi- unit spiking activity obtained at 10 ms timescale. LFP signals in the original dataset were extracted from the raw neural signals by band-pase filtering between 0.5 and 500 Hz and sampled at 2 kHz. From these L.PP signals, we computed L.PP power signals with window of size 256 ms (moved at 10 ms resolution) over 5 bands (0-4, 7-20, 70-115, 130-200 and 200-300 Hz), resulting in LFP power signals at 100 Hz. For the different timescale analyses, we downsampled LFP power signals to 20 Hz (50 ms timescale). The rest of the dataset generation details are the same as the previous dataset.
A.5.3 Behavior decoding
In our analyses, we took 2D cursor or manipulandum velocity in the x and y directions as the behavior variables for downstream decoding. For all methods, after we inferred latent factors for both training and test sets, we fitted a linear regression model from inferred latent factors of the training set to the corresponding behavior variables. Then, we used the same linear regression model to decode the behavior variables from the inferred latent factors in the test set. We quantified the behavior decoding accuracy by computing the CC between the true and reconstructed behavior variables across time and averaging over behavior dimensions.
When MRINE was trained with spiking activity and LFP signals with timescales of 10 ms and 50 ms, respectively, behavior decoding was performed at the 10 ms timescale for comparisons between MRINE and single-scale networks trained with spike channels. To provide a fair comparison between MRINE and single-scale networks trained with 50 ms LFP signals, inferred latent factors of MRINE were downsampled to 50 ms from 10 ms, and behavior was decoded at every 50 ms in these comparisons with LFP. For all baseline comparisons, behavior decoding with MRINE is performed in real-time (i.e., multiscale latent factors are inferred via cal-timmccausal Kalman filtering) unless otherwise stated.
A.5.4 Behavior decoding with same timescale signals
For both NHP grid reaching and NHP center-out reaching datasets, we also trained MRINE models with various combinations of 5, 10, and 20 channels of the same timescale spike and LFP signals. As shown in Fig. for both datasets, MRINE again successfully improved behavior decoding accuracies as LFP channels are fused with primary spike channels (Figs. 4a.c for NHP grid reaching dataset and Figs. Fe.g for NHP center-out reaching dataset), and when spike channels are added to primary LFP channels (Figs. 4p.d for NHP grid reaching dataset and Figs. [4 for NHP center-out reaching dataset). These results provide further evidence of MRINE's s information aggregation capabilities beyond its multiscale modeling.
Further, we performed the same baseline comparisons for all numbers of primary channels of the same timescale spike and LFP channels (similar to Table [ for different timescale signals). For this analysis, we did not perform imputation for mmPLRNN and MMGPVAE models since both spike and LFP channels have the same timescale, i.c., 10 ms. Note that we did not include MVAE in this analysis since it is not a dynamical model and it had the lowest performance among all methods in Table
As shown in Table 7 MRINE achieves the best behavior decoding performance both with its real-time and noncausal latent factor inference for the NHP grid reaching dataset across all information regimes (p < 10-5 5 n = 20, one-sided Wilcoxon signed-rank test). For the NHP center-out reaching dataset, MRINE achieves the best behavior decoding performance among all methods with noncausal latent factor inference (p < 0.006, n = 15, one-sided Wilcoxon signed-rank test). When MRINE performs real-time (causal) latent factor inference, MMGPVAE outperforms MRINE across low (5 channels)
31

---

C.
Spinn charnes
channets channes
Bakk channels
Spilio
(02
Sp chemchs
0.4
20 Spece dharnes
Sphic charnes
LFR channets
23
033
01
Bont-or dececin 355/
0.7
40
Improveeennin ootornt cobaviar doocding p00
0.0
85 2 18 40
30
0.6
0.4
0.3
Figure 4: Behavior decoding accuracies for the NHP grid reaching dataset (top row) and NHP center- out reaching dataset (bottom row) with same timescale spike and LFP signals. Figure conventions are the same as in Fig.
and high (20 channels) information regimes where MRINE is slightly better in medium (10 channels)
information regime, however, improvements are not statistically significant (p = 0.08 for 5 channels and MMGPVAE MRINE, p = 0.12 for 10 channels and MRINE > MMGPVAE, and p = 0.11 for 20 channels and MMGPVAE > MRINE, n 15, one-sided Wilcoxon signed-rank test).
Beyond the nonlinear multiscale modeling capabilities of MRINE which are crucial for neuroscien- tific application purposes and are not achieved by recent prior nonlinear benchmark methods (see performance differences of mmPLRNN and MMGPVAE between Table [ for different timescale results and Table 7 for same timescale results). these results suggest that MRINE achieves competi- tive performance through its encoder design and training objectives even when both modalities are recorded with the same timescale.
NHP grid reaching
NHP center-out reaching
Spike
10 Spike
2pkk
Spike
Spike
20 Spike
Method
L
LFP
20 LFP
L
10 LFP
20 LFP
MSID
0.452 0.015
0.483 + 0.015
0.544 + 0.016
0.467 + 0.019
0.548 + 0.020
0.596 + 0.022
miPLRNNN
0.455
0.012
0.478
0.011
0.533 + 0.012
0.530
0.022
0.556 + 0.024
0.591 + 0.027
MMGPVAE
0.424
0.012
0.511 # 0.014
0.579 + 0.010
0.558 = 0.022
0.624 + 0.022
0.670 : 0.021
MRINE
0.493 + 0.008
0.566 + 0.010
0.621 + 0.011
0.550 + 0.021
11628 + 0.022
0.663 + 0.020
MRINE noncausal
0.524 0.009
0.586 = 0.011
0.639 1 0.010
0.572 = 0.022
0.647 1 0.023 0.681 = 0.021
Table 7: Behavior decoding accuracies for the NHP grid reaching and center-out reaching datasets
with 5, 10, and 20 channels of same timescale (10 ms) spike and LFP signals for MSID, mmPLRNN. MMGPVAE, and MRINE (both with real-time and noncausal inference). The best-performing method in bold, the second best-performing method underlined, + represents SEM.
A.5.5 Behavior decoding with missing samples
In this analysis, we first trained all bascline methods with 20 spike and 20 LFP channels with different timescales (whose behavior decoding accuracies are shown in samples). To test the robustness of each method to missing samples, we randomly dropped samples in time during inference with fixed sample dropping probabilities for both modalities. Then, we inferred latent factors at all time-steps using only the available observations after sample dropping and performed behavior decoding as described above. Note that even though time-series observations were missing in time, behavior variables were available for all time-steps and were decoded at all time-steps. As shown in Fig. 5) MRINE outperformed all baseline methods for both datasets as it can leverage learned single-scale dynamics to account for missing samples within each data modality.
when there were no missing
32

---

a.
Behavior decoding accuracy
0.7
(CC. n 202
0.5
MRINE
0.4
MSID
MMGPVAE
mmPLRNN
0.4
0.6
0.8
0.2
0.4
0.6
0.8
Spike dropping probability
LFP dropping probability
d.
0.7
Behavior decoding accuracy
0.6
(CC. (15)
0.5
0.4
0.3
0
0.2
0.4 0.6
0.8
0
0.2
0.4
0.6
0.8
Spike dropping probability
LFP dropping probability
Figure 5: Behavior decoding accuracies all models for the NHP grid reaching dataset (top row) and NHP center-out reaching dataset (bottom row) when spike and LFP channels had missing samples and different timescales. Accuracies for models trained with 20 spike and 20 LFP channels. Sample dropping probability of LFPs was fixed at 0.2 while that of spikes was varied as shown on the x-axis. Lines represent mean and shaded areas represent SEM. b. Similar to a when sample dropping probability of spikes was fixed at 0.2 while that of LFPs was varied. c, d. Same as a, b but for NHP center-out reaching dataset.
A.5.6 Behavior decoding with higher channel counts
In our analyses. we mainly focused on channel counts up to 20 due several reasons. First, this low- to mid-inforration regime is quite important for testing multimodal aggregation capabilities, because the contribution from spikes and LFPs is more balanced and the added value combining modalities can be assessed more directly. Second, implantable devices such as chronic BCIs can experience signal loss over time due to various factors 84 85
such as scar tissue formation. Indeed, prior studies have shown that spiking activity from implanted electrodes can degrade faster than LFPs 86 89] Thus, multimodal information aggregation is especially helpful for BCIs in these low- to mid-information regimes because such aggregation can allow patients to continue using their implanted BCI systems for extended periods. Specifically. multimodal aggregation can improve BCI accuracy and robustness even in the face of signal loss by combining spikes and IFPs,, as also shown in prior studies [6. 28 90]. Finally, some recording devices can have lower channel count. For example, recent wireless BCI systems have explored lower channel count designs to reduce bandwidth and power demands, such as a neural system-on-chij using 16 channels with on-chip feature extraction [91].
Nevertheless, show that MRINE's information aggregation capabilities generalizes beyond low- to mid- information regimes, we also trained MRINE models on 30 spike and 30 LFP channels (30-30) as well as on spike and LFP channels (60-60) in the NHP grid reaching dataset, with the LFP and spikes having different timescales. For comparison, we also trained single-scale models, as well as MSID and MMGPVAE. this scenario, we performed behavior decoding from predicted firing rates (Az) as this yielded better performance, unlike in the analyses with up to 20 channels, where
33

---

decoding directly from latent factors was more accurate. As shown in Table 8 MRINE outperformed the multimodal baselines and single-scale models across both 30-30 and 60-60 channel regimes. indicating that its information aggregation capabilities generalize to high-channel count (information)
regimes.
30 Spike
60 Spike
Method
30 LFP
60 LFP
SS LFP
0.482 + 0.012
0,456 = 0.013
Spike
0.639 + 0.008
0.637 +: 0.006
MSID
0.573 + 0.012
0.597 + 0.0T0
MMGPVAE
0.516 + 0.006
0.568 1 0.006
MRINE
0.676 = 0.011
0.693 - 0.010
Table 8: Behavior decoding accuracies for the NHP grid reaching dataset with 30-30 and 60-60 channels of 10 ms spikes and 50 ms LFP for single-scale models (SS), MSID, MMGPVAE, and MRINE. The best-perforning method is in bold, the second best-performing method is underlined, represents SEM.
A.5.7 Behavior decoding performance with R metric
Throughout our analyses, we used CC as the main metric for comparing downstream behavior decoding performance across methods. In addition to this, we also computed the R metric for the 20-channel regime results presented in Table aas an example to show that MRINE's superior performance is not an artifact the reported metric. As shown in Table 9 MRINE again outperforms the baseline methods also when R metric is used.
NHP grid reaching
NHP center-out reaching
Method
CC
Rd
CC
Rd
MVAE
0.425 + 0.009
0.190 + 0.009
0.544 + 0.018
0.298 +: 0.022
MSID
0.519 + 0.012
0.273 + 0.013
0.561 0.020
0.343 1 0.030
mmPLRNN
0.540 1 0.011
0.302 0.013
0.538 = 0.032
0.294 = 0.036
MMGPVAE
0.479 0.017
0.334 0,012
0.601 + 0.021
0.351 + 0.029
MRINE
0.611 + 0.012
0.375 + 0.013
0.649 0.021
0.435 + 0.030
Table 9: Behavior decoding CC and R² accuracies for the NHP grid reaching and center-out reaching datasets with 20 channels of 10 ms spikes and 50 ms LFP for MVAE, MSID, mmPLRNN, MMGPVAE, and MRINE. The pers-rerformnng method is in bold, the second best-performing method is underlined, + represents SEM. CC results are the same in Table
A.5.8 Neural reconstruction
To evaluate each method's information aggregation capabilities beyond behavior decoding, we computed reconstruction accuracies of both modalities under various sample dropping probabilities in addition to having different timescales. To do that, the modality of interest was randomly dropped with varying sample dropping probabilities when the other modality was dropped with a fixed 0.2 probability (Fig. 6). Therefore, the reconstructions of the missing modality were generated by leveraging the learned modality-specific and multiscale dynamics. For all methods, neural reconstructions were obtained with non-causal smoothing. For mmPL.NNN and MMGPVAF, the missing timesteps for both spike and LFP signals were replaced by zeros, and the reconstruction metrics were computed between the true signals and model reconstructions of the missing timesteps. All models were trained with 20 spike and 20 LFP channels. For LFP signals modeled with Gaussian likelihood, we quantified the reconstruction accuracy by computing the CC between the reconstructed mean of the Gaussian likelihood distribution (μ((a,)) and the true observations across time.
For spike signals modeled with Poisson likelihood, reconstruction accuracy was quantified using the area under the curve (AUC) of the receiver operating characteristic (ROC) measure 24 We constructed the ROC by using the reconstructed firing rates, i.e., X(agr)) as the classification scores to determine whether a time-step contained spike or not 992 Both metries were averaged over observation dimensions.
34

---

NHP grid reaching
NHP center-out reaching
Spike
LFP
Avg.
Spike
LFP
Avg.
Total Avg.
Methods
MSID
2.6
3.0
2.8
2.8
3.2
3.0
2.9
1mmPLRNN
2.2
3.6
2.9
2.2
2.4
2.3
2.6
MMGPVAE
3.8
2.4
3.1
4.0
1.8
2.9
3.0
MRINE
1.4
1.0
1.2
1.0
2.6
1.8
1.5
Table 10: Neural reconstruction average ranks of each method for neural modalities and datasets. Average ranks for each neural modality and dataset (each cell) are computed by first ranking each method based on their reconstruction performances at a given sample dropping probability, and then averaging these ranks across the 5 sample dropping regimes in Fig. [ Also, compute the average rank across neural modalities within each dataset (denoted by Avg.) and across both datasets (denoted by Total Avg.). The best-performing method is in bold and the second best-performing method is underlined.
As shown in Table T0 and Fig. 6 MRINE achieved competitive performance compared baselines. To quantify each method's success in reconstructing neural modalities, we computed each method's average rank in terms of how well it reconstructed the neural modalities compared with other methods (Total Avg. in Table IO). To compute these average ranks in TTb we first ranked each method based on its reconstruction performance for each given sample dropping probability, neural modality, and dataset. We then averaged these ranks across all sample-dropping probabilities, neural modalities, and datasets for each method obtain the final average rank (Total Avg.). As shown in Fig. O there exists no method achieving the best performance across all sample dropping probabilities for both modalities and datasets. However, average ranks in Table TO show that MRINE overall achieves the best average performance (Total Avg.) compared to baseline methods, indicating its competitive and robust neural reconstruction performance compared to baseline methods. Therefore, in addition to MRINE's s superior performance in behavior decoding compared to baselines, these results show that MRINE is a valuable tool for neuroscientific studies that go beyond behavior decoding and study neural dynamics as well.
A.5.9 Visualizations of latent dynamics
To compute trial-averaged 3D PCA visualizations, for each algorithm, we first computed 3D PCA projections of latent factors, split them based on trial start and end indices, interpolated them to fixed length (due to variable-length trials), and then computed trial averages of PCA projections for each of 8 and 4 different reach directions for NHP grid reaching dataset and NHP center-out reaching dataset, respectively. As expected based on literature [93 94] all models recovered rotational neural population dynamics (see Fig. 7 for NHP grid reaching dataset and Fig. g for NHP center-out reaching dataset). Among all these algorithms, MRINE had the clearest rotations while each method revealed noisier trajectories for the NHP center-out reaching dataset, potentially due to consisting of smaller numbers of trials whose latent factors are averaged to obtain trial-averaged trajectories.
A.5.10 Visual stimuli dataset
To evaluate MRINE's generalizability beyond neural datasets during motor tasks, we trained MRINE models on high-dimensional (i.e., 800-D) visual stimuli dataset that contained neuropixel spiking activity and calcium imaging data sampled at 120 Hz and 30 Hz, respectively (67 68) Before training the models, we applied Gaussian smoothing on modalities with causal kernel with standard deviation of 8 ms (and treated both modalities with Gaussian observation models). In addition single-scale models trained on either modality, we also trained multisession CEBRA models using multimodal data. As the downstream task, we used the frame ID decoding task, as done in [17], where the goal is to predict the ID of the frame being shown to the subject, ranging from to 900 (30 Hz. and 30s movie). To align the latent factors of MRINE and CEBRA (extracted at 120 Hz following the faster modality) with the downstream target's frequency (30 Hz), we performed mean-pooling For downstream decoding using each method's latent factors, we used k-nearest neighbor and Bayes classifiers, and picked the one that achieved better performance. In addition to single-scale models (rows 3 and 4), CEBRA and MRINE, we also trained downstream decoders directly on neural activity (rows and 2). To quantify the prediction performance, we used the mean absolute error (MAE) between the true and predicted frame ID. As shown in Table II MRINE
35

---

b.
Spike reconstruction accuracy
0.85
LFP reconstrucitin aocuracy
0.90
(CC. 202)
0.80
(AUC,
0.70
0.70
0.85
0.60
0.60
0.55
0.40
0.2
0.4
0.6
0.8
0.2
Spike dropping probabillty
LFP dropping probability
C.
d.
reconstruction Bocuracy
0.75
LFP reconsiruction acouracy
0.95
0.85
(AUC. n= 15)
0.70
0..0
0.75
0.65
0.70
MRINE
0.65
MSID
0.60
0.60
MMGPVAE
0.55
mmPLRNN
0.55
0
0.2
0.4
0.6
0.8
0.50
Spike dropping probability
LFP dropping probability
Figure 6: Neural reconstruction accuracies of 20 spike and 20 LFP channels for the NHP grid reaching dataset (top row) and NHP center-out reaching dataset (bottom row) when spike and LFP channels had different timescales. The reconstruction accuracies of the spike channels. Sample dropping probability of spikes was varied as shown on the x-axis while LFP channels were dropped with 0.2 probability. Lines represent mean and shaded areas represent SEM. b. Similar when sample dropping probability LFPs was varied while spike channels were dropped with 0.2 probabibility. c, d Same but for NHP center-out reaching dataset.
successfully aggregated information across neuropixel spike and calcium imaging modalities of different timescales. Also, decoders trained MRINE's latent factors improved the frame ID prediction performance over decoders trained on single-scale models' latent factors and directly on neural modalities. In addition, MRINE outperformed the CEBRA baseline, potentially due to its explicit dynamical modeling.
Method
MAE
Calcium imaging
151.60
Neuropixel spike
67.76
SS Calcium imaging
22.66
SS Neuropixel spike
10,94
CEBRA
9.31
MRINE
5.09
Table Frame prediction performance the visual stimuli dataset. denotes single-scale model. The best performing method is in bold, the second best performing method is underlined.
A.6 Ablation Studies
A.6.1 Effect of Time-Dropout
To test the effectiveness of time-dropout, for the NHP grid reaching dataset, we performed an ablation study with the same setting used to generate Figs. 5p.b (see Section but we disabled time-dropout
(p: = 0). The remaining hyperparameters were in Table As shown in Fig. 9h. without time-
36

---

a.
MR NE
b.
20p
Cam n
1
Dum
d
MMSPIME
Dam
9
0
Figure 7: 3D PCA visualizations of trial-averes trial-averaged latent factors inferred for NHP grid reaching dataset by a) MRINE, b) MSID, c) mmPLRNN, and d) MMGPVAE.
a.
MRINE
MBID
Do
Din
d.
MYGPVAE
Plan
Dio
Dm 0
3D PCA visualizations of avealeave latent factors inferred for NHP center-out reaching dataset by a) MRINE, b) MSID, c) mmPLRNN, and d) MMGPVAE.
Figure 8:
dropout, the behavior decoding accuracies of MRINE decreased by 7.6% and 31.4% when 40% and 80% of spike samples were missing (in addition to 20% of LFP samples missing), whereas MRINE models trained with time-dropout experienced smaller performance drops of 5.4% and 20.4% the same missing samples settings (see Figs. 5)..). Similarly, MRINE models trained with lime-dropout were more robust to missing I.FP samples (Fig. I vs. Fig- (5)) but the performance drops were smaller due to spiking activity being the dominant modality for behavior decoding in this dataset. As
37

---

expected, the effect of time-dropout was more prominent in the high sample dropping probability regimes (i.e., more missing samples).
a.
b.
0.7-
aocuracy
Behavior decoding B
0.6
20
e0.5-
g
0.4
033
0.2
0.4 0.6
0
0.2
0.4
Spike dropping probability
LFP dropping probability
Figure 9: Behavior decoding accuracies in the NHP grid reaching dataset when time-dropout was disabled (pr = 0), and spike and LFP channels had both missing samples and different timescales. Figure conventions are the same as in Fig. 5
A.6.2 Effect of loss terms in the behavior decoding performance
To gain intuition on the improved performance with MRINE, for the NHP grid reaching dataset, we performed an ablation study on the effect of smoothness regularization terms Eq. TO and smoothed reconstruction term in Eq. 9 on behavior decoding performance. To achieve that, we trained MRINE
models with 20 channels of 10 ms spike and 20 channels of 10 ms LFP signals by removing Eq. y
and individual terms in Eq. 10 from the training objective in Eq. 12 Then, we performed behavior
decoding with these MRINE models as described in Section A.5.3
As shown in Table 12 both smoothness regularization terms in Eq. 10 and smoothed reconstruction term in Eq. 9 are important factors contributing to improved behavior decoding performance as MRINE model trained without Eq. and I0 (row 2) achieve worse performance than that of baseline methods in Table 7 We observed that applying smoothness regularization on I; an important contributing factor to improved performance (row 3 vs row 5) as well as smoothing reconstruction term (row 2 vs row 3). Even though smoothness regularizations on s; and y. may seem marginal when comparing results in rows 3 and 4, they play crucial role when combined with smoothness regularization on æ, (comparing row and row 5).
Behavior Decoding
Model
(CC)
MRINE
0.621 + 0.010
MRINE w/o
0.524 0.013
Eq. 9 and Eq. 10
MRINE w/o
0.565 0.012
Eq. 10
MRINE w/o
0.566 0.016
in Eq. 10
MRINE w/o
0.598 = 0.012
8, and y in Eq.
10
Table 12: Behavior decoding accuracies for the NHP grid reaching dataset with 20 channels of 10 ms spike and channels of 10 ms LFP signals for MRINE models trained without (w/o) loss terms denoted in the first column. The best-performing method is in bold, + represents SEM.
38

---

A.6.3 Effect of using different observation models
To better understand MRINE's performance compared to baseline methods, we trained MRINE models using Gaussian observation models for both modalities on the NHP grid reaching dataset. This approach allowed to evaluate whether the use of distinct observation models for spike and LFP modalities or multiscale modeling is a primary source MRINE's improvements. As shown in Table 13 modeling each modality with an appropriate observation model, i.e., Poisson for spikes and Gaussian for LFPs, significantly improves the performance as MRINE outperformed its variant where both modalities are modeled with the same (Gaussian) observation model (p < 0.0003, 20, one-sided Wilcoxon signed-rank test). Nonetheless, MRINE with the same observation model still outperformed the baseline methods shown in Table indicating that MRINE's improved performance is mainly caused by its multiscale modeling and other elements in the training objective (see Appendix A.6.2
Using the same observation model for both modalities can also enable a direct comparison to unimodal models, since the two modalities can simply be concatenated at the input level under shared Gaussian observation model. To test this, we trained LFADS models using concatenated same timescale (i.e., 10 ms) spike and LFP signals of the NHP grid reaching dataset (see Appendix [A. 2.9.9 for details), as well as CEBRA models (sce Appendix (A.28 for details). In this scenario, MRINE achieved a downstream decoding CC of 0.621 + 0.011, outperforming both LFADS (0.549 + 0.012), MHADD-multimodal (0.547 + 0.011), and CEBRA (0.433 + 0.002) models. These results further suggest that MRINE's superior performance stems from its multiscale architecture and training objective, rather than from simple input concatenation or observation model choices.
Behavior Decoding (CC)
Model
MRINE
0.611 = 0.012
MRINE w/
0.604 = 0.011
Same Observation Model
Table 13: Behavior decoding accuracies for the NHP grid reaching dataset with 20 channels of 10 ms spike and 20 channels of 50 ms LFP signals for MRINE trained with same and different observation models. + represents SEM.
A.6.4 Effect of multiscale encoder design
As discussed in Section 3.2, accounting for different sampling rates for neural signals is an important
consideration for MRINE'S s encoder design shown in Fig. []p. To achieve that, we learn modality- specific LDMs in MRINE's s encoder that can leverage within-modality state dynamics to account for missing samples whether due to timescale differences or missed measurements. Therefore, MRINE impute missing samples, such as zero-
can perform inference without relying on augmentations imputation as done in common practice [53 (7 that can yield suboptimal performance 52] Such suboptimal performance regimes can include either degraded performance due to misinformation presented by imputed signals or discarding the imputed signal almost completely especially in the high imputation regimes due to attempting reconstruct trivial imputations and focusing on only MRINE's performance degraded less than those of baseline methods when trained on different timescale signals due to these design choices.
one neural modality. As shown in Tables] D and |
To further investigate this, for the NHP grid reaching dataset, we trained MRINE models in a similar manner with baseline methods that do not account for training and inference on different timescale signals (i.c, mmPLRNN and MMGPVAE) where the missing LFP signals due to timescale difference are imputed by their global mean, i.c.. zeros due to z-scoring. In this setting. missing LFP timesteps are discarded (masked) in the training objective (as done for mmPLRNN and MMGPVAE) but those timesteps were not treated as missing samples during latent factor inference both during training and inference. In other words, we let MRINE process the misinformation presented by zero imputation similar to mmPLRNN and MMGPVAE due the recurrent nature of latent factor inference (even if their reconstuucions/ppecictioos are masked in the training objective).
In addition, we trained other versions of MRINE, mmPLRNN and MMGPVAE where imputed (missing) timesteps were not discarded in their training objectives (denoted by w/o Loss Masking in Table 14
39

---

Behavior Decoding
Model
(CC)
MRINE
0.611 + 0.012
MRINE w/ Zero Imputation
0.581 + 0.014
MRINE w/ Zero Imputation
0.523 + 0.013
and w/o Loss Masking
mmPLRNN
0.540 + 0.011
mmPLRNN w/o Loss Masking
0.498 0.009
MMGPVAE
0.479 + 0.017
MMGPVAE w/o Loss Masking
0,500 + 0.016
Table 14: Behavior decoding accuracies for the NHP grid reaching dataset with 20 channels of 10 ms spike and 20 channels of 10 ms zero-imputed LFP signals for MRINE, mmPLRNN, and MMGPVAE where zero-imputed LFP time-steps are either included or masked the training objective. The best-performing method is in bold, + represents SEM. MRINE, mmPLRNN, and MMGPVAE performances are taken from Table | for convenience.
As shown in Table the performance of MRINE models trained with zero-imputed LFP signals degraded compared to MRINE models trained with 50 ms LFP signals (row I vs rows 2 and 3). As expected, masking zero-imputed LFP time-steps in the loss function improved behavior decoding performance compared to the scenario where they are included in the loss function (row 2 vs row 3). However, removing zero-imputed LFP time-steps from the training objective still results in degraded performance for MRINE, showing the importance of multiscale encoder design (row 1 vs row 2).
Even though the performance of mmPLRNNN improved significantly when zero-imputed LFP time- steps were masked in the training objective (row vs row 5), allowing it to achieve performance comparable that in Table I unlike mmPLRNN, MMGPVAE achieved higher performance when imputed LFP signals were not masked in its training objective (row 7 vs. row 6). As discussed earlier, due operating in a high imputation regime (i.c., 4 out of every 5 LFP signals are imputed), is possible that including reconstruction of trivial imputed LFP signals may have led MMGPVAE to deprioritize LFP signals during latent factor inference, instead, focusing primarily on spike signals. This behavior suggests that fusing LFP signals with spiking signals can degrade MMGPVAE's performance, as including trivial LFP reconstructions in MMGPVAE's training objective results in higher performance. Therefore, carefully designed encoder is essential to process both modalities effectively and maximize downstream performance.
Overall, MRINE significantly outperformed all models when they were trained with zero-imputed different timescale signals, including its own variants (p < 0.02, n 20, one-sided Wilcoxon signed-rank test). These results show the importance of encoder design when modeling modalities with different timescales.
40


Reviews:
[{"review_id": "X0tEa4CXVt", "reviewer": "Reviewer_VG94", "summary": "The paper introduces a multi-modal encoder/decoder architecture that is specially designed to work with modalities that have different timescales (e.g. binned spikes, and LFP signals). The architecture is also robust to missing samples in either modality. Empirical evaluation is performed on synthetic, as well as real-world (offline) decoding datasets.", "questions": "- What does “primary” modality mean? From the description of MRINE, i don’t see any reason to consider one of the 2 modalities as primary. Does “primary” mean the modality with the fastest timescale? If this is the case, then how can Gaussian modality be primary?\n    \n- In Fig 2a and 3a, the points corresponding to number of LFP channels = 0 -- Are these results from single-modality and single-timescale versions of MRINE?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths:\n\n- Statistical significance tests have been performed for most results. This is highly appreciated!\n    \n- The description of the method is relatively clear.\n    \n- Excellent ablation studies presented in the Appendix.\n    \n\n## Weaknesses:\n\n- The method introduces a few hyperparameters. A description for how these were tuned would be beneficial for further use and adaption of the presented method.", "was_revised": true, "final_justification": "I think the paper is good. My questions were mostly conceptual and the authors have correctly answered them. I was confused between giving a 4 or a 5 rating. But leaving a review at a \"borderline\" recommendation doesn't seem good, so chose 5.", "rebuttal": "We thank the reviewer for their thoughtful feedback and for highlighting these important clarifications. Below, we address the reviewer’s comments:\n\n> Hyperparameters\n\nWe thank the reviewer for their comment. The HPs of MRINE, except time-dropout and smoothness regularization scales, are typical HPs in any deep-learning application, and we followed educated/common choices. For the smoothness regularization scales, the HP search on each dataset is performed using a random inner-training and inner-validation split from the training set of the first fold on the first available session over a small grid. We will make sure to include these details in the manuscript. \n\n> Meaning of primary\n\nWe apologize for the confusion caused by our use of the term “primary.” In our manuscript, we use “primary” to refer to the modality whose dimensionality is held fixed during scaling experiments, while the other modality is gradually added. Importantly, the designation of “primary” is symmetric and does not imply that one modality is inherently more important. For instance, in Figs. 3a and 3e, spike channels are treated as the “primary” modality because their dimensionality is fixed while the number of LFP channels increases. Conversely, in Figs. 3b and 3f, LFP is considered the “primary” modality. This framing was intended to demonstrate that performance improvements are bidirectional: adding either modality to the other enhances decoding performance.\n\n> About Figs 2a and 3a\n\nThe reviewer’s interpretation is correct. The points corresponding to 0 LFP channels (or 0 spike channels) represent the single-modality versions of MRINE. These serve as our unimodal baselines and reflect models trained using only one signal modality and its corresponding timescale."}, {"review_id": "qUjzXj7eEy", "reviewer": "Reviewer_FWMt", "summary": "The submission considers real-time bi-modal inference of a discrete and continuous neural modality evolving at different timescales with missing data. The authors combine several ideas in their work, including a multi-scale encoder and recursive variational inference based on Kalman filtering or smoothing. Empirical results on NHP indicate that the method improves the decoding and reconstruction performance compared to previous approaches.", "questions": "To clarify, (5) means that $a_t$ is deterministic given $s_t$ and $y_t$, so these factors have a point mass distribution. But then in (6) and (7), the modality-specific versions seem to be stochastic with Gaussian noise, and how they are related is not fully clear to me. As far as I understand $x_{t|t}^s$ is then just the mean of the filtering distribution of the latent factor, and does not capture the variance/covariances of the different factor dimensions? \n\nIs T’ subset T a necessary constraint, so that the discrete modality evolves faster?\n\nIn the model description 3.2, it is not clear where the last MLP from Fig.1 comes in?\n\nIt would be useful to analyse how previous work perform in the stochastic Lorenz attractor simulations. Does the relative performance of the suggested approach compared to previous work in this example depend on the specific task/evaluation, such as latent reconstruction accuracy, real-time k-step ahead data predictions? How sensitive are these evaluations with respect to the weighting hyperparameters in the loss function? \n\nWhy do the authors suggest the smoothness regularization term, instead of just using a KL penalty for sequential VAEs? Do the smoothness parameters adapt to the scale of each modality?", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "The combination of different regularisers, model architectures and inference models is indeed interesting. \n\nThe zero-imputation for missing values for the baseline methods seems a bit simple. It would be interesting to see how alternative approaches would peform using less naïve imputations, e.g. just a simple Kalman filter for the continuous modality. \n\nThe rationale why one needs the ‘multiscale embedding factors’ a are not fully clear to me. Why can I not just use x instead (and pre-condition with C, and incorporate the noise r into the noise w?).\n\nThe computational complexity of the method relative to the baselines should be made clearer. For example, I would assume that the method can be beneficial compared with the multi-modal GP approach which may require some approximations for many observations.\n\nI found the presentation to be a bit unclear at times. For example, how are the factors and parameters in (7) with ^s and ^y related to the counterparts in (1)? Does this mean that in the generative path (1), there is no independence assumption, while on the inference/encoder path, the two modalities are encoded independently until they are aggregated via an MLP at the final layer?", "was_revised": true, "final_justification": "The detailed response that has clarified most points (motivation for the multiscale embedding factor, details on the computational complexity of their approach in comparison to previous works, a new benchmark dataset, imputation approaches used in the baseline method). I will therefore increase my score to a weak accept.", "rebuttal": "We thank the reviewer for their questions and comments. Please see our answers below:\n\n> Meaning and rationale of ‘multiscale embedding factors’ $a_t$\n\nThe reviewer is correct that $a_t$ is deterministically computed from the input signals ($s_t$, $y_t$) via the multiscale encoder. Indeed, $a_t$ has a point mass conditioned on the inputs. This holds for both the multimodal case in Eq. (5) and the modality-specific versions in Eqs. (6) and (7). \n\nAfter obtaining $a_t^s$ and $a_t^y$ through the modality-specific MLPs in Fig. 1b, they serve as the observations for modality-specific LDMs. As such, we apply Kalman filtering on modality-specific LDMs to account for timescale differences and obtain $a_{t|t}^s$ and $a_{t|t}^y$. \n\nAlthough $a_t, a_t^s$, and $a_t^y$ are deterministically obtained from neural data, the latent factors $x_t, x_t^s$, and $x_t^y$ in multiscale and modality-specific LDMs are stochastic (conditioned on $a_t, a_t^s$, and $a_t^y$) due to these LDMs’ process and observation noise parameters, which are learned during training. We apply Kalman filtering to these LDMs, yielding $x_{t|t}$ as the posterior mean of the latent state. This filtered estimate is then linearly projected to obtain $a_{t|t} = C x_{t|t}$ (and similarly, $a_{t|t}^s = C_s x_{t|t}^s$ for spikes and $a_{t|t}^y = C_y x_{t|t}^y$ for LFP, which are later fused through the last MLP in Fig. 1b to form $a_t$). While $a_{t|t}$ is also a point estimate (mean of posterior), it is a denoised and temporally smoothed version of the raw $a_t$, benefiting from the learned temporal structure in the LDMs that enable Kalman filtering/denoising over time. Also, the Kalman filter computes this point estimate based on the trained process and observation noise covariances, which together dictate the Kalman gain and posterior covariance. In this way, although the posterior covariance is not explicitly used in our loss, the Kalman filtering process implicitly captures temporal uncertainty and regularity through the explicit use of process and observation noise covariances, which the posterior covariance is a function of.\n\nFurther, the reason for having $a_t$ in our modeling instead of using a single latent factor ($x_t$) is to form an LDM so that we can enable Kalman filtering. Specifically, here $a_t$ becomes the observations of the LDM, which relates it to $x_t$. Through Kalman filtering, we can effectively leverage the learned inherent dynamics for imputing intermediary/missing timesteps, and allow for both real-time recursive filtering and non-causal smoothing that provides further flexibility during inference time. If we directly connect $x_t$ to $y_t$ and $s_t$ without $a_t$, no LDM is formed and thus Kalman filtering cannot be utilized to achieve the above capabilities.\n\n> Imputation of missing samples\n\nWe definitely agree with the reviewer. This is exactly why we included MSID in our baselines, as it handles missing samples through Kalman filtering for continuous modality, and point process filtering for the discrete modality. We show that the MRINE outperforms MSID. For the deep-learning architectures, their model architectures do not support such a setting, so we applied zero-imputation as is common practice. \n\nWhile simple, this approach has been found effective in several domains. For instance, in a recent foundational modeling study on wearable behavioral data, zero-imputation of input signals outperformed learnable interpolation and tokenization-based alternatives designed to circumvent imputation entirely (Erturk et al., 2025).\n\n> Computational complexity\n\nWe thank the reviewer for their comment. We now report the per-timestep inference times for all methods to give an insight into the computational complexities. MRINE’s latent inference time per timestep is $1.82 \\pm 0.13$ ms with Intel i7-10700K processor, which is smaller than the timestep considered in this study (10 ms), and smaller than common brain-computer interface (BCI) timesteps (e.g., 50 ms in Stavisky et al., 2015), thus suggesting real-time BCI applicability. Also, the latent inference times per timestep are $1.83 \\pm 0.12$ ms for MMGPVAE, $0.03 \\pm 0.008$ ms for MSID, and $25.91 \\pm 3.44$ ms for mmPLRNN. The multi-modal GP approach had a similar inference time compared to MRINE, as it does not require time-recurrence but processes the whole data as a batch (which limits its real-time capabilities, as inferring latent factors requires using the future timesteps of neural data).\n\n> Independence assumption\n\nDuring the generation path, we make a conditional independence assumption that $p(y_t, s_t | a_t) = p(y_t | a_t) p(s_t | a_t)$. Multiscale embedding factors $a_t$’s are a function of their modality-specific counterparts, i.e., $a_t = h(a_t^s, a_t^y)$ where $h(\\cdot)$ is composed of modality-specific LDMs and the final fusion layer ($\\phi_m$). As our goal is to extract latent embeddings by fusing information across both modalities, as the reviewer stated, we do not make independence assumptions on the inference/encoder path similar to product-of-experts, or enforce disentanglement across modality-specific factors (factors being merged at the final MLP) through a loss term. We believe that incorporating such independence assumptions and disentanglement are some great future directions to investigate for a better performance.\n\n> $T$ and $T’$\n\nThis is a great question. We indeed assume that $T’$ is a subset of $T$, indicating that the discrete modality evolves faster, in line with the prior work (spikes have a faster sample rate than LFPs). However, our model formulation can be extended to scenarios where $T’$ is not a subset of $T$. In those scenarios, we would represent both signals in the greatest-common-divisor timescale with appropriate masks to indicate the availability of the signal samples. For instance, if spike and LFP signals were sampled at 20 and 50 ms, we could represent each signal at 10ms resolution during modeling, unmasking the spikes every other timestep and LFPs every 5 steps.\n\n> Last MLP in Fig. 1\n\nThe last MLP in Fig. 1b is the fusion module ($\\phi_m$) that fuses the $a_{t|t}^y$ and $a_{t|t}^s$ to form $a_t$, as we explained in L200.\n\n> Comparisons in other tasks\n\nThank you for the suggestion. The Lorenz simulations were designed primarily to validate MRINE’s information aggregation capabilities in a controlled setting, rather than to benchmark against prior work. We believe real neural datasets offer a more meaningful testbed, and accordingly, we evaluate all methods across two distinct datasets with different behavioral tasks and subjects, and across different timescale regimes, multiple information levels, and several sample-dropping scenarios. To further support MRINE’s generalization, we now include new results on a third dataset with another distinct task (see item 3 in our response to reviewer G325), where MRINE again outperforms CEBRA and single-scale models.\n\nAs shown in Appendix A.6, hyperparameters—particularly smoothness regularization—do impact performance, and removing them narrows the gap to baselines. This shows the benefit of these regularization components. Incorporating automated hyperparameter tuning (e.g., population-based training) is a promising direction for future work.\n\n> KL divergence term as in SAEs\n\nThank you for the insightful question. The smoothness regularization terms in our model are implemented as KL divergences between latent representations at consecutive timesteps. The corresponding weights are chosen to reflect the scale of each modality.\n\nDuring our development phase, we initially explored ELBO-based training, as in sequential VAEs (e.g., SAEs), but found it to be unstable and consistently underperforming in our setting. Our current objective provided more stable optimization and better decoding performance.\n\nTo further support this design choice, we now include comparisons with LFADS (an SAE variant) in item 2 of our response to reviewer BeqV, where MRINE demonstrates stronger performance.\n\nOverall, we hope that the above analyses, which provide further evidence on our model’s performance and generalization performance, address the reviewer’s concerns and encourage the reviewer to reconsider our work.\n\nReferences: \n\n- Erturk, E., Kamran, F., Abbaspourazad, S., Jewell, S., Sharma, H., Li, Y., Williamson, S., Foti, N. J., & Futoma, J. (2025). Beyond Sensor Data: Foundation Models of Behavioral Data from Wearables Improve Health Predictions. Forty-second International Conference on Machine Learning. \n- Stavisky, S. D., Kao, J. C., Nuyujukian, P., Ryu, S. I., & Shenoy, K. V. (2015). A high performing brain–machine interface driven by low-frequency local field potentials alone and together with spikes. Journal of Neural Engineering."}, {"review_id": "H6wN9QHHel", "reviewer": "Reviewer_G325", "summary": "The paper proposes a method based on dynamical systems for analyzing multi-modal, multi-scale neural data. The proposed method, unlike prior methods, accounts for differences in scales (sampling rates) in the modalities. It first encodes each modality separately and applies non-linear fusion to learn a shared representation.  From this representation, it reconstructs the data across modalities.  Through shared representation, it allows for missing data. Using both synthetic and real-world experiments, the paper demonstrates improved behavioral decoding over both single-modal model and multi-modal models, which do not account for timescale differences.", "questions": "Questions on architecture\n\n* On the experiments, including results with more common evaluation metrics such as R2, and considering realistic data settings (see weakness W1) could improve the paper. \n\n* In Figure 1b, the indexing for both modalities appears to be the same. Do you apply prediction for a common set to time indices, before concatenation? \n\n* If Fig 1b, And are the MLPs time invariant (the arrows indicate all inputs (across time) are fed to it at once)?  Maybe adding clarifying text could help. \n\n*  On the smoothness regularizers (Eqn. 10, 3rd term), the choice for fast/slow latent dimensions seems arbitrary. Any insights on the effect of this choice?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strength**\n \n\nS1: The paper is well-written, with a good literature review and a clear problem statement. The paper provides a good argument for why the non-causal methods fail in multi-scale settings, motivating a causal approach. \n\n \n\nS2: While other methods have proposed using deep learning based dynamical models for single-scale and also multi-scale linear LDS models, this paper considers a newer challenge (at least for deep dynamical models) of non-linearly fusing information across modalities, while still maintaining real-time inference, and handling missing data.  \n\n \n\nS3: The work conducted an interesting set of experiments to show the utility of fusing information across multiple modalities, while still accounting for their timescale differences of the modalities.  The ablation studies demonstrate the effectiveness of different components of the model in improving behavioral decoding in a multi-scale setting.  \n\nWeakness  \n\n \nW1: The number of channels used for the experiment (5-20) is rather small. This deviates from a realistic setting with recordings having 100s of simultaneously recorded neurons (channels). The plots shown in Figure 3 show the diminishing returns in about 10 channels. The current experiments are limited in assessing the impact of the method in more realistic settings.  \n\n\nW2: The evaluation metric for both latent reconstruction (section 4.1) and decoding accuracy (section 4.2) somewhat feels unjustified, since most prior approaches use R2 to directly measure the reconstruction error instead of correlation. How does the correlation inform us about the performance of the model?  Therefore, interpreting percentage improvement in this metric is hard. And plots comparing groundtruth and inferred latents/behavioral quantities could be helpful to see.  \n \n\nW3: The experiments appear to be limited to LFP and spikes, but cited prior methods like MMGPVAE explore wider ranges of multi-modal data, such as (temporal) video recordings of behavior, calcium imaging etc. Absence of other modalities beyond LFP and spikes limits the broader application of the method.", "was_revised": true, "final_justification": "The paper deals with a well-motivated problem. My main concern was with the experiment setup. The authors used a small number of channels in the experiments, which does not fit very well with realistic neural datasets. This limits the application of the method in real datasets. While the authors have provided some discussion during rebuttal, my concern is not fully addressed. Therefore, I maintain my original score.", "rebuttal": "We thank the reviewer for providing an excellent summary of our work and its strength. Below, we address their comments:\n\n> Small number of channels\n\nThis is a great insight, we thank the reviewer for raising that. There are two reasons why we focused on 20 channels in our original motor dataset. First, the number of channels it takes for motor decoding performance to saturate depends on the task complexity, and this number is typically around or less than 20 neurons for benchmark neural datasets, including the public data we had. For this reason, we went up to 20 neurons because, in that regime, performance was not saturated, so we could ask whether multimodal methods can aggregate information across modalities to improve performance. This saturation is also observed in the prior study that tested MSID. Second, while in laboratory neuroscience experiments, one can record more than 20 neurons, this is not the case for chronic BCIs because they are based on implantable devices with fewer channels, and further, because they lose the signal over time due to scar tissue formation. The main advantage of multimodal aggregation would be in low- to mid-information regimes for one modality, such that the other modality can help improve accuracy and robustness, as also shown in prior studies (Bansal et al., 2012; Ahmadipour et al., 2024; Hsieh et al., 2018). For example, prior studies show that spiking activity from implanted electrodes may degrade faster than LFPs (Stavisky et al., 2015), leading to a small number of spike channels. This is exactly the scenario where using the more robust LFP signals is the most advantageous, as it can provide significant robustness and performance improvements, as we show in the manuscript.\n\nNevertheless, to demonstrate that MRINE is not limited to low-dimensional inputs, we additionally evaluated it on a high-dimensional (800-channel) calcium imaging and spike dataset (see item 3 below). MRINE outperformed CEBRA and all other baselines on this task, supporting its scalability to more realistic, large-scale settings.\n\n> Decoding performance on R2:\n\nPer the reviewer’s request, we now provide the 20 Spike - 20 LFP results in Table 1 using R2 below. We show that the relative performance comparisons between MRINE and baseline methods are consistent across CC and R2.\n\nNHP grid reaching:\n| Model | Decoding R2 (mean $\\pm$ std) | \n| - | - | \n| MVAE | 0.190 $\\pm$ 0.009 | \n| MSID | 0.273 $\\pm$ 0.013 |\n| mmPLRNN | 0.302 $\\pm$ 0.013 |\n| MMGPVAE | 0.334 $\\pm$ 0.012 |\n| **MRINE** | **0.375 $\\pm$ 0.013** |\n\nNHP center-out reaching:\n| Model | Decoding R2 (mean $\\pm$ std) | \n| - | - | \n| MVAE | 0.298 $\\pm$ 0.022 | \n| MSID | 0.343 $\\pm$ 0.030 |\n| mmPLRNN | 0.294 $\\pm$ 0.036 |\n| MMGPVAE | 0.351 $\\pm$ 0.029 |\n| **MRINE** | **0.435 $\\pm$ 0.030** |\n\nWe will also include example decoded trajectories in our manuscript.\n\n> Broader application of the method\n\nWe thank the reviewer for their suggestion. We now tested MRINE on a completely new dataset that contained high-dimensional calcium imaging and spike data, where the downstream task is to predict the visual stimuli (ID of the frame) as done in CEBRA (Schneider et al. 2023). We compare the mean absolute errors (MAE) of frame ID prediction for MRINE, CEBRA, single-scale (SS) models trained on either calcium imaging or spike data, and for decoding directly from the data modalities using naive Bayes classifiers. MRINE outperforms all baselines:\n\n| Model | MAE|\n|-|-|\n|Calcium imaging| 151.60|\n|Spike|67.76|\n|SS calcium imaging|22.66|\n|SS spike|10.94|\n|CEBRA|9.31|\n|**MRINE**|**5.09**|\n\n> Indices in Fig. 1b\n\nWe used the same time index for latent representations of both modalities as modality-specific LDM of the slower timescale signals extracts latent representations in the timescale of the faster signal (i.e., spikes). We explain the procedure in L188-203.\n\n> MLP in Fig. 1b\n\nThe reviewer’s understanding is indeed correct, MLP modules are time invariant, and the dynamical information is leveraged through linear dynamical models. We will make sure to add text in the figure to make this clear.\n\n> Slow/fast scales on smoothness regularization\n\nThis is a great question. While we don’t report full results in the paper, we observed during development that applying smoothness regularization to half of the latent dimensions yielded equal or better performance than applying it to all dimensions. Indeed, this avoids overly strong regularization, which can bias the model toward overly smooth dynamics and suppress informative but irregular neural patterns. \n\nReferences: \n- Stavisky, S. D., Kao, J. C., Nuyujukian, P., Ryu, S. I., & Shenoy, K. V. (2015). A high performing brain–machine interface driven by low-frequency local field potentials alone and together with spikes. Journal of Neural Engineering.\n- Bansal, A. K., Truccolo, W., Vargas-Irwin, C. E., & Donoghue, J. P. (2012). Decoding 3D reach and grasp from hybrid signals in motor and premotor cortices: Spikes, multiunit activity, and local field potentials. Journal of Neurophysiology.\n- Ahmadipour, P., Sani, O. G., Pesaran, B., & Shanechi, M. M. (2024). Multimodal subspace identification for modeling discrete-continuous spiking and field potential population activity. Journal of Neural Engineering.\n- Hsieh, H.-L., Wong, Y. T., Pesaran, B., & Shanechi, M. M. (2018). Multiscale modeling and decoding algorithms for spike-field activity. Journal of Neural Engineering.\n- Schneider, S., Lee, J. H., & Mathis, M. W. (2023). Learnable latent embeddings for joint behavioural and neural analysis. Nature."}, {"review_id": "29hUHgBRjE", "reviewer": "Reviewer_BeqV", "summary": "The paper introduces MRINE, a method for real-time decoding of behavior from multimodal neural data with different timescales. By learning within-modality dynamics and fusing information nonlinearly, it improves decoding accuracy on simulations and NHP datasets compared to existing methods, showing potential for real-time BCI applications.", "questions": "1. Why did the authors only use 20 channels when the dataset provides two 96-channel Utah arrays? Later channels with lower R² can still carry meaningful information.\n\n2. Why is decoding accuracy evaluated only with CC? A CC around 0.6 is not ideal—what is the decoding R²?\n\n3. How is the real-time performance (latency, throughput) of the decoder evaluated in practice?\n\n4. Why does Table 7 use ranks, and what criteria were used for ranking?\n\n5. In Figure 6, what happens beyond 0.8? The trend has not yet plateaued at that point.\n\n6. How difficult is training this network in practice, and what ensures convergence during training?", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths\n- Good motivation: tackles multimodal fusion of LFP and spikes, missing data, and timescale mismatch.\n- Experiments are detailed and thorough.\n\nWeaknesses\n- Overall method novelty is limited: The MRINE method can be roughly understood as a combination of multimodal inputs (LFP and spike signals), nonlinear dynamics to extract latent variables, fusion of these latent variables, and a self-encoding structure to reconstruct neural signals. While this integration is useful, the overall novelty of the method is limited. The approach is largely based on existing techniques, such as multimodal fusion, nonlinear dynamic systems, and autoencoders, without introducing a fundamentally new concept. Additionally, many spike decoding methods bin spikes into time windows and model the firing rates as Gaussian distributions. This method is widely used and effectively captures the high-frequency nature of spike signals. It would be helpful to discuss the advantages and limitations of this approach and compare it with MRINE’s decoding strategy.\n- Uses too few channels and sessions: Although the paper claims the dataset contains 96 channels of neural data, the original dataset cited actually includes two 96-channel electrode arrays. However, the experiments in the paper only used 20 channels. This reduced number of channels is somewhat limiting for decoding performance, as even channels with low R² values may still carry valuable movement-related information. If the number of channels increases, a single modality signal (such as spike or LFP) might be sufficient to perform the motion decoding task, which makes the use and fusion of multimodal data seem less necessary. Furthermore, the number of sessions used in the study is limited, which raises concerns about the generalizability of the results. A larger dataset with more sessions would provide a more reliable assessment of the model's robustness and its ability to generalize across different experimental conditions.\n- Lacks comparison with some mainstream or latest decoding methods like LFADS, CEBRA, MINT, and BAND.\n- Does not consider non-stationarity of neural signals: Neural signals are inherently non-stationary, meaning their statistical properties can change over time. The relationship between neural signals and movement signals is likely time-dependent, and encoding functions may vary over time. However, MRINE assumes that the encoding model remains stationary, which overlooks the time-varying nature of neural data. This assumption could lead to unstable performance in decoding tasks, particularly in long-duration experiments or varying conditions. Incorporating mechanisms to adapt to the non-stationarity of neural signals would make the model more robust and applicable to real-time scenarios.\n- Real-time performance is not evaluated: Although the paper claims that MRINE is designed for real-time decoding, it does not provide a detailed evaluation of the model's decoding efficiency. Real-time performance, including latency, throughput, and computational demands, is critical for brain-computer interfaces (BCIs) and other real-time applications.", "was_revised": true, "final_justification": "The authors have added extensive additional experiments, which have resolved my concerns. I am satisfied with the revisions and will increase my score.", "rebuttal": "We thank the reviewer for acknowledging the strengths of our problem motivation and experimentation. Below, we address their comments.\n\n> Fundamentally new concepts\n\nThis is an important point to clarify. While the basic individual building blocks in MRINE (e.g., MLPs used as part of the encoder) exist, the novelty in MRINE is in how to use these basic blocks to develop a novel overall architecture, training method, and inference method that solves problems of importance to neuroscience and neurotechnology: real-time multimodal inference, different timescales, missing samples. **We clarify that this is similar to many important methods in ML that use existing basic blocks to come up with a new method, for example, all of our baselines are similar in this respect**. We emphasize that the main novelty of MRINE lies in designing a new technique for addressing problems in important applications in neuroscience, rather than proposing new individual deep learning components.\n\n> Gaussian modeling of spikes and additional baselines\n\nThis is a great point, and we agree with the reviewer that it’s a natural choice to model binned spikes through Gaussian distribution. Precisely for that reason, we have an ablation study in Appendix A.6.1 that shows that an MRINE variant, which modeled binned spikes through Gaussian distribution, achieved slightly worse but competitive performance compared to that of MRINE. Importantly, through this ablation study, we also show that MRINE’s improved performance does not simply result from modeling two modalities through different observation models (Gaussian vs. Poisson), but also from its model/inference design and training objective.\n\nAlso, the above choice would enable testing the performance of single-scale models of neural activity, such as LFADS, on multimodal signals by concatenating them at the input level. Following the reviewer’s recommendation, we now trained LFADS on same timescale spike and LFP signals of NHP grid-reaching dataset. In addition, we also trained CEBRA  models on these multimodal signals by following the CEBRA manuscript authors’ recommendation on multimodal training (i.e., treating each modality as a separate recording session and training a multi-session CEBRA model) (Schneider et al., 2023). As shown below, MRINE outperformed both methods in downstream behavior decoding. \n\nWe thank the reviewer for providing this opportunity to provide more evidence, in addition to our existing baseline results across 7 sessions from 2 different datasets on various settings such as different information regimes, different timescale scenarios, and several sample-dropping regimes.\n\n|Model | CC (mean $\\pm$ std)|\n|-|-|\n|**MRINE**|\t **0.621 $\\pm$ 0.011**|\n|LFADS|\t 0.547 $\\pm$ 0.011|\n|CEBRA| 0.433 $\\pm$ 0.002|\n\n> Using few channels\n\nWe thank the reviewer for raising this point. There are two reasons we focused on 20 channels in our original motor dataset. First, the number of channels it takes for motor decoding performance to saturate depends on the task complexity, and this number is typically around or less than 20 neurons for benchmark neural datasets, including the public data we had. For this reason, we went up to 20 neurons because, in that regime, performance was not saturated, so we could ask whether multimodal methods can aggregate information across modalities to improve performance. This saturation is also observed in the prior study that tested MSID. Second, while in laboratory neuroscience experiments, one can record more than 20 neurons, this is not the case for chronic BCIs because they are based on implantable devices with fewer channels, and further because they lose the signal over time due to scar tissue formation. The main advantage of multimodal aggregation would be in low- to mid-information regimes for one modality, such that the other modality can help improve accuracy and robustness, as also shown in prior studies (Bansal et al., 2012; Ahmadipour et al., 2024; Hsieh et al., 2018). For example, prior studies show that spiking activity from implanted electrodes may degrade faster than local field potentials (Stavisky et al., 2015), leading to a small number of spike channels. This is exactly the scenario where using the more robust LFP signals is the most advantageous, as it can provide significant robustness and performance improvements, as we show in the manuscript. \n\nNevertheless, also per reviewer G325’s request, we now compared MRINE to CEBRA in a 800-D calcium imaging and spike dataset, on which MRINE outperformed CEBRA in downstream visual stimuli decoding task (see item 3 in our response to reviewer G325). Overall, this new analysis shows that MRINE can well generalize to higher-dimensional datasets.\n\n> Real-time compatibility\n\nWe thank the reviewer for their comment. We agree with the reviewer that reporting the computational complexity of MRINE is important to ensure its real-time compatibility. Thus, we now report the per-timestep inference times for all methods to give an insight into the computational complexities. MRINE’s latent inference time per timestep is $1.82 \\pm 0.13$ ms with Intel i7-10700K processor, which is smaller than the timestep considered in this study (10 ms), and smaller than common brain-computer interfaces (BCI) timesteps (e.g., 50 ms in Stavisky et al., 2015), thus suggesting real-time BCI applicability. Also, the latent inference times per timestep are $1.83 \\pm 0.12$ ms for MMGPVAE, $0.03 \\pm 0.008$ ms for MSID, and $25.91 \\pm 3.44$ ms for mmPLRNN. While closed-loop real-time experiments with animals/humans in the loop constitute a major neurophysiological/experimental effort outside the scope of this modeling work, the above computational times show MRINE’s applicability for such real-time BCI experiments.\n\n> Stationarity of neural signals\n\nWe agree with the reviewer that accounting for non-stationarity is important for real-world decoding and discuss this as a limitation in the manuscript. Like many recent models—including MMGPVAE, mmPLRNN, LFADS, and CEBRA—MRINE assumes stationary encoding. Extending MRINE to handle distributional shifts or task-dependent changes is a valuable direction for future work, particularly for long-term or real-time applications. Indeed, developing calibration or alignment methods over time to address time variations is a complementary yet distinct direction of work (Farshchian et al., 2018; Karpowicz et al., 2025).\n\n> Decoding R2\n\nWe now provide the 20 Spike - 20 LFP results in Table 1 using R2 below. We show that the relative performance comparisons between MRINE and baseline methods are consistent across CC and R2.\n\n\nNHP grid reaching:\n| Model | Decoding R2 (mean $\\pm$ std) | \n| - | - | \n| MVAE | 0.190 $\\pm$ 0.009 | \n| MSID | 0.273 $\\pm$ 0.013 |\n| mmPLRNN | 0.302 $\\pm$ 0.013 |\n| MMGPVAE | 0.334 $\\pm$ 0.012 |\n| **MRINE** | **0.375 $\\pm$ 0.013** |\n\nNHP center-out reaching:\n| Model | Decoding R2 (mean $\\pm$ std) | \n| - | - | \n| MVAE | 0.298 $\\pm$ 0.022 | \n| MSID | 0.343 $\\pm$ 0.030 |\n| mmPLRNN | 0.294 $\\pm$ 0.036 |\n| MMGPVAE | 0.351 $\\pm$ 0.029 |\n| **MRINE** | **0.435 $\\pm$ 0.030** |\n\n> Using ranks in Table 7\n\nWe used ranks in Table 7 because, unlike Fig. 5 (behavior decoding), Fig. 6 (neural reconstruction) did not show a consistent top-performing model across all sample dropping regimes. To summarize overall performance, we ranked the four models (1 = best, 4 = worst) within each regime based on reconstruction accuracy, then averaged these ranks across all regimes. This provides a clearer aggregate comparison when no single model dominates across conditions, as described in Appendix A.4.6.\n\n> Regime beyond 0.8\n\nWe limited the sample dropping rate to a maximum of 0.8 because beyond that point, one modality is often almost entirely missing, making it difficult to meaningfully assess multimodal inference. At such extreme sparsity levels, the task shifts from inference to single-modality reconstruction, which is outside the scope of our current evaluation.\n\n> Difficulty of training the model\n\nTraining MRINE in practice is not difficult, since we observed that most of the default hyperparameters, such as in Tables 3,4,5 are the same, except for the smoothness regularization parameters for which we recommend a hyperparameter search over a small grid. We used the validation loss as the convergence criterion.\n \nOverall, we thank the reviewer for their careful evaluation and for highlighting several important points, which we have addressed in detail with additional analyses and clarifications. We hope that these responses and the new evidence provided will help the reviewer to reassess their overall evaluation of our work.\n\nReferences:\n\n- Stavisky et al. (2015). A high performing brain–machine interface driven by low-frequency local field potentials alone and together with spikes. Journal of Neural Engineering.\n- Bansal et al. (2012). Decoding 3D reach and grasp from hybrid signals in motor and premotor cortices: Spikes, multiunit activity, and local field potentials. Journal of Neurophysiology.\n- Ahmadipour et al. (2024). Multimodal subspace identification for modeling discrete-continuous spiking and field potential population activity. Journal of Neural Engineering.\n- Hsieh et al. (2018). Multiscale modeling and decoding algorithms for spike-field activity. Journal of Neural Engineering.\n- Farshchian et al. (2018). Adversarial Domain Adaptation for Stable Brain-Machine Interfaces. International Conference on Learning Representations. \n- Karpowicz et al. (2025). Stabilizing brain-computer interfaces through alignment of latent dynamics. Nature Communications.\n- Schneider et al. (2023). Learnable latent embeddings for joint behavioural and neural analysis. Nature."}, {"review_id": "3AGwO09B98", "reviewer": "Reviewer_aE6q", "summary": "In this paper, the authors introduces Multiscale Real-time Inference of Nonlinear Embeddings (MRINE) that is a novel learning framework designed for real-time, recursive decoding by dynamically and nonlinearly aggregating information from multiple neural time-series modalities.\n\nIt specifically addresses challenges posed by multiscale data, handling modalities with different timescales (e.g., fast spiking activity and slower LFP activity), varying probabilistic distributions (e.g., Poisson for spikes and Gaussian for LFP), and missing samples.\n\nKey to MRINE is its novel multiscale encoder design, which first learns within-modality dynamics using Linear Dynamical Models (LDMs) to account for timescale differences and predict missing data, then performs nonlinear information fusion via a fusion network. The training objective incorporates a multi-horizon k-step-ahead prediction loss, smoothed reconstruction, and smoothness regularization on latent dynamics to encourage learning meaningful representations and prevent overfitting. It also introduces a time-dropout technique to further enhance robustness to missing samples during training.\n\nEmpirically, MRINE has demonstrated its ability to improve real-time target decoding in both stochastic Lorenz attractor simulations and two independent nonhuman primate (NHP) brain datasets. It outperforms various linear and nonlinear multimodal benchmarks, showing particular strength when modalities have different timescales or missing samples, and also exhibits robustness to randomly missing samples during inference.", "questions": "- What is the computational complexity (in Big O notation) per time step, and what is the actual computation time required per time step in practice?\n- How effective is the model in terms of decoding performance when applied to reconstructed observations?", "limitations": "Yes.", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n*   The MRINE framework is designed for real-time recursive decoding from multiscale neural modalities. Its architecture, featuring a multiscale encoder with modality-specific Linear Dynamical Models (LDMs) and Kalman filtering, enables nonlinear information aggregation and intrinsically accounts for missing data, avoiding suboptimal zero-padding.\n*    MRINE consistently outperforms various linear and nonlinear multimodal benchmarks (MSID, mmPLRNN, MMGPVAE, MVAE) in both stochastic Lorenz attractor simulations and nonhuman primate (NHP) brain datasets. The model demonstrates strong robustness against randomly missing spike or LFP samples during inference, outperforming baselines. It also achieves competitive performance in neural reconstruction.\n*   Its real-time decoding capabilities and improved accuracy from multimodal fusion show potential for advanced brain-computer interfaces.\n\nWeaknesses\n* The model is based on the assumption of stationary linear dynamical systems, which may limit its capacity to effectively capture the complexity of neural dynamics.\n* Multiscale constraint: The framework requires that one time index be a subset of the other, which may restrict modeling flexibility.\n* Real-time performance: The paper lacks explicit discussion of computational efficiency or scalability. The term \"real-time\" appears to refer to the causal (filtering) nature of the inference procedure. To substantiate claims of real-time capability, the authors should demonstrate that the computational complexity per time step is constant with respect to the number of samples, and that the actual computation time per step does not exceed the duration of the time bin(s).\n* Evaluation methodology: The paper does not clearly specify whether and how k-step prediction was applied to the non-human primate (NHP) datasets. As a result, it is uncertain whether the proposed approach is effective at long-term forecasting or primarily suited for short-term tracking.\n* Comparison of latent factor decoding: The compared methods utilize different encoder and decoder architectures, typically involving nonlinear mappings between latent representations and observations. Therefore, evaluating performance based on linear decodability or Pearson correlation of latent factors may not provide a fair or meaningful comparison.\n* Related work: A technique similar to time-dropout has previously been explored in https://openreview.net/forum?id=Ln8ogihZ2S.\n\nTypos\n- 727:  MMGVAE => MMGPVAE\n- 857: Wilcoxon-signed => Wilcoxon signed-\n- 920: most clear => clearest\n- 981: missing sample => missing samples", "was_revised": true, "final_justification": "The authors have presented a strong and impactful contribution with MRINE, addressing critical challenges in decoding multiscale neural activity in real-time. The rebuttal effectively addressed most of the concerns raised, i.e.,\n  Real-time inference times: Authors provided concrete per-timestep inference times, confirming MRINE's practical applicability for Brain-Computer Interfaces (BCIs).\n Evaluation methodology: The justification for using linear regression for latent factor decoding was well-explained, aligning with standard practices for assessing representational quality.\n  Robustness to missing data: The added results on decoding from reconstructed inputs further highlight MRINE's versatility and strong handling of missing samples.\n  Time-dropout distinction: Clarification on time-dropout as a regularization technique was helpful.\nWhile acknowledged as a common limitation and future work, the assumption of time-invariant dynamics persists, which limits its direct applicability to highly non-stationary neural signals.\n  Multiscale constraint assumption that the discrete modality evolves faster is a current design constraint, though the authors state the model can be extended.", "rebuttal": "We thank the reviewer for their careful investigation of our work. Below, we address their comments.\n\n> Stationarity of neural signals\n\nWe agree with the reviewer that accounting for non-stationarity is important for real-world decoding and discuss this as a limitation in the manuscript. Like many recent models—including MMGPVAE, mmPLRNN, LFADS, and CEBRA—MRINE assumes stationary encoding. Extending MRINE to handle distributional shifts or task-dependent changes is a valuable direction for future work, particularly for long-term or real-time applications. Indeed, developing calibration or alignment methods over time to address time variations is a complementary yet distinct direction of work (Farshchian et al., 2018; Karpowicz et al., 2025).\n\n> Multiscale constraint\n\nWe indeed made an assumption on the multiscale nature of the signals by assuming that $T’$ is a subset of $T$, indicating that the discrete modality evolves faster, in line with the prior work. However, our model formulation can be extended to scenarios where $T’$ is not a subset of $T$. In those scenarios, we would represent both signals in the greatest-common-divisor timescale with appropriate masks to indicate the availability of the signal samples. For instance, if spike and LFP signals were sampled at 20 and 50 ms, we could represent each signal at 10 ms resolution during modeling, unmasking the spikes every other timestep and LFPs every 5 steps.\n\n> Real-time performance\n\nWe thank the reviewer for their comment. We now report the per-timestep inference times for all methods to give an insight into the computational complexities. MRINE’s latent inference time per timestep is $1.82 \\pm 0.13$ ms with Intel i7-10700K processor, which is smaller than the timestep considered in this study (10 ms), and smaller than common brain-computer interfaces (BCI) timesteps (e.g., 50 ms in Stavisky et al., 2015), thus suggesting real-time BCI applicability. Also, the latent inference times per timestep are $1.83 \\pm 0.12$ ms for MMGPVAE, $0.03 \\pm 0.008$ ms for MSID, and $25.91 \\pm 3.44$ ms for mmPLRNN. Overall, we can say that the time complexity of MRINE’s inference is $O(T)$, where $T$ is the the number of timesteps, due to MRINE’s recursive nature. Also, the complexity per time-step is constant due to MRINE’s recursive nature.\n\n> Evaluation methodology\n\nWe thank the reviewer for pointing this out. As detailed in Section 3.2 and L188-211, we perform k-step-ahead prediction of neural signals by first forward propagating the filtered multiscale latent factors ($x_{t|t}$) $k$-steps into the future through $x_{t+k|t} = A^k x_{t|t}$, and then obtain $k$-step-ahead predicted multisale embedding factors via $a_{t+k|t} = Cx_{t+k|t}$. These $k$-step-ahead predicted multiscale embedding factors are then passed through modality-specific decoder networks to obtain $k$-step-ahead prediction of input neural signals. \n\nThat said, MRINE is primarily designed to decode target signals more accurately through multimodal information aggregation from modalities with different temporal resolutions, with real-time inference support. Thus, we mainly focused on downstream target decoding performance rather than short-term or future-forecasting of neural signals. This choice also influenced our design choices. For instance, we preferred linear dynamics instead of nonlinear dynamics or time-varying dynamics. Through this design choice, MRINE can allow for handling missing samples through its internal dynamics and can leverage both real-time Kalman filtering and noncausal Kalman smoothing without requiring to train separate models for different inference modes, whereas prior models support only one of the two.\n\n> Comparison of latent factor decoding\n\nWe chose linear regression decoders to ensure a fair and standardized comparison of latent factor quality across models. This approach is widely used in both machine learning and neuroscience—for example, foundation models are often evaluated via linear probes, and in neuroscience, the Neural Latents Benchmark uses linear decoders to assess representational quality. For the specific example of Lorenz simulations, prior work also used similar linear metrics or alignment methods to quantify the alignment between inferred latent factors and true latent factors (attractors) through simple linear regressors or rotations (Sussillo et al., 2016; Gilpin, 2020).\n\n> Related work on time-dropout\n\nWe thank the reviewer for bringing this excellent work to our attention that we will cite and discuss appropriately in the revised manuscript. While both approaches address missing data, we believe MRINE's use of time-dropout is conceptually distinct. Specifically, we introduce artificial sample dropping during training to improve robustness to missing inputs—analogous to dropout regularization. In contrast, Dowling et al. propose a method that explicitly models missing data during both training and inference, similar in spirit to Kalman filtering (i.e., setting Kalman gain to zero vs. $\\alpha_t$ = 0). We will cite and discuss this work appropriately in the revised manuscript.\n\n> Decoding from reconstructed inputs\n\nWe thank the reviewer for their question. In our main analyses, we decoded behavior from the latent factors, following common practice. As requested, we now also evaluated decoding performance from the reconstructed observations for the NHP grid-reaching task (20 spike + 20 LFP channels). Using reconstructed spiking rates, MRINE achieved a decoding CC of $0.621 \\pm 0.011$, matching the performance from latent factors. As expected, decoding from reconstructed LFP signals yielded a lower CC of $0.490 \\pm 0.012$, likely due to its lower behaviorally relevant information and signal to noise ratio.\n\n> Typos\n\nWe sincerely thank the reviewer for bringing these typos to our attention. We will make sure to correct them.\n\nReferences:\n\n- Stavisky, S. D., Kao, J. C., Nuyujukian, P., Ryu, S. I., & Shenoy, K. V. (2015). A high performing brain–machine interface driven by low-frequency local field potentials alone and together with spikes. Journal of Neural Engineering.\n- Sussillo, D., Jozefowicz, R., Abbott, L. F., & Pandarinath, C. (2016). LFADS - Latent Factor Analysis via Dynamical Systems. arXiv.\n- Gilpin, W. (2020). Deep reconstruction of strange attractors from time series. Advances in Neural Information Processing Systems."}]

Author rebuttal:

Meta-review: In this manuscript titled "Dynamical modeling of nonlinear latent factors in multiscale neural activity with real-time inference", the authors proposed a multi-modal multi-timescale real-time decoding network architecture and loss function for brain-computer interface applications. The method consists of Linear Dynamical Models to account for temporal integration, timescale differences, and predict missing data, then performs nonlinear information fusion via a fusion network. It represents an original solution to an important engineering application that of some interest to the NeurIPS community. The conceptual novelty is in the combination of well-known building blocks and the experiments were somewhat unrealisitically used small number of channels, limiting its application in more realisitic scenario where array recordings are typically much higher dimension these days. However, the authors were able to provide very high-dimensional experiments during the rebuttal. I strongly encourage the authors to produce a significantly updated manuscript that includes the details of the additional experiments.
Decision: Accept (poster)