Across-animal odor decoding by probabilistic manifold alignment
Year: 2021

Abstract:
Identifying the common structure of neural dynamics across subjects is key for extracting unifying principles of brain computation and for many brain machine interface applications. Here, we propose a novel probabilistic approach for aligning stimulus-evoked responses from multiple animals in a common low dimensional manifold and use hierarchical inference to identify which stimulus drives neural activity in any given trial. Our probabilistic decoder is robust to a range of features of the neural responses and significantly outperforms existing neural alignment procedures. When applied to recordings from the mouse olfactory bulb, our approach reveals low-dimensional population dynamics that are odor specific and have consistent structure across animals. Thus, our decoder can be used for increasing the robustness and scalability of neural-based chemical detection.

heading: 1 Introduction
text: Collective network dynamics are the foundation of neural computation, from early sensory encoding [1], to working memory [2], decision making [3], or motor control [4]. Neural population activity often has low-dimensional structure that is qualitatively preserved across sessions [3], behavioral states [5], and even animals [6–8]. Aligning neural datasets from multiple recordings and animal subjects into a common latent space could provide a powerful tool for extracting unifying principles of brain computations. Nonetheless, progress is hampered by the lack of robust statistical tools for extracting shared neural population dynamics across datasets.
The alignment of low-dimensional neural manifolds is equally important in brain-computer interface (BCI) applications, where the assumption of a common low dimensional structure in neural responses is used to compensate for instability in neural recordings and ensure robust performance over time [9, 10]. This is not only restricted to motor control; in the sensory domain, a new generation of chemical detectors exploit the unsurpassed sensitivity of rodent olfactory receptors, but bypass the limitations of training animals to report odor identity by directly decoding it from neural responses [11–13]. Unfortunately, the applicability of this idea is limited by the need to learn the mapping between neural activity and chemical identity on an animal-by-animal basis, which involves costly data collection [14]. Statistical tools that make use of data from previous animals to quickly calibrate the decoder in a new animal could dramatically increase the practical use of such technologies.
The basic structure of the rodent olfactory bulb (OB) is well-understood [15]: a vast number of different olfactory receptors project to distinct OB subpopulations (glomeruli) whose spatial organization is largely preserved across animals [16, 17]. Odor identity is encoded in low-dimensional transients of their population dynamics [1, 18, 13], though how the glomerular topography translates into stimulus-dependent network dynamics is less clear. Nonetheless, given the anatomical and functional similarity of the rodent olfactory bulb, it is likely that the OB response to individual odors is itself
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
somewhat stereotyped across animals. With appropriate alignment, one could construct ‘universal’ odor decoders, applicable across animals.
Here we propose a novel probabilistic approach to aligning observations across neural recordings taken from different mice so as to extract commonalities in their circuit-level representations of different odors. Our estimation procedure, which uses our ‘aligned mixture of latent dynamical systems’ (amLDS) graphical model, learns independent mappings of different recordings into a common latent manifold.1 In this space, neural trajectories are similar across animals but distinct across odors which allows for robust odor identification. Moreover, our probabilistic decoder provides uncertainty estimates which can be used for optimal decision making. We used extended numerical simulations to evaluate the properties of the decoder in response to a range of biologically relevant manipulations, such as variations in measurement noise, dimensionality and across-trial variability in the latent dynamics, (mis)alignment of electrodes across animals, etc. We then tested our procedure on recordings from the OB of mice in response to several monomolecular odors. When fit to this data, our model confirmed that encoding of odors relies on transient dynamics in a shared low-dimensional manifold and that these dynamics are largely preserved across animals. Taking advantage of this reproducible structure, we were able to obtain remarkable improvements in odor identification performance relative to alternative approaches.

