Multi-modal Gaussian Process Variational Autoencoders for Neural and Behavioral Data
Authors: Rabia Gondur
Usama Bin Sikandar
Evan Schaffer
Mikio Christian Aoi
Stephen L Keeley
Venue: ICLR 2024 poster (2024)
Keywords: Gaussian Processes
Latent Variable Models
Variational Autoencoders
Neuroscience
Primary area: applications to neuroscience & cognitive science
TL;DR: A novel latent variable model that extracts smoothly evolving latent subspaces that are shared between or independent to distinct data modalities for use with modern systems neuroscience experiments.

Abstract: Characterizing the relationship between neural population activity and behavioral data is a central goal of neuroscience. While latent variable models (LVMs) are successful in describing high-dimensional data, they are typically only designed for a single type of data, making it difficult to identify structure shared across different experimental data modalities. Here, we address this shortcoming by proposing an unsupervised LVM which extracts shared and independent latents for distinct, simultaneously recorded experimental modalities. We do this by combining Gaussian Process Factor Analysis (GPFA), an interpretable LVM for neural spiking data with temporally smooth latent space, with Gaussian Process Variational Autoencoders (GP-VAEs), which similarly use a GP prior to characterize correlations in a latent space, but admit rich expressivity due to a deep neural network mapping to observations. We achieve interpretability in our model by partitioning latent variability into components that are either shared between or independent to each modality. We parameterize the latents of our model in the Fourier domain, and show improved latent identification using this approach over standard GP-VAE methods. We validate our model on simulated multi-modal data consisting of Poisson spike counts and MNIST images that scale and rotate smoothly over time. We show that the multi-modal GP-VAE (MM-GPVAE) is able to not only identify the shared and independent latent structure across modalities accurately, but provides good reconstructions of both images and neural rates on held-out trials. Finally, we demonstrate our framework on two real world multi-modal experimental settings: Drosophila whole-brain calcium imaging alongside tracked limb positions, and Manduca sexta spike train measurements from ten wing muscles as the animal tracks a visual stimulus.

Published as a conference paper at ICLR 2024
MULTI-MODAL GAUSSIAN PROCESS VARIATIONAL AU- TOENCODERS FOR NEURAL AND BEHAVIORAL DATA
Rabia Gondur"
Usama Bin Sikandar
Fordham University
Georgia Institute of Technology usama gatech. edu
rgondurê fordham. edu
Evan S. Schaffer
Mikio Aoi
Icahn School of Medicine at Mount Sinai evan. schafferhhrrrcce smm.s..memmmmmmmmmmmmmmmmmmm...m edu
University of California, San Diego maoi@uesd.cc
Stephen Keeley
Fordham University
skeeley1@fordhamy edu
ABSTRACT
Characterizing the relationship between neural population activity and behavioral data a central goal of neuroscience. While latent variable models (LVMs) are successful in describing high-dimensional time-series data, they are typically only designed for a single type of data, making it difficult to identify structure shared across different experimental data modalities. Here, we address this shortcoming by proposing an unsupervised LVM which extracts temporally evolving shared and independent latents for distinct, simultancously recorded experimental modal- ities. We do this by combining Gaussian Process Factor Analysis (GPFA), an interpretable LVM for neural spiking data with temporally smooth latent space, with Gaussian Process Variational Autoencoders (AP-VAEs). which similarly use GP prior characterize correlations in latent space, but admit rich expressivity due to deep neural network mapping to observations. We achieve interpretability in our model by partitioning latent variability into components that are either shared between or independent to each modality. We parameterize the latents of our model in the Fourier domain, and show improved latent identification using this approach over standard GP-VAB methods. We validate our model on simulated multi-modis data consisting of Poisson spike counts and MNIST images that scale and rotate smoothly over time. We show that the multi-modal GP-VA GP-VAE (MM-GPVAE) is able to not only identify the shared and independent latent structure across modalities accurately, but provides good reconstructions of both images and neural rates held-out trials. Finally, we demonstrate our framework on two real-world multi- modal experimental settings: Drosophila whole-brain caleium imaging alongside tracked limb positions, and Manduca sexta spike train measurements from ten wing muscles the animal tracks visual stimulus.
INTRODUCTION
Recent progress in experimental neuroscience has enabled researchers to record from a large number of neurons while animals perform naturalistic behaviors in sensory-rich environments. An important prerequisite to analyzing these data is to identify how the high-dimensional neural data is related to the corresponding behaviors and environmental settings. Traditionally, researchers often employ two-step approach, involving first dimensionality reduction on neural data followed by post-hoc investigation of the latent structure of neural recordings with respect to experimental variables of
interest such as stimuli or behavior (Wu et al. 2017 Zhao and Park 2017 Saxena and Cunningham
"Currently affiliated with Cold Spring Harhor Laboratory

---

Published as a conference paper at ICLR 2024
2019 Bahg et al. 2020 Pei et al. 2021). However, recent deep learning advancements allow for both experimental and behavioral variables to be part of a single autent-variate model (Schneider et al. 2023 Hurwitz et al. 2021; Sani et al.) 2021 Singh Alvarado et al. 2021). thus opening new avenue for unsupervised discoveries on the relationship between neural activity and behavior.
The existing approaches th jointly model neural activity and behavior are limited in that they often rely on a single latent space to describe data from both modalities (Schneider et al. [023]] Hurwitz] al. 202I), making it difficult for practitioners isolate latent features that are shared across and independent to neural activity or behavior. Approaches that do isolate modality specific and shared
latent structure either do so with no temporal structure (Gundersen et al. 2019 Kleinman et al.) 2024]
Singh Alvarado et al.) [021; Wu and Goodman 2018: [Shi et al. [019 Brenner et al. (0022), or use a
relatively inflexible linear dynamical system (Sani et al. al.] 2201]. J. Moreover, because of the deep neural
network mapping from latents to observations, many existing multi-modal approaches generate latent
spaces that are not obviously related to experimental variables of interest, and so additional model
features are often added to aid in interpretability and analysis in neuroscience settings (Schneider
ct al. 2023 Hurwitz et al., 2021; Zhou and Wei 2020).
However, LVMs developed for neural data often are able to uncover interpretable features using neural activity alone in an unsupervised fashion with minimal a priori assumptions. One example of these is Gaussian Process Factor Analysis (GPFA), a widely used LVM in neuroscience that finds smoothly evolving latent dynamics in neural population activity, and can illustrate different aspects of
neural processing (Yu et al. 2008 Duncker and Sahani 2018 Keeley al., 2020a] b
Balzani et al.
2022). GPFA constrains the latent space of neural activity through the use of a Gaussian Process (GP)
prior. GP priors have also been adapted to regularize the latent space of variational autoencoders
(GP-VAE or GP-prior VAE) with a variety of applications (Casale et al. 2018 Fortuin et al.] 2020
Ramchandran et al.] 2021] Yu et al., 2022). In each of these approaches, the GP prior provides a flexible constraint in the latent dimension. specifying correlations across auxiliary observations like viewing angle, lighting, or time. The GP prior often used for out-of-sample prediction in GP- VAE,, but in GPFA is frequently used to visualize latent structure in noisy neural population activity on trial-by-trial basis (Duncker and Sahani 2018 8 Zhao and Park, 2017: Keeley ct al.) (020a). Here, we wish to leverage the interpretability seen in unsupervised GPFA models with the power of GA-S to use with multi-modal time-series datasets in neuroscience.
We propose model for jointly observed neural and behavioral data that partitions shared and independent latent subspaces across data modalities while flexibly preserving temporal correlations using GP prior. We call this the multi-modal Gaussian Process variational autoencoder (MM- GPVAE). Our first innovation is to parameterize the latent space time-series in terms of a small number of Fourier frequencies, an approach that has been used before in the linear GPFA setting (Paciorek 2007 Aoi and Pillow 2017 Keeley et (202aa). We show that the Fourier representation dramatically improves latent identification for the standard GP-VAE. Our second innovation is to augment our new Fourier GP-VAE model to describe two data modalities simultaneously by mapping the latents linearly to neural data (like in GPFA), as well as nonlinearly to another experimental variable via a deep neural network (like GP-VAE). We leave the specific identity of this other experimental modality intentionally vague- it could be keypoints of limb positions as animals freely move or visual stimuli across time. Because this observation modality is characterized by deep neural network, our model can be adapted to any experimental variable of interest.
We validate our MM-GPVAE model on a simulated dataset of smoothly rotating and scaling MNIST digit alongside simulated Poisson neural activity, whereby the digit data and the spiking data share latent structure. We show that the MM-GPVAE model is able recover both shared and independent latent structures while simultaneously providing accurate reconstructions of data from both modalities. Lastly, we demonstrate the utility of our model by fitting the MM-GPVAE two real-world multi-modal neural datasets: 1) Drosophila (fly) whole-brain calcium imaging alongside tracked 16 limb positions, and 2) Manduca sexta (hawkmoth) spike train measurements from ten wing muscles alongside a continuously moving visual stimulus. In the former case, we are able to view Drosophila behavioral conditions across shared and independent subspaces and in the latter case, we show distinct time-varying latent components tracking muscle and stimuli movement in the experiment. By showing the MM-GPVAE in these two domains, we demonstrate that our model is adaptable to range of diverse experimental preparations in systems neuroscience.
2

---

