A Bayesian hierarchical latent variable model for spike train data analysis

A common approach to analyzing spike train data in stimulus-response experiments is to estimate spike rates relative to the stimulus onset. These experiments typically involve collecting measurements from the same subject across several trials and sessions, and performing the same experiment across multiple subjects. Accounting for these sources of variability is important to accurately estimate the population response. However, current spike rate estimation methods only allow for estimating trial-level, session-level or group-level spike rates, neglecting the variability between levels. We develop a Bayesian hierarchical latent variable model for estimating population and subject-level spike rate functions from multi-subject, multi-session and multi-trial spike train data. We model the contribution of a stimulus to spiking dynamics through a stochastic ordinary differential equation model, where a stimulus has a cumulative effect on latent dynamics that decays exponentially upon its offset. This effect is governed by parameters that can be specified as prior distributions. Thus, the model can capture a broad range of dynamics including sustained inhibition or excitation and transient inhibition or excitation. Our approach provides goodness-of-fit assessment and a Monte Carlo algorithm for computing confidence intervals for population and subject-level spike rate functions, and for performing uncertainty quantification on model parameters. We assess the accuracy of our model in simulation. We apply our model to data from previously published stimulus-response experiments to characterize the effect of: 1. thalamic stimulation on cortical spiking in non-human primates under propofol-induced unconsciousness; 2. a repetitive whisker stimulus on simultaneous recordings from individual thalamic barreloids in the rat somatosensory whisker/barrel system; 3. photoinhibition of the anterior lateral motor cortex during a pole location discrimination task. By using a hierarchical Bayesian framework, our model pools information across levels and can estimate population-level spike rates that accurately describe data across all levels of the hierarchy.