heading: 2 Shared odor encoding dynamics across animals
text: We start by mathematically formalizing the idea of shared stimulus-specific neural dynamics with animal-specific measurements (Fig. 1A). The amLDS graphical model assumes d-dimensional latent variables zt, which evolve according to stochastic linear dynamics with odor-specific parameters (Fig. 1B):
zt = Akzt−1 + bk,t + wt, (1)
where t = {2 . . . T} indexes time within trial; Ak and bk denote parameters specifying the dynamics in response to odor k, with independent noise, wt ∼ N (0,Qk). The prior for the initial conditions z1 is normal, with zero mean and covariance Q0.
The shared odor-specific dynamics are mapped into measured neural responses via an animal-specific linear Gaussian observation model, shared across all odor conditions
x (m) t = Cmzt + vt, (2)
with i.i.d. noise vt ∼ N (0,Rm); Cm and Rm are observation model parameters specific to animal m. We will use Nm to denote the size of observations for animal m. Note that the animal-specific observation model means that it is not strictly required for the neural measurements to be tightly matched. For instance, different animals may have different numbers of measurement channels (e.g. due to some electrodes becoming damaged). Moreover, the model accounts for variability in the SNR of different electrodes (across channels and devices). Assuming that the measurement process
1Code available: github.com/pedroherrerovidal/amLDS.
is consistent in terms of the properties of the electrodes and their alignment on the olfactory bulb surface, one can additionally define a prior that encourages Cs to be similar across animals. To keep the procedure general, we only consider the unconstrained version here.
Parameter learning
The dataset D combines measurements from M animals, jointly covering all K stimulus conditions. For each animal-stimulus pair, (m, k), we have Im,k trials.2 The goal of the parameter estimation procedure is to derive maximum likelihood estimates of all model parameters, jointly denoted by θ = {A1:K ,b1:K ,Q0:K ,C1:M ,R1:M}. We use expectation maximization (EM) to optimize the parameters, by introducing the variational distribution q(z) and using it to lower-bound the marginal likelihood [19]:3
log P(D|θ) = ∑ i,m,k ∫ z log P(x|θ)q(z)dz
= ∑ i,m,k ∫ z q(z) (log P(z,x|θ)− log P(z|x, θ)) dz = ∑ i,m,k ∫ z q(z) log P(x, z|θ) q(z) dz− ∑ i,m,k ∫ z q(z) log P(z|x, θ) q(z) dz
 = L(q, θ) +
