Venue: ICLR 2021 Conference

Reviews:

final_score: rating: 5: Marginally below acceptance threshold
review_content: This submission extends a previous nonlinear data analysis method ( Latent Factor Analysis of Dynamical Systems, LFADS) to deal with the calcium imaging data (CaLFADS). 

Quality: Developing appropriate methods to handle the calcium imaging data is an important question as calcium imaging becomes popular. This submission represents one such attempt. Although some of the results are promising, unfortunately the work is rather incremental. Furthermore, the robustness of the method (as well as the comparisons to alternative methods)  under more realistic noise regime is not carefully tested. Also, the improvement over simple baseline model is small.
Clarity: The writing is overall good, although the paper would benefit by more clearly describing and motivating the assumptions/approximations of the algorithm ( see later for a few more detailed comments).
Originality: The method for handling the calcium traces is novel, although it is a relatively straightforward modification of LFADS.
Significance:  As there is only small improvement over simpler, alternative methods, I feel the significance is relatively low. Also, no new insights were shown by applying these methods compared to simple alternatives.
Pros: 
1. CaLFADS is a holistic approach to infer the latent dynamics from fluorescence traces and avoids the two-step procedure which is nice. 
2. The method can obtain measures of the uncertainty- an issue that was ignored in most of the previous work on calcium data analysis.
3. The method shows some improvement over a few alternative methods when fitting the data in certain regimes.
Cons:
1. The work is rather incremental. The improvement over the alternative models is small.  This could be seen in Table 1, Fig S2, S3 and Fig. 3 . In Table 1, the margin between the different modes is small. I worry that, in most scientific applications, the small improvement of CaLFADS over Gaussian-LFADS and OASIS-LFADS is unlikely to matter. Take the Lorenz Attractor for example. If the 3-d state-space reconstruction is plotted, the basic structure of the Lorenz Attractor will be recovered for all three models based on other results shown there. It appears that no additional insights are revealed by CaLFADS.
2. The robustness of the method regarding parameter mis-specification and realistic measurement noise is not evaluated. The simulated calcium traces shown in the figures suggest that authors were focusing a high SNR regime- a regime that is very rare practically for large-scale calcium imaging. For example, real fluorescence traces when imaging hundreds of neurons are much noisier than the traces shown in Figure S2. For this to be more convincing for real applications, I would strongly suggest the authors to increase the noise when simulating the fluorescence traces, and test and compare the different models. I am concern that the results may change dramatically in the more realistic noise regime.
3. It is unclear whether OASIS-LFADS is worse than CaLFADS  is simply due that OASIS-LFADS method did not carefully model the noise characteristics of the deconvovled spikes from calcium. For the CaLFADS, the authors mentioned that they used a continuous relaxation to transient amplitudes inspired by a zero-inflated gamma model by Wei et al 2019. While whether this additional ingredient helped with the performance of CaLFADS is unclear, I wonder if the authors had similarly incorporated this zero-inflated gamma model to the OASIS-LFADS method to make it a more fair comparison. 
4. Another potential concern is that, practically, it might be difficult to use as the proposed method is complicated and seems to require a fair amount of hyper-parameter tuning. If methods like PCA could already reveal the key insight from a dataset (e.g., Fig. 3), one may question how much value this complicated nonlinear method adds.

More specific comments:
* How well does the Gaussian-LFADS method perform on the real data (i.e., mouse primary visual cortex in Sect. 3.3)? Is it significantly worse than CaLFADS? I am puzzled why this is not reported given Gaussian-LFADS  is one of the main alternative models used in the simulated data. Also, it seems that CaLFADS only slightly outperforms PCA (which is a very basic baseline) for the real data, which again is concerning. 
* The parameter of the AR(1) is inferred from a separate model (Friedrich et al 2017). It is a bit disappointing that the model could not infer these parameters directly, as one would think that modeling the underlying dynamics of the computation would further constrained the AR(1) model and help with the deconvolution process (I had considered this as one potential advantage of this holistic approach, but unfortunately it didn’t quite work out). How much do the results depend on the accuracy of these parameters? This is relevant, because practically Friedrich et al 2017 may not perform particularly well in recovering the AR(1) parameters when perform large-scale imaging. Thus it would be useful to check the robustness of the proposed method.
* Related to the last point- when fitting the simulated data, do the authors use the ground truth parameters from the AR(1) model? If that’s the case, it would be useful to test if the model is robust to the model mis-specification, and check if the CaLFADS model still output alternatives in this regime.

