Estimating a Separably-Markov Random Field (SMuRF) from Binary Observations

09/27/2017 ∙ by Yingzhuo Zhang, et al. ∙ 0

A fundamental problem in neuroscience is to characterize the dynamics of spiking from the neurons in a circuit that is involved in learning about a stimulus or a contingency. A key limitation of current methods to analyze neural spiking data is the need to collapse neural activity over time or trials, which may cause the loss of information pertinent to understanding the function of a neuron or circuit. We introduce a new method that can determine not only the trial-to-trial dynamics that accompany the learning of a contingency by a neuron, but also the latency of this learning with respect to the onset of a conditioned stimulus. The backbone of the method is a separable two-dimensional (2D) random field (RF) model of neural spike rasters, in which the joint conditional intensity function of a neuron over time and trials depends on two latent Markovian state sequences that evolve separately but in parallel. Classical tools to estimate state-space models cannot be applied readily to our 2D separable RF model. We develop efficient statistical and computational tools to estimate the parameters of the separable 2D RF model. We apply these to data collected from neurons in the pre-frontal cortex (PFC) in an experiment designed to characterize the neural underpinnings of the associative learning of fear in mice. Overall, the separable 2D RF model provides a detailed, interpretable, characterization of the dynamics of neural spiking that accompany the learning of a contingency.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 25

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

A fundamental problem in the analysis of electrophysiological data from neuroscience experiments is to determine the trial, and time within said trial, when a neuron or circuit first exhibits a conditioned response to a stimulus. This is a challenging problem because neural spike rasters resulting from such experiments can exhibit variability both within a given trial and across trials (Czanner ., 2008). Fear conditioning experiments (Allsop ., 2014) are a prime example of a scenario when this situation arises: a neutral stimulus, present across all trials of an experiment, gives rises to stereotypical within-trial spiking dynamics, while the associated aversive stimulus leads to changes in spiking dynamics across a subset of the trials.

State-of-the-art methods for analyzing neural spike rasters fall primarily within two classes. The most pervasive class of such methods neglect the inherent two-dimensional nature of neural spike rasters by aggregating the raster data either across time or trials, and subsequently applying techniques applicable to one-dimensional signals (Smith  Brown, 2003; Zammit-Mangion ., 2012; Yuan ., 2012; Scott  Pillow, 2012). In contrast to these one-dimensional methods, two-dimensional methods model both the within and cross-trial dynamics of neural spiking (Czanner ., 2008; Rad  Paninski, 2010). Within the class of one-dimensional methods, the past decade has seen a growing interest in approaches based on state-space models of neural spiking activity. These approaches treat neural spiking data as realizations of a stochastic point process whose conditional intensity function obeys a stochastic smoothness constraint in the form of a Markov process followed by a nonlinearity. The main challenge is to estimate the parameters of the model, and various solutions have been proposed towards this end (Smith  Brown, 2003; Zammit-Mangion ., 2012; Yuan ., 2012; Scott  Pillow, 2012). The main drawback of one-dimensional approaches applied to the analysis of neural spike rasters is the need, preceding analysis, for aggregation across one of the dimensions. Among one-dimensional methods, non-parametric methods based on rank tests (e.g. Wilcoxon rank sum test) have been the most popular, primarily due to their ease of application. In addition to the need to collapse neural activity of time or trials, two common pitfalls of non-parametric methods are their reliance on large sample assumptions to justify comparing neural spiking rates, and the need to correct for multiple comparisons. For instance, tests that rely on estimates of the neural spiking rate based on empirical averages are hard to justify when it is of interest to characterize the dynamics of neural spiking at the millisecond time scale. Consider a neural spike raster for which it is of interest to assess differences in instantaneous spiking rates between distinct time/trial pair. At the millisecond time scale, there would only be one observation per time/trial pair, violating the large sample assumptions that such non-parametric methods rely upon. To the best of our knowledge, the work of (Czanner ., 2008) remains the most successful attempt to characterize simultaneously the within and cross-trial dynamics of neural spiking. This approach uses a state-space model of the cross-trial dynamics, in conjunction with a parametric

model of the within-trial dynamics. The use of a parametric model for the within-trial dynamics is convenient because it enables the estimation of the model parameters by Expectation-Maximization (EM), using a combination of point-process filtering and smoothing in the E-step (to fill-in the missing cross-trial effect), and an M-step for the within-trial parameters that resembles a GLM 

(Truccolo ., 2005). The main drawbacks of this approach are, on the one hand, the high-dimensionality of the state-space model that captures the cross-trial dynamics, and on the other hand the lack of a simple interpretation, as in the one-dimensional models (Smith  Brown, 2003; Zammit-Mangion ., 2012; Yuan ., 2012), for the state sequence. Lastly, a two-dimensional approach based on Gaussian processes was proposed in (Rad  Paninski, 2010). One advantage of this approach, which is based on a Gaussian process prior of the neural spiking rate surface, is its ability to model the interaction between the two dimensions through the use of a two dimensional kernel. As is common with kernel methods, it does not scale well to multiple dimensions.