∑ i,m,k KL(q(z)||P(z|x, θ)).
Due to the Markovian structure of the dynamics, the joint probability of the latent and observations has a relatively simple joint Gaussian structure. The E-step optimizes L w.r.t. q, by setting q(z) = P(z|x, θ); this is computed by traditional Kalman smoothing using the appropriate trial-specific parameters (Suppl. S1).
The M-step optimizes L w.r.t. to θ. The cost function decomposes into the sum of terms that only depend on the subject-specific observation parameters (Cm,Rm) and the condition-specific latent dynamic parameters (Ak,Qk,Q0,bk,t), respectively, and can be optimized separately. Taking the derivative of the loss and setting it to zero results in parameter updates:
bnewk,1 = 1
Ik ∑ i E[z1] (3)
Qnew0 = 1
Ik ∑ i ( E[z1z1>] − E[z1]E[z1]> ) , (4)
where the indicator i covers all trials of the same stimulus type, across all animals Ik = ∑
m Im,k; expectations are taken under posterior for trial i (explicit indexing omitted for brevity).
The same set of trials are used for updating the latent dynamics, and the input drive:
Anewk = ( Ik∑ i=1 T∑ t=2 ( E[ztzt−1>]− bk,tE[zt−1>] ))( Ik∑ k=1 T∑ t=2 E[zt−1zt−1>] )−1 (5)
bnewk,t = 1
Ik Ik∑ i=1 (E[zt]−Anewk E[zt−1]) . (6)
Finally, the updated noise covariance is computed as
Qnewk = 1
Ik(T − 1) Ik∑ i=1 T∑ t=2 ( E[ztzt>]−Anewk E[zt−1zt>]− E[ztzt−1>]Anewk >
+ Anewk E[zt−1zt−1>]Anewk > − E[zt]bnewk,t > − bnewk,t E[zt]> + Anewk E[zt−1]bnewk,t > +bnewk,t E[zt−1]>Anewk > + bnewk,t bnewk,t > ) . (7)
2Note that it is not necessary that all stimuli are presented to all animals, nor that the number of stimulus presentations is identical across animals and conditions.
3To keep notation simple, we drop the explicit indexing by time, trial, stimulus type, and animal of x and z.
The updates for the observation model are the standard Kalman ones, with the distinction that the index i iterates across all trials and stimulus conditions for a particular animal (total number Im = ∑ k Im,k):
Cnewm = ( Im∑ i=1 T∑ t=1 xtE[zt] )( Im∑ i=1 T∑ t=1 E[ztzt>] )−1 (8)
Rnewm = 1
ImT ∑ i,t ( xtxt > −Cnewm E[zt]xt> − xtE[zt]>Cnewm > + Cnewm E[ztzt>]Cnewm > ) (9)
All the required posterior moments are obtained during the E step via Kalman smoothing as E[zt] = µt|T , E[ztzt>] = Σt|T − µt|Tµt|T> (see Suppl. S1 for further details).
As an important practical side note, the quality of the estimated parameter depends critically on a good initialization of θ. We achieve this by fitting a simpler factor analysis model (FA), which ignores all temporal dependencies. This directly yields initial conditions for Cm and Rm. The inputs bk,t are initialized to the average of the corresponding inferred latent variable zt, using all trials from stimulus condition k, while parameters Ak and Q0:K are initialized heuristically by linear regression using the FA-inferred posterior mean estimates of the latents.
We used Bayesian model comparison to determine the dimensionality of the latent space from data (evaluated on a separate validation set). As a more intuitive measure of goodness-of-fit, we also estimated the reconstruction error when predicting unobserved measurement dimensions (also known as leave-neuron-out error), commonly used for assessing the quality of models of neural population dynamics [20, 21]. More specifically, we use trials from the validation dataset, infer the posterior over the latent trajectory given data measurements in all electrodes except the j-th by Kalman smoothing P (z1:T |x¬j , θ), which is used to predict the values of the response of the jth neuron, x̂jt = Cjmµ (¬j) t|T , where Cjm is the jth row of Cm. We report reconstruction error as the squared distance between the reconstructed and measured activity, 〈‖xjt − x̂ j t‖22〉j,i, where ‖ · ‖2 is the L2 norm and 〈·〉j,i is the mean over trials and choices for j.
Odor decoding as hierarchical inference
At test time, one needs to infer the identity of the stimulus given the observations x1:T in animal M , and the estimated parameters θ, obtained by Bayes rule, P(k|x1:T , θ) ∝ P(x1:T |k, θ)P(k). This requires estimating the marginal likelihood for each possible stimulus condition, which can be computed iteratively during the Kalman filtering procedure. More specifically, using the chain rule and then the Markov structure of the stimulus-specific LDS, we have
P(x1:T |k, θ) = T∏
t=1
P(xt|x1:t−1, k, θ) = T∏
t=1
∫ P(xt, zt|x1:t−1, k, θ)dzt.
Taking the logarithm, rearranging the terms and taking advantage of the fact that the posterior marginals are normal with parameters obtained by Kalman filtering (see Suppl. Info. S1), P(zt|x1:t−1, k, θ) = N (µt|t−1,Σt|t−1), yields:
log P(x1:T |k, θ) = T∑
t=1
log ∫ P(xt|zt,x1:t−1, θk)P(zt|x1:t−1, k, θ)dzt
= T∑ t=1 logN ( xt; Cmµt|t−1,CmΣt|t−1Cm > + Rm )
Lastly, here we assume a uniform prior over possible stimuli, to match the statistics of our OB dataset. Nonetheless, for real world applications the prior should reflect natural stimulus statistics.

