CaLFADS: latent factor analysis of dynamical systems in calcium imaging data
Authors: Luke Yuri Prince
Shahab Bakhtiari
Colleen J Gillon
Blake Aaron Richards
Venue:  (2021)
Keywords: latent variable modelling
lfads
neuroscience
variational autoencoders
dynamical systems
calcium imaging
neural data analysis

Abstract: Dynamic latent variable modelling has been a hugely powerful tool in understanding how spiking activity in populations of neurons can perform computations necessary for adaptive behaviour. The success of such approaches has been enabled by the ability to construct models derived with the characterization of spiking activity as point-processes since spiking dynamics occur on a much faster time-scale than the computational dynamics being inferred. Other experimental techniques, such as calcium imaging, pose a problem for latent variable modelling of computational dynamics, since the time-scales of calcium dynamics and computational dynamics overlap. As such, the success of dynamic latent variable modelling in calcium imaging data rests on being able to disentangle the contribution of these two sources of variation. Here we extend recent advances using variational autoencoders to analyze neural data, by incorporating a ladder architecture that can infer a hierarchy of dynamical systems. Using built-in inductive biases for calcium dynamics, we can capture calcium flux as well as underlying dynamics of neural computation. First, we demonstrate with synthetic calcium data that we can correctly infer an underlying Lorenz attractor at the same time as calcium dynamics. Next, we show that we can infer appropriate rotational dynamics in spiking data from macaque motor cortex after it has been converted into calcium fluorescence data via a calcium dynamics model. Finally, we show that our method applied to real calcium imaging data from primary visual cortex in mice allows us to infer latent factors that carry salient sensory information about unexpected stimuli. These results demonstrate that variational ladder autoencoders are a promising approach for inferring hierarchical dynamics in experimental settings where the measured variable has its own slow dynamics, such as calcium imaging data, thereby providing the neuroscience community with a new analysis tool for a wider array of data modalities.

Under review as a conference paper at ICLR 2021
CALFADS: LATENT FACTOR ANALYSIS OF DYNAMI- CAL SYSTEMS IN CALCIUM IMAGING DATA
Anonymous authors Paper under double-blind review
ABSTRACT
Dynamic latent variable modelling has provided powerful tool for understand- ing how populations of neurons compute. For spiking data, such latent variable modelling can treat the data as set of point-prccesss, due to the fact that spik- ing dynamics occur on much faster timescale than the computational dynamics being inferred. contrast, for other experimental techniques the slow dynam- ics governing the observed data are similar in timescale to the computational dy- namics that researchers want to infer. An example of this is in calcium imaging data, where calcium dynamics can have time timescales of milliseconds. As such, the successful application of dynamic latent variable modelling to modalities like calcium imaging data will rest the ability to dis- entangle the higher- and shallower-level dynamical systems' contributions to the data. To-date, no techniques have been developed to directly achieve this. Here we solve this problem by extending recent advances using sequential variational autoencoders for dynamic latent variable modelling of neural data. We solve the problem of disentangling higher- and shallower- dynamics by incorporating a ladder architecture that can infer hierarchy of dynamical systems. Using some built-in inductive biases for calcium dynamics, we show that we can disentan- gle calcium flux from the underlying dynamics of neural computation. First, we demonstrate with synthetic calcium data that we can correctly disentangle an un- derlying Lorenz attractor from calcium dynamics. Next, we show that we can infer appropriate rotational dynamics in spiking data from macaque motor cortex after it has been converted into calcium fluorescence data via a calcium dynam- ics model. Finally, we show that our method applied real calcium imaging data from primary visual cortex in mice allows us to infer latent factors that carry salient sensory information about unexpected stimuli. These results demonstrate that variational ladder autoencoders are a promising approach for inferring hier- archical dynamics in experimental settings where the measured variable has its own slow dynamics, such as calcium imaging data. Our new, open source tool, thereby provides the neuroscience community with the ability to apply dynamic latent variable modelling to a wider array of data modalities.
the order of hundreds
INTRODUCTION
Dynamic latent variable modelling has been a hugely successful approach to understanding the func- tion of neural circuits. For example, it has been used to uncover previously unknown mechanisms for computation in the motor cortex (Churchland et al. 2012 Pandarinath et al. 2018). somatosen- sory cortex (Wei al. 2019b). and hippocampus (Chaudhuri et al., 2019). However, the success of this approach is largely limited to datasets where the observed variables have dynamics whose timescales are much faster than the dynamics of the underlying computations. This the case, for example, with spiking data, where the dynamics governing the generation of individual spikes much faster than the dynamics of computation across the circuit. This makes it possible to charac- terise the observed data, e.g. the spiking data, set of point-processes that be used directly for inferring latent variables.
However, many datasets in the life sciences are generated by hierarchy of dynamical systems, wherein the shallower-leve dynamical systems that directly generate the observed data have tempo- ral dynamics whose timescales overlap with that of the deeper-level dynamical system to be inferred.

---

Under review as a conference paper at ICLR 2021
A clear example of this is in-vivo calcium imaging data, which is widely used in neuroscience. Thanks to advances in imaging technology and genetically encoded calcium indicators, calcium imaging enables monitoring of the activity of large populations of genetically targeted neurons in awake hehaving animals (Yang & Yuste 2017 Lin & Schnitzer, 2016). However, calcium imaging introduces an additional layer of relatively slow dynamical system between the computations oc- curring in the brain and the measurements that neuroscientists make. This problem is outlined in Figure DA in which calcium fluorescence observations, 2, depend on the state of calcium flux, 21, which is governed by a shallower-level dynamical system with temporal dynamics on the order of hundreds of milliseconds. These dynamics are driven, in part, by perturbations due to spikes, U1, which are themselves governed by computational state, 22, with similar timescale in its dynamics to the calcium flux (and which itself can be perturbed independently by unknown inputs, u2 that may also have slow dynamics). Due to the overlap in timescales, it is impossible to identify immediately which components of the calcium fluorescence data are driven by the dynamics of calcium flux, and which are driven by the deeper-level latent dynamics of neural computation. Ideally, neuroscientists would have a method for inferring both the shallower-level calcium dynamics and the deeper-level computational dynamics in order to uncover the hierarchical dynamical systems that generated their data. Such tool would be significantly benefit the systems neuroscience community.
Currently, these problems are treated separately. For situations where the observed data can be treated as a point process, we have good techniques for inferring the deeper-level dynamics. For example, recent applications of sequential variational autoencoders have seen great success in in- ferring underlying computations from extracellular spiking data (Pandarinath et al. 2018)) This technique, named Latent Factor Analysis of Dynamical Systems (LFADS), has improved neuro- scientists' ability to infer underlying neural computations from spiking data, e.g. it has been used to identify latent rotational reaching dynamics and decode reaching behaviour of macaques and humans with higher fidelity than other techniques.
Although LFADS has significantly advanced our ability to analyze neural data in the form of spike trains, it does not address the problem highlighted above for calcium imaging, wherein the calcium dynamics introduce an additional shallower-level dynamical system whose timescale overlaps with the timescale of neural computation. Theoretically, this problem could be solved independently by first inferring spikes from calcium data, whether by deconvolution (e.g., OASIS) (Friedrich et al.) 2017), variational inference (e.g.. DeepSpike) (Speiser et al. 2017), dynamic programming (Deneux ct al. 2016), or any other method (Berens et al.) 2018 Pachitariu et al., 2018, Evans et al.) 2019 Wei et al.) 2020), and then applying LFADS. However, this approach treats each calcium trace as a completely independent variable when inferring calcium dynamics. This ignores correlations in population activity that inform the separation of calcium dynamics (which are independent of popu- lation activity) from computational dynamics (which are not independent of population activity). If this separation is sub-optimal, then inference of the deeper-level system will be impaired.
Here, we address this problem by extending LFADS with a variational ladder autoencoder architec- ture (Zhao et al.) 2017) that folds the calcium dynamics inference into the larger inference problem (Fig. [B). Our system, CaLFADS, incorporates inductive biases for calcium dynamics and, thanks to the ladder architecture, is able to infer the deeper-level dynamical system better than an approach that treats inference of calcium dynamics and deeper-level dynamics as separate problems. For further discussion on how CaLFADS relates to previous work, see Appendix B
First, we show using synthetic data that we are able to reconstruct ground-truth latent dynamics from synthetic calcium traces. Next, we apply CaLFADS to spiking data from macaque motor cortex that has been converted into calcium fluorescence traces using a calcium dynamics model. We show that we are able to recover rotational dynamics from this "calcified" data just as LFADS identifies rotational dynamics from spiking data. Finally, we show using real 2-photon calcium imaging data from mouse primary visual cortex that CaLFADD can identify deeper-level latent factors that carry information about unexpected visual stimuli. Altogether, our work shows the benefits of incorpo- rating the calcium dynamics inference procedure into the larger inference problem. It also provides neuroscientists with new, open-source tool for analyzing calcium imaging data in order identify deeper-level dynamics. Given the importance of calcium imaging to modern systems neuroscience, we believe that CaLFADS will be a very useful analysis tool for the community. Furthermore, our CaLFADS could be adapted to other neuroscience data modalities, such as fMRI data which, like calcium imaging, comprises both fast deeper-level latent brain dynamics and slower shallower-level blood oxygen measurement dynamics (Ollinger et al. 2001). Beyond this, believe that our ap-
2

---

Under review as a conference paper at ICLR 2021
proach may be applicable to other life sciences domains, for example insurance claim modelling, where better identifying the relatively slow dynamics underlying the observed variables, such as claim submission distributions over time, could improve forecasting accuracy (Schmidt 2014).
A
④ E)
fi(21, u1 (22))
2 f2(22, (62)
Ca2+
wummmmmn
Ca 2+
X
h2
h2
Inference network
Q(z|x)
o²
μ
o²
μ
U2121
Generative network
92.t
U1,0
Gaussian Likelihood
Poisson
Likelihood
&
91.6
Figure A) Hierarchy of dynamical systems (Top). schema of calcium and Lorenz dynamics (Bottom). B) Schema of our hierarchical model. Latent dynamics model in blue, calcium dynamics model red.
MODEL
CaLFADS extends the variational ladder autoencoder (VLAE) approach (Zhao al., 2017 with RNNs (Fig (B).. which enables inference of hierarchical latent dynamical systems The full di- rected acyclic graph for CaLFADS can found in Figure SI The VLAE approach inferring hierarchies removes the sequential dependence of each latent variable, and instead loads parts of flat latent code with more abstract features using more expressive neural networks (see appendix C for more background on VLAEs). For CaLFADS, this means we separate the shallower-level infer- ence and the deeper-level inference into parallel pathways for inference and generation ( orange arrow), and subsequently, recombines these levels in the generative network [13, pink arrow). This approach has two major advantages: 1) solves the problem of latent variable collapse in which all latent features are loaded onto the lowest-order variables in the hierarchy. 2) It retains very sim- ple formulation for the variationallevidencee lower bound (ELBO) that is easily reparameterisable
Sønderby et al. 2016 Zhao et al. 2017).
'Supporting code for this submission be found https: 4open wcrnnnscceeccoe en repository//3e56001-33 -9dc-- a55980c2al