We propose a two-dimensional (2D) random field (RF) model of neural spike rasters–termed Separably-Markov Random Field (SMuRF)–in which the joint conditional intensity function of a neuron over time and trials depends on two latent Markovian state sequences that evolve separately but in parallel. Conventional methods for estimating state-space models from binary observations (Smith  Brown, 2003; Zammit-Mangion ., 2012; Yuan ., 2012) are not applicable to SMuRF. We derive a Monte Carlo Expectation-Maximization algorithm to maximize the marginal likelihood of observed data under the SMuRF model. In the E-step, we leverage the Polya-Gamma (Polson ., 2013)

representation of Bernoulli random variables to generate samples from the joint posterior distribution of the state sequences by Gibbs sampling. A similar strategy was adopted in 

(Scott  Pillow, 2012) for a one-dimensional state-space model. The sampler uses a highly efficient forward-filtering backward-sampling algorithm for which the forward step can be implemented exactly

and elegantly as Kalman filter, while the backward step uses Bayes’ rule to correct the filter samples. The SMuRF model obviates the need for aggregation across either time or trials, and yields a

low-dimensional 2D characterization of neural spike rasters that is interpretable in the sense that the posterior of the two state sequences capture the variability within and across trials respectively. Moreover, being model-based, the SMuRF model, unlike non-parametric methods, yields a characterization of the joint posterior (over all trials and time within a trial) distribution of the instantaneous rate of spiking, thus allowing us to precisely determine the dynamics of neural spiking that accompany the learning of a contingency. To demonstrate this, we apply the model to data collected from neurons in the pre-frontal cortex (PFC) in an experiment designed to characterize the neural underpinnings of the associative learning of fear in mice. We find that the trial at which the cortical neurons begin to exhibit a conditioned response to the auditory conditioned stimulus is robust across cells, occurring to trials into the conditioning period. We also find that the time with respect to conditioned stimulus onset when we observe a significant change in neural spiking compared to baseline activity varies significantly from cell to cell, occurring between to ms after conditioned stimulus onset. These findings are likely reflective of the variability in synaptic strength and connectivity that accompany learning, as well as the location of the neurons in the population.

The rest of our treatment begins in Section 2 where we motivate the SMuRF model, define it and introduce our notation. In Section 3, we present the Monte-Carlo EM algorithm for parameter estimation in the SMuRF model, as well as our process for inferring the dynamics of neural spiking that accompany the learning of a contingency by a neuron. The reader may find derivations relevant to this section in the Appendix. We present an application to the cortical data in Section 4, and conclude in Section 5.

2 Notation and SMuRF model

We begin this section with a continuous-time point-process formalism of a neural spike raster, characterized by a trial-dependent conditional intensity function (CIF). Then, we introduce the SMuRF model, a model for the discrete-time version of the CIF.

2.1 Continuous-time point-process observation model

We consider an experiment that consists of successive trials. During each trial, we record the activity of a neuronal spiking unit. We assume, without loss of generality, that the duration of the observation interval during each trial is . For trial , , let the sequence correspond to the times of occurrence of events from the neuronal unit, that is to say the times when the membrane potential of the neuron crosses a given threshold. We assume that is the realization in of a stochastic point-process with counting process , where is the indicator function in of . A point-process is fully characterized by its CIF. Let denote the trial-dependent CIF of defined as

(1)

where is the history of the point process up to time .

We denote by , the discrete-time process obtained by sampling at a resolution of , . Let denote the discrete-time, trial-dependent, CIF of the neuron.

2.2 Separably-Markov Random Field (SMuRF) model of within and cross-trial neural spiking dynamics

Definition 1

Let be a collection of random variables. We say that this collection is a separable random field if s.t. unique s.t. .
If in addition and are Markov processes, we say that is a separably-Markov random field or “SMuRF”.

A 2D random field is a collection of random variables indexed over a subset of

. We call this collection a separable field if there exists latent random vectors

and (each indexed over a subset of ) such that are independent conditioned on and and only a function of the outer product between and . If, in addition, and are Markov, we say that the field is a SMuRF. Intuitively, a separable random field is a random field that admits a stochastic rank-one decomposition.

We propose the following SMuRF model of the discrete-time, trial-dependent, CIF of a neuronal spiking unit

(2)

By construction, this is a SMuRF of the trial-dependent CIF of a neuron. and are indicator functions of presence of cue. To provide some intuition, if we assume is small, then the SMuRF model approximates the trial-dependent CIF as , that is, as the product of a within-trial component in units of Hz (spikes/s) and a unitless quantity . For a given trial , represents the excess spiking rate above what can be expected from the within-trial component at that trial, which we call the cross-trial component of the CIF. The within and cross-trial components from the SMuRF model are functions of two independent state sequences, and , that evolve smoothly according to a first-order stochastic difference equation. The parameters , , , , and , which govern the smoothness of and , must be estimated from the raster data.
Remark 1: We note that, in its generality, our model does not assume that . In our model, . The approximation holds for a neuron with small neural spiking rate (Truccolo ., 2005).

Figure 1

shows a graphical representation of the SMuRF model as a Bayesian network. It is not mathematically possible to rewrite the state equations from the SMuRF model in standard state-space form

without increasing significantly the dimension of the state space. We give a sketch of an argument as to why in the Appendix. Therefore, in an Expectation-Maximization (EM) algorithm for parameter estimation, one cannot simply apply classical (approximate) binary filtering and smoothing in the E-step (Smith  Brown, 2003). We derive a Monte-Carlo Expectation-Maximization algorithm to maximize the likelihood of observed data under the SMuRF model, with respect to the parameter vector .