Clarification questions: 
* “The approximation results from the use of an L1 norm on the spikes,…a range that is poorly handled by the log-gamma function. ” This is a bit confusing. What “data” were referred to here? Inferred “spikes” from calcium traces or the spike train data from electrophysiology?
* The variables in Fig1 Panel B are not explained in the figure caption.
The following paper is relevant as they also proposed a holistic approach to infer the underlying neural structure directly from the raw fluorescence trace.   Aitchison, L., Russell, L., Packer, A. M., Yan, J., Castonguay, P., Hausser, M., & Turaga, S. C. (2017). Model-based Bayesian inference of neural activity and connectivity from all-optical interrogation of a neural circuit. In Advances in Neural Information Processing Systems (pp. 3486-3495).
*The paper by Hernandez et al (2020) is also related to this paper. They proposed methods to model the neural dynamics based on the ephys and calcium data.    Daniel Hernandez et al., Nonlinear Evolution via Spatially-Dependent Linear Dynamics for Electrophysiology and Calcium Data, the NBDT journal, 2020.
* “it would ignore the population dynamics that could be informative for separating out the calcium dynamics from the higher-order computational dynamics. ”  It is unclear what this means. Please clarify the idea here, in particular what the high-order computational dynamics means. There also seems be an issue of timescale here, which the paper didn’t touch on. If the neural dynamics is much slower than the calcium dynamics, it does not seem to matter to much in terms of whether calcium dynamics is carefully modeled. 
* The last section before Section 2.1, “The hidden state of this GRU is…to provide an approximation to spike counts”. This needs to better explained. 

%%%% after the rebuttal:

I would like to thank for the authors for their effort to address my concerns. The manuscript is now improved, and I am raising my score to 5.  After considering other discussions/comments, overall I still think the manuscript is below the threshold. 

review_title: An extension of a nonlinear method for analyzing the neural dynamics of calcium imaging data 

final_score: rating: 7: Good paper, accept
review_content: ##########################################################################

Summary: 

Traditionally, the problem of disentangling neural activity from calcium dynamics *and* inferring latent variables from neural activity have been treated separately; the neural activity is first inferred from calcium traces (using a variety of different approaches) and then data analysis is performed on the inferred neural activity. Instead, the authors tackle both of these problems **jointly** by learning a hierarchy of dynamical systems where one is responsible for modeling the calcium dynamics which is driven by another latent system that is responsible for modeling the low-D latent dynamics underlying the population activity; the approach is titled CaLFADS which is a novel extension to both LFADS and variational ladder autoencoder (VLAE) to sequential data.

##########################################################################

Pros: 

I thought the proposed approach was super cool! To the best of my knowledge, this is the first time I have seen a hierarchical VAE applied to sequential data which is very exciting and very non-trivial. The empirical results presented on both synthetic and real data are also great as the authors demonstrate the ability to disentangle the neural factor from the calcium dynamics.

##########################################################################

Cons: 

My biggest con is that the paper is not clear enough. As I mentioned previously, this work is a non-trivial extension to not only LFADS but **also** to hierarchical VAEs. Sadly, not enough information was given explaining the details. First, the authors should have at least provided a brief introduction to variational ladder autoencoders. 

Second, the mathematical form of the generative model does not appear anywhere in the manuscript. Instead, it is explained using figure 1 and the accompanying paragraphs (the first two paragraphs in section 2) but in my opinion, these explanations are opaque without equations accompanying them. For instance, in the original LFADS paper all the variables were a deterministic transformation of the latent variables (the initial condition and the inferred inputs); this allows the system to be completely characterized by these two groups of latent variables. From reading the two paragraphs in section 2, I thought this was also the case for CaLFADS but once I looked at the algorithm table in the appendix I found that $g_{1,t}\sim \mathcal{N}(\mu_{g_{1,t}}, \sigma^2_{g_{1,t}})$! (Note that this stochastic mapping prevents the system to be completely characterized by just the initial conditions and the inferred inputs as now there is uncertainty over $g_{1, t}$ which effectively makes $g_{1, t}$ a latent variable.)