heading: 3 Numerical experiments on simulated data
text: To validate the parameter estimation and decoding procedure, we constructed artificial datasets with the same statistical regularities assumed in the graphical model. More specifically, we defined a
shared low dimensional manifold (d = 3) and embedded K = 50 latent trajectories evolving over T = 41 time steps (Fig.2A). The stimulus-dependent inputs bk were constructed using a common template with stimulus-specific amplitude (individual dimensions scaled by a value drawn from N (1; 0.0004)) and rotation (evenly spaced over 170 degrees). The other latent dynamics parameters were set randomly: matrices Ak have diagonal entries drawn from N (0.4; 0.01) and off-diagonal drawn from N (0; 0.04) and the noise covariances Qk and Q0 are diagonal with variances drawn from N (0.55; 0.0025). To model across-animal variability in the observation model we started from a prototype projection matrix C∗, which we then morph into animal specific parameters Cm by corrupting the individual latent axes with independent additive noise N (0;α2), and rescaling their norm to a noisy value around 1.0 (gaussian noise with variance 0.03). Parameter α controls the degree of alignment of measurements across animals, with default value 0.1. Finally, the observation noise covariance is set to be diagonal with variances set to the absolute value of independent draws from N (0; 0.25) (Fig.2B). This setup results in latent dynamics qualitatively similar to odor responses, and has the advantage of explicitly parametrizing the similarity across animals and stimuli.
As a first test of the estimator we simulated 5 animals, with Im,k = 50 trials for each stimulus condition, which is on the same order of magnitude to the amount of data one may be able to collect in experiments. amLDS has low computational and memory demands; on a 2.9GHz CPU it takes 58 minutes for parameter learning (with 12500 trials; average over 5 runs) and 4 seconds to infer
stimulus class in a given trial (averaged over 2500 trials). In this setup, model comparison based on reconstruction error was able to correctly discover the dimensionality of the latent dynamics (Fig.2D) and reconstruct the underlying trajectories (Fig.2C), up to the expected degeneracies in the latent space scale and alignment [19]. Moreover, the decoding of stimulus identity of unseen test trials (20 per condition) revealed that hierarchical inference is able to correctly identify the true stimulus out of the 50 possible classes (Fig.2E), where we measure decoding accuracy as the fraction of test trials for which the MAP estimate of the stimulus class recovers the ground truth stimulus identity.
Qualitatively similar results can be obtained for a wide range of model choices. Performance is consistent across a wide range of noise levels in the observations (Fig.2F, top), which account for differences in recording quality. The same robustness can be measured for noise in the latent dynamics, due to potential changes in encoding owed to masking odors, variability in sniffing, etc (Q ∼ N (µ, 0.0025), µ ∈ [0.1, 0.9]; Fig.2F, bottom). The estimation is also relatively robust to different degrees of variability of the projections Cm across animals (Fig.2G, left), which may arise for example due to variability in the alignment of the electrodes on the OB surface. We also varied the dimensionality of the latent dynamics (Fig.2G, middle) and the number of stimulus classes (Fig.2G, right) and found that our decoding procedure shows robust performance for a wide range of settings. Given the relatively modest data requirements (50 trials per subject, and condition), these observations suggest that our approach is robust enough to be applicable to experimentally relevant data regimes.
To better quantify decoding performance as a function of the amount of data that needs to be collected in a new animal, we fixed the amount of data collected from 4 ‘source’ animals (K = 50, Im,k = 50) and varied the amount of new data collected in a fifth ‘target’ animal, for which we evaluated decoding accuracy. We found that in this setup amLDS yields high performance with a minimal amount of additional data (Fig.2F, red; one measurement per stimulus condition). Reaching the same decoding performance with data collected in the target animal alone would require many more trials, because the dynamics of all K stimuli need to be learned de novo (Fig.2F, blue). In the case of amLDS however, (effectively) one only needs to align the new measurements to the manifolds estimated based on data from previous animals.4 In fact, observing a few trials for a single stimulus is enough to achieve reasonable performance in the target animal (Fig.2F, green).
We compared the decoding performance of our method to three state-of-the-art methods: 1) across condition alignment by canonical correlation analysis (CCA) [22, 23] 2) multi-set CCA (mCCA) [24] and 3) Procrustes alignment with FA-based dimensionality reduction (FA+Procrustes) [9, 6]. In general, for all these approaches the latent dimensions are determined from data that combines all stimulus conditions, without explicit knowledge of the stimulus identity. Decoders are trained based on the resulting representation using stimulus labels.
CCA takes pairs of measurements from two animals and projects them into a shared low dimensional manifold so as to maximize the correlation between them. Since there is no meaningful way of pairing single trials, CCA requires averaging responses across trials within condition and pairing these average responses by stimulus type. Moreover, unlike our approach, both procedures treat measurements as independent over time, ignoring temporal structure in the neural responses. Since simple CCA can only be used for pairs of animals, we increased the number of trials in the single source animal to match the other methods. Additionally, we also used multi-set CCA, a generalization of CCA for more than two animals.
The alternative FA+Procrustes alignment is similar to CCA in that it requires across-trial averages and only considers similarity between pairs of measures; however, it can be applied multiple times to align animals 2−M to the prototype axes given by animal 1 (the target animal data was used for FA and Procrustes parameter estimates in this case). More formally, we use orthogonal Procrustes [25, 9] where given a set of prototype zk∗1:T and misaligned Oz k 1:T trajectory pairs finds a rotation between
the two, by minimizing the average misalignment after applying a rotation O, argminO ∑
k ‖zk∗1:T − Ozk1:T ‖2F , under the constraint of orthogonality, OO> = I; where ‖ · ‖F is the Frobenius norm. The stimulus-specific latent trajectories zk1:T are obtained via a dimensionality reduction step, performed independently in each animal by factor analysis. As there is no standard way for establishing the shared latent dimensionality based on data for (m)CCA and FA+Procrustes, we set it to the ground truth d in all cases. Having established a linear map between observations and a common latent
4While the estimation procedure still uses EM based on the full data, the fractional contribution of the new measurements to the latent parameter updates is comparatively small, and becomes negligible as the number of sources increases.
space by either procedure, individual trials across conditions were all projected in this space. We used a linear SVM with L2 regularization to decode stimulus identity from the latent trajectories [13, 26]. All controls were constructed using the library scikit-learn [27]. We found that although all three approaches performed above chance, none came close to matching the across animal decoding performance of amLDS.
Robustness to fluctuations in concentration, odor mixtures
Odor decoding in the wild needs to overcome further challenges that are not explicitly captured in the amLDS graphical model, such as fluctuations in concentration that rescale the amplitude and speed of odor specific responses [28]. Moreover, odorants are rarely present in isolation and have to be detected when presented in mixtures [29]. We used further numerical simulations to test how the estimator would respond in such scenarios.
We model across-trial variability in odor concentration as changes in the amplitude of the input drive bt,k [30], with all other parameters the same as before. The parameter learning used data from the highest concentration and the decoder was applied to test trials with lower concentrations (Fig.3A). We found that the decoding accuracy degrades gracefully as the amplitude of responses decreases, remaining above chance for a wide range of concentrations.
It remains unclear how mixtures are encoded in OB responses [31], but a reasonable first approximation is to use a linear interpolation of sample trajectories for the individual odors present [32, 33]. For simplicity, we used K = 2 and further tied the parameters A and Q of the latent dynamics. We constructed two component mixtures which interpolate between stimuli labeled A and B (Fig.3B, left). Using a model trained on individual odors, we estimated decoding performance on mixture test trials while varying the ratio between the two components (Fig.3B, middle). The posterior over the two odors appropriately reflects the inherent ambiguity in the observations (Fig.3B, right). Overall, these results suggest that the decoding procedure is robust to naturalistic odor variations.