Figure 1: Representation of the SMuRF model as a Bayesian network. The SMuRF model approximates the trial-dependent CIF the product of a within-trial component in units of Hz (spikes/s) and a unitless quantity. The within and cross-trial components are functions of two independent state sequences, and

, that evolve smoothly, each according to a first-order stochastic difference equation. Observations from the raster are Bernoulli random variables whose probability of of occurrence is a nonlinear function of the sum of the two sequences.

3 Parameter Estimation in the SMuRF by Maximum Likelihood

3.1 Maximum Likelihood Estimation by Expectation-Maximization

Let , , and . The goal is to maximize, with respect to , the likelihood of the SMuRF model

(3)

This is a challenging problem because of the high-dimensional integral that must be carried out in Equation 3. We propose to maximize the likelihood by EM.
Remark 2

: For the moment, we treat

and as missing data; in the sequel, we will augment the model with additional missing data that will simplify the EM algorithm.
Given a candidate solution , EM (Dempster ., 1977) maximizes by building a sequence of successive approximations of (the so-called E-step) such that maximizing these approximations, which in general is simpler than directly maximizing , is guaranteed to not decrease . That is, each iteration of EM generates a new candidate solution such that . By iterating this process, EM generates a sequence of iterates that, under regularity conditions, converge to a local optimum of  (Dempster ., 1977).

In the context of the SMuRF model, the key challenge of EM is to compute defined as

(4)

the expected value of the complete-data likelihood with respect to the joint posterior distribution of the missing data conditioned on the observed data and the candidate solution . This expectation is not tractable, i.e. it cannot be computed in closed-form. The intractability stems not only from the lack of conjugacy between the Bernoulli observation model and our Gaussian priors–also an issue for one-dimensional models (Smith  Brown, 2003)

–but also because, as mentioned previously, the SMuRF model cannot be reduced to a standard state-space model. We propose to approximate the required expectations using Markov-Chain Monte-Carlo (MCMC) samples from

. In particular, we will use Gibbs sampling (Casella  George, 1992), a Monte-Carlo technique, to generate samples from a distribution by sampling from its so called full conditionals (conditional distribution of one variable given all others), thus generating a Markov chain that, under regularity conditions, can be shown to converge to a sample from the desired distribution. Gibbs sampling is attractive in cases where sampling from the full-conditionals is simple. However, it is prone to the drawbacks of MCMC methods, such as poor mixing and slow convergence, particularly if one is not careful in selecting the full-conditionals from which to generate samples from. Two observations are in order, that will lead to the derivation of an elegant block Gibbs sampler with attractive properties

  • Conditioned on

    , the joint distribution,

    , of and is equivalent to the joint distribution from a one-dimensional state-space model with binary observations (Smith  Brown, 2003). By symmetry, this is also true for . This readily motivates a block Gibbs sampler that alternates between sampling from and . This leaves us with one challenge: how to obtain samples from the posterior distribution of the state in a one-dimensional state-space model with Bernoulli (more generally binomial) observations?

  • We introduce a new collection of i.i.d., Polya-Gamma distributed 

    (Polson ., 2013) random variables , such that sampling from is equivalent to sampling from the posterior of the state in a linear Gaussian state-space model (we will prove this in the Appendix) using a forward-filtering backward-sampling algorithm (Frühwirth-Schnatter, 1994). Moreover, it has been shown that the Gibbs sampler based on this Polya-Gamma augmentation scheme (Choi  Hobert, 2013)

    is uniformly ergodic and possesses superior mixing properties to alternate data-augmentation scheme for logit-based models 

    (Polson ., 2013; Choi  Hobert, 2013). The intuition behind the introduction of the Polya-Gamma random variables is the following: they are missing data that, if we could observe, would make the Bernoulli observations Gaussian. Stated otherwise, the Polya-Gamma random variables are scale variables in a Gaussian scale mixture (Andrews  Mallows, 1974) representation of Bernoulli random variables.

Remark 3: The random vector in the preceding bullet point is the vector additional missing data alluded to in Remark 2.

Together, these two observations form the basis of an efficient block-Gibbs sampler we use for maximum-likelihood estimation of the parameters from the SMuRF model by Monte-Carlo EM, also referred to as empirical Bayes (Casella, 2001). We introduce the basic ideas behind PG augmentation and its utility in Bayesian estimation for logit-based models. In the Appendix, we provide detailed derivations for the PG sampler adapted to the SMuRF model.

3.2 Polya-Gamma augmentation and sampling in one dimension

Let and and suppose that, conditioned on , is Bernoulli with mean , i.e.

(5)

We begin with a definition of Polya-Gamma (PG) random variables, followed by a PG augmentation scheme for the Bernoulli/binomial likelihood. We will see that the augmentation scheme leads to an attractive form for the posterior of given the observation and the augmented variable. Finally, we will see that the posterior of the augmented variable itself follows a PG distribution. Our treatment follows closely that of (Choi  Hobert, 2013).

Definition of Polya-Gamma random variables: Let be a sequence of i.i.d. exponential random variable with parameter equal to 1. The random variable

(6)

follows a PG(1,0) distribution, where