Last, I am still confused about the functional form of the likelihood relating fluorescence to the latent calcium dynamics. In algorithm table 1 line 28 it says  $x_{t} = \textrm{DyeKinetics}(x_{t-1}, z_{1, t})$ but  $\textrm{DyeKinetics}$ is not defined. Towards the end of the second paragraph of section 2, it says “While it is possible to transform calcium dynamics into fluorescence using nonlinear models of fluorescence indicator kinetics (Speiser et al., 2017; Deneux et al., 2016), we found that using the AR(1) process to model the transformation of spike counts to fluorescence was sufficient for accurate reconstruction of the latent space, even in synthetic data with nonlinear calcium transient generation.” Does this mean that $p(x_t \vert z_{1,t}) = \mathcal{N}(Az_{1, t}, \sigma^2)$?

##########################################################################

Questions/Comments:

1) In section 2.1 the authors state that an extra term is added to the loss function that enforces the inferred spike counts, $u_{1, t}$ to follow the dynamics of the higher-order system $z_{2,t}$ through the rate function, $\lambda_t$ i.e. $\mathbb{E}[\log p(u_{1, t} \vert z_{2,t})]$. This seems to be at odds with the VLAE approach but would be enforced if a traditional hierarchical VAE was utilized. Could the authors elucidate why a VLAE with this added regularizer was chosen as opposed to a traditional hierarchical VAE? If a VLAE doesn't work as well without this regularizer it seems like this is a signal that a traditional hierarchical VAE should be used.

2) Just like LFADS, the control inputs are a function of the previous latent factors, thus the variational distribution should depend on the previous control inputs and the initial condition i.e. $Q(u_{2, t} \vert g_{2, 0}, u_{2, < t})$.

3) it would be nice if the authors could show that CaLFADS could do prediction as it would showcase that the model is not only good at inferring the latent factor but is also capable of generating data. One easy way to do this is given N trajectories of length T, 80% of the trajectories are used for training, and then to estimate the other 20% the model generates forward. While I know that the model is evaluated on a held-out test, the initial condition for the system is determined by the full trajectory. Thus, it would be interesting to see how the model works when only a partial length of the trajectory is given.

4) At the bottom of page 16 above equation 10, $gamma$ -> $\gamma$.

5) Many of the figures are not vector graphics.
##########################################################################

To reiterate, I am a big fan of this work and I think it is a very cool contribution. Sadly, more work needs to be done on the structure of the manuscript, specifically on details of the proposed approach. I suggest the authors cut down the intro and use the extra page to really flesh out the details. I will be more than happy to raise my score if the authors do this! 

##########################################################################
Post-rebuttal: I have increased my score to a 7. I have laid out my reasoning in response to the author's comments.

review_title: The authors introduce a novel extension of hierarchical VAEs to sequential data where in this work they focus on data obtained from calcium imaging. While the proposed approach is great and the results are very promising, the structure of the paper needs work as many details are not explained well or missing.

final_score: rating: 5: Marginally below acceptance threshold
review_content: This paper, as I understand it, extends a very nice method for Latent Factor Analysis of Dynamical Systems (LFADS) from spike train data to calcium imaging data.  I am familiar with the LFADS paper, as well as much of the literature on calcium imaging for spike inference, and latent dynamical systems more generally.  I am in no way an expert on anything deep learning, nor nonlinear dynamics such as Lorenz attractors.  I had trouble building an intuition about the model of neural dynamics, the approach, and the results, possibly due to the presentation style, which I was not used to, and I explain below.

