Dynamical modeling of nonlinear latent factors in multiscale neural activity with real-time inference

We thank the reviewer for their careful investigation of our work. Below, we address their comments.

Stationarity of neural signals

We 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).

Multiscale constraint

We 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.

Real-time performance

We 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.

Evaluation methodology

We 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.

That 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.

Comparison of latent factor decoding

We 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).

Related work on time-dropout

We 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.

Decoding from reconstructed inputs

We 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.

Typos

We sincerely thank the reviewer for bringing these typos to our attention. We will make sure to correct them.

References:

• 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.
• Sussillo, D., Jozefowicz, R., Abbott, L. F., & Pandarinath, C. (2016). LFADS - Latent Factor Analysis via Dynamical Systems. arXiv.
• Gilpin, W. (2020). Deep reconstruction of strange attractors from time series. Advances in Neural Information Processing Systems.

We thank the reviewer for acknowledging the strengths of our problem motivation and experimentation. Below, we address their comments.

Fundamentally new concepts

This 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.

Gaussian modeling of spikes and additional baselines

This 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.

Also, 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.

We 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.

Using few channels

We 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.

Nevertheless, 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.

Real-time compatibility

We 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.

Stationarity of neural signals

We 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).

Decoding R2

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.

NHP grid reaching:

NHP center-out reaching:

Using ranks in Table 7

We 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.

Regime beyond 0.8

We 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.

Difficulty of training the model

Training 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.

Overall, 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.

References:

• 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.
• 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.
• Ahmadipour et al. (2024). Multimodal subspace identification for modeling discrete-continuous spiking and field potential population activity. Journal of Neural Engineering.
• Hsieh et al. (2018). Multiscale modeling and decoding algorithms for spike-field activity. Journal of Neural Engineering.
• Farshchian et al. (2018). Adversarial Domain Adaptation for Stable Brain-Machine Interfaces. International Conference on Learning Representations.
• Karpowicz et al. (2025). Stabilizing brain-computer interfaces through alignment of latent dynamics. Nature Communications.
• Schneider et al. (2023). Learnable latent embeddings for joint behavioural and neural analysis. Nature.

Number of channels

We thank the reviewer for raising these great points. To demonstrate MRINE’s performance beyond 20-20 channel regime per reviewer’s request, we now trained MRINE, single-scale spike and LFP models, MSID, and MMGPVAE in 30-30 and 60-60 spike and LFP channel regimes.

As shown in the table above, MRINE outperformed all single-scale and multimodal baselines in these high-channel regimes. This result shows that MRINE is also advantageous in the high-channel count regime in both the spike-LFP motor dataset as well as the new calcium/neuropixel dataset. We sincerely thank the reviewer for their great suggestion.

That being said, below we provide reasons why performance comparisons in the 20-channel regime are also valuable for future BCI research:

1. Viable recording channels drop over time post-implantation: While we agree that some chronic BCI systems deploy 96-channel arrays, the number of reliably recording channels post-implantation is often much lower. For example, the BrainGate study (Harris and Tyler, 2013) reported that the number of active spike units can drop to approximately 24 units over time. Similarly, Barrese et al. (2013) identified several electrode failure modes and found that the percentage of viable recording channels can fall to as low as 5–24%. These observations reflect real-world constraints that show the importance of lower-channel regimes.

2. Spike signals degrade faster than LFPs: Beyond hardware degradation above, signal stability differs between modalities. Prior work has shown that LFP signals exhibit significantly greater stability over time compared to spikes, with spikes degrading faster (Flint et al., 2016; Wang et al., 2014; Heldman and Moran, 2020; Sharma et al., 2015). This can result in a substantial reduction in usable spike channels over time. In such scenarios, robust multimodal methods like MRINE can be particularly valuable by aggregating information from a limited number of spike channels with the more abundant, stable LFP signals (please note the performance difference between single-scale LFP model and MRINE). Indeed, multimodal methods can make a big difference to implanted patients in these low-channel regimes: they allow these patients to keep using their implanted BCI for longer.