denotes equality in distribution. The moment generating function of

is

(7)

An expression for its density , expressed as an infinite sum, can be found in (Choi  Hobert, 2013) and (Polson ., 2013). The PG(1,) random variable is obtained by exponential tiling of the density of a PG(1,0) random variable. Letting denote the density of a PG(1,) random variable,

(8)

PG augmentation preserves the Bernoulli likelihood: Following the treatment of (Choi  Hobert, 2013), conditioned on , let be a PG(1,) random variable. Further suppose that, conditioned on , and are independent. Then

(9)

Integrating out , we see that the augmentation scheme does not alter . One may then ask, what is the utility of the augmentation scheme? The answer lies in the following identity, discussed in detail in (Choi  Hobert, 2013), and which is the key ideal behind PG augmentation

(10)

where , and indicates that we are dropping terms independent of . Equation 10 states that, given and a logit model, a Bernoulli random variable is, up to a constant independent of , a scale mixture of Gaussian (Andrews  Mallows, 1974)

, i.e. a Gaussian random variable with random variance

, where follows a PG distribution (Polson ., 2013; Choi  Hobert, 2013). If we assume , then (Choi  Hobert, 2013)

(11)

Implications of augmentation on and :

(12)

where we make use of Equation 11. If is Gaussian, then is Gaussian and available in closed-form! (Appendix).

(13)

i.e. . Together, Equations 12 and 13 form the basis of a uniformly ergodic (Choi  Hobert, 2013) Gibbs sampler to obtain sample from .

3.3 Block Gibbs sampler for PG-augmented SMuRF model

Consider the following version of the SMuRF model with PG augmentation:

(14)

We can apply the basic results from the previous subsection to derive the following result (proof in Appendix):

Theorem 2

Suppose , , and come from the PG-augmented SMuRF model (equation ), then is equivalent in distribution to the following linear-Gaussian state-space model

(15)

Following the discussion from the previous subsection, it is not hard to see that such a result would hold. The proof of this result is in the appendix, as well as the derivation of an elegant forward-filtering backward-sampling algorithm (Frühwirth-Schnatter, 1994) for drawing samples from . By symmetry, it is not hard to see that a similar result holds for .

Block Gibbs sampling from PG-augmented SMuRF model:

The E-step of the Monte-Carlo EM algorithm consists in sampling from by drawing from using a block Gibbs sampler that uses the following full-conditionals

  • , which according to the theorem above is equivalent to the posterior distribution of the state sequence in a linear-Gaussian state-space model.

  • , which obeys properties similar to the previous full-conditional (by symmetry).

  •  (Polson ., 2013).

In the Appendix, we detail how we initialize the algorithm and monitor convergence.

In practice, we found that estimating and is difficult. We hypothesize that including those parameters yields an unwieldy likelihood function. In the results we report, we assume , and focus on estimating a simple model with two parameters and . The assumption gives the random walk priors more freedom, thus allowing us to be capture the variability of the within and cross-trial processes. We have run simulations, not reported here, that show that the joint estimation of , , and is stable and that our EM algorithm converges. This demonstrates the ability of the SMuRF model (Equation (14)) to incorporate exogenous input stimuli.

Assuming , in the M-step, the update equations for the parameters , , and follow standard formulas (Smith  Brown, 2003)

(16)
(17)
(18)
(19)

where we set , and we approximate the expectations with respect to and using Gibbs samples from the E-step.

3.4 Assessment of within-trial and cross-trial spiking dynamics

Bayesian estimation of the SMuRF model (Equation (14)) enables us to infer detailed changes in neural dynamics, in particular to extract the within-trial and cross-trial components of the neural spiking dynamics that accompany the learning of a contingency by a neuron. This is because, following estimation, inference in the SMuRF model yields the joint posterior distribution of the instantaneous spiking rate of a neuron as a function of trials, and time within a trial, conditioned on the observed data. We can use this posterior distribution, in turn, to assess instantaneous changes in neural spiking dynamics, and without the need to correct for multiple comparisons as with non-parametric methods.

In what follows, we let denote the posterior distribution of and , given the raster data and the maximum likelihood estimate of . In what follows, it is understood that we use Gibbs samples from to obtain an empirical estimate of the distribution.

Posterior distribution of the joint CIF over time and trials: We can use these posterior samples to approximate the posterior distribution, at , of any quantities of interest. Indeed, it is well known from basic probability that if is a sample from , then is a sample from . In particular, if the instantaneous spiking rate of a neuron a time and trial , we can use the Gibbs samples to approximate the joint posterior distribution of given and .

Let be the random variable that represents the a posteriori instantaneous spiking rate of the neuron at time trial and time within that trial. The superscript ‘p’ highlights the conditioning on the data and , and the fact that this quantity is a function of distributed according to .

Within-trial effect: We define the within-trial effect as the a posteriori instantaneous spiking rate at time , average over all trials

(20)

It is important to note that the averaging is performed after characterization of the joint CIF as a function of time and trials, which is not the same as first aggregating the data across trials and applying one of the one-dimensional methods for analyzing neural data (Smith  Brown, 2003; Zammit-Mangion ., 2012; Yuan ., 2012). In practice, every Gibbs sample pair leads to a scalar quantity

(21)