---

Under review as a conference paper at ICLR 2021
CaLFADS infers the dynamics and unknown events (e.g. spikes) that determine the calcium fluores- cence signal (Fig BB, red network nodes), and the underlying computational dynamical system that controls these events (Fig BB. blue network nodes), through successive encoding and decoding lay- ers in which latent variable distributions are modelled in separated 'rungs' of the ladder architecture. This VLAE architecture enables CaI.FADS to disentangle the shallower-level calcium dynamics and the deeper-level computational dynamics. Algorithm I provides details on the forward pass of CaL- FADS. Additionally, we provide complete description of the model in Appendix D Briefly, layers of Gated Recurrent Units (GRUs) create a deterministic hierarchical embedding of the data in the inference network (Fig [ above the dashed line). These embeddings are then concatenated with time-shifted, layer-specific deterministic projections of the latent state. This concatenated set of activities is then linearly transformed into the parameters of a factorised Gaussian posterior distribu- tion for the latent variables in the generative network [ B, below the dashed line). This process provides the inference network with information about the history of the latent states via feedback connections BB. purple feedback arrows).
As in any variational autoencoder (VAE) system, latent states are sampled from the posterior distri- butions determined by the inference network. For the inference of the computational dynamical sys- tem, these latent states represent the initial conditions (92.4)) and unknown inputs (W22.) to a learned dynamical system modelled discretely by GRU in the generative network (g2,z). The hidden state of this GRU is then linearly transformed into a set of potentially low-dimensional time-varying la- tent factors (ft). Concurrently, the latent states for the calcium dynamics (2t) are concatenated with the dynamic factors from the layer above (ft, and undergo a linear transformation followed by peetified-epponentie non-linearity f(x) ReLU (exp(a)) 1) to provide an approximation to spike counts (#1,e). This non-linearity was chosen as a continuous relaxation to discrete calcium transient amplitudes inspired by the static distribution of amplitudes extracted via OASIS (Wei et al., (019a)). In our case, we model approximate spike counts using log-normal distribution with truncation at O. The approximate spike counts are then used as inputs to an order autoregressive (AR(1)) process model of calcium dynamics (g1,c). Since we found that the loss function was extremely sensitive to changes in the AR(1) model parameters, these parameters were fixed using values empirically derived from the DASISSdeconvvolution algorithm (Friedrich et al., 2017). While it is possible to transform calcium dynamics into fluorescence using nonlinear models of fluorescence indicator ki- netics (Speiser et al. 2017 Deneux 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. Model hyperparameters for different datasets are shown in Sa
2.1 LOSS FUNCTION AND TRAINING
One of the advantages of using VLAES is that because the stochastic nodes operate in parallel, the ELBO formulation is the same as that for standard VAEs despite the hierarchical latent space (Zhao ct al. 2017 As such, the CaLFADS cost function remains very similar to that of LFADS. It is described in Algorithm 2 One key addition to the LFADS loss function is that because the data is fluorescence data, which includes normal emissions noise, and the inferred shallower-level factors calcium concentrations, the likelihood function of the data given the estimated calcium concen- tration, ( is modelled as a Gaussian distribution Ie ~ ~N o²), where o² is determined by OASIS. Another key addition is that although spike counts W1,e are not modelled discretely, their probability distribution given the inferred spike rates, P(u (e)), is modelled as an approximate Poisson process such that ((( & U1,t log(A,) (Ae + u1,e). The approximation results from the use of an L1 norm on the spikes, instead of the ogg-gamma term that would normally be used in a Poisson process, which was necessary because most of the data sits within range that is poorly handled by the log-gamma function. The use of this approximate Poisson process pro- vides an additional regularization term in the cost function, which acts as an inductive bias to ensure that the inferred spike trains in the shallower-level dynamical system are constrained to follow the dynamics of the leeeer-level dynamical system's is spike rates Aç. Per the VLAE framework, the shallower-level and deeper-level latent variables are treated as independent, with hierarchy imposed in the deterministic paths, i.e., Q( I2, U2, (z) = (2) (=20))0000==))(((uuaa)
Parameters of CaLFADS were optimised with ADAM, with an initial learning rate of 0.01. Learning rates were decreased by a factor of 0.95 when plateaus in training error were detected. As in LFADS training, KI. and 1.2 terms in the cost function were 'warmed up", i.e., had a weight (eE(12) between
4

---

Under review as a conference paper at ICLR 2021
O and 1 aapllied that gradually increased (Bowman et al. 2016 Sonderby et al. 2016). Warm-up for the deeper parameters (1 2, blue modules in Figure was delayed until warm-up for shallower parameters (1 1, red modules in Figure I completed. We found in preliminary tests that this delayed warm-up sequence was necessary for good hierarchical inference compared to simultaneous warm-up. As with aggressive inference network training methods (He al.) 2019), this prevents issues with training of the deeper-level dynamical system parameters converging more slowly than the shallower-level dynamical system parameters in the initial stages of training.
3 RESULTS
3.1 LORENZ ATTRACTOR SYNTHETIC DATA
3.1.1 DATA DESCRIPTION
As an initial test of CaLFADS, we examined its ability to infer the hierarchical dynamical systems of synthetic calcium fluorescence data generated from a simple model of spiking neurons. We examined synthetic data generated from network with a Lorenz attractor embedded in its dynamics. This tests our ability recover ground-truth deeper-level dynamical system from the data, and has been used as a benchmark by many others (Zhao & Park 2017 2019 She& Wu 2019: Koppe et al. 2019).
Pandarinath et al. 2018 Clark et al.,
In this system, network dynamics were used to generate a set of spike rates in simulated neurons. We then used model of calcium dynamics and emissions noise to transform the spikes into synthetic fluorescence data. We tested CaLFADS on data generated from two models of synthetic calcium dynamics: a simple linear model equivalent to an AR(I) process, and a nonlinear model in which an AR(I) process is transformed by Hill function nonlinearity. We then added white noise to the resulting traces. This procedure is described in more detail in Appendix H
3.1.2 MODEL COMPARISON
We measured the performance of CaLFADS by computing R2 goodness-of-fi with the embedded Lorenz attractor state variables. We also measured the ability of CaLFADS to reconstruct other ground-truth variables in the synthetic data (spike counts, spike rates, and fluorescence traces) using R2, but these were not critical for assessing model performance.
CaLFADS performance was compared against two baselines. First, we use LFADS with a Gaus- sian likelihood observation model to account for fluorescence (Gaussian-LFAADS) Second, we con- sider the situation where spike counts are first estimated separately using the OASIS deconvolution algorithm (Friedrich et al. 2017), a robust, popular, and computationally inexpensive method of spike extraction widely used by systems neuroscientists (Pachitariu et al. 2018 [Evas et al.) 2019). LFADS is then used to infer the deeper-level dynamical system and reconstruct the spike rates from the deconvolved spike trains. We refer to this approach as OASIS+LFADS throughout.
Table compares R. goodness-of-ftt in reconstructing ground-truth dynamic latent variables in held-out validation data. Table [ shows CaLFADS is able to reconstruct the Lorenz attractor state for synthetic data with either linear or nonlinear synthetic calcium models (see Table D caption for significance test results). Indeed, the inclusion of nonlinearities in calcium dynamics does not affect CaLFADS ability to reconstruct the Lorenz attractor state as much it does with OA- SIS+LFADS. CaLFADS is also better able to reconstruct other network state variables compared to OASIS+LFADS as can also be seen in Table However, it is worth noting that high accuracy in reconstructing other network state variables does not appear to be necessary for accurate reconstruc- tion of the underlying Lorenz attractor. This indicates that CaLFADS is better suited to handling uncertainty due to spike-timing and is capable of separating sources of slow dynamics to recon- struct the embedded latent space. Example outputs of CaLFADS and OASIS+LFADS are shown in Appendix 1.1 and 1.2
Note that CaLFADS reconstructs the Lorenz dynamics almost as well as LEADS does when ap- plied to the spiking data. Of course, it is to be expected that LFADS applied to true spiking data performs better than CaLFADS, since there is an additional source of observation noise from the generation of fluorescence transients. But, the fact
get very close to the same level of perfor-
5

---

Under review as a conference paper at ICLR 2021
Table 1: Comparison of model performance (R² mean+sem) on synthetic Lorenz datasets generated with 15 different seeds. A hyphen indicates that the variable cannot be compared, as the model does not infer it. The top row is italicised as the performance of LFADS on this task is considered the upper limit, having no additional observation noise from fluorescence. Results of paired t-tests for testing significance of differences in performance between OASIS+LFADS and CaLFADS within seeds for a) linear calcium model: Lorenz state 114 2.85, p 0.013; Spike Counts 414 20.67, p < 0.001; Rates 114 = 4.81, P < 0.001; Fluorescence h14 = 16.2, p < 0.001. b) nonlinear calcium model: Lorenz state f14 3.29, p 0.005; Spike Counts t14 17.16, p 0.001;
Rates f14 6.35, < 0.001; Fluorescence 14 11.16, < 0.001.
Calcium
Model
Lorenz state
Spike Counts
Firing Rates
Fluorescence
NA
LFADS
0.975 + 0.001
0.956 + 0.003
Gaussian-LFADS
0.818 + 0.005
0.701 + 0,004
Linear
OASIS+LFADS
0,960 + 0.001
0.580 + 0.006
0.912 + 0.002
0.807 + 0,005
CaLFADS
0.965 = 0.001
0.817
0.007
0.937 + 0.002
0.909
0.003
Caussian-LFADS
0,808 = 0.005
0.650 = 0.006
Nonlinear
OASIS+LFADS
0.877 = 0.003
0.351 = 0.004
0.682 = 0.005
0.843 - 0.004
CaLFFAD
0.902 + 0.006
0.561 0.010
0.756 + 0.008
0.932 0.005
mance indicates that CaLFADS is effective at performing the same inference calcium data that LFADS performs on spiking data. We speculate that the uncertainty over ARI process parameters is overcome in CaLFADS by constraining the reconstructed spike counts with the latent dynamics to something which cannot
make the inference of the observed dynamics process "population-awart" be done when using LFADS applied to spike counts obtained by OASIS deconvolution. This may mitigate errors in spike count reconstruction that occur by deconvolution due to the absence of infor- mation about whole-network dynamics. Thus, CaLFADS can use population-level dynamics when conducting the inference of network variables, and this helps it to separate out the shallower-level dynamics more accurately.
3.2 ROTATIONAL DYNAMICS IN MONKEY MOTOR CORTEX
Next, we wanted to test CaLFADS on a real neural dataset to which LFADS has previously been ap- plied, and successfully used to uncover meaningful latent dynamics. Specifically, it has been shown previously that rotational dynamics underlie neuronal responses in monkey and human motor cortex during reaching behaviour (Churchland et al. 2012). LFADS has been successful in uncovering these known rotational dynamics for single-trial spikes recorded from primary motor cortex (M1) and dorsal premotor cortex (PMd) in macaques (Pandarinath et al. 20188) To test whether CaLFADS could do the same, we took the original spiking data and converted into semi-synthetic calcium traces, using the same model of calcium dynamics that was used for synthetic data generation in the previous sections 2A). We used data from the monkey clectrophysiology dataset previously deseribed, along with the reaching task, by Churchland et al. (2012), which was kindly provided to us by the original authors. Briefly, monkeys were trained to reach a target under 108 different reach conditions while multi-electrode recordings were made in MI and PMd. Reaches started from a specified location on a screen, and monkeys were rewarded for correctly reaching toward target while avoiding on-screen obstacles (Fig 2B).
First, to replicate the previous LFADS results, we applied jPCA (Churchland et al. 2012) to condition condition-averaged firing rates (trial averages for each of the 108 reach conditions), as well as single- trial firing rates inferred from spike data using LFADS. FFAA is a dimensionality reduction technique that finds orthogonal projections capturing rotational dynamics that explain variability in firing rates. The original LFADS paper showed that rotational dynamics explained a large amount of the vari- ance in firing rates. As in the original papers, we identified both condition-averager (Fig (Fig top) and single-trial (Fig - bottom) rotational dynamics from firing rates inferred by LFADS, which explained a large amount of the variance (R² = 0.81), thereby replicating the original findings. Then, we evaluated
CaLFADS model in uncovering the rotational dynamics of both condition- averaged and single trials from firing rates. As shown in Fig our model also successfully uncovers rotational dynamics from calcium traces, explaining a large amount of variance (R² 0.78) for both

---

Under review as a conference paper at ICLR 2021
A
C
LFADS (spikes)
D
CaLFADS (fluor)
E
LFADS (fluor)
Condition-an -oatiioo-averaeed averaged
projection anto JPC; (a.0.)
Ca
projection onto JPC,
Single-trial
(
B
projection anto JPC, (a.u.)
projection onto JPC, (a.u.)
projection onto JPC, (e.u.)
Figure 2: A) Schematic for the process of converting spikes to calcium traces, and B) for the reaching task. C) Rotational dynamics inferred from condition-averaged (Top) and single-trial spikes (Bot- tom) using LFADS. D) Rotational dynamics inferred from condition-as condition-averaged aveerged (Top) and single-trial calcium traces (Bottom) using CaLFADS. E) Rotational dynamics inferred from condition-averagee (Top) and single-trial calcium traces (Bottom) using OASIS+LFADS Traces are coloured based on their initial state values along jPC1 (from green red for increasingly larger values)
condition-averaun (Fig 2D top) and single trials (Fig 2D bottom). These results demonstrate that CaLFADS capable of identifying latent dynamics real neural data, similar to LFADS, even when the spiking data is transformed by calcium dynamics and emissions noise. We also found that the LASIS*LFADS approach was not as successful as CaLFADS in uncovering rotational dynam- ics, explaining only half the variance (R2 0.53, Fig 2B.. This result corroborates the importance of the hierarchical modelling of calcium and neuronal dynamics in CaLFADS.
3.3 UNEXPECTED STIMULUS DETECTION IN MOUSE PRIMARY VISUAL. CORTEX
Finally, we wanted to test CaLFADS on an entirely new calcium imaging dataset, to determine whether CaLFADS can infer dynamic computational factors that carry information about relevant features of the outside world. We chose to analyse data from mouse primary visual cortex (VisP), region widely studied in systems neuroscience research where calcium imaging is standard tool. There is evidence that the visual cortex of mammals performs a predictive function, anticipating expected stimuli (Dayan al. 1995 Rao & Ballard [1999] Fiser et al. [016 Leinweber et al. 2017). As a result, unexpected stimuli can induce perturbations in network dynamics (Fiser et al. 2016
Leinweber al. 2017). Thus, wanted to determine whether CaLFFADS could infer dynamic computational factors that carry information distinguishing unexpected from expected visual stimuli.
To this end, trained CaLFADS calcium imaging data from the mouse visual cortex collected by the Allen Institute for Brain Science (for a detailed description, see Appendix J.I While awake behaving mice were presented with visual stimuli a screen (Fig AA), calcium fluorescence re- sponses in cortical layer 2/3 of VisP were recorded using 2-photon microscopy (Fig B)). The mice were familiarised with sequences of stimulus frames that followed simple probabilistic rules. Briefly, each sequence consisted of four randomised Gabor pattern frames followed by grey screen (A, B, C. D, grey) with orientations drawn for each trial from the same distribution. After familiarization with this sequence, the expected D frame was replaced in 10% of trials by a Gabor pattern (E) with orthogonal orientations (Fig 3C) CaLFADS inferred smooth latent factors that did not show signs of over-fitting noise (Appendix J.2. Fig S7 To determine if salient information about these unexpected stimulus features was present in the inferred factors from CaLFADS, we trained a non- linear decoder to classify whether the inferred factors came from trial with expected D frame or an unexpected E frame (for details, see Appendix J.3). We also trained the non-lineat decoder on factors inferred by OASIS+LFADS. and fluorescence traces reduced in dimensionality with prin-

---

Under review as a conference paper at ICLR 2021
A
B
C
D
0.9
Average recall
B
0.7
D or
Fhao
Gaussion
LFADS
OASIS LFADS CaLFADS
gr
Figure 3: A) Schematic of 2-photon calcium imaging recording setup. Mice are head-fixed on running wheel under the imaging objective, and visual stimuli are projected onto screen to the side, B) 2-photon calcium imaging recording plane, C) Schematic example of visual stimuli for expected AA-B-C-D-grey and unexpected (A-B-C-E-grey) trials. D) Average recall (mean = sem) of expected vs. unexpected trials across non-linear decoders trained on principal components of fluorescence traces (0.815 + 0.011, green), Gaussian-LFADS factors (0.500 + 0.000, blue), LASIS+LFADS factors (0.810 + 0.005, purple), and CaLFADS factors (0.871 + 0.004, maroon). All pairs but one (fluor. vs. OASIS+LFADS factors, marked with n.s.) were significantly different at < 0.05 in 2-tailed independent t-tests with Bonferroni correction for multiple comparisons.
cipal components analysis. To compare the models fairly, despite the imbalance between expected and unexpected trial frequencies, we measured average recall. We find that using the CaLFADS factors leads to the highest performance on stimulus -trial identity decoding, with the average re- call being significantly higher than when OASIS+ FFASS factors are used (Fig 3D). This indicates that CaL.AADS infers latent dynamics from real visual cortex calcium imaging data that reflect the predictability of visual stimuli, corroborating the ability of CaLFADS to infer meaningful dynamic factors in real data.
4 DISCUSSION
In this paper we presented CaLFADS, hierarchical recurrent variational autoencoder model ca- pable of reconstructing latent computational dynamics. We confirmed CaLFADS' ability to recon- struct known underlying dynamics using synthetic datasets where ground-truth was known. We also showed that CaLFADS is able to infer sensorimotor dynamics from real neural recordings. This indicates that CaLFADS a promising method for analyzing calcium imaging data in neuroscience.
There are two key advantages of CaLFADS over the use of deconvolution of calcium traces fol- lowed by application of LFADS. The first is that the latent dynamics, and the latent spike counts. The second is that CaLFADS performs better than OASIS+LFADS when nonlinearities in caleium dynamics are present, is the case in real experi- mental data where there are many sources of calcium influx related to spikes.
obtain measures of the uncertainty in both
We designed CaL sequential variational ladder autoencoders, in which different layers of recurrent neural networks can be used infer and generate different layers of hierarchical dynamical system. In this sense, this class of models is modular and composable, meaning it could be readily adapted to other do- mains. For example, it should be possible to replace LFADS as the deeper-level dynamical system model with any other differentiable model. Likewise the same is true for our AR1 based model of calcium dynamics: we could replace the calcium dynamics model with other models as required by the experimental set-up. Furthermore, there is no need to stop at two layers in the hierarchy: the brain is comprised of many interconnected recurrent neural networks sending long range signals to one another, and it should be possible to add in additional modules to capture additional hierarchies and dynamics within the brain.
fit into much broader class of artificial neural network models, namely
In summary, CaLFADS is new, open source tool for inferring hierarchical dynamics from calcium imaging data that also has great potential for being modified and applied to other data modalities. This could be of real benefit to the thousands of neuroscience laboratories around the world con- ducting calcium imaging experiments.
8