Published as a conference paper at ICLR 2024
2 THE GAUSSIAN PROCESS VARIATIONAL AUTOENCODER
The Gaussian Process variational autoencoder (GPVAE) uses high-dimensional data (c.g. images) accompanied by auxiliary information, like viewing angle, lighting, object identity, or observation time. This auxiliary information provides the indices in the GP prior latent representation, specifying correlations across the latent space and allowing for out-of -ample predictions new auxiliary values (Casale et al. 2018 Fortuin et al., 2020 Ramchandran et al. 2021). However, in our case we consider continuously observed time-series data in experimental neuroscience experiments, our auxiliary information here is evenly sampled time-bins. Because of this, we can leverage the advantages of the Fourier-domain GP representation (Aoi and Pillow 2017 Keeley al., 2020b Paciorek] (007).
Formally, consider smoothly-varying image data across timepoints, represented by the pixels-by-time matrix Y N=*2 For latent variable z(t) E each latent zp(t) (f E (1,2...7)) .')) evolves according to a Gaussian process, (t) ~ GP(O, ka(-:-)), with covariance kernel kg. The time-by- time covariance matrix of each (p( then given by the Gram matrix K₁ corresponding to kg. We use squared exponential (RBF) kernel for K governed by marginal variance and length scale {p. E} with an additive diagonal term a to help with inference (Yu et al. 2009 Casale et al.) 2018). The likelihood of the image data is Gaussian with mean given by the latent values at any timepoint t passed through deep neural network (-(() with parameters ei and whose covariance is T
The GP- is learned using standard VAE amortized inference (Kingma and Welling 2013), where the parameters H&o and o of a variational distribution (s are given as neural network functions the observed data Y. parameterized by o. Here, the evidence lower bound (ELBO) maximized with respect to the variational parameters Φ, model parameters tb, and GP hyperparameters. In this work, set to fixed value of 1c-2 for all experiments except for the final data analysis example, where it is set to value of 1e-4. While the ELBO may be expressed in variety of ways, we will follow the form that includes the variational entropy term. For details about this approach for the standard GP-VAE, see Casale et al. (018).
2.1 FOURIER-DOMAII REPRESENTATION OF THE GP-VAE
We consider a version of the GP-VAE whose auxiliary variables are a Fourier frequency representation of the time domain, opposed to timepoints sampled parameterize frequency representation of the latent variables that is Gaussian distributed according to 3,p ~ (0, Ko) where the original RBF GP covariance matrix K may be diagonalized by K = BKB Here, B is the orthonormal discrete Fourier transform matrix and ž, is single Fourier-domain latent vector whose length equal to F, the total number of Fourier frequencies representing the latent. The model prior can now be written as
regular lattice. This allows to
2 0 a)) - II (=, K. +aIF)
Where Ž represents a P x F matrix of the frequency representation of the P latent variables the frequencies and K2 is an covariance matrix. The model likelihood retains the same form as the standard GP-VAE,
p PY v. 0²)
II 1=1
((
where z, represents the th time index of Z, P T matrix where each row the time-domain
representation of the nd latent, given by 2p Bž,.
We define the variational posterior so that it factorizes over the Fourier frequencies. The variational distribution can be written
q0 7 Y
II N (z) Pso (Y'). diag (o% (Y"))).
(3)
where i indexes a trial of times-series images and w indexes the Fourier frequencies (w e 1, .... F). In contrast to the standard GP-VAr (Casale et al. 2018), the Fourier representation
3

---

Published as a conference paper at ICLR 2024
enoodor
map to spaco
variational freeçe
time-domain
docodor
kstent
Free Mep
DNN
[4]
(d)
Fourler GPVEE
Poure-GPPAE (pruned) Fourer
50
Fourer
GPVEE
GPYME
6
6
6
OPVAE
G
300
VAE
6
20
40
20
40
20
Timebine
Timebine
Tiieebin
Figure 1: (a) schematic for the Fourier domain GP-VAE. All images at all timepoints for a given trial are encoded via a deep neural network into variational parameters of a pruned Fourier representation of the latent space. This Fourier representation is then mapped back into the time domain before being passed through a decoder network to give the image reconstruction at each timepoint. (b) Image reconstructions of the standard VAE, GP-VAE and Fourier domain GP-VAE. (c) Estimated latent for each model alongside the true underlying latent angle. (d) Mean squared error (MSE) of estimated latents and true latents for 60 held-out trials. Error bars indicate standard error.
requires that, for each trial, images at all timepoints are mapped to single P x F-dimensional Fourier representation. We accomplish this in two steps first, we use a deep neural network for each image at each time point yi. which will result in 7 total network embeddings for a single trial, each of dimension P. We follow that with single linear layer each for the mean and variance the variational distribution, la()), laa (-), that will map from the number of timepoints to the number of Fourier frequencies. 1/, laz PxT + RPF A schematic for the Fourier domain variational and generative architecture displayed in Figure
There are number of advantages to representing the GP- VAE latents and variational parameters in the Fourier domain. For one, the diagonal representation of K means we avoid costly matrix inversion when evaluating the GP prior (Aoi and Pillow 2017 Paciorek] 2007 Hensman et al. 2017] Wilson and Adams 2013 Loper et al.) 2021). Secondly, we can prune the high frequencies in the Fourier domain, effectively sparsilying the variational parameters, and hence frequencies (F < T), while enforcing smooth latent representation. Lastly, the Fourier representation of the variational parameters mean we can retain both temporal correlations as well as the advantages of using a mean-field approximation for the variational posterior qo. These advantages have been seen in simpler GP models. For more information, see (Aoi and Pillow 2017 Keeley et al., (020a)..
The Fourier-represeented GP-VAL GAF dramatically improves the ability of the GA-VAE learn true underlying smoothly evolving latent for high-dimension data in a non-linear model. We demonstrate this on a simulated example using an MNIST digit that rotates by time-varying angle given by draw from a GP with an RBF kernel. We fit our Fourier GP-VAE model as well as the standard VAE and standard GP- VAE models to these trials. In each case, we are looking to recover the true underlying generative latent angle and hence we assert a one-dimensional latent dimensionality for the model. For more information on training and testing, see appendix.
Because the Fourier variational distribution preserves temporal correlations and prunes the high- frequency components of the latents z, the Fourier GP-VAE model shows much smoother underlying latent representations on held-out trials (Fig (c)) though each model retains the ability accurately reconstruct the images through the network mapping from the latent space (Fig [b)). When the latent space is mapped through an affine transform to align latents held-out trials to the true latent space,
the Fourier domain GP- VAE does much better uncovering the true generative latent angle (Fig (c).
(d)).
THE MULTI-MODAL. GAUSSIAN PROCESS VARIATIONAL. AUTOENCODER
We now focus on extending the GP-VAE to model data of two modalities simultaneously. The examples we will emphasize here involve neural activity alongside some other experimental variable
4

---

Published as a conference paper at ICLR 2024
(a)
(b)
I...F
w
linear layer (GPFA)
8, a
W
DNN
spikes
linear layer
IFFT
linearity
estimated rates
P
DNN
z
DNN
P
latents
behavior
estimated behavior
a²
(2)
P. P. Y )
Figure 2: (a) Graphical model of the multimodal GP-VAN PAE. A set of Fourier frequencies describe the Fourier representation of shared and independent latents across modality with GP prior over each latent. Latents transformed to the time domain and combined to generate data for each modality. (b) Schematic the MM-GPVAE.
such as naturalistic movement or high-dimensional stimuli. The observations of our model are two distinet data modalities represented by matrices YA IR and Y' t E (1.2... T) denotes time-bin indices, i denotes trials, and N and M denote the dimension of the observations for each modality. We assume, before, that all data generated by smoothly varying low-dimensional latent variables that are now either independent to or shared hetween data modalities. As before, latents initially parameterized in a pruned Fourier representation before being mapped to the time domain. The latents are then partitioned by a loadings matrix W which linearly combines the shared latent representation with the independent latents.
(d) RM Here,
x TB
WA
WsI W
WB
0
AA
zs
+d
xB_
Here, zs refers to latents that are shared between modalities, while 2,1 and æB denote modality- specific latent variables and d represents an additive offset. The outputs after this linear mixing result mpoeaiit--specifi embeddings which we call X,A and XB. Here, modality A, the "neural modality". passed through a pointwise nonlinearity to enforce non-negativity of Poisson rates, and modality B, the "behavioral modality", is passed through deep neural network. The likelihood of the data given the embedding is:
p(yAlæA) P(f(xA)), p(yal=)) (Z ININ,
(5)
Where the function (-) represents decoder neural network for the behavioral data and we use the exponential nonlinearity for f. Learning is performed as before the mean and variance of mean-field Gaussian variational distribution is given by the data passed through a neural network with parameters 6. Here, our lower bound is similar to that found in Casale et al. (2018) with an additional term for the neural data modality.
Poisson Likelihond (Neural Rates)
Gaussian I Likelihood (ather modality)
ELBO- E Egngg
Σ (lg +
log W (yn gh ,
I
(6)
GP Prior
Entropy
log
T. T)
(qs)
A graphical depiction of the generative model of the MM-GPVAE is shown in Figure 2 a). The hyperparameters of the model are the latent-specific kernel parameters f that provide the GP co- variance structure well entire schematic of the MM-GPVAE model, including the generative and variational distributions, is depicted in Figure b).
constant additive diagonal offset (fixed during inference). The
5

---

Published as a conference paper at ICLR 2024
(a)
(b)
(c)
image-one onllynnl lalent
estimated
only (GPFA)
image only (GP-VAE)
(scale)
3
3
3
3
GPVAE
Image reconstruction
neural reconstruction
shared latent estimate
true
3
(angle)
shared
lalent
?
4400
1.5
MSE
4000
neura-ooly
latent
0.5-
neurall
2
3600
-0.5
0.5
20 limebins
40
60
0
20
40
60
Model
Model
Model
timebins
Figure 3: (a) True and estimated latents for the MM-GPVAE trained on simulated neural spiking data as well as a smoothly scaling and rotating MNIST digit. (b) Estimated neural rates on an example trial for example neurons well as 4 example reconstructed images of the MNIST digit at different angles and scales. (c) (left) Reconstruction accuracy from the image data trained on the images alone (GPVAE) compared to training with both modalities simultancously (MM-GPVA)). (middle) Accuracy of estimated neural rates (left) trained on neural activity alone (Poisson GPFA) compared to MM-GPVAE. (right). Accuracy of shared latent estimated from the MM-GPVAE compared to single-modality model variants. Error bars are standard error.
4 EXPERIMENTS
4.1 SIMULATED DATA
We first assess the performance the MM-GPVAE model using a simulated dataset two modalities: smoothly rotating and scaling MNIST digit, and simultaneously accompanied Poisson spike counts from 100 neurons. A total three latents were drawn from an underlying GP kernel to generate single-dimensional shared, neural, and image subspaces. There were a total of 300 trials, each trial consisting of 60 timebins. The data was split into 80% for training and 20% for testing. For this simulated example, one latent represents an interpretable modulation of the image as directly affects the scaling of the MNIST digit. So as the latent values change along the trial the image scales smaller and larger continuously. Another latent, corresponding to the independent component of the neural modality, provides one component of the log rates of the Poisson spiking data. The final latent reflects the shared variability between both modalities. Here, the latent is again interpretable with respect to the image in that this latent corresponds to the angle of the smoothly rotating MNIST digit. This shared latent also linearly combines with the neural-only subspace to provide the log rates of the spiking Poisson population. The MM-GPVAE trained simultaneously on both images and spikes can successfully recover the dimenssoaa latent structure across these two modalities. Figure 3)) a) shows the underlying true latent in held-out trial alongside the estimated latent structure extracted from the MM-GPVAE. As before, the held-out latent space is scaled to align with the true latent space as closely as possible, the MM-GPVAE model latent space is invariant to scaling transform. In addition to accurately extracting the latent structure, the MM-GPVAE also has the ability to accurately reconstruct both neural rates and images across the time series from the latent space. Four example images and three example neural rates are shown in Figure 3 (b) For more examples, see Figure g and
The MM-GPVAE is also able to exploit information across modalities to better model the observed data. The left side of Figure shows the mean squared error (MSE) across pixels for reconstructed images in held-out trials. By leveraging latent information from the spiking modality, the MM- GPVAE is able to reconstruct the images better than single-modality GP-VAE trained on the images alone (Casale et al. 2018). This improvement is additionally dependent the number of neurons analyzed (Figl 12). Similarly, the held-out neural rates are better estimated when the image data is included, though the effect here is more modest. Finally, we also show that the shared latent variable can be more accurately identified when data from both modalities is used (c)). Here we assess the accuracy of the shared latent estimate test trials when use data from each modality alone