Performing this computation over all Gibbs samples and times leads to a joint empirical distribution for the within-trial effect .

Cross-trial effect: We define the cross-trial effect as the a posteriori excess instantaneous spiking at trial and time (above the within-trial effect effect ) averaged across all times

(22)

In practice, every Gibbs sample pair leads to a scalar quantity

(23)

Performing this computation over all Gibbs samples and trials of interest, , leads to a joint empirical distribution for the cross-trial effect .
Remark 4: The following paragraph explains the meaning of in the context of an associative learning experiment.

3.5 Assessment of neural spiking dynamics across time and trials

Consider an associative learning (conditioning) experiment characterized by the pairing of a conditioned stimulus (e.g. auditory) to an aversive stimulus (e.g. a shock). Let be the number of conditioning trials and the length of the habituation period. Gibbs samples from the SMuRF model (Equation (14)) paramaterized by let us approximate the a posteriori probability that the spiking rate at a given point (Point C in Figure 2) during one of the conditioning trials (trials through in this example) is bigger than the baseline spiking rate at that trial (Region A in Figure 2) and the average spiking rate at the same time during the habituation period (Region B in Figure 2). This yields a probabilistic description of the the intricate dynamics of neural spiking that accompany the learning of the contingency by a neuron. Let

Event U (24)
Event V (25)

For a given pair s.t. , this probability is

(26)
(27)

where the second line approximates the probability of the event of interest using its frequency of occurrence in the posterior samples. As we demonstrate in the following section, we thus obtain an detailed characterization of the dynamics of neural spiking that accompany learning.

Figure 2: Regions defined to quantify changes in neural spiking dynamics in an associative learning experiment. The SMuRF model lets us approximate the a posteriori probability that the spiking rate at a given point C during one of the conditioning trials (trials through in this example) is bigger than the baseline spiking rate at that trial (Region A) and the average spiking rate at the same time during the habituation period (Region B). This yields a probabilistic description of the the intricate dynamics of neural spiking that accompany the learning of the contingency in the experiment by a neuron.

In the following section, we use simulated and real data examples to demonstrate the utility of the SMuRF model (Equation (14)) for the characterization of detailed neural spiking dynamics.

4 Applications

4.1 Simulation studies

We simulated neural spike raster data from a neuron that exhibits a conditioned response to the conditioned stimulus (Figure 4) in an associative learning experiment. The experiment consists of trials, each of which lasts s. The conditioned stimulus becomes active s into a trial, while the aversive stimulus becomes active after trial . We obtain the simulated data by dividing the raster into two pre-defined regions as shown in Figure 3. Region A consists of all trials before trial , along with the period from all trials before the conditioned stimulus is presented. We assume that the rate of spiking of the neuron is Hz. Region B consists of the period from trials following trial after the conditioned stimulus is presented. The rate of spiking of the neuron in this region is Hz.

Figure 3: Set up used to simulate neural spike raster data from a neuron that exhibits a conditioned response to the conditioned stimulus in an associative learning experiment. The experiment consists of trials, each of which lasts s. The conditioned stimulus becomes active s into a trial, while the aversive stimulus becomes active from the trial onwards. An idealized neuron that exhibits a conditioned response would exhibit two distinct regions of activity. Region A consists of all trials preceding trial , along with the period from all trials before the conditioned stimulus is presented. We assume that the rate of spiking of the neuron is Hz. Region B consists of the period from trials following trial after the conditioned stimulus is presented. The rate of spiking of the neuron in this region is Hz.

We applied the SMuRF model (Equation (14)) to the analysis of this simulated neural spike raster. Figure 4(a) shows that, during conditioning, learning is accompanied by a doubling of the spiking rate above the within-trial spiking rate of the neuron. Indeed, the left hand panel of the figure shows the cross-trial effect which, following conditioning, increases above its average initial value of to . Figure 4(b) provides a more detailed characterization of the neural spiking dynamics. With probability close to , the spiking rate at a given time/trial pair–following the conditioned stimulus and during conditioning (Figure 2 C)–is bigger than the average rate at the same trial (Figure 2 A) and the average rate at the same time (Figure 2 B). We conclude that, with high probability, the simulated neuron exhibits a conditioned response to the conditioned stimulus.

Figure 4: (a) Simulated neural spike raster, along with estimated within and cross trial effects from the SMuRF model. The horizontal red line indicates the beginning of conditioning. The vertical green line indicates the onset of the conditioned stimulus. The left panel of the figure shows the cross-trial effect which, following conditioning, increases above its average initial value of to . (b) Empirical probability that spiking rate at a given trial and time is bigger than the average rate at the same trial and the average rate during habituation at the same time. With probability close to , the spiking rate at a given time/trial pair–following the conditioned stimulus and during conditioning (Figure 2 C)–is bigger than the average rate at the same trial (Figure 2 A) and the average rate at the same time (Figure 2 B)

4.2 Neural dynamics during associative learning of fear

Basic experimental paradigm:

The ability to learn through observation is a powerful means of learning about aversive stimuli without direct experience of said stimuli. We use the SMuRF model (Equation (14)) to analyze data from a fear conditioning paradigm designed to elucidate the nature of the circuits that facilitate the associative learning of fear. The experimental paradigm is described in detail in (Allsop ., 2017). Briefly, an observer mouse observes a demonstrator receive conditioned stimulus-shock pairings through a perforated transparent divider. The experiment consists of to trials, divided into two phases. During the first trials of the experiment, termed the habituation period, both the observer and the demonstrator simply hear an auditory conditioned stimulus. From the trial onwards, the auditory conditioned stimulus is followed by the delivery of a shock to the demonstrator. The data are recorded from the pre-frontal cortex (PFC) of the observer mouse.

Results:

Figure 5(a) shows the within and cross-trial effects estimated using the SMuRF model applied to a cortical neuron from the experiment described above. The estimated within-trial (bottom) and cross-trial (left) components indicate significant changes respectively in response to the conditioned stimulus and to conditioning. By definition (Equation 22), the cross-trial effects takes into account the increase in spiking rate due to the presentation of the conditioned stimulus. The bottom panel suggests that this neuron exhibits a delayed response to the conditioned stimulus, beginning at ms following conditioned stimulus presentation. Accounting for this increase in within-trial spiking rate due to the conditioned stimulus, the left panel shows a multiplicative increase in spiking rate due to conditioning from an average initial value of (indicative of suppression, as can be seen through the sparseness of the raster during trials through ) to a peak average value of at trial . This increase, however, does not persist as in the case of the simulated data (Figure 4), suggesting that conditioning is accompanied by intricate dynamics in neural modulation.

Figure 5: (a) Raster, along with SMuRF within and cross-trial components, from a cortical neuron that exhibits a conditioned response to the conditioned stimulus. The horizontal red line indicates the beginning of conditioning. The vertical green line indicates conditioned stimulus onset. The bottom panel suggests that this neuron exhibits a delayed response to the conditioned stimulus, beginning at ms following conditioned stimulus presentation. Accounting for this increase in within-trial spiking rate due to the conditioned stimulus, the left panel shows a multiplicative increase in spiking rate due to conditioning from an average initial value of (indicative of suppression, as can be seen through the sparseness of the raster during trials through ) to a peak average value of at trial . This increase, however, does not persist as in the case of the simulated data (Figure 4).
Figure 5: (b) Empirical probability that spiking rate at a given trial in and time is bigger than the average rate at the same trial and the average rate during habituation at the same time (refer to Figure 2). This panel indicates that this neuron exhibit a delayed conditioning to the conditioned stimulus (beginning ms following conditioned stimulus presentation) and that the extent of the condition is highest first between trials and and then between trials and . Panel (a), and panel (b) in particular, suggest that conditioning is accompanied by intricate dynamics in neural modulation.

Figure 5 (b) provides a more detailed characterization of the neural spiking dynamics of this neuron. The figure shows the evolution, as a function of time and trials, of the probability that the spiking rate at a given time/trial pair (Figure 2 C) is bigger than the average rate at the same trial (Figure 2 A) and the average rate at the same time (Figure 2 B). The figure indicates that this neuron exhibit a delayed conditioning to the conditioned stimulus (beginning ms following conditioned stimulus presentation) and that the extent of the condition is highest first between trials and and then between trials and .

Figure 6 shows an application of the SMuRF model (Equation (14)) to a cortical neuron that does not exhibit a conditioned response to the conditioned stimulus. The bottom of panel (a) indicates no significant increase in the within-trial spiking rate in response to the conditioned stimulus, while the left panel shows that the cross-trial effect remains constant throughout the experiment an average value of . This indicates that conditioning does not result in a significant increase in spiking rate. Panel (b) corroborates these findings: for all points C following the conditioned stimulus and during conditioning, there is a small probability that the instantaneous spiking rate is significantly different from the average spiking rates in Regions A and B.

Figure 6: (a) Raster, along with SMuRF within and cross-trial components, from a cortical neuron that does not exhibit a conditioned response to the conditioned stimulus. The horizontal red line indicates the beginning of conditioning. The vertical green line indicates the onset of the conditioned stimulus. This neuron does not exhibit a conditioned response to the conditioned stimulus. The bottom of the panel indicates no significant increase in the within-trial spiking rate in response to the conditioned stimulus, while the left panel shows that the cross-trial effect remains constant throughout the experiment an average value of . This indicates that conditioning does not result in a significant increase in spiking rate.
Figure 6: (b) Empirical probability that spiking rate at a given trial and time is bigger than the average rate at the same trial and the average rate during habituation at the same time This panel corroborates the observations from panel (a): for all points C following the conditioned stimulus and during conditioning, there is a small probability that the instantaneous spiking rate is significantly different from the average spiking rates in Regions A and B.

Figures 10 and 11 show results for two additional cortical neurons that exhibit a transient conditioned response to the conditioned stimulus.

Using SMuRF inference to determine a neuron’s learning time and trial

The power of the Bayesian approach, and the SMuRF model (Equation (14)) in particular, lies in the fact that it lets us approximate the a posteriori probability that the spiking rate at a given point (point C in Figure 2) during one of the conditioning trials (trials through in this example) is bigger than the baseline spiking rate at that trial (Region A in Figure 2) and the average spiking rate at the same time during the habituation period (Region B in Figure 2) (Equation 27). This yields an instantaneous probabilistic quantification of the extent of learning for any given time and trial pair.