---

Under review as a conference paper at ICLR 2021
REFERENCES
Laurence Aitchison, Lloyd Russell, Adam M Packer, Jinyao Yan, Philippe Castonguay, Michael Hausser, and Srinivas Turaga. Model-based bayesian inference of neural activity and connec- tivity from all-optical interrogation of neural circuit. In L. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30, pp. 1486-3495 Curran Associates, Inc., 2017.
A. Anonymous. Identifying the role of top-down feedback connections in cortical learning and inference (AIBS openscope grant proposal), 2018.
Philipp Berens, Jeremy Freeman, Thomas Deneux. Nikolay Chenkov, Thomas McColgan, Artur Speiser, Jakob H. Macke, Srinivas C. Turaga, Patrick Mineault, Peter Rupprecht, Stephan Ger- hard, Rainer W. Friedrich, Johannes Friedrich, Liam Paninski, Marius Pachitariu, Kenneth D. Harris, Ben Bolte, Timothy A. Machado, Dario Ringach, Jasmine Stone, Luke E. Rogerson, Nico- las J. Sofroniew, Jacob Reimer, Emmanouil Froudarakis, Thomas Euler, Miroslav Román Rosón, Lucas Theis, Andreas S. Tolias, and Matthias Bethge. Community- benchmarking improves spike rate inference from wo-phooon calcium imaging data. PLOS Computational Biology, 14 (5):1-13, 05 2018. doi: 10. 1371 1ft.///oorall.pcbb........ 106.m.11....
1006157.
Samuel Bowman, Luke Vilnis, Oriol Vinyals, Andrew Dai, Rafal Jozefowicz, and Samy Bengio. Generating sentences from a continuous space. In Proceedings of The 20th SIGNLL Conference on Computational Natural Language Learning, pp. 10-21, Berlin, Germany, August 2016. Asso- ciation for Computational Linguistics. doi: 10.18653/v1// 1/tp:////1111-11 1002.
Rishidev Chaudhuri, Berk Gerçek, Biraj Pandey, Adrien Peyrache, and lla Fiete. The intrinsic attrac- tor manifold and population dynamics of a canonical cognitive circuit across waking and sleep. Nature Neuroscience, 1219):1512-1520 September 2019. doi: 10.1033//41593-009-
Number: 9 Publisher: Nature Publishing Group.
T. W. Chen, T. J. Wardill, Y. Sun, S. R. Pulver, S. L. Renninger, A. Baohan, E. R. Schreiter, R. A. Kerr, M. B. Orger, Jayaraman, L. L. Looger, K. Svoboda, and D. S. Kim. Ultrasensitive fluorescent proteins for imaging neuronal activity. Nature, 49957458) 2013.
Mark M. Churchland, John P. Cunningham, Matthew T. Kaufman, Justin Foster, Paul Nuyu- jukian, Stephen I. Ryu, and Krishna V. Shenoy. Neural population dynamics during reaching. Nature, 587(7405)551-56, July 2012. doi: 1038////aa Number: 7405 Publisher: Na- ture Publishing Group.
David G. Clark, Jesse Livezey, and Kristofer E. Bouchard. Unsupervised discovery of tempo- ral structure in noisy data with dynamical components analysis. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' Alché-Buc, Emily B. Fox, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 32: Annual Conference on Neural Informa- tion Processing Systems 2019, NeurIPS 2019, 8- 14 December 2019, Vancouver, BC, Canada, pp. 14267-14278, 2019.
Peter Dayan, Geoffrey E. Hinton, Radford M. Neal, and Richard S. Zemel. The Helmholtz machine. Neural Comput., ((5):889 904, September 1995. doi: 10.1162/neco.19955
Saskia EJ de Vries, Jerome A Lecoq, Michael A Buice, Peter A Groblewski, Gabriel K Ocker, Michael Oliver, David Feng, Nicholas Cain, Peter Ledochowitsch, Daniel Millman, et al. A large-scale standardized physiological survey reveals functional organization of the mouse visual cortex. Nature Neuroscience, 23(1):138-1) 2020.
Thomas Deneux, Attila Kaszas, Gergely Szalay, Gergely Katona, Tamás Lakner, Amiram Grinvald, Balázs Rózsa, and Ivo Vanzetta. Accurate spike estimation from noisy calcium signals for ultrafast three-dimensional imaging large neuronal populations in vivo. Nature Communications, 7(1): 12190, July 2016. doi: 10.1038/ncomms/22
Mathew H. Evans, Rasmus S. Petersen, and Mark D. Humphries. On the use of calcium deconvolu- tion algorithms in practical contexts. bioRxiv, 2019. doi: 10. 10.110111100 1011110//71133 137.