---

Published as a conference paper at ICLR 2024
(A)
lmm position tracking
(b)
neural data
behavioral data
(9)
eetimated
true
0.8-
Vanabilty
0.4
Inb posiions
shared
Ind
shared
(d)
behavioral subspace
neural subepacee
sharod subspace
1.5
10
0.5
neural aetiii
PC2
0
0
stil bending
-1.5
nunning
-4
0
0
-4
0
0.2 (seconds)
0.4
PC1
PC1
PC1
Figure 4: (a) Limb position tracking Drosophila (b) (top) Contribution of the variability in the data across trials of the neural (left) and behavioral (right) modalities due to shared and independent subspaces (c) (top) Limb position estimates and true values and (boltom) six randomly selected calcium trace estimates and true values for given trial. (d) Visualization of average latent value across time in the neural, behavioral and independent subspaces for 3 behavioral categories.
(akin to standard GP-VAE and GPFA) compared to both modalities simultaneously (MM-GPVAE). As expected, data from both modalities better identifies shared structure across the data types. This additionally does not depend the structure of nature of the encoding distributions, (Fig 10 and is markedly improved by representing the MM-GPVAE with pruned Fourier frequencies (Fig 11
4.2 APPLICATION TO FLY EXPERIMENTAL DATA
Next, we look to evaluate the MM-GPVAE on a real-world whole-brain calcium imaging from adult, behaving Drosophila (Schaffer et al.] 2021). We isolate 1000 calcium traces recorded using SCAPE microscopy (Voleti et 20.9) from an animal while performs a variety of distinct behaviors fixed on spherical treadmill. These 1000 calcium traces are dispersed widely and uniformly across the central fly brain and were randomly picked from the dataset use in our model. Alongside the neural measurements, eight 2-D limb positions are extracted from a recorded video using the software tool DeepLabCut (Mathis et al.) (018). These simultaneously recorded calcium traces and limb position measurements are split into 318 trials of 35 time-bins sampled at 70hz. Each trial has one of 5 corresponding behavioral labels (still, running, front grooming, back grooming and abdomen bending) determined via a semi-superviied approach from the tracked limb position measurements (Whiteway et al.) (221)). Importantly, these behavioral labels are not used when fitting the MM-GPVAE. For additional information on how we isolate this behavioral data and neural traces, the appendix. We consider 7-dimensional latent behavioral subspace and 26-dimensional neural subspace, with 5 of these dimensions being shared across modalities. These choices were made through initial exploration of the model and examining cross-validated model performance (see appendix), though the results we show are robust to wide range of dimensionality choices for each of the subspaces.
dataset. Here, we consider a
The MM-GPVAE is able to successfully reconstruct both behavioral trajectories and calcium traces in held-out trials for these data. The top of Figure 4c) shows the true and decoded 16 limb position measurements with 3 highlighted for clarity. The bottom part of Figure shows 6 randomly selected caleium traces alongside their model reconstructions. In each case, the capture the temporal trends in this multi-modal dataset. To analyze the shared and independent latent subspaces of these data we consider in that subspace calculated each trial. We additionally color-code the trial according the behavioral label. In these shared and independent latent subspaces, we find that the "still" behavioral state is well separated from many of the other behavioral conditions in the neural and shared subspaces. Here, show two others (bending and running) as an illustration Figg(d) (For visualizing all 5 behaviors in the latent space, Figure 19. Lastly, illustrate how much each subspace contributes to the
can roughly
projection each latent space, and plot the mean latent value
7

---

Published as a conference paper at ICLR 2024
100
<
10
ViGual 30muuuu
moddo IDCOrGllation
601
stared
noural
exparieertal selegg
nesuion ID
10
Purpa Pusces
0.4
(192)
(So0))
Figure 5: (a) Experimental set-up for Manduca sexta. Spikes from 10 muscle groups are recorded as an animal tracks a Hz moving flower stimulus (b) top: Spikes from 3 example motor neurons middle: estimated Poisson rates bottom: The neural latent along torque measurement from the hawkmoth. The 22Hz modulation reflects wing-flapping (c) Visual stimulus and reconstructions from MM-GPVAE (d) Weights of hawkmoth spike decoder for neural-only and shared latents. (e) A one-dimensional latent from the visual- visual-modality subspace closely tracks the stimulus position. (f) The shared latent between modalities plotted alongside the torque measurement.
overall variability of neural and behavioral data. Figure shows the variance of the neural (left) and behavioral (right) data calculated for each trial where the reconstruction of the data is generated either from the shared subspace, or either of the independent subspaces. We find that the shared subspace for each modality contributes more to the overall variability in the data than either of the independent subspaces, suggesting large fraction of shared variability between the data modalities.
4.3 APPLICATION TO MOTH EXPERIMENTAL DATA
Next, we evaluate the MM-GPVAE on a dataset of a hawkmoth (Manduca sexta) tracking a moving flower. The hawkmoth is an agile flier that is able to closely follow swiftly moving targets while hovering in midair, making it model organism for the study of sensorimotor control (Sprayberry and Daniel 2007). The modalities in our dataset consist of a time series of images of visual stimulus temporally paired with electromyography (EMG) signals recorded during 20-second experimental trials. Each image is an event-based reconstruction (inspired by insect visual system) of white floral target moving laterally during trial as viewed in the hawkmoth's 's frame of reference (Sikandar al. 2023) aa)). The floral target moves sinusoidally at Hz. against a black background. Each pixel in a video frame can occupy polarity of -1, O or +1 based on whether its luminosity has respectively decreased, unchanged or increased as compared to the previous frame in the temporal sequence. Temporally paired with the floral stimulus, EMG recordings are precisely measured spike trains from the 10 major flight muscles of the hawkmoth responding to the floral target while flying on tether and flapping at about 22 Hz. Altogether, the 10-muscle recordings form near-complete motor program that controls its flight during target tracking response (Putney et al. 2019).
We fit the MM-GPVAE model these data using a 2-dimensional latent space to describe the moving visual stimulus and a 2-dimensional latent space to describe the Poisson rates from the 10 motor units. Of these, one latent dimension is shared between modalities. Here, small number of dimensions are used because only a few dimensions were needed to be able to accurately reconstruct the data of each modality (for more information about dimensionality selection, see appendix).
We find that the MM-GPVAE is able to accurately reconstruct both the neural rates (Fig 5 (() and visual stimulus (Fig c)) from its low dimensional latent space. We see that the neural rates are strongly modulated by the wing-flap oscillation at approximately 22 Hz. This is also seen in the neural latent space identified by the MM-GPVAE. At the bottom of Figure 5)b) we see that the 1 1-dimensional latent space closely aligns with the torque measurement along the z-axis of the moth motion, a measurement previously identified to be closely associated with wing beating (Putney et al., 2019). The MM-GPVAE is also able to reconstruct the visual stimulus data, seen ((c) When viewing the latent space corresponding to the visual stimulus modality, we find that the value of the smoothly evolving latent variable closely tracks the position of the stimulus in dimension (the only dimension along which it moves, see Fig e()). Finally, the shared latent dimension shows some modulation with wing-flapping but exhibits temporal structure not obviously related to the flapping
8

---

Published as a conference paper at ICLR 2024
frequency or stimulus motion. This may correspond to longer- motor behavior such as the variation in how strongly the hawkmoth is responding to track the stimulus, which could depend on the degree of its attention to the stimulus and its state of motivation. Importantly, we find that including a shared dimension better reconstructs the data than not including one, suggesting that find that if
there is component of neural dynamics shared with stimulus dynamics. Specifically, we remove the shared latent we see approximately 70% reduction the reconstruction accuracy of the image modality and about 5% reduction in the reconstruction of the rates (Fig 2T.. Because the neural rates have strong wing-flapping modulation, we expect the slower-moving dynamics to account for only a small amount of the variability here. Lastly, to demonstrate the interpretability advantage of linear decoder layer used for the neural likelihood, we visualize how each of the 10 motor neurons loaded onto either the neural-only dynamics or the dynamics shared with the moving 5 d). For closer look individual neurons modulations by shared and independent latent components, see Figure21[b)
stimulus Figure
5 CONCLUSION
In this work, we have introduced the multi-modal Gaussian Process variational autoencoder (MM- GPVAE) to identify temporally evolving latent variables for jointly recorded neural activity and behavioral or stimulus measurements. We first parameterize the model in the Fourier domain to better extract identifiable temporal structure from high-dimensional time-series data. We then combine a linear decoder for neural data with a deep neural network decoder for behavioral or stimuli data to enhance the applicability of the MM-GPVAE across experimental settings while retaining the ability to identify neural tuning with respect to shared and independent temporally evolving latents. We show that our multi modal GPVAE can accurately recover latent structure across data modalities, which we demonstrate using a smoothly rotating and scaling MNIST digit alongside simulated neural spike trains. We then show that the MM-GPVAL can be flexibly adapted to multiple weal-world multi-modal settings. First, we analyze calcium imaging traces and tracked limb positions of fruit flies (Drosophila) exhibiting various distinct behaviors. Then, we apply the model to spike trains from hawkmoth flight muscles as well as the oscillating visual stimulus the animal tracks during flight.
Comparison to multi-modal approaches in neuroscience: Our work contrasts in important ways with two similar accetll-developed multi-modal latent variable models for time series data in neuroscience. One closely related model, Targeted Neural Dynamical Modeling (TNDM) (Hurwitz al. 2021), ), was specifically designed to non-linearly separate the neural latent representation into behaviorally relevant and irrelevant subspaces using a latent space governed by a recurrent neural network (RNN). However, this work was evaluated on a dataset with a relatively simple behavioral paradigm, and fails to accurately model the behavioral modality when it is higher dimensional or more complex (Fig 6). Another related multi-modal time-series neuroscience model, Preferential Subspace Identification (PSID) (Sani et al. 2021). linearly separates the neural latent representation into behaviorally relevant and irrelevant subspaces. However, PSID may be limited in that uses latent dynamics are governed by linear state-space model, and, similar to TNDM, is restricted in capturing complex behavioral modalities well due a linear decoder to behavioral data (Fig (). For a more thorough evaluation and discussion to these approaches, see the appendix.
Choice of prior: In this work we focus on the use of GP priors to describe latent dynamics, though other temporal dependencies are of course possible. Though TNDM and PSID have certain limitations when compared with the MM-GPVAE (Fig 6). these are not necessarily due to their latent temporal prescriptions. These models could be modified with more flexible decoding schemes to better isolate the role of the different latent dynamics in characterizing shared and independent subspaces compared to the MM-GPVAE. In addition, other choices such as switching linear or non-linear latent dynamics models are powerful in single modality settings (Hernandez et al.) 2018 Linderman et al.) 2017] Glaser et al.) (0200) and could also be adapted to multi-modal settings. Our choice of a GP prior specifically is well suited to experimental paradigms where the latent dynamics are hypothesized to be smooth. Other approaches might be more appropriate when this smoothness assumption is unlikely to hold. Thus, development of neuroscience and evaluating their strengths and limitations alongside the MM-GPVAE remains an important avenue of future work in neuroscience.
versions of different latent dynamical models in