heading: 4 Across-animal decoding in the rodent olfactory bulb
text: We tested our model on neural recordings from a 64-site grid-electrode stereotaxically implanted over the dorsal part of the olfactory bulb in five mice [13]. We simultaneously recorded a pressure readout of the animal’s sniff cycle. The animals were presented with five different monomolecular odorants (methyl valerate, ethyl tiglate, benzaldehyde, hexanal and background air), delivered in a temporally precise window (20-100 trials each). We extracted glomerular responses in a time window of 210ms from inhalation onset. The neural responses were preprocessed by removing the high frequency component (100Hz low-pass filter) and electrical noise (60Hz notch filter). Single electrode signals with peaks exceeding 2mV in a 5 second period were considered damaged and removed from the study. The electrode-specific signals were extracted by subtracting the instantaneous mean across all viable electrode sites.
We split the animals into 4 sources, and one target. All data from the source animals and 10% of the trials from the target animal were used for parameter learning and hyperparameter estimation. The
remaining data from the target animal was used to assess odor decoding performance. When fitting the model to the neural recordings we found that the inferred latent trajectories were separable across stimuli and similar across animals (Fig. 4A). Some across-animal variability in the responses to the same odor would be expected in the model for a limited number of trials – due to the independent instantiation of the noise in the latent space, but it may also reflect unaccounted for sources of across-animal variability. The model comparison revealed a relatively small latent dimensionality in the data (Fig. 4B; we used d = 7 for all subsequent analyses). The model accurately captured variability in observed electrodes in the target animals, despite the very limited target data (2-6 trials per stimulus; Fig. 4C). The inferred map from the neural measurements into the latent space was similar across animals (Fig. 4D; see also Suppl. Fig. S1), likely reflecting the stereotypy of OB encoding and a precise alignment of the electrodes on the bulb surface.
We evaluated the decoding performance of amLDS on a target animal which revealed high accuracy MAP inference based decoding (Fig. 4E), with generally sharp single trial posteriors (Fig. 4F). We measured decoding accuracy for all possible choices of target animal and compared it to alternative estimators (Fig. 4G). The CCA, mCCA and FA+Procrustes estimation procedures were similar to those described for the simulated data, with the same latent dimensionality inferred for amLDS. We used across trial averages within each stimulus condition from the source animals and the 10% fraction of the target animal data for dimensionality reduction and across-animal alignment. We used the same data to train a linear SVM classifier based on the latent trajectories, which we then used to predict odor identity for the withheld target data. For CCA [22] we used single source animals and the training target data for the alignment and SVM training, averaging decoding performance on test data across all possible sources. Our probabilistic decoder significantly outperformed all three alternative alignment methods. For all possible choices of target, amLDS also systematically outperformed probabilistic decoding based on data from the target animal alone, confirming that pooling together data across animals improves data efficiency and overall decoding performance. While the single animal performance is close to that of amLDS given enough data, 82%± 8% vs. 76%± 3% when training with 50% of the target animal’s dataset (Fig. 4G), the difference between the two becomes dramatic in the low data regime used for amLDS (10% of target data used for training), with single animal performance going down to 60%± 8%. This reinforces the idea that across-animal decoding is key to limiting the amount of data that needs to be collected for calibrating any new animal at deployment time.