---

Under review as a conference paper at ICLR 2021
A. Fiser, D. Mahringer, H. K. Oyibo, A. V. Petersen, M. Leinweber, and G. B. Keller. Experience- dependent spatial expectations in mouse visual cortex. Nature Neuroscience, 19(12):1658-1664 2016.
Johannes Friedrich, Pengcheng Zhou, and Liam Paninski. Fast online deconvolution of calcium imaging data. PLOS Computational Biology, 13(3)=e1004423, March 2017.
Karol Gregor, Ivo Danihelka, Alex Graves, Danilo Rezende, and Daan Wierstra. Draw: A recurrent neural network for image generation. volume 37 of Proceedings of Machine Learning Research, pp. 1462-1471, Lille, France, 07-09 Jul 2015. PMLR.
Junxian He, Daniel Spokoyny, Graham Neubig, and Taylor Berg- Kirkpatrick. Lagging inference net- works and posterior collapse in variational autoencoders. In International Conference on Learning Representations, 2019.
Daniel Hernandez, Antonio Khalil Moretti, Ziqiang Wei, Shreya Saxena, John Cunningham, and Liam Paninski. Nonlinear evolution via spatially deeeadenn linear dynamics for electrophysiol- ogy and calcium data, 2020.
J. Homann, S. A. Koay, A. M. Glidden, D. W. Tank, and M. J. Berry. Predictive coding of novel versus familiar stimuli the primary visual cortex. bioRxiv, pp. 17, 2017. doi: 10. 10.1101/197608.
Elke Kirschbaum, Manuel Haußmann, Steffen Wolf, Hannah Sonntag, Justus Schneider, Sheha- beldin Elzoheiry, Oliver Kann, Daniel Durstewitz, and Fred A Hamprecht. LemoNADe: Learned motif and neuronal assembly detection in calcium imaging videos. In International Conference on Learning Representations, 2019.
Georgia Koppe, Hazem Toutounji, Peter Kirsch, Stefanie Lis, and Daniel Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fmri. PLOS Computational Biology, 15(8):1-35, 08 2019, doi: 10. 1joun/jorllllccbb 1007263.
Marcus Leinweber, Daniel R. Ward, Jan M. Sobezak, Alexander Attinger, and Georg B. Keller. A Sensorimotor Circuit in Mouse Cortex for Visual Flow Predictions. Neuron, 95(6):1400-1432.e. September 2017. doi: 10.1016/j 200....neroon.21..8888 20..008..066.. Publisher: Elsevier
Michael Z. Lin and Mark Schnitzer. Genetically encoded indicators of neuronal activity. Nature neuroscience, 1999),1142, 2016.
JM Ollinger, M Corbetta, and GL Shulman. Separating processes within a trial in event-related functional mri. Neuroimage, 13(1):218 2298, 2001.
Marius Pachitariu, Carsen Stringer, and Kenneth D. Harris. Robustness of Spike Deconvolution for Neuronal Calcium Imaging. The Journal of Neuroscience, 38(37):7966-7885 September 2018.
Chethan Pandarinath, Daniel O'Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman, Stephen I. 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.
Rajesh P. N. Rao and Dana H. Ballard. Predictive coding in the visual cortex: a functional interpre- tation of some extra-classical receptive-field effects. Nature Neuroscience, 2(1):79-87, January 1999. doi: 10.1038/4580. Number: Publisher: Nature Publishing Group.
Thorsten Schmidt. Catastrophe insurance modeled by shot-noise processes. Risks, 2(1):3-24, 2014.
Qi She and Anqi Wu. Neural dynamics discovery via gaussian process recurrent neural networks. In Amir Globerson and Ricardo Silva (eds.), Proceedings of the Thirty-Fifth Conference on Un- certainty in Artificial Intelligence, UAI 2019, Tel Aviy, Israel, July 22-25 2019, pp. 159. AUAI Press, 2019.
Casper Kaae Sønderby, Tapani Raiko, Lars Maaløe, Søren Kaae Sønderby, and Ole Winther. Ladder Variational Autoencoders. In D. D. Lee, M. Sugiyama, U. V. Luxburg. I. Guyon, and R. Gamett (eds.), Advances in Neural Information Processing Systems 29, pp. 3738-3746 Curran Asso- ciates, Inc., 2016.
10