---

Published as a conference paper at ICLR 2024
REPRODUCIBILITY STATEMENT
All models this manuscript were trained end-to-end PyTorch using the ADAM optimizer. Training done Macbook Pro with Apple MI max chip and all evaluations took less than an hour to fit. All encoder and decoder neural networks were standard feedforward neural networks with ELU activation functions. Additional details on neural network architectures for the experiments can be found in the appendix. An implementation of MM-GPVAE can be found at: GitHub Repository for MM-GPVAE
ETHICS STATEMENT
Our work provides general model for exploratory data analysis for multi-modal time-series data. Though we develop the model with neuroscience experiments in mind, in principle the model could be adapted to other settings where the assumption of smoothly evolving latent dynamics holds across distinct data modalities. We do not anticipate any potential negative societal impacts of our work.
10

---

Published as a conference paper at ICLR 2024
REFERENCES
Aoi, M. and Pillow, J. W. (2017). Scalable bayesian inference for high-dimensional neural receptive fields. bioRxiv, page 212217.
Bahg, G., Evans, D. G., Galdo, M., and Turner, B. M. (2020). Gaussian process linking functions for mind, brain, and behavior. Proceedings of the National Academy of Sciences, 117(47):29398- 29406.
Balzani, E., Noel, P., Herrero- Vidal, P., Angelaki, D. E., and Savin, C. (2022). probabilis- tic framework for task-aligned intra-and inter-area neural manifold estimation. arXiv preprint arXiv:2209.02816.
Brenner, M., Koppe, G., and Durstewitz, (2022). Multimodal teacher foreing for reconstructing nonlinear dynamical systems. arXiv preprint arXiv:2212.07899
Casale, F. P., Dalca, A., Saglietti, L., Listgarten, J., and Fusi. N. (2018), Gaussian process prior variational autoencoders. Advances in neural information processing systems, 31.
Duncker, L. and Sahani, M. (2018). Temporal alignment and latent gaussian process factor inference in population spike trains. In Advances in Neural Information Processing Systems, pages 10445- 10455.
Fortuin, V., Baranchuk, D., Rätsch, G., and Mandt, S. (2020). Gp-vae: Deep probabilistic time series imputation. In International Conference on Artificial Intelligence and Statistics, pages 1651-1661. PMLR.
Glaser, J., Whiteway, M., Cunningham, J. P., Paninski, L.., and Linderman, S. (2020). Recurrent switching dynamical systems models for multiple interacting neural populations. Advances in neural information processing systems, 33:14867-14878,
Gundersen, G., Dumitrascu, B., Ash, J. T., and Engelhardt, B. E. (2019). End-to-end training of deep probabilistic cca on paired biomedical observations. In Uncertainty in artificial intelligence.
Hensman, J., Durrande, N., and Solin, A. (2017). Variational fourier features for gaussian processes. The Journal of Machine Learning Research, 18(1):5537-5588.
Hernandez, D., Moretti, K., Wei, Z., Saxena, S., Cunningham, J., and Paninski, L. (2018). Nonlinear evolution via spatially-deeendent linear dynamics for electrophysiology and calcium data. arXiv preprint arXiv: 1811 1.0249..
Hurwitz, C., Srivastava, A., Xu, K., Jude, J., Perich, M., Miller, L., and Hennig, M. (2021). Targeted neural dynamical modeling. Advances in Neural Information Processing Systems, 34:29379-29399
Keeley, S., Aoi, M., Yu, Y., Smith, S., and Pillow, J. W. (2020a). Identifying signal and noise structure in neural population activity with gaussian process factor models. Advances in neural information processing systems, 33: 15:13395-18800
Keeley, S., Zoltowski, D., Yu, Y., Smith, S., and Pillow, J. (2020b). Efficient non-conjugate gaussian process factor models for spike count data using polynomial approximations. In International Conference Machine Learning, pages 5177-5186. PMLR.
Kingma, D. P. and Welling, M. (2013). Auto-encoding variational bayes. 112.v:1312.6114
arXiv preprint
Kleinman, M., Achille, A., Soatto, S., and Kao, (2024). Gacs-korner common information variational autoencoder. Advances in Neural Information Processing Systems, 36.
Linderman, S., Johnson, M., Miller, A., Adams, R., Blei, D., and Paninski, L. (2017). Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial Intelligence and Statistics, pages 914-922, PMLR.
Loper, J., Blei, D., Cunningham, J. P., and Paninski, L. (2021). A general linear-time inference method for gaussian processes on one dimension. The Journal of Machine Learning Research, 22(1):05500-106)
11

---

Published as a conference paper at ICLR 2024
Mathis, A., Mamidanna, P., Cury, K. M., Abe, T., Murthy, V. N., Mathis, M. W., and Bethge, M. (2018). Deeplabeut: markerless pose estimation user-defined body parts with deep learning. Nature neuroscience, 21(9):1281-1289.
Paciorek, C. J. (2007). Bayesian smoothing with gaussian processes using fourier basis functions in the spectralgp package. Journal of statistical software, 19(2):nihpa22751.
Pei, F., Ye, J., Zoltowski, D., Wu, A., Chowdhury, R. H., Sohn, H., O' O'Doherty, J. E., Shenoy, K. V., Kaufman, M. T., Churchland, M., et al. (2021). Neural latents benchmark 21: evaluating latent variable models of neural population activity. arXiv preprint arXiv:2109.0444
Putney, J.. Conn, R., and Sponberg, S. (2019). Precise timing is ubiquitous, consistent, and coordinated across a comprehensive, spike-resolved flight motor program. Proceedings the National Academy of Sciences, 116652):2951-26690 -69660.
Ramchandran, S., Tikhonov, G., Kujanpää, K., Koskinen, M., and Lähdesmäki, H. (2021). Longitudi- nal variational autoencoder. In International Conference on Artificial Intelligence and Statistics, pages 3898-3906. PMLR.
Sani, G., Abbaspourazad. H.. Wong, Y. T., Pesaran, B., and Shanechi, M. M. (2021). Modeling behaviorally relevant neural dynamics enabled by preferential subspace identification. Nature Neuroscience, 24(1): 24(1):140-149.
Saxena, S. and Cunningham, (2019). Towards the neural population doctrine. Current opinion in neurobiology, 55:103-111.
Schaffer, E. S., Mishra, N., Whiteway, M. R., Li, W., Vancura, M. B., Freedman, J., Patel, K. B., Voleti, V., Paninski, L., Hillman, E. M., et al. (2021). Flygenvectors: the spatial and temporal structure of neural activity across the fly brain. bioRxiv, pages 2021-09.
Schneider, S., Lee, H., and Mathis, M. W. (2023), Learnable latent embeddings for joint behavioural and neural analysis. Nature, 617(7960):360-33) 368.
Shi, Y., Paige, B., Torr, P., et al. (2019). Variational coattuur--f-expperre autoencoders for multi-modal deep generative models. Advances in neural information processing systems, 32.
Sikandar, U. B., Choi, H., Putney, J., Yang, H.,, Ferrari, S., and Sponberg, S. (2023). Predicting moduuaty precisely-timed spikes across a coordinated and comprehensive motor program. In 2023 International joins conference on neural networks (IJCNN). IEEE.
Singh Alvarado, J., Goffinet, J., Michael, V., Liberti III, W., Hatfield, J., Gardner, T., Pearson, J., and Mooney, (2021). Neural dynamics underlying birdsong practice and performance. Nature, 599(7886):635-63,
Sprayberry, J. D. and Daniel, T. L. (2007). Flower tracking in hawkmoths: behavior and energetics. Journal of Experimental Biology, 210(1):37-45.
Voleti, V., Patel, K. B., Li, W., Perez Campos, C., Bharadwaj, S., Yu, H., Ford, C., Casper, M. J., Yan, R. W., Liang, W., et al. (2019). Real-time volumetric microscopy of in vivo dynamics and large-scale samples with scape 2.0. Nature methods, 16110):1054-166
Whiteway, M. R., Schaffer, E. S., Wu, A., Buchanan, E. K., Onder, O. F., Mishra, N., and Paninski, L. (2021). Semi-supervised sequence modeling for improved behavioral segmentation. bioRxiv, pages 2021-06.
Wilson, A. and Adams, R. (2013), Gaussian process kernels for pattern discovery and extrapolation. In International conference on machine learning, pages 1067-1075. PMLR.
Wu, A., Roy, N. A., Keeley, S., and Pillow, (2017). Gaussian process based nonlinear latent structure discovery in multivariate spike train data. Advances in neural information processing systems, 30.
Wu, M. and Goodman, N. (2018). Multimodal generative models for scalable learning. Advances neural information processing systems, 31.
12

---

Published as a conference paper at ICLR 2024
Yu, B., Cunningham, J., Santhanam, G., Ryu, S., Shenoy, K., and Sahani, M. (2009). Gaussian- process factor analysis for low-d single-tri trial analysis of neural population activity. In Frontiers in Systems Neuroscience. Conference Abstract: Computational and systems neuroscience.
Yu, B. M.. Cunningham, J. P., Santhanam, G., Ryu, S., Shenoy, K. V., and Sahani, M. (2008). Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. Advances neural information processing systems, 21.
Yu, C., Soulat, H., Burgess, N., and Sahani, M. (2022). Amortised inference structured generative models with explaining away. Advances in Neural Information Processing Systems.
Zhao, Y. and Park, I. M. (2017). Variational latent gaussian process for recovering single-trial dynamics from population spike trains. Neural computation, 2965):129--136.
Zhou, D. and Wei, X.-X. (2020). Learning identifiable and interpretable latent models of high- dimensional neural activity using pi-vae. Advances Neural Information Processing Systems, 33:7234-7247.
13

---