heading: 5 Discussion
text: Animals are unparalleled in their ability to sense odors, which makes them an invaluable resource for a wide range of applications from security to medicine [34–36]. BCI approaches to odor detection attempt to bypass costly behavioral training by reading out odor information directly from the olfactory system, but they remain data limited [14]. Here we proposed a novel probabilistic procedure for mapping neural responses to odors into a shared manifold. This representation allows for robust odor decoding in a new animal with minimal amounts of further data collection. Calibrating the decoder to a new animal can be achieved with a small subset of stimuli, limiting the demand for rare/expensive chemicals. For instance, in our numerical experiments the decoder could extract odor identity in a new animal with as few as 2 trials for each stimulus, or 50 trials of a single stimulus. Moreover, we also found the estimation procedure to be robust to modeled fluctuations in odor concentration and the presence of background odors. These benefits were confirmed in data from mouse OB, where we found that decoding quality improved dramatically compared to previously-proposed alignment procedures, particularly when using a limited amount of data.
Our across-animal probabilistic decoding procedure brings neural chemical detection closer to practical applicability. Not only does it improve scalability and robustness, but it also significantly reduces data collection requirements. Furthermore, most odor detection applications have asymmetric costs and need to negotiate complex trade-offs between sensitivity and accuracy. For instance, when detecting Clostridium difficile [37], the medical costs of missing a positive diagnosis needs to be balanced against the financial costs of unnecessarily quarantining someone due to a false alarm. Our probabilistic decoder provides an explicit report of uncertainty associated with possible outcomes, which is critical for optimal decision making in these contexts.
Alternative procedures for neural datasets alignment are usually not explicitly probabilistic.5 Instead, they rely on a sequence of independent steps: first a (usually off-the-shelf) dimensionality reduction procedure. Dimensionality reduction is followed by across-condition data alignment within the latent manifold, using approaches such as CCA (and its multivariate generalizations [39, 40]), or Procrustes alignment [6, 41, 9]. Lastly, a common decoder is trained on the aligned across-condition data [22, 10, 9]. These procedures are popular as a replacement for costly and inconvenient BCI re-calibration in prosthetics [42], but their application to across-animal alignment remains rare [7, 6]. The multi-step procedure, although formally suboptimal, allows for some degree of flexibility in the individual component steps. For instance, the map into the latent manifold could be nonlinear, e.g. using Isomap [7, 43]. There also exist several nonlinear alignment procedures, such as kernel CCA or Distance Covariance Analysis, although such extensions have yet to translate into BCI practice [44]. In contrast, the hierarchical model proposed here is jointly estimated and linear in both the map and the latent dynamics. The linearity assumption here is not just a simplifying mathematical assumption, but it reflects domain-specific knowledge about the underlying neurophysiology of the system.
Neural activity has complex temporal dynamics and structured across-trial variability [45]. Previous approaches to aligning population activity across conditions via CCA or Procrustes alignment fail to capture either of these key statistical regularities. First, they treat the measurements across time as independent. Second, they intrinsically rely on paired measures, and since the correspondence of individual moments in time across trials in different animals is essentially arbitrary, they are restricted to across-trial average responses that may obscure important features of the underlying computation [46, 47]. It is also less clear how to determine the latent dimensionality of the data using these procedures. In contrast, amLDS extracts common low-dimensional dynamics using single trials, exploiting known features of the underlying computation. Moreover, the animal-specific observation model provides a principled framework for determining data dimensionality and for handling different sources of experimental variability (e.g. due to grid misalignment or missing data from defective electrodes). The probabilistic formulation also lends itself to a variety of generalizations. One could easily convert from Gaussian to Poisson observation noise for modeling spiking data instead of
5Some exceptions exist, such as probabilistic CCA [38], but they are not commonly used in practice.
the raw electric signals in our data [48, 49], or replace the linear dynamics with shared nonlinear equivalents [50, 51]. It should also be possible to extend the prior over trajectories to allow for some degree of heterogeneity across animals – for instance by allowing the parameters of the latent dynamics to slightly vary across animals, or to fall into a few distinct ‘dynamics classes’ (formally, an LDS mixture) [8]. Hence, amLDS can be thought of as yet another building block in the expanding statistical toolkit for understanding neural computation through the lens of dynamics in low-dimensional manifolds.
The representational stereotypy in the rodent olfactory bulb and the low-dimensional structure of its odor-induced dynamics are well documented [16–18, 13] make this system an ideal test case for out model. Nonetheless, we expect the same statistical tools to prove useful in other systems. The shared latent dynamics are best thought as reflecting a common computational solution, not a literal neuron-to-neuron, or subpopulation-to-subpopulation, match (in contrast to the hierarchical switching linear dynamical system approach of Ref. [8]). Hence, our model assumptions remain reasonable in circuits where representations are learned, rather than genetically encoded, as long as the dynamical system underlying the computation is low dimensional and the task strongly restricts the topology of the solution [52] (e.g. attractor dynamics for evidence integration [3]). We hope that our results inspire the search for common dynamic signatures in other brain areas and model systems.
Author contributions. CS and PHV conceived the original idea, developed the model and wrote the manuscript. PHV performed the numerical simulations and the data analysis. CS supervised the project. PHV collected the neural data in DR lab, previously published in Shor et al. 2022 [13].
Declaration of competing interest. DR is one of the founders of a company that commercializes related bio-electronic technology [13].
Acknowledgements. PHV is supported by training grant R90DA043849 (NIH). CS is supported by NIMH award 1R01MH125571-01, NSF award 1922658 and a Google research faculty award. DR was supported by DARPA BAA 15-35. We thank Caroline Haimerl, Edoardo Balzani and Erez Shor for their constructive feedback.