3. Some devices can have lower channel count: Recent developments in wireless BCI systems have explored lower channel count designs to reduce bandwidth and power demands. For example, Uran et al. (2022) introduced a neural SoC system using 16 channels with on-chip feature extraction, illustrating that such configurations remain relevant in certain application contexts.

Drop beyond 0.8

Per reviewer’s request, for the NHP grid-reaching dataset and for MRINE, MMGPVAE, and MSID, we computed the neural reconstructions in Fig. 6 also when the modality of interest is completely missing but reconstructed from the other modality. In this scenario, reconstruction quality significantly decreased for all models, for instance, MRINE, MSID and MMGPVAE achieved 0.592, 0.584, and 0.571 AUC for spikes and 0.244, 0.244, and 0.041 CC for LFP signals. We believe that this is expected since neither of these models’ training objectives is cross-modal reconstruction, but rather on self-prediction or self-reconstruction by aggregating information. To achieve competitive cross-modal reconstruction, the training objective should be modified accordingly, and we will include this limitation (that is not only inherent to MRINE) in the discussion paragraph.

References:

• Harris and Tyler (2013). Biological, Mechanical, and Technological Considerations Affecting the Longevity of Intracortical Electrode Recordings. Critical Reviews in Biomedical Engineering.

• Barrese et al. (2013). Failure mode analysis of silicon-based intracortical microelectrode arrays in non-human primates. JNE.

• Flint et al. (2016). Long-Term Stability of Motor Cortical Activity: Implications for Brain Machine Interfaces and Optimal Feedback Control. JNE.

• Sharma et al. (2015). Time Stability and Coherence Analysis of Multiunit, Single-Unit and Local Field Potential Neuronal Signals in Chronically Implanted Brain Electrodes. Bioelectronic Medicine.

• Heldman and Moran (2020). Chapter 20—Local field potentials for BCI control. In N. F. Ramsey & J. del R. Millán (Eds.), Handbook of Clinical Neurology (Vol. 168, pp. 279–288). Elsevier.

• Wang et al. (2014). Long-term decoding stability of local field potentials from silicon arrays in primate motor cortex during a 2D center out task. JNE.

• Uran et al. (2022). A 16-Channel Neural Recording System-on-Chip With CHT Feature Extraction Processor in 65-nm CMOS. IEEE Journal of Solid-State Circuits.

We thank the reviewer for providing an excellent summary of our work and its strength. Below, we address their comments:

Small number of channels

This 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.

Nevertheless, 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.

Decoding performance on R2:

Per 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.

NHP grid reaching:

NHP center-out reaching:

We will also include example decoded trajectories in our manuscript.

Broader application of the method

We 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:

Indices in Fig. 1b

We 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.

MLP in Fig. 1b

The 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.

Slow/fast scales on smoothness regularization

This 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.

References:

• 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.
• 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.
• 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.
• 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.
• Schneider, S., Lee, J. H., & Mathis, M. W. (2023). Learnable latent embeddings for joint behavioural and neural analysis. Nature.

I thank the authors for answering my questions and doing additional experiments on calcium data.

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.

Per the reviewer’s request, we now provide the 20 Spike - 20 LFP results in Table 1 using R2 below.

I thank the authors for being candid about these limitations. However, I still think the choice of dataset may not be appropriate for your model and does not reflect the realistic settings. Therefore, I will keep my score.

We thank the reviewer for their questions and comments. Please see our answers below:

Meaning and rationale of ‘multiscale embedding factors’ $a_t$

The 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).

After 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$.

Although $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.

Further, 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.

Imputation of missing samples

We 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.

While 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).

Computational complexity

We 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).

Independence assumption

During 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.

$T$ and $T’$

This 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.

Last MLP in Fig. 1

The 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.

Comparisons in other tasks

Thank 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.

As 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.

KL divergence term as in SAEs

Thank 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.

During 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.

To 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.

Overall, 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.

References:

• 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.
• 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.

We thank the reviewer for their thoughtful feedback and for highlighting these important clarifications. Below, we address the reviewer’s comments:

Hyperparameters

We 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.

Meaning of primary

We 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.

About Figs 2a and 3a

The 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.