Panel (b) of Figures 5610 and 11 provide a detailed characterizations of the dynamics of learning and its extent for all times following the onset of the conditioned stimulus, all conditioning trials.

Here, we provide some guidance for practitioners to summarize the results of our inference to a single learning time/trial pair. We would like to stress, however, that the power of our methods lies in the detail provided by panel (b) of Figures 5610 and 11. Since the SMuRF model enables us to compute a empirical probability that the spiking rate at a given time/trial pair (Figure 2, Point C) is bigger than the average rate at the same trial (Figure 2, Region A) and the average rate at the same time (Figure 2, Region B), we can identify learning time and learning trial for each neuron by finding the first time after cue and after conditioning that this probability exceeds a certain threshold. Table  1 reports the learning time and trial computed using a threshold of of 95%. Note that the learning time is computed with respect to the onset of the conditioned stimulus (time ms).

Neuron in Figure 5 Neuron in Figure 6 Neuron in Figure 10 Neuron in Figure 11
Learning time (ms) 617 1316 202 20
Learning trial 16 34 18 18
Table 1: Learning trial and time computed for the cortical neurons analyzed. The learning times and trials reported are consistent with the detailed inference provided by the respective Figures for these neurons. The cortical unit from Figure 5, for instance, shows a delayed response, significant ms after conditioned stimulus onset and at trial . The cortical unit from Figure 6 only exhibits a significant change in neural spiking ms following the conditioned stimulus and at trial . This is consistent with our previous observation from Figure 6 that this neuron does not exhibit a conditioned response to the stimulus.

The learning times and trials reported in Table 1 are consistent with the detailed inference provided by the respective Figures for these neurons. Indeed, the cortical unit from Figure 5 shows a delayed response, significant ms after conditioned stimulus onset and at trial . The cortical unit from Figure 6 only exhibits a significant change in neural spiking ms following the conditioned stimulus and at trial . This is consistent with our previous observation from Figure 6 that this neuron does not exhibit a conditioned response to the stimulus.

In the Appendix, we perform a simulation that demonstrates the ability of SMuRF inference to identify learning time and trial when learning of a contingency is accompanied by sustained changes in neural spiking following conditioned stimulus onset and during conditioning. We also demonstrate through simulation that the SMuRF model is robust to the presence of error trials.

4.3 Application of SMuRF model to a non-separable example

We demonstrate the limitations of the separability assumption in the SMuRF model (Equation (14)) by applying it to the neural spike raster data from (Czanner ., 2008). We briefly describe the experiment here and refer the reader to (Wirth ., 2003) for a more detailed description. Panel (a) of Figure 7 shows neural spiking activity from a hippocampal neuron recorded during an experiment designed for a location-scene association learning task. The same scene was shown to a Macaque monkey across 55 trials, and each trial lasted 1700 ms. The first 300 ms of every trial is fixation period, and the scene is presented to the monkey from 300 to 500 ms. A delay period takes place from 800 to 1500 ms, followed by a response period from 1500 to 1700 ms.

The data from the experiment are shown in the center of panel (a) from Figure 7. The raster suggests that the time and trial-dependent CIF of this neuron is not separable. this Intuitively, this can be seen from the fact that the region in which there are significant changes in neural spiking does not follow the rectangular form from Figure 3. Nevertheless, the CIF could be well approximated by a separable model. We apply the SMuRF model to these data to uncover some of its limitations in non-separable settings. The bottom panel of Figure 7(a) shows the estimate of the within-trial effect from the SMuRF model, while the left panel shows the cross-trial effect. These two figures indicate that the SMuRF model is able to capture within and cross-trial dynamic changes in the spiking activity of the neuron. Figure  7(b) shows the estimate of the a-posteriori mean instantaneous spiking rate (in Hz) of the neuron at time trial and time within that trial. This figure shows that, while the SMuRF model is able to characterize the detailed changes in spiking dynamics, it does not fully capture the non-separable nature of the raster data.

Remark 5: Unlike for the cortical neurons, this experiment does not have a conditioning period. That’s why, it does not make sense to generate plots such as Figure 5(b).

Figure 7: Application of the SMuRF model to neural spiking activity from a hippocampal neuron recorded during an experiment designed for a a location-scene association learning task.
Figure 7: The green lines represent the end of the fixation, scene presentation and delay periods respectively. (a) Neural spike raster from the hippocampal neuron. The bottom panel shows the estimate of the within-trial effect from the SMuRF model, while the left panel shows the cross-trial effect. (b) Estimate of the a-posteriori instantaneous spiking rate (in Hz) of the hippocampal neuron at time trial and time within that trial. This figure shows that, while the SMuRF model is able to characterize the detailed changes in spiking dynamics of this neuron during the task, it fails to capture the non-separable nature of the raster data.

5 Conclusion

We proposed a 2D separably-Markov random field (SMuRF) for the analysis of neural spike rasters that obviates the need to aggregate data across time or trials, as in classical one-dimensional methods (Smith  Brown, 2003; Zammit-Mangion ., 2012; Yuan ., 2012), while retaining their interpretability. The SMuRF model approximates the trial-dependent conditional intensify function (CIF) of a neuron as the product of a within-trial component, in units of Hz (spikes/s), and a unitless quantity, which we call the cross-trial effect, that represents the excess spiking rate above what can be expected from the within-trial component at that trial. One key advantage of our 2D model-based approach over non-parametric methods stems from the fact that it yields a characterization of the joint posterior (over all trials and times within a trial) distribution of the instantaneous rate of spiking of as a function of both time and trials given the data. This not only obviates the need to correct for multiple comparisons, but also enables us to compare the instantaneous rate of any two trial time pairs at the millisecond resolution, where non-parametric methods break down because the sample size is 1.