1. I am very comfortable with equations, including explicit generative models of the calcium, and inference/encoder models as well.  The LFADS paper, for example, explained both quite clearly.  In this paper, however, I could not find the relevant equations, other than those for the synthetic data generator. I struggled to understand therefore, precisely what was the implicit model, which was informally described in the 2nd paragraph of the paper.  Similarly, the encoding model was described informally in section 2, and with Figure 1, and Table S1 and Figure S1 helped, but I really wished there were equations.   The LFADS paper has equations that I understand, so I guessed one possibility was that this generative model generalized that by adding another random variable "on top" of the spikes to characterize the calcium activity. Then I would have guessed that the encoder/inference models would be modified similarly, but I could not tell.  Algorithm 1 helped for the forward method, but even there some of the notation was not clear (e.g., what precisely does 'Dense' mean?).

2. The discussion around lower-order and higher-order dynamics confused me a bit as well, as I could not find formal definitions of those terms.  I guessed that lower-order dynamics referred essentially to first-order markov chains, and second-order dependencies counted as higher-order, but that logic didn't quite work throughout.  If parsing the lower vs higher order dependencies is a key desiderata of this method, it would have been nice to explicitly see the approach doing that in the results. 

3. The results in Table 1 confused me a bit.  Why are three significant digits shown for R^2, there were 20 different seeds, so precision down to 1/1000 does not seem possible.  Also, I lack intuition as to what the different true estimated latent states look like.  Do I care about 0.005 better R^2? I would think not.  t-tests are not obviously valid tests here either.  

4. Figure 2 was very nice and clear and helpful.  It clearly shows that on the synthetic data, CaLFADS does something closer (in some qualitative, though not quantiative) sense, than OASIS+LFADS.  I did wonder, however, OASIS is a fully online algorithm, whereas LFADS and CaLFADS are batch algorithms.  So, if I understand correctly, OASIS is it a large disadvantage because it cannot use data after each spike to help with the inference.  This is particularly important for fast dynamics, because the post-spike effect can be much more informative than the pre-spike information.  So, I would think that a more fair approach would be to use a batch, rather than online, comparison. 

5. Figure 3 was clear enough about recall.  I wonder, from the experimentalist perspective, if they can get 0.8 simply by operating on the raw data, and can get 0.82 by running this complicated model, would they?  I guess the answer depends on how difficult and time-consuming it is for them and their computer to run this, but no information is provided about that which I can find.  I looked through the code briefly, and it seems maybe easy enough, but I didn't see a tutorial for using it on my data.  Also, one of the notebooks somewhat un-anonymizes you because it includes direct paths on your computer, rather than relative paths :)

To provide some concrete recommendations:

1. I'd like more equations, the more the merrier for me, could be in the appendix.
2. I'd like figures showing the data, both real and synthetic, as well as the estimated latent states, so we can understand what the summary of the errors means. 
3. I'd like more explicit clarity on any changes you made to LFADS other than adding a layer on top, what precisely those changes were, what precisely motivated those changes, and whether they were worth it.

In conclusion, I think this paper adds a relatively incremental addition to LFADS, which could potentially be interesting to neuroscientists, but I failed to clearly understand the ingenuity of the approach, or the impact of the empirical results. 

===============================================================
==== updates after rebuttal period =================================

1. The authors added more equations, I now understand the methodological contributions required substantial effort.  The VAE/etc. work is outside my area of expertise.
2. The authors added some figures showing both real and synthetic data. I've looked through them, and the statistics carefully.  Personally, reading through all of this, I would use PCA+OASIS rather than this new technique.  I understand R^2 is higher for this technique, but looking at the data, without actual confirmation of any of the inferences. There exists data with joint calcium imaging and ephys one can use to calibrate and evaluate methods, those are the best data to use for such purposes.  In the absence of those data, we are guessing.  My guess is that the OASIS+LFADS output looked better than the CaLFADS output. R^2 is a funny metric, especially when we are trying to compare spike trains, because translations utterly break R^2, unless the output is smooth, in which case, it does not.  It looked to me that CaLFADS had a smoother output than OASIS+LFADS, which could make R^2 higher, but is not actually what I typically want when analyzing calcium imaging data. Indeed, if we were ok with smoothing, we could simply operate on the calcium imaging data itself, without worrying about spike train inference.
3. I got more clarity on the differences between LFADS and CaLFADS, though again, this is outside my area of expertise. 
review_title: Potentially very interesting, but I lack clarity on certain key issues