Published as a conference paper at ICLR 2024
APPENDIX
CONTRASTING THE MM-GPVAE WITH EXISTING APPROACHES
We evaluate our MM-GPVAE model alongside both Targeted Neural Dynamical Modeling (TNDM) (Hurwitz et al. 2021) and Preferential Subspace Identification (PSID) (Sani et al. 2021 using our synthetic data of a rotating and scaling MNIST digit '3' and corresponding Poisson neural rates. In Figure 6 show the reconstruction of both the behavioral modality, the MNIST digit "3", and Poisson neural rates for each model. In addition to these, we also provide the error of the reconstruction performance on both data modalities held-out trials. We find that the MM-GPVAE is able to better reconstruct data from both modalities in this simulated example than the competing approaches.
(a)
mm-
3
3
3
(b)
est
true
gpvae
7.5
tndm
mm-
gpvae
2.5
psid
3
20
7.5
(c)
tndm
neurra rate recon
image recon
mse (competing
mse (competing model)
3000
30000
model)
7.5
1500
psid tndim
15000
psid tndm
psid
2.5
1500
3000
0
15000 30000
20
40
80
20
40
60
mse (mm-gpvae)
(mm-gpvae)
time
time
Figure 6: (a) Reconstruction of the scaling/rotating MNIST digit '3' with MM-GPVAE, TNDM and
PSID. (b) Reconstruction of neural rates with MM-GPVAE, TNDM, and PSID. (c) MSE for neural rate reconstruction (left), and MSE for image reconstruction (right). Here, each dot indicates one trial msee from MM-GPVAE vs competing model. The majority of the trials errors fall above the unity line for both models, indicating overall better reconstruction with MM-GPVAE.
Next, to compare these models in real-world setting, implement PSID, TNDM, and MM-GPVAE on data from a simpler neuroscience experiment where models are capable recovering both behavioral trajectories well neural rates. We focus on simpler behavioral task because these competing approaches are not well-suited to reconstruct complex behavioral modalities due to a lack of neural network decoder. We use the primate reaching data from Hurwitz et al.) (2021) for these evaluations. Here, linear decoder as well as our deep neural network can reconstruct both modalities well. The visualizations of the latent spaces for these data can be seen in Figure 7 While both TNDM and PSID can provide only two subspaces in which only one can separate the 8 reaching positions (relevant subspace). MM-GPVAE can provide visualization for 3 subspaces, 2 modality specific, and shared, and we can see a clear separation with both the behavioral-only subspace as well as shared subspace.
14

---

Published as a conference paper at ICLR 2024
(a)
psid
(b)
tndm
relevant
non-relevant
relevant
non-relevant
2
0.2
PC2
PC2
0
-0.5
-0.2
0.5
-0.8
5
2
PC1
PC1
PC1
PC1
(c)
mm-gppae
neural subspace
behavioral subspace
shared subspace
0.4
0.2
PC2
PC2
-0.4
-0.2
-0.5
0.5
-0.2
0.2
0.6
-2
PC1
PC1
PC1
Figure 7: (a) Separation of all 8 reaching directions in the relevant subspace using PSID. (b) Separation of all 8 reaching directions in the relevant subspace using TNDM. (c) Separation all reaching directions in behavior-only and shared subspaces. Dots here indicate the mean latent value across the entire trial. The neural subspace shows no behavioral separation in the latent-space whereas the behavioral and shared subspaces show strong behavioral separation. This result closely parallels with (Hurwitz et al. 2021 Sani et al. 2021). which isolates behaviorally relevant and irrelevant neural subspaces. However, in contrast (Hurwitz et al., 2021 Sani et al. (021) the MM-GPVAE isolates a distinct shared subspace as well as both neural and behavioral independent subspaces from a raw unsupervised partitioning of both the behavioral and neural data.
15

---

Published as a conference paper at ICLR 2024
(8)
astimated
3
3
3
3
16
neural
40
10
20
Image-only
latent
0.B
0.5
neural-nnj -oura-only
1.5
1.5
latent
0.5
0.5
-0.5
-0.5
shared lalent
50
50
M
50
100
-50
500
550
1000
100
20
40
60
20
40
60
0
20
40
60
0
20
40
60
20
40
60
timebins
Figure 8: (a) Additional estimated images and neural rate trajectories of the MNIST digit '3' as well as independent and shared latent estimates.
ADDITIONAL EVALUATIONS OF MM-GPVAE ON SIMULATED DATA
To provide a more complete picture of the ability of the MM-GPVAE to both reconstruct simulated behavioral and neural data, and accurately recover the true underlying latent trajectories, we show additional performance evaluations here in simulated setting. Figure g shows an example of 24 reconstructed from the evaluation shown in Figure 3 of the manuscript. Again, here we reconstruct scaling and rotating MNIST digit '3' alongside 100 neural spike trains, where one latent dimension shared across modalities. Figure alss shows latent reconstructions and 5 example neural rates in 5 held out trials. Using the same set-up, we additionally run the MM-GPVAE with MNIST digit '2' and show an example of 24 reconstructed images as well as example neural rates (of 100 neurons) on 5 held-out trials in Figure 9
In addition to these, we extend Figure 3(a) from our manuscript with all the latent trajectories for our simulated data comparisons. Because
unimodal (GPFA and GP-VAF) comparisons, there is ambiguity as to which latent is "shared", we simply select the latent trajectory which most closely matches the true underlying latent variable for each comparison. Panels (b) and (c) in Figure [0 show the reconstruction mean squared error of the MM-GPVAE compared to the unimodal models GP-VAE and GPFA. Here, we include version of each of the unimodal models with two different encoder representations, one which represents what the GPFA and MM-GPVAE would typically have access simply the data of their own modality (i.e. neural data for GPFA, image/behavioral data for GP-VAE) and second where the encoding distribution is function of data from both modalities, matching the encoder exactly with that of the MM-GVPAE, In each these cases, we find that the encoding representation has little effect model performance. We visualize both the average mean mean-squared error on held-out data (c) as well as a per-trial mcan-squared error scatter- plot (b) to provide better sense of performance across all models and encoders. We note though that the extent to which MM-GPVAE more accurately identifies latent structure than the unimodal variants is sensitive the nature of the simulated data and can vary across model fits. For (b) and (c) specifically, we average model performance across 4 random seeds to report overall reconstruction accuracy. We also found through model exploration that GPFA specifically often does comparatively well in identifying latent structure when neural rates are high or have large variance. This is likely true because identifying latents within linear map trivial under the standard GPFA model, and so latent accuracy can be very good here if the neural data is sufficiently robust. This means that the metrics reported for the identification of latent structure in Figures 10 and 3 are somewhat dependent on the statistics of the neural rates used to generate and fit the model.
16

---

Published as a conference paper at ICLR 2024
estimated
true
R
A
neural
rate
20
image- only
1.75
latent
0.4
0.75
-.2
-028
neural-only latent
1.5
0.5
0.2
0.5
0.4
-0.5
-0.2
00.5
-0.4
100
40
shared
lstent
80
60
20
50
60
40
80
20
40
60
20 timebins
40
80
40
60
Figure 9: (a) Example of MNIST digit '2' with reconstructed digits and neural rates well as independent and shared latent estimation.
(x)
mmngggva gDvao
gpta true
modalttly
both modalities
image latent
moe (gpla)
mne sapeen
neura/ only
(mm-gpvae)
(mm govag)
meg
shared
0.5
mm-ggovee
grvad
griaad
Sola
gula
20
40
60
20
40
60
20
40
60
time
9me
Time
Figure 10: (a) Extension of Fig 3(a) from the manuscript with best latent trajectories, showing the best matched latent from GPFA and GP-VAE to the true shared and independent latents. (b) Identification acccuracy of the best performing latents estimating the shared latent for GPFA (left) and GP- VAE (right). Each dot indicates a trial. (c) Average MSE for recovery of shared latent for the best possible latents from GP- EPE, GP- VAE (encoding both modalities) GPFA, and GPFA (encoding both modalities).
ASSESSING THE IMPORTANCE PRUNING THE FOURIER FREQUENCIES
Figure demonstrates the necessity of pruning the Fourier frequencies when evaluating the MM- GPVAE on our simulated data. Here, we present standard and pruned Fourier MM-GPVAE model fits (Fig a)) and show their ability in reconstructing the true generative smooth latent across each latent subspace (Fig b)). Without constraining the number of Fourier features for these data, the latent identification fails, akin to what is seen when we ignore the Fourier representation entirely. We note here that this feature is function of the length of the trial, and the failure is not as stark for shorter-length trials. For these trials, our pruned condition the length scale of the GP to a value of 10 (i.e. the hyperparameter l can go smaller than 10).
17

---

Published as a conference paper at ICLR 2024
gpvae
gmvve (pruned)
(o)
image only latent
neural only latent
shared latent
(b)
image only latent
neurel only latent
shared ladent
0.8
0.50
1.8
1.50
0.5
mes
0.8
0.40
1.4
me
m98
0.4
0.30
1.2
0.2
0.20
model
model
model
20
40
60
20
40
60
0
20
time
timp
time
Figure 11: (a) Accuracy of identification of all three latents for the pruned and standard version of the Fourier MM-GPVAE, (b) Example latent trajectories for each subspace. The time-domain formulation fails to identify the true underlying latents accurately, similar to what is seen in Fig
VARYING THE NUMBER OF NEURONS
To demonstrate how the recovered latent of the shared and independent subspaces might trade-off with the amount of data practitioner may have for given modality, consider fits of our MM-GPVAE in our simulated setting with varying number of neurons. As expected, increasing the number of neurons increases the ability of the MM-GPVAE identify the neural latent, well as its ability to accurately reconstruct neural data. Importantly, because structure is shared across modalities, increasing the number of neurons also increases the ability of the MM-GPVAE identify shared latent structure, as well as shows some improvement in the MM-GF AM-GPVAE to generate image data. Increasing the number of neurons has no effect on identifying the image latent structure. (See Fig 12
image only latent
shared latent
neural only Istent
(b)
neural reconstruction
image reconstruction
4.6
0.26
1.3
0.8
4.2
4000
0.6
ege
0.20
3800
3.4
0.7
0.2
3
3200
100
100
5
100
100
100
number of neurons
number of neurons
number of neurons
number of neurone
number of neurons
Figure 12: (a) Latent identification accuracy varying number of neural rates are available the model. (b) image and rate reconstruction accuracy as a varying number of neural rates are available the model.
ROBUSTNESS MEASUREMENT ERROR
To assess how the MM-GPVAE might perform the presence of noisy imperfect behavioral measurements, consider a simulated experiment where a randomly selected 20% of the images are greyed on either the top or bottom half. We the MM-GPVAE jointly to the image data of this form well uncorrupted neural rates of the same form as in the main manuscript. We find that the MM-GPVAE able to exploit smoothness, rendering well-generated images even during these occluded trials (Fig 13 a)). As expected, find that standard MM-VAE (without GP prior over time) is not able to reconstruct the true image as well, but does a better job reconstructing the occluded image, treats each frame independently (Fig (3). We also show that this occlusion procedure has no bearing on the ability of the MM-GPVAE or MM-VAE to reconstruct the neural rates (Fig c)).
18