We applied the SMuRF model to data collected from neurons in the pre-frontal cortex (PFC) in an experiment designed to characterize the neural underpinnings of the associative learning of fear in mice. We found that, as a group, the recorded cortical neurons exhibit a conditioned response to the auditory conditioned stimulus, occurring to trials into conditioning. We also found intricate and varied dynamics of the extent to which the cortical neurons exhibit a conditioned response (e.g. delays, short-term conditioning). This is likely reflective of the variability in synaptic strength, connectivity and location of the neurons in the population.

In future work, we plan to investigate non-separable random field models of neural spike rasters, such as Markov random fields (Besag, 1974) (MRFs). Compared to the SMuRF model, MRFs are 2D models for which the dimensionality of the putative state-space is as large as the dimensionality of the raster, suggesting that MRFs may provide a more detailed characterizations of neural spike rasters. Indeed, the SMuRF model makes the strong assumption that the neural spiking dynamics are decomposable into two time scales, with the additional simplifying assumption that there is only one component per time scale. This simplifying assumption is motivated by one-dimensional state-space models of neural data (Smith  Brown, 2003) in which a neuron’s time-dependent CIF is only a function of one hidden state sequence. We will investigate the inclusion of additional components in future work. We also plan to investigate analogues of the SMuRF model for population level data. MRFs, multi-component and population-level SMuRF models, naturally lead to model selection problems, and to the investigation of tools, based on sequential Monte-Carlo methods (Chopin ., 2013) (aka particle filters), to compare state-space models of neural spike rasters (such as one-dimensional models (Smith  Brown, 2003), the SMuRF model, and MRFs). The development of such tools for model comparisons is, in our opinion, the ultimate measure of the ability of different models to capture the intricate dynamics present in neural spike rasters. Lastly, as previously mentioned, the SMuRF model can be interpreted as a two-dimensional Gaussian process prior on the neural spiking rate surface (Rad  Paninski, 2010), with a separable kernel that is the Kronecker product of kernels from Gauss-Markov processes (one process for each dimension). The choice of kernels in the SMuRF model leads to the very efficient algorithms for estimation and inference derived in this article. Moreover, these algorithms scale well to more than two dimensions unlike classical kernel methods. We plan to explore this connection to Gaussian process inference in future work.

Acknowledgements

We would like to thank Dr. Anne C. Smith for her generous feedback on this manuscript and extensive discussions regarding the SMuRF model. Demba thanks the Alfred P. Sloan Foundation. K.M.T. is a New York Stem Cell Foundation–Robertson Investigator and McKnight Scholar and this work was supported by funding from the JPB Foundation, the PIIF and PIIF Engineering Award, PNDRF, JFDP, Alfred P Sloan Foundation, New York Stem Cell Foundation, McKnight Foundation, R01-MH102441-01 (NIMH), RF1-AG047661-01 (NIA), R01-AA023305-01 (NIAAA) and NIH Director’s New Innovator Award DP2-DK-102256-01 (NIDDK).

Appendix

The SMuRF model cannot be converted easily to a standard state-space model

We focus on the simple case when , and

Let , , . The index is obtained by “unstacking” the raster trials and serializing them.

The question we ask is whether the state equations from the SMuRF model can be turned into ones of the form

(28)

where . Let and denote the first and second components of respectively. We ask that because of the two dimensions present in the SMuRF model. Allowing the dimensionality of to increase up to would allow a representation of the form of Equation 28. However, this would become a very high-dimensional, unwieldy state-space model.

Intuitively, this cannot be done for the following reason: the dimensionality of the latent states in the SMuRF model is , while the dimensionality of the state sequence in Equation 28 is . For there to be an equivalence, the sequence must necessarily be redundant, i.e. some of the states must be copies of previous states. Storing these copies, would necessarily mean having to increase the dimensionality of the state space!

Let . The quantity gives the trial index corresponding time index . The within-trial index corresponding to index is then obtained by substracting from .

Note, for instance, that and . In general, for some if and only if s.t. for some integer , where we assume without loss of generality that . That is, two different indices and share the same within-trial component if and only if they are apart by an integer multiple of . Stated otherwise, the first component of exhibit circular symmetry! Therefore, for Equation 28 to hold, must equal , which is not possible because and are not apart by an integer multiple of !

The argument above shows that, in order to write the SMuRF state equations in the form of Equation 28, one would need to augment the state to dimension , which would lead to a very high dimensional standard state-space model, thus increasing the complexity of performing inference.

Derivation of Gibbs sampler for PG-augmented SMuRF model

We first derive Theorem 2, which leads to the forward-filter backward-sampling algorithm from the full-conditionals for and in the Gibbs sampler.

Since is drawn from a PG distribution, we can write the log pdf of as,

(29)

The complete data likelihood of the SMuRF model is,

(30)
(31)

The log of the complete data likelihood is therefore,

(32)
(33)
(34)