---

Under review as a conference paper at ICLR 2021
Artur Speiser, Jinyao Yan, Evan W Archer, Lars Buesing, Srinivas C Turaga, and Jakob H. Macke. Fast amortized inference of neural activity from calcium imaging data with variational autoen- coders. In I. Guyon, U. V. Luxburg, S. Bengio. H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems 30, pp. 4024-4034. Cur- Associates, Inc., 2017.
Xue-Xin Wei, Ding Zhou, Andres Grosmark, Zaki Ajabi, Fraser Sparks, Pengcheng Zhou, Mark Brandon, Attila Losonczy, and Liam Paninski. A zero-inflated gamma model for post- deconvolved calcium imaging traces. bioRxiv, 2019a. doi: 10. 101//6666765.
Ziqiang Wei, Hidehiko Inagaki, Nuo Li, Karel Svoboda, and Shaul Druckmann. An orderly single- trial organization of population dynamics in premotor cortex predicts behavioral variability. Na- ture Communications, 10(1):216, Jamuary 2019b. doi: 10.1038//41467-01-088/8/ 01-10-111161866666.1 Number: Publisher: Nature Publishing Group.
Ziqiang Wei, Bei-Jung Lin, Ten Chen, Kayvon Daie, Karel Svoboda, and Shaul Druckmann. A comparison of neuronal population dynamics measured with calcium imaging and electrophys- iology. PLOS Computational Biology, 16(9):1-29, 09 2020. doi: 10.137 L/jurnall.pcbi. 1008198.
Weijian Yang and Rafael Yuste. In vivo imaging of neural activity. Nature methods, 14(4):349, 2017.
Shengjia Zhao, Jiaming Song, and Stefano Ermon. Learning Hierarchical Features from Deep Gen- erative Models. In International Conference on Machine Learning, pp. 4091-4099, July 2017.
Yuan Zhao and Il Memming Park. Variational Latent Gaussian Process for Recovering Single- Trial Dynamics from Population Spike Trains. Neural Computation, 1965):1293-1316, March 2017.

---

Under review as a conference paper at ICLR 2021
A
VARIABLE GLOSSARY
Table S1: Glossary variables, with descriptions and dimensionality across datasets
Dimensions
Variable
Description
Lorenz
RNN
M1/PMd
V1
2t
Fluorescence signal
30
50
202
96
Calcium dynamics embedding encoder state
128
64
200
128
h2,
Computational dynamics embedding encoder state
64
128
100
128
Cia
Calcium dynamics embedding controller state
128
64
200
128
C2
Computational dynamics embedding controller state Calcium dynamics embedding state
128
64
g1.t
100
64
200
128
22,1
Computational dynamics embedding state
64
200
100
200
U1.1
Approximate spike counts
30
50
202
96
U2,t
Computational dynamics perturbations
g1.t
Calcium dynamics
30
50
202
96
ft
Dynamic factors
20
40
32
dt
Reconstructed fluorescence signal
30
50
202
96
12

---

