A key problem in neuroscience, and life sciences more generally, is that data is generated by a hierarchy of dynamical systems. One example of this is in \textit{in-vivo} calcium imaging data, where data is generated by a lower-order dynamical system governing calcium flux in neurons, which itself is driven by a higher-order dynamical system of neural computation. Ideally, life scientists would be able to infer the dynamics of both the lower-order systems and the higher-order systems, but this is difficult in high-dimensional regimes. A recent approach using sequential variational auto-encoders demonstrated it was possible to learn the latent dynamics of a single dynamical system for computations during reaching behaviour in the brain, using spiking data modelled as a Poisson process. Here we extend this approach using a ladder method to infer a hierarchy of dynamical systems, allowing us to capture calcium dynamics as well as neural computation. In this approach, spiking events drive lower-order calcium dynamics, and are themselves controlled by a higher-order latent dynamical system. We generate synthetic data by generating firing rates, sampling spike trains, and converting spike trains to fluorescence transients, from two dynamical systems that have been used as key benchmarks in recent literature: a Lorenz attractor, and a chaotic recurrent neural network. We show that our model is better able to reconstruct Lorenz dynamics from fluorescence data than competing methods. However, though our model can reconstruct underlying spike rates and calcium transients from the chaotic neural network well, it does not perform as well at reconstructing firing rates as basic techniques for inferring spikes from calcium data. These demonstrate that VLAEs are a promising approach for modelling hierarchical dynamical systems data in the life sciences, but that inferring the dynamics of lower-order systems can potentially be better achieved with simpler methods. Many datasets in the life sciences are generated by a hierarchy of dynamical systems, wherein lower-order dynamical systems that directly generate the data are driven by higher-order dynamical systems that are not observable. This problem is outlined in figure 1A, in which noisy observations x depend on the state z 1 of a low-order dynamical system that is perturbed by inputs u 1. The state of this dynamical system is also coupled to the state z 2 of a higher-order dynamical system, which can be perturbed independently by inputs u 2. One example of such a system in in-vivo two-photon calcium imaging from neuroscience. Calcium imaging provides systems neuroscientists with the ability to observe the activity of hundreds of neurons simultaneously during behavioural experiments. Such experiments have allowed neuroscientists to ask questions about the underlying computations and algorithms that neural circuits are implementing in perception, decision-making, memory, and many other processes. Such experiments can be characterized as observing a hierarchical dynamical system (Fig 1B) in which measurable calcium fluorescence is primarily determined by dynamics based on voltage-gated calcium channels and calcium binding to fluorescence dyes, and the rate of fluorescence transients controlled by the underlying computation. Recent applications of sequential variational autoencoders to neural data analysis has seen great success in inferring underlying computations in populations of cells in macaque and human motor cortex. By characterizing neural computation as low-dimensional dynamic factors in a non-hierarchical dynamical systems, showed that these dynamic factors trained to generate the inhomogeneous intensity functions explaining the rate of spikes assumed to follow a Poisson process. Crucially, these low-dimensional factors could also decode reaching behaviour of macaques and humans with much higher fidelity than any other dimensionality reduction method. Although this is a significant advance in our ability to analyze neural data in the form of spikes trains, two-photon calcium imaging poses the additional problem of identifying latent spike trains in fluorescence traces. This problem has been independently addressed in a number of different ways, including deconvolution and variational inference. If we continue to model the frequency of events as being generated by a Poisson process, this can be seen as hierarchy of dynamical systems (Fig 1A), in which low dimensional dynamics generate spike probabilities that in turn drive fluctuations in biophysical dynamics of calcium activity (Fig 1B. Here we propose a method that extends LFADS to accommodate calcium activity using this hierarchical dynamical systems approach, in which we can infer both the latent dynamics and the latent spike trains from the observed calcium fluorescence signal. The model is a variational ladder autoencoder (VLAE) with recurrent neural networks (RNNs) that supports uncovering latent dynamical systems (Fig 1C, full directed acyclic graph in Fig A1). It can be seen as a unification of two recent applications of variational autoencoders (VAEs) in neuroscience: 1) Latent Factor Analysis for Dynamical Systems (LFADS) and 2) DeepSpike, a VAE approach to inferring spike counts from calcium imaging data . We choose the VLAE approach since it has been shown to avoid the problem of trying to reconstruct the data solely from the lower-order features by separating latent variables via divergent deterministic pathways in the encoder, and convergent deterministic pathways in the generator (Sønderby et al., 2016; . Model hyperparameters are shown in table A1. The inferred dynamical system underlying the frequency of calcium events in the data is identical to that of LFADS (Fig 1C, blue modules). The prior distribution of initial conditions g 0 and external inputs u t are modelled as Gaussian distributions P (g 0) = N (µ g0, σ 2 g0), and P (u t) = N (µ ut, σ 2 ut). The underlying dynamical systemġ = G(g t, u t) is modelled by a Gated Recurrent Unit (GRU) taking the initial hidden state g 0 and inputs u t. Low dimensional factors f t are calculated as a linear transformation of the generator hidden state f t = W f ac g t. These factors are used to reconstruct the Poisson process intensity function with a fully connected layer and exponential non-linearity Inferred spike counts s t are generated by sampling z t from Gaussian distributions P (z t) = N (µ zt, σ 2 zt) and projecting these through an affine transformation and non-linearity along with the factors from the deeper layer, i.e., We assume a simple model of calcium dynamics:ẏ = −y t /τ y + α y s t + β y where the parameters τ y, α y, β y are measured from the data, however it is a topic for future research to fit the calcium dynamics simultaneously. In our synthetic data, these are valued at 0.4 s, 1, and 0 respectively. The value of τ y is chosen as it is the known decay time constant of GCamP6, a commonly used calcium fluorescence indicator used in calcium imaging experiments. The calcium fluorescence signal x t is mapped onto the parameters of the variational posterior distributions Q(z t |x), Q(g 0 |x), Q(u t |x). These distributions are all modelled as Gaussians, with the mean and standard deviations parameterized by a stack of bidirectional GRUs, con. The final hidden states of E 1 gen and E 2 gen are mapped onto the parameters of Q(z 0 |x) and Q(g 0 |x) respectively with fully connected layers. The hidden states E 1 con and E 2 con are concatenated and passed as inputs to single direction GRUs C 1 and C 2. The hidden states of C 1 and C 2 are concatenated at each time step t with s t−1 and f t−1. Subsequently these concatenated activities are mapped onto the parameters of Q(z t |x) and Q(u t |x) with fully connected layers. One of the advantages of using VLAEs is that the evidence lower bound (ELBO) formulation is the same as for VAEs despite the hierarchical latent space. As such, our cost function remains very similar to that of LFADS. g0,u∼Q(z,g0,u|x) [log(s|λ)] Where W gg is the weight matrix in the GRU that take g as input. The likelihood function P (x t |y t) is modelled as a Gaussian distribution x t ∼ N (y t, σ 2 y), where σ 2 y is learned. Although s t is not discrete, P (s t |λ t) is treated as an approximate Poisson process s t ∼ P oisson(λ t) = s λt t exp(−λ t)/Γ(s t + 1). This is an inductive bias to ensure that the inferred spike trains in the lower-order dynamic system is constrained to follow the dynamics of the higher-order dynamic system. The latent variables are treated as independent, with hierarchy imposed in the deterministic paths, i.e., Q(z, g 0, u|x) = Q(z|x)Q(g 0 |x)Q(u|x) Parameters of our model were optimized with ADAM, with an initial learning rate of 0.01, which decreased by a factor of 0.95 whenever plateaus in training error were detected. As in LFADS training, KL and L2 terms in the cost function were'warmed up', i.e., had a scaling factor being 0 and 1 applied that gradually increased. Warm-up for the deeper parameters (blue modules in Figure 1) was delayed until warm-up for shallower parameters was completed (red modules in Figure 1). The model was tested on synthetic data with Lorenz dynamics embedded in the frequency of calcium fluorescence transients, as described by , where generated spikes were convolved with an exponential kernel with a time constant of 0.4 ms, and white noise added to the ing traces. We measure the performance of the model in three ways: 1) uncovering the underlying Lorenz dynamics, 2) reconstructing the rates of calcium transients an inhomogeneous Poisson intensity functions, 3) reconstructing the spike counts contributing to increases in the calcium fluorescence signal. The model was compared against a ground-truth where the spike counts are known, and LFADS is used to reconstruct the latent dynamics and intensity function, and against a model where spike counts are extracted using a deconvolution algorithm before using LFADS to reconstruct the rates and intensity function (OASIS + LFADS). It was also tested against a model that used a 1-D convolution of the intensity function to reconstruct either the first two (Gaussian-LFADS) or four (Edgeworth-LFADS) time-varying moments of fluorescence, as used previously in estimating the intensity functions of filtered Poisson processes in neuroscience . Figure 2 shows examples of performance of our model in reconstructing the fluorescence traces (Fig 2A), Poisson intensity functions (Fig 2B), spikes (Fig 2C) and Lorenz dynamics (Fig 2D). Visually, the model provides very close fit to the fluorescence traces, intensity functions, and Lorenz dynamics. The model also captures spike-timing, although these spike trains appear smoothed. Table 1 compares the R 2 goodness-of-fit on reconstructing held-out validation data with ground-truth latent dynamic structure. Of all approaches, our model easily performs best at reconstructing fluorescence traces. It should be noted that the reason the Gaussian and Edgeworth LFADS models perform so poorly at reconstructing fluorescence is that it only attempts to reconstruct the time-varying moments of fluorescence as a 1-D convolution of the intensity function, which stems from the approximation of exponential-kernel shot noise processes . Additionally, our model almost performs as well as LFADS in reconstructing the Lorenz dynamics. It is to be expected that LFADS performs better than our method, since there is an additional source of observation noise in our synthetic dataset generating fluorescence transients from spikes. Notably, our model does not perform as well as the deconvolution method OASIS in reconstructing spike trains, however this does not impact the ability of our model to reconstruct the latent dynamics. In fact, constraining the reconstructed spikes by the latent dynamics may mitigate any errors in spike train reconstruction that occur by deconvolution, since the deconvolution algorithm may erroneously drop spikes during high rates, whereas our model should be less likely to do so. It will be necessary to assess this possibility further. It should be noted that the deconvolution algorithm performs much better at reconstructing spike trains in our synthetic dataset than in real datasets where ground-truth spiking is known . To our knowledge, there are no known dual recordings of population 2-photon calcium imaging with ground-truth electrophysiology in a subpopulation of neurons in-vivo during behaviour driven by hypothesized low-dimensional dynamics that we would be able to validate this with. Nevertheless, since the relationship between calcium dynamics and somatic spiking is highly non-linear, especially in dendrites, it remains to be seen how useful it is to faithfully reproduce unseen spike trains in calcium fluorescence activity. The model was then tested on synthetic data generated from a 50 cell recurrent neural network with chaotic dynamics, and an external perturbation at a random time, as described in , with parameters adjusted to make the data more representative of firing rates and time scales observed in calcium imaging experiments. Spike trains were transformed into fluorescence signals using the same procedure as with the Lorenz system dataset. (Fig 3A), spike trains (Fig 3B), and intensity functions (Fig 3C). Visually, the fluorescence traces and spike trains have been reconstructed reasonably well, whereas reconstructions of the the intensity functions are highly noisy. The latent dynamics represented by the factors (Fig 3D) and external inputs (Fig 3E) also show a lot of noise. This appears to be due to difficulty in inferring external inputs, in which it is not clear whether the timing of external perturbation has been accurately inferred, despite a slight transient at roughly the time of the perturbation in the example. Table 2 compares the R 2 goodness-of-fit on reconstructing held-out validation data, which demonstrates that our model performs very well in reconstructing the fluorescence signal and does reasonably well at reconstructing spikes. However, in this more complex benchmark where the deeper dynamic system is perturbed by external input, the pre-processing of fluorescence transients with the OASIS deconvolution algorithm reconstructs firing rates far better than our method. We present a hierarchical recurrent variational autoencoder model capable of reconstructing latent dynamics, latent spike trains, and calcium fluorescence traces in a benchmark synthetic dataset. Of the four methods tested, our model is the only one capable of reconstructing all three. Furthermore, our model performed best in reconstructing latent dynamics in our synthetic dataset We will need to assess our model on further synthetic benchmark data to assess the validity of our approach. Since our model is trained end-to-end, it should be possible to extend to reconstructing raw 2-photon imaging videos, which could enable us to train models to uncover latent dynamics from arbitrarily shaped neuronal structures. This would of great use to neuroscientists who are largely restricted to techniques that extract fluorescence traces from regions of interest with somatic shapes, whereas the morphological diversity of dendrites is much greater. An additional advantage of using our hierarchical model is that we can obtain measures of the uncertainty in both the latent dynamics, and the latent spike trains. The correlation in uncertainty between layers of this hierarchy may be what allows superior inference of the latent dynamics, despite less accurate reconstructions of the spike trains than OASIS, which provides no measure of uncertainty. We hope to improve our model to better capture the relationships between layers of this hierarchy in future. We describe a use-case in neuroscience (2-photon calcium imaging data) for which this model may be very useful. However, we are keen to investigate the general case of hierarchical dynamical systems and their utility in uncovering structure in datasets outside this domain. A APPENDIX
We extend a successful recurrent variational autoencoder for dynamic systems to model an instance of dynamic systems hierarchy in neuroscience using the ladder method.