Review summaries:
Recent advances in recording techniques allow us to record from a large population of neurons - this high dimensional activity can typically be described by low dimensional latent factors. The authors leverage this low dimensional structure to capture the shared dynamical structure across animals in neural recordings from the olfactory bulb as animals are exposed to various odor stimuli. They fit a linear dynamical system (LDS) model with latent parameters and trajectories estimated for each odor that are shared across animals. These latent trajectories are subsequently mapped onto the neural recordings by employing animal specific observation parameters. The LDS model parameters are optimized using expectation maximization (EM). In order to decode the odor identities, the marginal likelihood for all stimulus conditions is computed. They validate the model on simulated data as well as neural recordings from the olfactory bulb. 

The report provides a method for decoding stimulus identity (odor) using data from multiple animals. The framework is straightforward and the performance appears respectable. I anticipate that this approach could be of use in neuroscience labs and could help with developing models and BCI. Given its general applicability, it would strengthen the impact if its use was demonstrated on spiking data and/or other types of stimuli or behaviors. This is a well written paper with helpful illustrations. I expect it would support interesting discussions.

This submission is a very well done application of the probabilistic inference to decoding of odour identity from neural recordings in mice. Interestingly, the latent code is assumed to have temporal dynamics shared across animals. The main benefits (high accuracy, data savings,  robustness to incomplete data) are shown in simulation and neural recordings. The method outperforms the existing (non-probabilistic) neural alignment procedures.  

The authors develop a method for aligning the neural activity of multiple animals. Specifically, they address the challenge of odor detection, where neural activity from multiple animals is recorded in response to many stimuli (odors), and then the odor needs to be decoded from a given animal’s neural activity. They use a hierarchical and probabilistic approach, in which there is a shared low-dimensional space of neural activity that is shared across animals, with different low->high dimensional mappings for each animal. The authors demonstrate the utility of the approach in simulations and on experimental neuroscience data.