Under review as a conference paper at ICLR 2021
B RELATED WORK AND CONTRIBUTION
This work brings together two important strands of computational neuroscience research that have been of interest to the machine learning community, namely probabílistic latent factor models of neural dynamics and inference of spikes in calcium fluorescence traces. We solve these two prob- lems as hierarchical inference problem using variational ladder autoencoders. We give brief introduction to variational ladder autoencoders in Appendix C
LFADS power latent factor model of neural dynamics proposed by Pandarinath et al. (2018).
It is an extension of the sequential variational autoencoder DRAW (Gregor et al.) (015) applied to spiking data that adds an additional RNN to the encoder architecture to infer unknown initial conditions to a generator RNN. This enables interpretation of the generative model as nonlinear dynamical system that generates observed spikes. Low dimensional projections of the generator hidden state (latent factors) were found to be capable of decoding single-trial reaching movements with very high accuracy. This was a hugely influential and important contribution to the field of neuroscience it demonstrated means to find latent dynamical systems that solved computational tasks required for executing behaviour: However, LFADS could not be applied to neural data where the spiking was itself unknown, in two-photon calcium imaging.
SeepSpike by Speiser et al. (017) is a variational autoencoder adapted from MLSpike, a proba- bilistic model of calcium fluorescence recordings (Deneux et al. 2016)) It was a contribution the SpikeFinder competition (Berens et al. 2018) in which participants were given the task of inferring an unknown spike train from caleium fluorescence recordings. Unlike in most two-photon record- ings, these data also included ground-trutt whole-cell electrophysiology. i.c., ground-truth spike trains. While other techniques had higher performance, it was the first attempt at solving this prob- lem using variational inference. One issue with this model is that it used a Bernoulli distribution to model the latent spike train. Since this is not reparameterisable, Speiser et al. (2017) needed to use a high variance Monte Carlo objective to train the model. The use of this objective means it would not scale well to 2-photon measurements from many cells. Furthermore, such a model does not capture how neurons spike multiple times within the time-bin bin of a typical two-photon imaging experiment.
We note some similarities between our generative model and Aitchison et al. 2017). A key differ- ence here is that we try to find nonlinear latent dynamics, whereas their model comprises a linear dynamical system in the latent space. Furthermore, theirs is fully Bayesian approach, whereas CaLFADS has many more deterministic components.
We also note the capabilities of VIND (Hernandez et al. 2020 in inferring non-linear latent dynam- ics in wide-field calcium imaging data, which has insufficient spatial resolution to observe individual cells. Hernandez et al. (020) propose variational inference technique in which latent dynamics are modelled by a locally linear dynamical system and find informative latent trajectories in the recorded smoothed population activity. We note however that this approach is not used to infer latent dynam- ics in two-photon recordings, where it would be necessary to disentangle calcium dynamics from population dynamics.
Finally, the variational autoencoder LeMoNADe proposed by Kirschbaum et al. (2019) similarly proposes an end-to-end system of inferring latent activities from calcium imaging. In contrast to CaLFADS, they are attempting to infer the presence of neuronal ensembles, sets of cells grouped together by their activities. While this is a related goal to inferring latent dynamical system, there is requirement for the latent factors to be as interpretable as in LeMoNADe. Furthermore, the neuronal ensembles identified by LcMoNADe were not demonstrated to help solve relevant behavioural tasks. However, LeMoNADe were the first to demonstrate how neuronal activities could be modelled by a continuous relaxation to the Bernoulii distribution using the Gumbel-Softmax trick.
CaLFADS ties together this related work. We propose a hierarchical variational autoencoder capable of inferring latent dynamic factors in calcium imaging data by separating shallower-level calcium dynamics and deeper-level computational dynamics. We propose a novel continuous relaxation to spike counts that allows an arbitrary number of spikes to be counted within the typically low sampling rate of calcium fluorescence recordings (10-30 Hz), We also propose a novel KL annealing technique for allowing hierarchical latent dynamics to be inferred.
13

---

Under review as a conference paper at ICLR 2021
C INTRODUCTION TO VARIATIONAL LADDER AUTOENCODERS
Variational Ladder Autoencoders were proposed as a way to ensure that higher- and lower-level features could be stored in different parts of the latent state. They propose a ladder architecture in which more expressive networks sit deeper in the hierarchy and no hierarchy is assumed among latent variables directly. For layers in the hierarchy
P(zi,... zT.) II P(z)
(1)
((21 fi(hi) h 1-
(2)
where ho = and and fi are neural networks with increasing degrees of expressivity as / increases. This encourages more abstract features to be represented certain portions of the latent state.
This is in contrast traditional hierarchical autoencoders that assume Markovian structure in the latent variables
2(x,z)= P(x|z) II P P(z
The key insight from (Zhao al. 2017) is that under idcalised assumptions for sampling and optimization the ELBO, the ELBO for hierarchical VAEs is optimised by minimizing only (2 P(z, [x)). This means that traditional hierarchical VAEs will load all representa- tions onto the lowest level of the hierarchy.
In contrast, by splitting the latent code by abstractness and parameterising the factorised posterior distribution with increasingly expressive functions, the VLAE can be trained with the simple ELBO formulation used in non oirrhieaalllcal VAEs.
14

---

Under review as a conference paper at ICLR 2021
D CALFADS PROBABILISTIC MODEL
Here we outline the assumed probabilistic model for CaLFADS. In doing so, we will briefly review the model of LFADS (Pandarinath al., 2018) and the extensions that make it possible to infer dynamic latent variables from calcium fluorescence traces.
D.I LFADS GENERATIVE MODEL
LFADS assumes that neural activity x is generated from low dimensional non-linear dynamical system with unknown initial conditions subjected unknown control inputs. The state of the low- dimensional non non-linear dynamical system is modelled by the hidden state of Gated Recurrent Unit receiving control imputs ue. The prior for the initial conditions go is assumed to be normal with mean Mgo and variance of External inputs are assumed to have an autoregressive process prior with process mean Mux process variance of and time constant Tz,
P(g+) NN ogv))
(5)
P(u1) = -N o.H)
(6)
P(uelue-1) Nu - (u - -ul), 02 -
(7)
gt = GRU Gg (1 ut) The hidden states of the GRU are then projected to now-dimensional set of factors ft
(8)
fr W gt. + b/
In Pandarinath et al. (2018), the neural activity being modelled was single-trial spikes extracted from electrophysiolgical data macaque motor cortex. LFADS assumes simple observation model for spiking data where the spike count observed in given time-bin di follows an independent Poisson process parameterised by the instantancous firing rate A, To obtain the instantaneous firing rate, low-dimensional factors ft are linearly projected to the dimensionality of the data and non non-linearity enforeing positivity applied, in this case the exponential function
de = expex + -b²) P(x+|A+)= = Poisson(Ae)
(10)
(11)
D.2 CALCIUM FLUORESCENCE GENERATIVE MODEL
Before describing how the low dimensional dynamic factors are used to generate calcium fluores- cence data, first describe probabilistic model for calcium fluorescence transients. We assume that calcium fluorescence transients are generated by an autoregressive process e occluded by white emissions noise and perturbed by unknown spike train St-
Probabilistic models of unknown spike counts in calcium fluorescence data typically assume a Bemoulli distribution, some continuous approximation thereof (Speiser al. 2017). This is somewhat problematic given that fluorescence is sampled at rate of 10-30 Hz, leaving time-bins large enough observe more than one spike per time-bin. Recently, Wei et al. (00 9a) demonstrated that spike counts inferred by the deconvolution algorithm OASIS (Friedrich et al. 2017) could be closely modelled by a zero-inflated gamma model. Similar this, we modelled spike counts using a similar continuous approximation
P(st) = Ru (LogN(Looooollllll o2,) -1
(12)
This approximation is easily reparameterisable since sample from lognormal distribution by transforming samples from a normal distribution with an exponential function. The ReLU function and subtraction of I from samples ensures that many counts inflated gamma model which is a mixture distribution, it would require modelling a hierarchical discrete-contimuous latent variable which is not simply reparameterisable for backpropagation.
In contrast, since the zero-
15

---

Under review as a conference paper at ICLR 2021
Spike counts then used to perturb autoregressive process modelling the slow decay of calcium fluorescence after a spike. A simple autoregressive model for calcium fluorescence is the AR(1) process. With discretised time, this can be modelled with Euler updates
Ct ct- 3) ast
(13)
where 3 = At with At the size of time bin and Te the decay of calcium fluorescence, and with AL Te. Or is constant gain term controlling the effect of spikes on fluorescence.
A nonlinear function modelling dye kinetics d(-)) then be applied. For these experiments have not done so, but it would be a relatively simple extension. As such, the observation model for fluorescence F becomes
= ((((c))
(14)
where in our case, d(ct) = Ct.
D.3 CALFADS GENERATIVE MODEL
To tie LFADS and our observation model for calcium fluorescence together, we will adjust the notation slightly to make the hierarchical positioning of latent variables clearer. The variables of the deeper-level model of LFADS are denoted by the subscript shallower-level model of observed calcium fluorescence are denoted by the subscript / 1. For consistency and notation simplicity, denote variables that are the state of dynamical system by Si. and inputs to the dynamical system as ur,t-
2, whereas the variables of the
We observe neural activity in the form of single-trial calcium fluorescence traces. We assume that fluorescence follows a simple autoregressive process gi=1 perturbed by an unknown set of spike counts u==1 and subject to additive white emissions noise. We also assume that spike counts u1=1 are influenced by the state of a nonlinear dynamical system with unknown initial conditions g1-2,t=0 and unknown inputs 11-2-
As in LEADS, the nonlinear dynamical system is modelled by a GRU. The hidden state of the GRU 2t= is projected low-dimensional set of factors ft
P(g..0) = N (P o (200)
(15)
P(u2,1) N (Muuu ou2) -1 P(u2,(u22e-1) N Nmn, +e Two (U2,t-1
(16)
-1
paz), ob (1 Tuz)
(17)
g2 (g. 12,,)) f = W g2,t + bd
(18)
(19)
For spike counts to be informed by the dynamic factors f. sample spike counts according
Zt P(z.) (2)
(20)
U1.t ReLU (exp( W [ze fú] + b*) -
(21)
using the continuous approximation to spike counts described in section D.2
Although this factorization separates the generation of approximate spikes counts from the dynamic latent factors, for the formulation of the loss function we treat U1,t if it is generated from homogencous Poisson process with intensity At exp(W>Ft + b^),
P( Axe) R 1 Poisson(Xe)
(22)
We then use these approximated spikes counts to perturb the state gt of an AR(1) process
g1t t B1 B) + at11,t
(23)
The likelihood of the observed fluorescence in given time bin x is then modelled normal random variable with mean dependent on the state of the AR(1) process.
((xlggı,)) (Nd(g..,))
(24)
where again d(-) is optional function modelling dye kinetics, In our case d(g1,t) gi,
16

---

Under review as a conference paper at ICLR 2021
D.4 CALFADS INFERENCE NETWORK
The inference network of CaLFADS parameterises the approximate posterior distributions Qgg.o(a)) (2)2 and Q=1(²). As per the ladder architecture (Zhao et al., 2017) these poste-
rior distributions take the form
Q(z[x) qa (x)(
(25)
Q(u2,( = (x(p((x))
(26)
where p()) g(-), r(-) are neural networks parameterised by Φ, 0, v respectively. The combina- tion of latent space factorization and parameter sharing ensures that dynamics can be disentangled effectively.
As in LFADS, the approximate posterior distribution over the initial conditions to the nonlinear dynamical system (02,0) is parameterised separately using bidirectional Gated Recurrent Unit
Qg2,,0 = (Ng22.0 g2.0)
(27)
Mgz.o
W
92.0 + b "22.0
(28)
02.0
exp(
2z.0 000 b b220))
(29)
where h92.0 is the final hidden state of a bidirectional Gated Recurrent unit
92.0
[h h
fund
hg hg o2.b.b
(30)
h BIG 90.0 xi)
(31)
GRU GRu
.h411
x()
(32)
where h h.2.0./wd fwdd is the state of a GRU running forward sequentially over the input, and h (2.,, bund the state of a GRU running backward sequentially over the data.
The functions f-). ge(-), )), (-) are modelled by pairs of recurrent neural networks. As LFADS, we use bidirectional GRUs followed by unidirectional controller GRUs receiving representations from the generative network. The bidirectional GRU is defined before
h BiGRU (h", Xx)
(33)
h (U
(34)
From this point the inference model splits into the levels of the hierarchy. The hidden state h is passed through the a controller GRU which takes feedback from the generated samples
C1, GRU (hf u==1)
(35)
C22 GRU [h f--1))
(36)
The states of the controller GRU are then linearly transformed onto the parameters of normal distribution
Q(2((x) = Nu o2,) Q(u,,(x)= = NNpr,,
(37)
(38)
Nz. -W + btz
(39)
oz, = expp WD w C1.1 bT))
(40)
Pus. = W + Huz
(41)
Juz. = expl W W + bba2)
(42)
17

---

Under review as a conference paper at ICLR 2021
D.5 CALFADS LOSS FUNCTION
To construct the cost function of CaLFADS, we start from the evidence lower bound (ELBO)
P(x) 2 E2, -e..a....0-0((a............ [lo P(æ, un 11/,, A)]
(43)
-Qkl(Q(2, u2, (2.02))))] 42: 92.0))
=
where x denotes the observed data and & denotes the reconstructed data. With our model, this factorises as
CELBO Ez 00 [log P(x[*)]
(44)
1 www.ga..-Q(2..u.ee. log P(ui
DKL (K 82,0 [x] P (2..)))))))
Since "1, continuous random variable that we are treating as if it came from an independent, homogeneous Poisson process with parameter At, we approximate the log-likelihood with
log P(u E log P(u
(45)
~ Σ un log A - -- = -aapike
(46)
This approximation very similar to the true log-likelihood of Poisson distribution. The key difference is that the log-factorial term over the observations is replaced with the identity over the observations, i.e., log u1.e! & U1 In deep learning applications, the log-factorial is often dropped entirely from the Poisson log-likelihood (see e.g., the Pytorch implementation and default setting It is sometimes replaced with shifted log-gamma function log k! s log mation via Stirling's formula log k! 2 log 2 Since Stirling's formula blows up near zero, and since log-gan tog-gamma function has a local minimum between zero and one we decided against using these.
[[ + 1) or an approxi-
However, we observe that this term helps to enforce sparsity in the observations. Since this is desirable for inferred spike counts, we replace this term with what essentially the L1 norm u1,6 since by construction i is forced to be positive.
The remaining terms of the ELBO were not approximated in any special way. As such we breakdown LELBO follows,
Crecon -E, - flo P(*[a]]
(47)
KL1 DKL (Q()[=)]((()))
(48)
Kl DRI. DKI,(QQua, ((l
(49)
The final cost function of CaLFADS becomes
LCaLPADS Crecon + Lepikes + WKL, CKL, + WKL,
(50)
where WKL, and WKL, are KL warmup terms that are typically used in training variational au- toencoders. These terms slowly transition the model from a traditional autoencoder to a variational autoencoder, which helps to prevent pathological behaviour where the KL terms are minimised too early in training. This typically results in the model then generating trivial solutions. In our case, use a staggered KL warm-up which WKL, linearly increases from to over epochs 0-100, then WKL, linearly increases from 0 to over epochs 100-200. This allows the model to converge to solution that overfits spikes to the fluorescence data first, then gradually regularises this proposed set of spikes by ensuring they well explained by the low-dimensional latent factors.
htttps:/pyttrcch.. http://///.ooooo/oo////////o//tr
stabl e/ httee//eeetteettorr..ieee nn. Pai ssonNI....oss html
18

---

Under review as a conference paper at ICLR 2021
E
CALFADS DIRECTED ACYCLIC GRAPH
A
A2
Ap-1
AT
fo
f1
fr-1
fr
+
A
1
->
g²
g²;
gº
g2
00²0
A
u²
1
222
22
C2
T-1
X2
ho2
fwid h2
Fwid hu
fwid
fwid
fwd
hu2
h2
huz
h
XT-1
hd2
bwd
XT
ho2
g
92
g
g
ha2 T+1
u²
u²
1
Z1
Z2
ZT-1
Zr
4
2
47
2
0'T1
0²
1
C²
fwid
fwd
Foed
fwid
fwid
fwd
hul
hul
hul
hu
hul
h
X1
XT.1
XT
Figure S1: Directed acyclic graph for hierarchical model. Solid arrows denote deterministic map- pings, open arrows denote sampling steps.
19

---

Under review as a conference paper at ICLR 2021
F CALFADS FORWARD METHOD PSEUDOCODE
Algorithm 1: Forward method for CaLFADS
Forward:
Input Fluorescence (dF/F) traces with dimensions T&N Output: Denoised dF/F é, Inferred spike train u1,t Inferred dynamic factors fz.
[N/vdd 1 h hawd 01100 2t=1:1
for t < to do
for 1 < 1 to 2 do
h 11t
GRU (hg fwed [hffdd Eluud)
h bued 11, GRU 2 1: 1,T -.T--1) 1)
had
fwdd h
end
end
Dense ( (12,7444 12200
hurd)
pgnat=o
10 02 exp(0.5Densed fwd hbwd 12)
11 92 N (o o 20)
fi=0 (gn
13 for t &- to do
14
GRU(CI, C₂ G
[12, 122
fwd bud 11)
15
[h
fiod pubbrd
22, f(-1l)
Pug.1 = (ense(c2,))
16
17
#2. = exp( Deep(Denee(C...
2
U2 ((uu )
19
22.t GRU (9,1-12 U2211)
29
ft = Dense(g2.1
A exp(0.5Dense(fi)
21
MaL Dense(cr.t)
22
a exp(0.5Dense(cu.))
24
11, :
41, ReLU (eeneeeeeDDes((essee fol) 1))
25
I11 AR1 (1
27
&t Dineiinlics(g..
28
end
20

---

Under review as a conference paper at ICLR 2021
G CALFADS LOSS FUNCTION PSEUDOCODE
Algorithm 2: Loss function for CaLFADS
CaLFADS Loss Function;
Input Fluorescence (dF/F) traces æ, reconstructed dF/F æ, Intensity function X, Inferred
spikes u1, Posterior distribution (' 93. aagas May 2zy
Output: Loss
Lrecon (z-((())) llog P(a
Lapike = -2 (x) [log(ui|\]]
LKL, wi DKI. (Q(2)=))((())
LKI u2 DKI. (Q u2 (x) P ((gz. W2))
C = Lrecon + Lypike + LKL, + LKL.
return C
21

---

Under review as a conference paper at ICLR 2021
H SYNTHETIC DATA GENERATION
As discussed the main text, synthetic calcium fluorescence data was generated by embedding a Lorenz Attractor in the dynamics of spiking neural network. Two models of calcium dynamics were then used to transform spikes to calcium transients, with additive emissions noise added to the resulting traces.
The Lorenz Attractor a nonlinear dynamical system with 3 states 21, 22, Dy commonly used to study chaotic dynamies. The dynamical system is defined by,
dz
0 - -x1)
(51)
dt
dzz
I (p -23) - 222
(52)
di
dr3
D1 - B333-
(53)
dt
(54)
This system was parameterised in typical chaotic regime with a 10, 8/3, p- = These states were then normalised to have mean of zero and a range of [-1,
The state of this system x was then randomly projected onto the firing rate a population of = 30
neurons,
A exp Wx; + b)
(55)
with W NOO, 1 I), and VN using a Poisson distribution
Spike counts in time-bin for neuron 8 were then sampled
$i, ~ Poisson(Nr,AAt)
(56)
where At is the width of a time-bin. For our simulations, we used At = 0.1s.
Calcium fluorescence traces were then modelled one of two ways. The first method modelled calcium concentration in neuron as an exponentially decaying variable with time constant T
0.3 and perturbations by spikes
de;
(57)
dt
Fluorescence was then modelled by adding white emissions noise with standard deviation ayp 0.2.
F t e~ c()
(58)
The second method for generating synthetic calcium traces added an intermediary step between calcium influx and observed fluorescence. As previously, spikes were integrated with a slow-varying exponentially decaying variable e as described in Equation 57] A hill function was used to capture the nonlinear binding kinetics of calcium to the indicator dye,
d,
(59)
with hill coefficient 2, and 7 0.0001. These parameters were chosen based on parameter fits in Deneux et al. (2016). Finally white noise was added to these resulting traces to provide synthetic fluorescence traces,
F. = 4 c0..00
(60)
22

---

Under review as a conference paper at ICLR 2021
1 CALFADS EXAMPLE OUTPUTS AND SCATTER PLOTS ON SYNTHETIC BENCHMARK DATASETS
1.1 LORENZ ATTRACTOR: LINEAR CALCIUM MODEL
CaLFADS Example Output
A
B
E
F
Fluorescence
Rates
(Hz)
12
40
Spike Rate
OFF
Truth
oo
Reconstruction
12
Reconstruction 40
C
H
Spike Counts
Lorenz State
oW Spike Counts
Factors
Truth
Truth
-1
(s)
10
Time
10
Reconstruction
Reconstruction
LFAIS Example Output
A
B
E
Fluorescence
F
Rates
0. Spike Rate (H2)
12
OFFF
Truth
Truth
-1
Reconstruction
Reconstruction
C
D
G
H
Spike Counts
Lorenz State
Ow Spike Counts
1
Factors
he
Truth
Truth
Reconstruction
Reconstruction
Time
10
Time (s)
10
Figure S2(Top), S3 Example traces of A) Observed Fluorescence, B) Spike Rates, C) Spike Counts, D) Dynamic Factors. Red: Ground-truth, Blue: Reconstructed. Scatter plots compar- ing reconstructed and ground-truth of E) Observed Fluorescence, F), Spike Rates, G) Spike Counts, H) Lorenz State
Figure S2 shows examples the performance of CaLFFAS reconstructing the fluorescence traces (Fig A2), spike rates (Fig S2BB, spike counts S2C) and Lorenz attractor states (Fig S2D) when using a linear model of synthetic calcium dynamics. Visually, it can seen that the model achieves very close fit to the fluorescence traces, spike rates, and Lorenz dynamics. The model also captures spike-timing, with these spike trains appearing smoothed.
For comparison, Figure S3 show examples of OASIS+LFADS reconstructions. A key difference is that the spike counts are not smooth as spike counts are discretised from the deconvolved transient amplitudes, and there is greater variance in reconstruction quality high spike rates shown by
23

---

Under review as a conference paper at ICLR 2021
the scatter plot in Figure S3F. Nevertheless, OASIS+LFADS also provides a good fit fluorescence traces, spike rates, spike counts and lorenz dynamics.
1.2
LORENZ ATTRACTOR: NONLINEAR CALCIUM MODEL
CaLFADS Example Output
A
B
E
F
Fluorescence
Rates
Spike Rate (Hz)
60
40
Truth
OFFF
OF
-1
Reconstruction
60
0 Reconstructtoon 40
Spike Counts C
D
G
Spike Counts
H
Lorenz State
Factors
Truth
Truth
Time 10
10
Time (s)
Reconstruction
Reconstruction
LFSIS Example Output
A
Fluorescence
F
Rates
Spike Rate (H2)
60
dF/F
Truth
Truth
8
000
Reconstruction
Reconstruction
C
G
Spike Counts
Lorenz State
Spike Counts
Factors
Truth
Truth
000
10
9
Reconstruction
1
Time (s)
Time (s)
Reconstruction
Figure S4(Top), (5(Bottom): Example traces of A) Observed Fluorescence, B) Spike Rates, C) Spike Counts, D) Dynamic Factors. Red: Ground-truth, Blue: Reconstructed, Scatter plots compar- ing reconstructed and ground-truti of E) Observed Fluorescence, F), Spike Rates, G) Spike Counts, H) Lorenz State
Figure S4 shows examples of the performance of CaLFADS in reconstructing the fluorescence traces (4A) spike rates (Fig S4B), spike counts (Fig S4C) and Lorenz attractor states (Fig (4D) when using a nonlinear model to generated synthetic calcium traces. As with the linear calcium model, it can be that the model achieves close reconstruction of to the Lorenz dynamics. However, the reconstructions of other network variables have deteriorated to a much greater extent. similar trend can be seen in Figure $5 for FFSIS.
24

---

Under review as a conference paper at ICLR 2021
J
MOUSE VISP DATASET, DECODING NETWORK AND EXAMPLE FACTORS
CaLFADS example output
A
Fluorescence
B
Spikes
C
Factors
A
Actual Reconstructed
Time
Time
Time
Time
OASIS*LFADS example output
A Fluorescence
B Spikes
C
Factors
Activity
Activitty
Actuvty
Actvvty
(6)
Time
Time
Time
samples of actual and reconstructed fluorescence traces (S6A, left column) and 6 samples of fluorescence traces (S6A right column and S7A), B) 6 samples of inferred spike count traces corresponding to the fluorescence traces (right column, for S6A), C) 12 samples of inferred latent factor traces.
Figure S6(Top), S7(Btttom):
J.1 MOUSE VISP DATASET
The 2-photon calcium imaging data used here is part of the OpenScope project dataset (Anonymous) 2018), collected by the Allen Institute for Brain Science (AIBS) in Seattle, WA. All animal proce- dures were approved by the Institutional Animal Care and Use Committee at the AIBS. Neuronal activity was recorded head-fixed, awake Cux2-Cre mice on running wheel (de Vries et al. 2020)
(Fig 3A) expressing the calcium indicator GCamP6f (Chen et al. 2013) in layer 2/3 pyramidal cells of VisP BB) The mice were first habituated to repeating, expected stimulus over days, after which unexpected trials were introduced. The stimuli were adapted from (Homann al. 2017). In each expected trial, 4 consecutive sets of 30 Gabor patches appeared in sequence for 300 ms each, followed by 300 of grey screen (A, B, C, D, grey). For each set, the locations and sizes of the
25