final_score: rating: 4: Ok but not good enough - rejection
review_content: ------------- UPDATE FOLLOWING DISCUSSION -------------

The authors were responsive and suggested a series of updates to address the concerns raised here. I find the inclusion of a theory section valuable. Overall, though, my stance did not change: I’d like to encourage the authors to pursue their idea, but I don’t think the paper is ready for acceptance yet. I’d recommend focusing on two aspects:
(i) Simulations and numerical results are especially confusing. Please consider a major revision in light of the comments raised by different reviewers.
(ii) I find the application to fMRI interesting and potentially impactful. If the authors can include such analysis, I think it’d improve the paper significantly. Being a completely different experiment, it would also alleviate the concerns on numerics. (significance, effect size, stability)

------------- END OF UPDATE -------------

Summarizing the dynamics observed via high-dimensional recordings of neuronal activity is an important problem in neuroscience. Imaging the calcium dynamics of neurons via calcium binding fluorescent molecules, rather than recording the spiking electrophysiological activities, enables higher throughputs and has become popular. This, however, creates the additional step of inferring the underlying spikes from the smooth calcium activity recordings. This paper addresses the problem of inferring the latent dynamics directly from calcium activity recordings by jointly optimizing for both steps with a sequential ladder variational autoencoding architecture.  

The paper demonstrates (somewhat small) improvements over existing methods, and it may be useful for neuroscientists. However, I do not support the acceptance of this submission due to the following major concerns:  

--- Limited methodological novelty: The approach seems to be a direct application of a particular sequential hierarchical autoencoder to this particular calcium-imaging problem. There is limited methodological discussion. In this sense, the manuscript may be a better fit for a neuroscience journal.  
--- Focusing on a single domain: The experiments cover a narrow domain for an application-focused paper. The authors do not provide example use cases beyond calcium imaging in neuroscience.  
--- Comparison with the baseline method (OASIS+LFADS): OASIS was designed for fast, online spike inference from calcium imaging data.  
     a) How much time is needed for OASIS+LFADS vs CaLFADS? If OASIS+LFADS is faster, could more complicated settings/architectures be used for OASIS/LFADS to improve performance? (e.g, lag window size in OASIS, recurrent network size in LFADS)  
     b) The OASIS paper mentions different working modes, where only some of the possible parameters are optimized. How was OASIS run in this manuscript?  
--- Exposition: The paper focuses on neuroscience-oriented readers only, and may not appeal much to the more general ML crowd. For instance, I thought it would help a lot to display the formal probabilistic model using equations.  

Other concerns:  
--- The paper refers to the different dynamics within the hierarchy as higher-/lower-order dynamical systems. The order of a dynamical system has an established definition in the literature already. Unless the authors cite previous such usage, I suggest changing this potentially confusing terminology.  
--- C. Pseudocode: (i) T\times N, (ii) lines 5-6, right hand side doesn’t depend on l. (iii) Variables aren’t properly defined. e.g., What is E_{l,t}^{bwd}?  
--- What does 14 represent in t_{14} in Table 1?
review_title: Application of sequential hierarchical autoencoders to calcium imaging in neuroscience

Meta-review: Inferring latent trajectory from noisy Ca time series is an important and timely problem and the current study shows some progress in the inference problem. Although the proposed model has some originality, there are remaining issues rendering the manuscript not ready for publication yet. Reviewers raised issues on readability, lack of details, statistics of experiments, Ca time constant estimation, effect size, and lack of comparison. Through an extensive discussion and revisions, I'm happy to see the manuscript was greatly improved in readability and additional statistics were provided. However, the Gaussian-LFADS' performance at exactly chance raises red flags, effect size is small, and the significance of the scientific findings remain weak. The model is presented as a variational ladder autoencoder system with 2 layers, but the shallow latent representation is tied to the continuous approximation of the point process likelihood. Hence, I view the model as an extension of LFADS rather than a flat-hierarchical VAE.

Overall, this paper has a potential of becoming a solid contribution for statistical neuroscience, once the above shortcomings are addressed.
Decision: Reject