---

Published as a conference paper at ICLR 2024
(a)
(b)
true
occluded
mm- gpvae
3
3
3
3
3
image reconstruction
image reconstuction
5400
4800
5000
mse
4400
mm-- vae
9
4600
4000
4200
3600
(c)
model
model
oocluded
3
3
3.5
rate reconstruction
mm-
7
gpvse
mse
3
3
3
3
GG66
2.5
2
model
Figure 13: (a) Image reconstructions of MMMVAE and MM-GPVAE trained on occluded threes alongside the true and occluded threes. (b) (left) MSE MM-GPVAE and MM-VAE for the true, not occluded data, (right) MSE of MM-GPVAE and MM-VAE for the occluded data. (c) MSE of rate reconstructions.
NEURAL NET ARCHITECTURES
For the MM-GPVAE, both the simulated and real-world multi-modal evaluations used similar neural architectures. However, there were some modifications of the nodes/layers that were unique to each evaluation. This was necessary as the latent dimensionality was different across our different evaluations, and certain behavioral reconstructions, especially those that were higher-dimensional, required richer neural network parameterizations. The schematic for the MM-GPVAE neural network architecture in our simulated example can be found in Figure [4] and the schematics for the neural network architectures for the real-world multi-modk datasets Figure T6 (hawkmoth). For our evaluation on the dataset used (Hurwitz et al. 2021], whose results above, the architecture can be found in Figure 17 Note that for all experiments, each modality is encoded to its own set of variational means and variances (transformed into the Fourier domain). The encoded means and variances representing the shared latents are then summed give the encoded shared latents means and variances. Across all evaluations, we parameterized our latents in the Fourier domain and converted back to the time domain before decoding.
be found in Figure 15 (fly) and
input layer (60, 3024)
Input layer (60, 100)
latent space (12, 3)
LL ELU
LL ELU
hidden layer (60,784)
IFFT
JIFFT
hidden layer (60,80)
(c)
(d)
LL ELU
hidden layer (60, 2)
hidden layer (60,512)
LL ELU
hidden layer (60,256)
ELU
hidden layer (60,40)
hidden layer (60,256)
LL
hidden layer (60,
LL+ ELU
layer (60,20)
LL ELU
LL ELU
hidden layer (60,256)
hidden layer (60,256)
LL ELU
hidden layer (60,512)
mean (60, 2)
log varlance (60, 2)
mean (60,
log variance (60,
output layer
hidden layer (60,784)
LL
(60,100)
map to frequency space
latent space (12,
map to frequency
output layer (60,1296)
LL ELU
space
Figure 14: MM-GPVAE architecture for simulated data. (a) Encoder network of the MNIST digit. (b) Encoder network of the neural information. (c) Decoder network of the MNIST digit. (d) Decoder network of the neural information.
19

---

Published as a conference paper at ICLR 2024
inpul layer (35.16) layer (35, 12)
input layer (35,1000)
latent space (26, 28)
LL ELU
LL ELU
hidden layer (35, 100)
IFFT hidden layer (35, 7)
JIFT
(d)
ELU
LL ELU
hidden layer (35, 56)
hidden layer (35, 50)
hidden layer
hidden layer (35, 7)
(35. 26)
mean (35,
log variance (35,
mean (35, 26)
log vartancel (35,
LL ELU
LL
hidden layer (35, 56)
LL ELU
output layer
frequency space
latent space (26, 28)
map frequency space
output layer (36, 16)
(35,1000)
Figure 15: Neural network architecture for evaluation on fly dataset. (a) Encoder network for behavior. (b) Encoder network for neural data. (c) Decoder network for behavior. (d) Decoder network for neural data.
input layer (3324, 3024
input layer (3324, 10)
latent space (480, 4)
LL ELU
hidden layer (3324, 784)
LL ELU
IFFT
JIFFT
hidden layer (3324.80)
(c)
(d)
ELU
hidden layer (3324, 2)
hidden layer (3324,
LL+ ELU
hidden layer (3324,
LL ELU
hidden layer (3324,40)
hidden layer (3324,256)
hidden layer (3324, 3)
LL ELU
hidden layer (3324,20)
LL ELU
LL ELU
hidden layer (3324, 256)
hidden layer (3324,256)
LL
LL ELU
hidden layer (3324,512)
mean (3324,
log variance (3324, 2)
mean (3324,
log variance (3324, 3)
output layer
hidden layer (3324,784)
LL ELU
(3324, 10)
map to frequency
Iatent space (480,
map to frequency space
output layer (3324,3024)
LL ELU
space
Figure 16: Neural network architecture for evaluation on hawkmoth dataset. (a) Encoder network for behavior. (b) Encoder network for neural data. (c) Decoder network for behavior. (d) Decoder network for neural data.
(8)
(b)
input layer (73, 70)
latent space (16, 64)
Input layer (73. 2)
LL ELU
hidden layer (73, 60)
IFFT
IIFT
LL
LL ELU
layer (73, 40)
hidden layer (73, 8)
hidden layer (73, 48)
hidden layer (73, 48)
hidden layer (73, 20)
LL ELU
LL ELU
mean (73, 48)
log variance (73, 48)
mean (73, 48)
log variance (73,
oulput layer (73,
bulput layer (73, 70)
map to frequency space
space (16, 64)
map frequency
Figure 17: Neural network architecture for Monkey reaching data, evaluated above. (a) Encoder network for the behavior. (b) Encoder network for the neural activity. (c) Decoder network for the behavior. (d) Decoder network for the neural activity.
20

---

Published as a conference paper at ICLR 2024
ADDITIONAL INFORMATION ON FLY EXPERIMENTAL DATA
Pre-processing: We isolated 1000 raw calcium traces from (Schafter al., 202I) for use with our model for the evaluations on the Drosophila dataset. For this set-up, we consider variant of the MM-GPVAE where we have removed the non-linearity for the neural data modality, and instead consider Gaussian observations with an additional parameter controlling the observation variance for the fluorescence traces (akin to the original formulation GPFA (Yu et al.) 2009)). This dataset also contained x.y positions from 8 tracked limbs positions (Mathis et al. 2018)) The data came as one continuous recording, and at every time point there was a behavioral categorization probability for of 6 distinct behaviors (undefined, still, running, front grooming. back grooming, abdomen bending) determined via the algorithm outlined in (Whiteway et al.] (021). To split this continuous recording into trials, found segments of the recording where, for 35 continuous samples, single behavior was estimated at 2 %60 probability. This generated 318 total trials where each trial was of one of possible behaviors. There was no section of the recording where there was an 'undefined" behavior for 35 continuous samples, so there are no undefined behavioral trials in our analysis.
Additional evaluation: Though we show the time- comeccouse of the behavioral reconstruction in the main manuscript, we show a 2-d reconstruction on 2 example trials in Figure 18 Here we plot true x.y positions for trails for each of the 8 limb positions alongside our model's reconstruction. We can see here that our model also captures the spatial information in the fimb-positing data. We also show for comparison the 2-d depiction of the shared and independent latent representation all trials the dataset (Figure 19), with all five of the behaviors labelled. You can here that the "still" behavior separates well from all other behaviors in the shared and neural subspaces.
(a)
True behavior
True behsvvirr
y-positoo (nomalized)
y-posiiion (nomalized) 05
0.5
position
position
(normalized]
(normalized)
Reconstructed behaviar
Reconstructed behavior
y-poetion (normalzee) O
(normalizee)
y-postion
-0.8
x x-pasition (normalized)
x-position (normaizee]
Figure 18: (a) Reconstruction of 8 fly limb positions in 2 held-out trials. Here, we see the MM- GPVAE is able to reconstruct the spatial information of the behavioral modality the 8 tracked limb positions.
(8)
behavioral subspace
neural subspace
shared subspace
Still
10
Back
grooming
08
Front
grooming
-10
Bending
-12
-5
-2
Running
PC1
6
PC1
PC1
Figure 19: (a) Separation of all 5 fly behaviors in the behavior-only, neural-only, and shared subspaces. Similar to what is seen in the main manuscript, the 'still' behavior is well separated from trials with other behavioral labels in the neural and shared latent subspaces.
21

---

Published as a conference paper at ICLR 2024
ADDITIONAL INFORMATION ON HAWKMOTH EXPERIMENTAL DATA
Pre-processing: For our hawkmoth data, the original synthetic visual stimuli were sampled at 125 Hz and the neural and torque recordings were sampled at 10K Hz (Sprayberry and Daniel 2007 Sikandar et al., 2023 Putney et al., 2019). To prepare these data for evaluation with the MM-GPVAE we first downsampled neural information and torque measurements to 1K Hz by binning the spike counts and averaging the torque measurements at this temporal resolution. To align these measurements with the visual stimuli, we upsampled the images 8-fold. The entire recording was 20 seconds long and was split into 6 evenly-divided trials. Here, it is worth noting that the trials are much longer (~ 3300 timebins) than in the other experiments in this work.
The visual stimuli contained 3024 pixels and there were ten total hawkmoth motor neurons. We set the ceueel-indeenndent subspace to I-dimensional, the images-independed subspace to 1-dimensional, and an additional I dimension for the shared subspace. To encourage slow-evolving smooth latents in the shared and image subspaces, and faster-evolving neural latents, we initialized the length scale parameters for each latent dimension to different values. The length scale latents, 150 for the shared latent, and 300 for the image latent. This biased the model fits to capture slower dynamics in the image subspace and the faster dynamics the neural subspace.
10 for the neural
Additional evaluation; The neural rates are well-captured by the MM-GPVAE for all 10 hawkmoth motor neurons, which each show strong periodicity. Figure 20 shows the strong fast-oscillation spike rates captured by the MM-GPVAE for all ten hawkmoth neurons alongside the recorded spikes for 1.5 second period. Similar periodicity exists across the entire 20 second recording.
To emphasize the added interpretability we get from the linear layer in the neural likelihood and the role of the shared latent in generating the data from each modality, in Figure 21 we show the reconstructions for both modalities using one no shared latent variables (Fig (a)) as well as 2 example neurons' tuning to the neural and shared latent. Here we can see that neuron 5 is more positively modulated by the shared latent while neuron 9 is negatively modulated by the shared latent.
(a)
10
spikes
5-
estimated rates
0.3-
0.2-
0.1-
0.2
0.6
1.4
time (sec)
Figure 20: (a) True spikes and estimated rates of all 10 hawkmoth neurons over 1500 time points.
ADDITIONAL MODEL FITTING INFORMATION
All data had a 80-20 split for training and testing respectively. The Fourier frequency pruning was set to the minimum length scale of 10, 10. 3, and 16 for GP- VAE (simulated), MM-GPVAE (simulated), MM-GPVAE (fly), and MM-GPVAE (moth) respectively. GP length scales parameters were initialized to a value of 30 for all except for the hawkmoth evaluations (where initial values are indicated above), and jointly optimized with the ELBO. The covariance parameter value of le-2, 1e-2, 1e-3, and le-4 for GP- VAE (simulated), MM-GPVAE (simulated), MM-GPVAE (fly), and MM-GPVAE (moth) respectively. We additionally initialized the offsets d of the neural modality to the average log-rate of the neural data.
set at fixed
Latent dimensionality selection: In the of the simulated examples in this work, the dimension- ality was always set to the true generative dimensionality. For the real neural dataset examples in the main paper, the dimensionality chosen using distinct approach for each modality. For the
22

---

Published as a conference paper at ICLR 2024
VisuAl comuuuu
mottll rooonstructio
model recomstuuctin Ino shaned)
model inoondtrutioo momoool reconatrucion (8)
(10 anud)
shaaed Toadirga
neural loading
edimiti
1206
(D0c)
C.4
022
(a) Stimulus reconstructions using one shared latent vs no shared latent variables (left). Neural rate reconstructions using one shared latent vs no shared latent (right). (b) Fast moving neural latent (top) along with slow moving shared latent (bottom). The loadings matrix (right) display the tuning of each neuron to our latents, here we highlight neuron 5 which is positively modulated by the shared latent and neuron which is negatively modulated modulated by the shared latent.
Figure 21:
neural dataset, we use PCA and select the number of dimensions that capture greater than 97% of the variance. For the other modality (limb position or visual image) the dimensionality was chosen by increasing the number of dimensions in unimodal Fourier GP-VAE until cross-validated perfor- mance did not show significant improvement. Then, we identify the appropriate number of shared dimensions by fitting the full MM-GPVAE model to both modalities and systematically increasing the number of dimensions shared between datasets until predictive performance no longer improved. Across both datasets, we always saw that adding some shared dimensions increased overall predictive performance compared to none, suggesting shared structure across the data modalities for these experiments.
We also note that the dimensionality selection procedure was used to approximate shared and independent dimensionality best possible, but across CV folds the variance in cross validated log- ligllikelhoodd was high. However, small changes to the number of shared and independent dimensions did change latent representations or decoding performance much in our models.
For the evaluation where we compare MM-GPVAE to TNDM and PSID, we choose approximately the same dimensionality in the original TNDM paper (Hurwitz et al. (221) readers can easily compare fits across models.
CHOICE OF GAUSSIAN OBSERVATION VARIANCE INITIALIZATION. a²
We found that how we initially set the value of x² could effect the performance of the MM-GPVAE, especially regarding the reconstruction of the behavioral data. Although we learned the r² value jointly with the other parameters during optimization, the value of o² tended to vary minimally from its initial value. Therefore, our setting of o² tended to act reconstruction penalty, balancing the ability of the MM-GPVAE to prefer to reconstruct either the behavioral or neural data. Such a scaling term in the ELBO reconstruction has been seen in other models (Hurwitz et al., 2021), and in the original GP- VAE this parameter was chosen carefully through a specific cross- validation approach (Casale al. 2018). Here, simply set a² to value proportional to the dimensionality of the behavioral modality, which tended to balance the neural and behavioral terms in the ELBO, and generate good reconstructions for each modality. The initial values for x² were 1000, 100, 1e-6, and 1 for GP-VAE (simulated), MM-GPVAE (simulated), MM-GPVAE (fly), and MM-GPVAE (moth) respectively.
DERIVATION OF THE EVIDENCE LOWER BOUND (ELBO)
Here show the derivation of the evidence lower bound used for the MM-GPVAE. For clarity, will start by deriving the ELBO for the MM-GPVAE in the time domain, and ignore the Fourier representation in this derivation. The derivation trivially applies to the Fourier representation as the mapping between the two domains is linear. We set the GP prior parameters 0 = {a, f). W to be the loadings matrix and offsets from equation 4 in the manuscript, T to be the observed timepoints,
well,
23

---

Published as a conference paper at ICLR 2024
and Z to be the collection of all the shared and independent latents Z {za, Zs. 2bb.
logp (Y
Y
T, T. W.
o² 0)
p Y Y
Z. Φ. W. d2) p( T,O)
log
Irp 2 Y Bd
go ( YA. YB)
p Y
Z. Φ. W. o²) p(Z T,0)
log
9u(2 Y . dZ
( Y Y
EN [ogp(Y")
B Φ 0² 25 kbb + logg( Y
W, 2x, 2a) + log Z [T,0]
log stp z Y) Y)dZ
E. Ezzago
Σ log. (yB 90 (BB) ; I2IN) + Σ (og(P(yal((aA))
+ logg(( T, 0) + Hgs))
The Fourier domain representation of the ELBO only requires sampling over variational parameters in the Fourier space, but only changes the expression of the GP prior term p( Z)).
GP prior: The expectation of the GP prior term be expressed in the Fourier domain
Eg B2
pl Z w.w)
Σ E zmag
log.V (2e10) O0, (Kllw))
32
log( K Klpow + a)) + (Kplura + +a)
Zzzgg
Σ log( [ + a)) + piae (Y) (i (Y) p.ui
2
(K + a))
where the double sum is due to the variational distribution q being a mean field Gaussian, and here indexes latents and w indexes Fourier frequencies.
Neural likelihood: Since the estimated log-rates of the neural data is a linear transform of the shared and neural latent variables, we can also evaluate the expectation of the neural-modality likelihood term in closed form. Recall that
X W = W ŽB
where X is the matrix of embeddings, and Ž is the P x F matrix of Fourier-domain latent variables. We may therefore write the embedding for measurement at time f as
Titt = :.: Zb:
where w; is the 7th row of W and bJ is the 4th column of B If we let ip.. E Ez Fz and correspondingly let jú be the x F matrix of ip 's then Ez [1i,e] = w. jub.. We may equivalently write Zi,t as
Fs,t vec(I;,)) = vec(w, Zb.)
vec(b, 2 TTu,,
(w) @ b, )²
where & = vee( TT). We may therefore we may derive the variance of 24, as
Var Ergq [1i.]]
Var Zwgo [w.] b. )=]
Eg [[ww]
T
b. )=) E 2~ [w. T @ b, )=²]
E2zng [race[[w, @ b,) & b, )&8 1 - (w. jáb²)²
Trace[(w, 3 (:) 2 b. (V + vec(@)vec((i) (w. j2b))
-ww," @ 6. (w; @ br),
24

---

Published as a conference paper at ICLR 2024
where V is the diagonal posterior covariance of & whose elements the encoded Fourier variational variances, 0 (Y). Therefore, we observe that under the variational posterior Zi,e| Y ~ (mms,, v(a)) where m,, = w; μb, and V/, (w ! 2 b V (w, @ b.).
We note that for X t follows a long-normal distribution, meaning that, for Fi,t N (miii vi.))
then [[A,,] e m This allows us to specify the posterior expectation of the Poisson likelihood in closed form. Namely,
R Rz-as [og (og((vull(()x))) i ~as Miit log deit + Xo.t log yi.e!
yb, 2 a~gg [o6,]] E Ez-an [r,ee + const 2
Vssemmee + emi t01, + const
25


Reviews:
[{"review_id": "W1fYLDWrJz", "reviewer": "Reviewer_Kdv4", "summary": "The authors proposed a multi-modal Fourier-domain Gaussian Process variational autoencoder (MM-GPVAE), aiming at jointly modeling observed neural and behavioral data in neuroscience applications with both shared and independent latent subspaces. MM-GPVAE combines the ideas in previous GP-VAE, GP Factor Analysis (GPFA), and spectral representation. The authors have also presented simulations and two neuroscience case studies showing the effectiveness of MM-GPVAE modeling neural measurement / spiking data with either visual stimulus or position data.", "questions": "1. How sensitive is the performance with respect to the setting of several hyperparameters, including $F$ and the dimension of shared and independent latent subspace, etc.? \n\n2. How optimization was done with FFT/iFFT involved in MM-GPVAE?", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. This submission exploits the multi-modal latent variable models based on GP-VAE, using Poisson and normal likelihoods for spiking and image data respectively, which can help modeling both shared and independent factors across modalities in neuroscience. \n\n2. Experiments on simulated and real-world datasets with visualized results showing the potential in jointly analyzing neuroscience data. \n\n3. The presentation is clear with articulated motivations of the presented mm-GPVAE. Source code is also provided.", "weaknesses": "1. Methodological contributions may be limited. Many papers have done similar things, such as extending the latent space into frequency domain [1,3] and using the Poisson model for spiking data [1,2]. The method proposed in this paper is combination of existing methods, not fundamentally different. \n\n2. The experimental results can be more comprehensive: (1) More baselines, especially similar methods such as [1,4] need to be compared; (2) The authors stated one main advantage of the Fourier domain is that high frequencies can be pruned and thus we can sparsify the variational parameters. However, no results are provided about this point; (3) An ablation study between the models with Frequency prior and with other priors (e.g., a temporal form) should be added; (4) More discussions on the claimed interpretability should be provided. \n\n[1] Keeley S, Aoi M, Yu Y, et al. Identifying signal and noise structure in neural population activity with Gaussian process factor models. Advances in neural information processing systems, 2020, 33: 13795-13805.\n\n[2] Duncker L, Sahani M. Temporal alignment and latent Gaussian process factor inference in population spike trains. Advances in neural information processing systems, 2018, 31.\n\n[3] Hensman J, Durrande N, Solin A. Variational Fourier Features for Gaussian Processes. J. Mach. Learn. Res., 2017, 18(1): 5537-5588.\n\n[4] Pearce M. The Gaussian Process prior VAE for interpretable latent dynamics from pixels. Symposium on advances in approximate Bayesian inference. PMLR, 2020: 1-12.", "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": "", "extra_scores": {}}, {"review_id": "fQpBHs5b2y", "reviewer": "Reviewer_YQRh", "summary": "This paper proposes a multi-modal gaussian process VAE model for neuroscientific data. Due to the time-varying nature of neuronal data, the authors use gaussian priors for Fourier frequencies of the data. Thus the network features FFT and IFFT layers. Results are shown on synthetic data as well as on calcium signaling data from Drosophila and Moth EMG data.", "questions": "The synthetic data does not make sense to me, what exactly is the connection between the digits and the \"spikes\"?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The use of Fourier frequencies as latent variables can have some advantages in learning neural features that correlate with motor or other behaviors. The authors also separate shared from independent variability in the neural data which can be useful in understanding what part of behavior could actually be predicted from specific neuronal measurements.", "weaknesses": "-The synthetic data seems weak and unnecessary. Why not use a simulator like NEST or augment one and create \"behavioral\" data with known connections to this. MNIST images have a different structure because they are a two dimensional image. \n-The authors lack comparison to anything but an ablation of their own model, and that too on the above artificial dataset. \n-At least the authors should compare to RNN/transformers/ODE or other sequential models for this data\n-The authors seem to criticize DNNs for transforming the data too much such that the latent variables are \"not obviously related\" to external variables, and yet they too use DNNs. \n-There is no theory or study of identifiability of the dynamics or conditions/noise under which dynamics are recovered.", "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": "", "extra_scores": {}}, {"review_id": "VnRm5wl0Hl", "reviewer": "Reviewer_bxWt", "summary": "This paper introduces a method, MM-GP-VAE, for modeling multimodal data, presented in the context of analyzing neural and behavioral data.  The balance that the proposed method strikes is retaining the interpretability and disentanglement enjoyed by simpler models, but retaining the expressivity and ease-of-use of deep methods.  Low-dimensional latent states describe each channel individually, and a shared low-dimensional state describes the shared covariance.  The recurrent model is implemented as a Fourier-domain GP-VAE, to encourage smoothness in the latent state and increase expressivity.  The model is benchmarked on a synthetic MNIST-like experiment, and is then applied to neural data.  The model appears to perform well against some sensible baselines.", "questions": "I don’t have any further questions from those outlined in Weaknesses.\n\n**Note:** If the authors can allay my concerns, then I fully intend to increase my review score.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "The paper is incredibly well presented.  Figures are prepared exceptionally well.  The prose is clear, and presents a self-contained introduction to all of the necessary techniques and considerations. \n\nThe method itself is gracefully simple.  Although there is not a huge methodological contribution, the correct components parts are assembled and deployed in a way that is very experimentally useful.  I could see this methodology having real uptake within the community, and engendering multiple follow-up works.", "weaknesses": "I do have several queries/concerns however:\n\n- **a. Fixed time horizon**:  The use of an MLP to convert the per-timestep embeddings into per-sequence Fourier coefficients means that you can only consider fixed-length sequences.  This seems to me to be a real limitation, since often neural/behavioral data – especially naturalistic behavior – is not of a fixed length.  This could be remedied by using an RNN or neural process in place of the MLP, so this is not catastrophic as far as I can tell.  However, I at least expect to see this noted as a limitation of the method, and, preferably, substitute in an RNN or neural process for the MLP in one of the examples, just to concretely demonstrate that this is not a fundamental limitation. \n\n- **b. Hidden hyperparameters and scaling issues**:  Is there a problem if the losses/likelihoods from the channels are “unbalanced”?  E.g. if the behavioral data is 1080p video footage, and you have say 5 EEG channels, then a model with limited capacity may just ignore the EEG data.  This is not mentioned anywhere.  I think this can be hacked by including a $\\lambda$ multiplier on the first term of (6) or raising one of the loss terms to some power (under some sensible regularization), trading off the losses incurred by each channel and making sure the model pays attention to all the data.  I am not 100% sure about this though.  Please can the authors comment. \n\n- **c.  Missing experiments**:  There are a couple of experiments/baselines that I think should be added. \n  - Firstly, in Figure 3, I'd like to see a model that uses the data independently to estimate the latent states and reconstruction.  It seems unfair to compare multimodal methods to methods that use just one channel.  I’m not 100% sure what this would look like, but an acceptable baseline would be averaging the predictions of image-only and neuron-only models (co-trained with this loss).  At least then all models have access to the same data, and it is your novel structure that is increasing the performance.\n  - Secondly, I would like to see an experiment sweeping over the number of observed neurons in the MNIST experiment.  If you have just one neuron, then performance of MM-GP-VAE should be basically equivalent to GP-VAE.  If you have 1,000,000 neurons, then you should have near-perfect latent imputations (for a sufficiently large model), which can be attributed solely to the neural module.  This should be a relatively easy experiment to add and is a good sanity check.\n  - Finally, and similarly to above, i’d like to see an experiment where the image is occluded (half of the image is randomly blacked out).  This (a) simulates the irregularity that is often present in neural/behavioral data (e.g. keypoint detection failed for some mice in some frames), and (b) would allow us to inspect the long-range “inference” capacity of the model, as opposed to a nearly-supervised reconstruction task.  \n\n  Again, these should be reasonably easy experiments to run.  I’d expect to see all of these experiments included in a final version (unless the authors can convince me otherwise).\n\n- **d.  Slightly lacking analysis**:  This is not a deal-breaker for me, but the analysis of the inferred latents is somewhat lacking.  I’d like to see some more incisive analysis of what the individual and shared features pull out of the data – are there shared latent states that indicate “speed”, or is this confined to the individual behavioral latent?  Could we decode a stimulus type from the continuous latent states?  How does decoding accuracy from each of the three different $z$ terms differ? etc.  I think this sort of analysis is the point of training and deploying models like this, and so I was disappointed to not see any attempt at such an analysis.  This would just help drive home the benefits of the method.  \n\n\n### Minor weaknesses / typographical errors:\n\n1.  Page 3:  why are $\\mu_{\\psi}$ and $\\sigma_{\\psi}^2$ indexed by $\\psi$?  These are variational posteriors and are a function of the data; whereas $\\psi$ are static model parameters. \n\n2.  Use \\citet{} for textual citations (e.g. “GP-VAE, see (Casale et al., 2018).” ->   “GP-VAE, see Casale et al. (2018).”)\n\n3.  The discussion of existing work is incredibly limited (basically two citations).  There is a plethora of work out there tackling computational ethology/neural data analysis/interpretable methods.  This notable weakens the paper in my opinion, because it paints a bit of an incomplete picture of the field, and actually obfuscates why this method is so appealing!  I expect to see a much more thorough literature review in any final version.\n\n4.  Text in Figure 5 is illegible.  \n\n5.  Only proper nouns should be capitalized (c.f. Pg 2 “Gaussian Process” -> “Gaussian process”), and all proper nouns should be capitalized (c.f. Pg 7 “figure 4(c)”).\n\n6.  Figure 1(a):  Is there are sampling step to obtain $\\tilde{\\mu}$ and $\\tilde{\\sigma}^2$?  This sample step should be added, because right now it looks like a deterministic map.  \n\n7.  I think “truncate” is more standard than “prune” for omitting higher-frequency Fourier terms.\n\n8.  I find the use of “A” and “B” very confusing – the fact that A is Behaviour, and B is Neural?  I’m not sure what better terms are.  I would suggest B for Behavioural – and then maybe A for neural?  Or A for (what is currently referred to as) behavioral, but be consistent (sometimes you call it “other”) and refer to it as Auxiliary or Alternative data, and then B is “Brain” data or something.  \n\n9.  The weakest section in terms of writing is Section 3.  The prose in there could do with some tightening.  (It’s not terrible, but it’s not as polished as the rest of the text).\n\n10.  Use backticks for quotes (e.g. ‘behavioral modality’ -> ``behavioral modality’’).", "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": "", "extra_scores": {}}, {"review_id": "kz2CHhKk0f", "reviewer": "Reviewer_y53H", "summary": "The paper proposes a variational auto encoder framework to jointly model neural activities and behavior data. Their contribution is to add Gaussian Process Prior and parameterize the latent in Fourier domain.  Their goal is to distinguish the separate or joint latents from multiple modalities, and assume a linear combination of those latents. The work is evaluated on one simulated dataset (MNIST + spike train) and two real-world datasets from different animal species (drosophila, Manduca).", "questions": "1. Is nonlinearity needed in eqn (4)?\n2. How to decide dimensions for separate and joint latent?", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "1. The work is well-motivated,  jointly train latent variable models from multi-modal neural data is an important question.\n2. The contribution to use Gaussian Process prior to capture correlations in latent space, and parameterize the latent in Fourier space are novel contributions.\n3. The work is demonstrated on multiple datasets to improve its soundness, including one synthetic dataset and two real-world datasets from different animal species .", "weaknesses": "Methods:\n1. The assumption of linearity between latents in eqn (4) seems over-simplified, other nonlinearity variants could be further investigated. \n2. The construction of simulated dataset is confusing, it is not clear to see how the rotation angle of MNIST image is inherently related to spiking distribution. And it is puzzling that how it is tightly linked with real-world datasets. \n3. The paper proposed parameterize the latents in Fourier space to constrain them more distinguishable, and attempt to prune the high-frequency components. While it didn't provide solid ablation studies and compared to the original temporal space on real datasets, and it is problematic to claim no information preserved in high-frequency domains. \n4. The authors didn't provide clear explanation or guidance about important hyperparameters of latent dimensions. \n\nEvaluations:\n1. No quantitative measurements on real datasets only with a few qualitative samples, and not compared with other baselines.\n2. In Fig 4(d), different behavior states are not clearly separable in these subspaces.\n3. Blurriness and mismatch of reconstruction outputs and latents in the results.\n\nWriting:\n1. The paper is not well organized, missing quantitative evaluations and comparisons with baselines, a few important contents are pointed to supplementary, while it is not clearly described there either.", "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": "", "extra_scores": {}}, {"review_id": "aLAXaPr5m8", "reviewer": "Reviewer_9t6f", "summary": "This paper proposes a model for multi-modal data based on a latent space using both shared and mode-specific latent subspaces and Gaussian Process priors. These GP priors are intended to capture smooth variations in the latent space, akin to Gaussian Process Factor Analysis (GPFA) for neural time series data. Another key ingredient is the use of a Fourier basis in the latent space, which allows for a natural low-frequency regularization on latent dynamics. Experiments include both simple synthetic data setups and applications to two real data sets in which neurophysiological recordings are combined with behavior.", "questions": "- As noted above, the GP method would seem to need data to be broken into snippets of reasonable size. For distinct trials, this is pretty obvious, but what if the time series are very long (e.g., natural behavior)? Does it work to simply snip the behavior at random? Does the model handle this? Do all snippets need to be the same size?\n- Is it possible to put some shrinkage priors on the $W$ matrices in (4) so that unneeded dimensions are pruned away? Can the model distinguish between \"true\" and \"false\" shared latents? Is there some notion of parsimony here?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- Use of GP priors to smooth latent trajectories is a well-matched strategy for smoothing highly noisy neural data.\n- Moving to the Fourier basis is a clever idea that allows for a very natural regularization. \n- The idea of shared and specific latent subspaces has been done before but adds to the interpretability of the model.\n- Several tricks are combined here that render the notoriously inefficient GP inference to be performed efficiently.\n- Good set of experiments on real neural data.", "weaknesses": "- The synthetic experiments tend to focus on smooth, continuous variation in \"incidental\" properties (e.g., size and angle of a digit) that may not be a good approximation in many real data sets.\n- The GP model may not work well in cases where data do not have a natural trial structure or produce very long behavioral bouts.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)