---

gr Under review as a conference paper at ICLR 2021
Gabor patches were held constant within session. However, within each trial, each Gabor patch's orientation was sampled from a von Mises distribution, with trial mean sampled from {0, 45, 90, 135} and standard deviation 0.25%. In expected trials, occurring 10% of the time, the D set was replaced with a distinct set E, with its own locations and sizes. addition, the Gabor patch orien- tations for E sets were sampled from (Fig Processed fluorescence (dF/F) traces were extracted from the calcium imaging recordings for the putative neurons identified, as described dd Vries al. 2020).
Mises distribution with mean shifted positively by 90º
J.2 MOUSE VIsP DATA: EXAMPLE OUTPUTS
Figures S6 and S7 show samples of actual fluorescence traces (A), the corresponding inferred spike count traces (B), and inferred latent factors (C) generated as the output the CaLFADS and OA- SIS+LFADS models, respectively. Figure S additionally shows examples of reconstructed and actual fluorescence traces (A, left column). As can be seen the figures, the inferred latent factors of CaLFADS are smoother and less noisy compared to those of the OASIS+LFADS. This is impor- tant as the quality of the inferred factors is crucial for decoding purposes, such as the classification task described in Appendix J.3
J.3 MOUSE VIsP DECODING NETWORK AND EXAMPLE FACTORS
A
C
DIE
Factors
:
Time (s)
1.5
Figure S7: Schema of the non-linear model trained to decode expected vs. unexpected stimulus trials. Example latent factors inferred by CaLFADS are plotted for expected (blue) and unex- pected (red) trial. Each stimulus frame (A, B, C, D/E, grey) is labelled, and its onset and offset in the trial marked with dotted lines. Latent factors passed through GRU and linear modules (green), followed by sigmoid decision function.
Figure S7 shows schematic of the non-linear model used to decode from the latent factors whether a trial contained an expected or unexpected frame (decoding performances reported in Fig 3D). In short, the non-linear decoder comprised of a single-layer GRU, followed by linear layer and sigmoid decision function. We show an example set of latent factors (Figure [S7 left) inferred in expected (blue) and unexpected (red) trials. From these, one can see that CaLFADS was able to infer smooth latent factors from the calcium imaging data that did not show signs of over-fitting to noise.
26


Reviews:
[{"review_id": "pCfJ7xISc8", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "------------- UPDATE FOLLOWING DISCUSSION -------------\n\nThe 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:\n(i) Simulations and numerical results are especially confusing. Please consider a major revision in light of the comments raised by different reviewers.\n(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)\n\n------------- END OF UPDATE -------------\n\nSummarizing 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.  \n\nThe 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:  \n\n--- 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.  \n--- 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.  \n--- Comparison with the baseline method (OASIS+LFADS): OASIS was designed for fast, online spike inference from calcium imaging data.  \n     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)  \n     b) The OASIS paper mentions different working modes, where only some of the possible parameters are optimized. How was OASIS run in this manuscript?  \n--- 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.  \n\nOther concerns:  \n--- 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.  \n--- 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}?  \n--- What does 14 represent in t_{14} in Table 1?", "extra_scores": {}, "extra_text": {}}, {"review_id": "XoSscj2joHl", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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.\n\n1. 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?).\n\n2. 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. \n\n3. 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.  \n\n4. 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. \n\n5. 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 :)\n\nTo provide some concrete recommendations:\n\n1. I'd like more equations, the more the merrier for me, could be in the appendix.\n2. 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. \n3. 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.\n\nIn 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. \n\n===============================================================\n==== updates after rebuttal period =================================\n\n1. The authors added more equations, I now understand the methodological contributions required substantial effort.  The VAE/etc. work is outside my area of expertise.\n2. 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.\n3. I got more clarity on the differences between LFADS and CaLFADS, though again, this is outside my area of expertise. ", "extra_scores": {}, "extra_text": {}}, {"review_id": "QczkeV0Tu0M", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This submission extends a previous nonlinear data analysis method ( Latent Factor Analysis of Dynamical Systems, LFADS) to deal with the calcium imaging data (CaLFADS). \n\nQuality: 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.\nClarity: 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).\nOriginality: The method for handling the calcium traces is novel, although it is a relatively straightforward modification of LFADS.\nSignificance:  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.\nPros: \n1. CaLFADS is a holistic approach to infer the latent dynamics from fluorescence traces and avoids the two-step procedure which is nice. \n2. The method can obtain measures of the uncertainty- an issue that was ignored in most of the previous work on calcium data analysis.\n3. The method shows some improvement over a few alternative methods when fitting the data in certain regimes.\nCons:\n1. 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.\n2. 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.\n3. 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. \n4. 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.\n\nMore specific comments:\n* 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. \n* 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.\n* 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.\n\nClarification questions: \n* “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?\n* The variables in Fig1 Panel B are not explained in the figure caption.\nThe 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).\n*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.\n* “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. \n* 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. \n\n%%%% after the rebuttal:\n\nI 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. \n", "extra_scores": {}, "extra_text": {}}, {"review_id": "1c2x6RaAdo6", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "##########################################################################\n\nSummary: \n\nTraditionally, 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.\n\n##########################################################################\n\nPros: \n\nI 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.\n\n##########################################################################\n\nCons: \n\nMy 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. \n\nSecond, 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.)\n\nLast, 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)$?\n\n##########################################################################\n\nQuestions/Comments:\n\n1) 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.\n\n2) 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})$.\n\n3) 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.\n\n4) At the bottom of page 16 above equation 10, $gamma$ -> $\\gamma$.\n\n5) Many of the figures are not vector graphics.\n##########################################################################\n\nTo 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! \n\n##########################################################################\nPost-rebuttal: I have increased my score to a 7. I have laid out my reasoning in response to the author's comments.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

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