arXiv:1709.09723v1 [stat.ME] 27 Sep 2017

1

Estimating a Separably-Markov Random Field
(SMuRF) from Binary Observations
Yingzhuo Zhang1 , Noa Malem-Shinitski2 , Stephen A Allsop3 , Kay
Tye3 and Demba Ba1
1 Harvard University, John A. Paulson School of Engineering and Applied Sciences.
2 Technische Universität Berlin.
3 Massachusetts Institute of Technology, Department of Brain and Cognitive Sciences.

Keywords: Spike rasters, Dynamics, Random field, Bayesian Estimation, associative
learning.

Abstract
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.

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 et al., 2008). Fear conditioning experi-

2

ments (Allsop et al., 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 et al., 2012; Yuan et al., 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 et al., 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 et al., 2012; Yuan et al., 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

3

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 nonparametric methods rely upon. To the best of our knowledge, the work of (Czanner et
al., 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 et al., 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 et al., 2012;
Yuan et al., 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,

4

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
et al., 2012; Yuan et al., 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
et al., 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

5

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 3 to 4 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 20 to 600 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.

6

2.1

Continuous-time point-process observation model

We consider an experiment that consists of R 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 (0, T ]. For trial r,
r = 1, · · · , R, let the sequence 0 < tr,s < · · · < tr,Sr < T correspond to the times of
occurrence of events from the neuronal unit, that is to say the times when the membrane
r
is the realpotential of the neuron crosses a given threshold. We assume that {tr,s }Ss=1

ization in (0, T ] of a stochastic point-process with counting process Nr (t) =

Rt
0

dNr (u),

r
where dNr (t) is the indicator function in (0, T ] of {tr,s }Ss=1
. A point-process is fully

characterized by its CIF. Let λr (t|Ht ) denote the trial-dependent CIF of dNr (t) defined
as
P [Nr (t + ∆) − Nr (t) = 1|Ht ]
,
∆→0
∆

λr (t|Ht ) = lim

(1)

where Ht is the history of the point process up to time t.
We denote by {∆Nk,r }K,R
k=1,r=1 , the discrete-time process obtained by sampling dNr (t)
at a resolution of ∆, K =

T 
. Let {λk,r }K,R
k=1,r=1 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

K×R
Definition 1 Let {yk,r }K,R
be a collection of random variables. We say
k=1,r=1 ∈ R

that this collection is a separable random field if ∃ x ∈ RK , z ∈ RR s.t. ∀k, r ∃ unique
(xk , zr ) ∈ x × z s.t. yk,r |(x, z) ∼ f (xk , zr ).
7

If in addition x and z are Markov processes, we say that {yk,r }K,R
k=1,r=1 is a separablyMarkov random field or “SMuRF”.
K×R
A 2D random field {yk,r }K,R
is a collection of random variables ink=1,r=1 ∈ R

dexed over a subset of N+ × N+ . We call this collection a separable field if there exists
latent random vectors x and z (each indexed over a subset of N+ ) such that {yk,r }K,R
k=1,r=1
are independent conditioned on x and z and only a function of the outer product between x and z. If, in addition, x and z 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
{λk,r }K,R
k=1,r=1 of a neuronal spiking unit




xk = ρx xk−1 + αx ux,k + k , k ∼ N (0, σ2 )







 zr = ρz zr−1 + αz uz,k + δr , δr ∼ N (0, σδ2 )

(2)


λ ∆


log 1−λk,rk,r ∆ = xk + zr







 ∆Nk,r |xk , zr ∼ Bernoulli(λk,r ∆)
By construction, this is a SMuRF of the trial-dependent CIF of a neuron. ux,k and
uz,k are indicator functions of presence of cue. To provide some intuition, if we assume xk + zr is small, then the SMuRF model approximates the trial-dependent CIF as
λk,r ∆ ≈ ezr · exk , that is, as the product of a within-trial component exk in units of Hz
(spikes/s) and a unitless quantity ezr . For a given trial r, ezr 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

8

from the SMuRF model are functions of two independent state sequences, (xk )K
k=1 and
(zr )R
r=1 , that evolve smoothly according to a first-order stochastic difference equation.
The parameters ρx , αx , σ2 , ρz , αz and σδ2 , which govern the smoothness of (xk )K
k=1 and
(zr )R
r=1 , must be estimated from the raster data.
Remark 1: We note that, in its generality, our model does not assume that λk,r ∆ =
ezr · exk . In our model, λk,r ∆ =

exk +zr
.
1+exk +zr

The approximation λk,r ∆ ≈ ezr · exk holds

for a neuron with small neural spiking rate (Truccolo et al., 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 θ = (ρx , αx , σ2 , ρz , αz , σδ2 ).

9

∆NK,R

zR

∆N1,R

∆N2,R

···

∆Nk,R

···

···

···

···

···

···

···

zr

∆N1,r

∆N2,r

···

∆Nk,r

···

∆NK,r

···

···

···

···

···

···

···

z1

∆N1,1

∆N2,1

···

∆Nk,1

···

∆NK,1

x1

x2

···

xk

···

xK

···

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 compoR
nents are functions of two independent state sequences, (xk )K
k=1 and (zr )r=1 , 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.

10

3

Parameter Estimation in the SMuRF by Maximum
Likelihood

3.1

Maximum Likelihood Estimation by Expectation-Maximization

Let x = (x1 , · · · , xK )T , z = (z1 , · · · , zR )T , and ∆N = {∆Nk,r }K,R
k=1,r=1 . The goal is to
maximize, with respect to θ, the likelihood L(θ|∆N) of the SMuRF model
Z
L(θ|∆N) = log p(∆N; θ) = log

p(∆N, x, z; θ)dxdz.

(3)

x,z

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 x and z 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 et al., 1977) maximizes L(θ|∆N) by
building a sequence of successive approximations Q(θ|θ(`) ) of L(θ|∆N) (the so-called
E-step) such that maximizing these approximations, which in general is simpler than
directly maximizing L(θ|∆N), is guaranteed to not decrease L(θ|∆N). That is, each
iteration of EM generates a new candidate solution θ(`+1) such that L(θ(`+1) |∆N) ≥
L(θ(`) |∆N). By iterating this process, EM generates a sequence of iterates {θ(`) }∞
`=1
that, under regularity conditions, converge to a local optimum of L(θ|∆N) (Dempster
et al., 1977).
In the context of the SMuRF model, the key challenge of EM is to compute Q(θ|θ(`) )
defined as


Q(θ|θ(`) ) = Ex,z log p(∆N, x, z; θ)|∆N, θ(`) ,
11

(4)

the expected value of the complete-data likelihood with respect to the joint posterior
distribution of the missing data (x, z) conditioned on the observed data ∆N 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
p(x, z|∆N; θ(`) ). 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 z, the joint distribution, p(∆N, x|z; θ), of x and ∆N 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
p(∆N, z|x; θ). This readily motivates a block Gibbs sampler that alternates
between sampling from x|∆N; θ and z|∆N; θ. This leaves us with one chal12

lenge: 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 et al.,
2013) random variables w = {wk,r }K,R
k=1,r=1 , such that sampling from x|∆N, w; θ
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 (Frü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 et al., 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 PolyaGamma random variables are scale variables in a Gaussian scale mixture (Andrews & Mallows, 1974) representation of Bernoulli random variables.
Remark 3: The random vector w 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 logitbased models. In the Appendix, we provide detailed derivations for the PG sampler
13

adapted to the SMuRF model.

3.2

Polya-Gamma augmentation and sampling in one dimension

Let ∆N ∈ {0, 1} and X ∈ R and suppose that, conditioned on X = x, ∆N is Bernoulli
with mean

ex
,
1+ex

i.e.
p(∆N |x) =

(ex )∆N
1 + ex

(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 x given the observation ∆N 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 {Em }∞
m=1 be a sequence of i.i.d.
exponential random variable with parameter equal to 1. The random variable
d

W =

∞
Em
2 X
π 2 m=1 (2m − 1)2

(6)

d

follows a PG(1,0) distribution, where = denotes equality in distribution. The moment
generating function of W is
−tW

E[e

] = cosh

−1

√ 
t
.
2

(7)

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

14

the density of a PG(1,c) random variable,
pW (w|c) = cosh

c
2

e−

c2 w
2

pW (w).

(8)

PG augmentation preserves the Bernoulli likelihood: Following the treatment of (Choi
& Hobert, 2013), conditioned on X = x, let W be a PG(1,|x|) random variable. Further
suppose that, conditioned on X = x, ∆N and W are independent. Then
p(∆N, w|x) = p(∆N |x)pW (w|x).

(9)

Integrating out W , we see that the augmentation scheme does not alter p(∆N |x). 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
(ỹ−x)2
(ex )∆N
− 21 1/w
·
PG(1,
|x|)
∝
e
∝N
p(∆N |x)pW (w|x) =
1 + ex

where ỹ =

y− 12
w



1
ỹ; x,
w


,

(10)

, and ∝ indicates that we are dropping terms independent of x. Equa-

tion 10 states that, given X = x and a logit model, a Bernoulli random variable is, up to
a constant independent of x, a scale mixture of Gaussian (Andrews & Mallows, 1974),
i.e. a Gaussian random variable with random variance 1/w, where W = w follows a
PG distribution (Polson et al., 2013; Choi & Hobert, 2013). If we assume X ∼ pX (x),
then (Choi & Hobert, 2013)


1
p(∆N, x, w) = p(∆N |x, w)pW (w|x)pX (x) ∝ N ỹ; x,
pX (x).
w

(11)

Implications of augmentation on p(x|∆N, w) and p(w|∆N, x):
p(∆N, x, w)
p(x|∆N, w) =
∝ p(∆N |x, w)pW (w|x)pX (x) ∝ N
p(∆N, w)



1
ỹ; x,
w


pX (x),
(12)

15

where we make use of Equation 11. If X is Gaussian, then p(x|∆N, w) is Gaussian and
available in closed-form! (Appendix).

pW (w|∆N, x) =

p(∆N, x, w)
p(∆N |x)pW (w|x)pX (x)
=R
= pW (w|x),
p(∆N, x)
p(∆N
|x)p
(w|x)p
(x)
W
X
w

(13)

i.e. p(w|∆N, x) = pW (w|x) = PG(1, |x|). Together, Equations 12 and 13 form the
basis of a uniformly ergodic (Choi & Hobert, 2013) Gibbs sampler to obtain sample
from p(x, w|∆N ).

3.3

Block Gibbs sampler for PG-augmented SMuRF model

Consider the following version of the SMuRF model with PG augmentation:




xk = ρx xk−1 + αx ux,k + k , k ∼ N (0, σ2 )








zr = ρz zr−1 + αz uz,k + δr , δr ∼ N (0, σδ2 )




x +zr

e k
λk,r ∆ = 1+e
xk +zr







∆Nk,r |xk , zr ∼ Bernoulli(λk,r ∆)







 wk,r |xk , zr ∼ PG(1, |xr + zr |), k = 1, · · · , K; r = 1, · · · , R.

(14)

We can apply the basic results from the previous subsection to derive the following
result (proof in Appendix):
Theorem 2 Suppose w, x, z and ∆N come from the PG-augmented SMuRF model
(equation ), then p(∆N, x|w, z; θ) is equivalent in distribution to the following linearGaussian state-space model




xk = ρx xk−1 + αx ux,k + k , k ∼ N (0, σ2 )





∆Ñk,r = xk + zr + ṽk,r , ṽk,r ∼ N (0, wk,r −1 ), i.i.d. , r = 1, · · · , R




1


 ∆Ñk,r = ∆Nk,r − 2 .
wk,r
16

(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 (Frühwirth-Schnatter, 1994)
for drawing samples from p(∆N, x|w, z; θ). By symmetry, it is not hard to see that a
similar result holds for p(∆N, z|w, x; θ).

Block Gibbs sampling from PG-augmented SMuRF model: The E-step of the MonteCarlo EM algorithm consists in sampling from p(x, z|∆N; θ(`) ) by drawing from p(x, z, w|∆N ; θ(`) )
using a block Gibbs sampler that uses the following full-conditionals
• p(x|∆N, w, z; θ(`) ), which according to the theorem above is equivalent to the
posterior distribution of the state sequence in a linear-Gaussian state-space model.
• p(z|∆N, w, x; θ(`) ), which obeys properties similar to the previous full-conditional
(by symmetry).
• p(wk,r |x, z) = p(wk,r |xk , zr ) = PG(1, |xk +zr |), k = 1, · · · , K, r = 1, · · · , R (Polson et al., 2013).
In the Appendix, we detail how we initialize the algorithm and monitor convergence.
In practice, we found that estimating ρx and ρz is difficult. We hypothesize that
including those parameters yields an unwieldy likelihood function. In the results we
report, we assume ρx = ρz = 1, αx = αz = 0 and focus on estimating a simple
model with two parameters σ2 and σδ2 . The assumption ρx = ρz = 1 gives the random
walk priors more freedom, thus allowing us to be capture the variability of the within
17

and cross-trial processes. We have run simulations, not reported here, that show that the
joint estimation of σδ2 , σ2 , αx and αz is stable and that our EM algorithm converges. This
demonstrates the ability of the SMuRF model (Equation (14)) to incorporate exogenous
input stimuli.
Assuming ρx = ρz = 1, in the M-step, the update equations for the parameters σ2 , σδ2 ,
αx and αz follow standard formulas (Smith & Brown, 2003)
PK
αx(`+1) =

k=1

PR
αz(`+1) =

r=1



Ex xk − xk−1 |∆N, θ(`) ux,k
,
PK 2
k=1 ux,k


Ez zr − zr−1 |∆N, θ(`) uz,r
,
PR 2
r=1 uz,r

#
K
1 X
= Ex
(xk − xk−1 − αx(`+1) uz,k )2 |∆N, θ(`) ,
K k=1
"
#
R
1X
= Ez
(zr − zr−1 − αz(`+1) uz,k )2 |∆N, θ(`) ,
R r=1

(16)
(17)

"

σ2(`+1)
2(`+1)

σδ

(18)

(19)

where we set x0 = z0 = 0, and we approximate the expectations with respect to
p(x|∆N, θ(`) ) and p(z|∆N, θ(`) ) 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 dynam-

18

ics, and without the need to correct for multiple comparisons as with non-parametric
methods.

In what follows, we let p(x, z|∆N ; θ̂M L ) denote the posterior distribution of x and
z, given the raster data ∆N and the maximum likelihood estimate θ̂M L of θ. In what
follows, it is understood that we use Gibbs samples (xi , zi )ni=1 from p(x, z|∆N ; θ̂M L )
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 θ̂M L , of any quantities of
interest. Indeed, it is well known from basic probability that if (xi , zi ) is a sample from
p(x, z|∆N ; θ̂M L ), then f (xi , zi ) is a sample from p(f (x, z)|∆N ; θ̂M L ). In particular, if
the instantaneous spiking rate of a neuron a time k and trial r λk,r ∆ =

exk +zr
,
1+exk +zr

we can

use the Gibbs samples to approximate the joint posterior distribution of {λk,r ∆}K,R
k=1,r=1
given ∆N and θ̂M L .

Let {λpk,r ∆}K,R
k=1,r=1 be the random variable that represents the a posteriori instantaneous
spiking rate of the neuron at time trial r and time k within that trial. The superscript ‘p’
highlights the conditioning on the data ∆N and θ̂M L , and the fact that this quantity is a
function of (x, z) distributed according to p(x, z|∆N ; θ̂M L ).
Within-trial effect: We define the within-trial effect as the a posteriori instantaneous
spiking rate at time k, average over all trials
R

eWT
k =

1X p
λ (xk , zr ), k = 1, · · · , K.
R r=1 k,r

(20)

It is important to note that the averaging is performed after characterization of the joint
19

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 et al., 2012; Yuan et al., 2012). In
practice, every Gibbs sample pair (xi , zi ), i = 1, · · · , n leads to a scalar quantity
R

1X
=
λk,r (xi,k , zi,r ), k = 1, · · · , K.
R r=1

êWT
i,k

(21)

Performing this computation over all Gibbs samples and times k = 1, · · · , K leads to a
T K
joint empirical distribution for the within-trial effect {eW
k }k=1 .

Cross-trial effect: We define the cross-trial effect as the a posteriori excess instantaT
neous spiking at trial r and time k (above the within-trial effect effect eW
k ) averaged

across all times k
eCT
r

p
K
1 X λk,r (xk , zr )
, r = 1, · · · , R.
=
T
K k=1
eW
k

(22)

In practice, every Gibbs sample pair (xi , zi ), i = 1, · · · , n leads to a scalar quantity
êCT
i,r

K
1 X λk,r (xi,k , zi,r )
=
, r = 1, · · · , R.
T
K k=1
êW
i,k

(23)

Performing this computation over all Gibbs samples and Rc trials of interest, r =
R − Rc + 1, · · · , R, leads to a joint empirical distribution for the cross-trial effect
R
{eCT
r }r=R−Rc .

Remark 4: The following paragraph explains the meaning of Rc 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
20

Rc be the number of conditioning trials and Kh the length of the habituation period.
Gibbs samples from the SMuRF model (Equation (14)) paramaterized by θ̂M L 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 16 through 45 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


Average rate in Region A 

z
}|
{






R
c


X
1
p
p
λ (xk , zm )
Event U = λk,r (xk , zr ) >


Rc m=1 k,m










Average rate in Region B


z
}|
{






K
h


X
1
p
p
Event V = λk,r (xk , zr ) >
λ (xs , zr )


Kh s=1 s,r









(24)

(25)

For a given pair (k, r) s.t. k ≥ Kh , r ≥ Rc , this probability is
P [Event U ∩ Event V]

(26)

n

1X
≈
I
n i=1 {λpk,r (xi,k ,zi,r )> R1c

R
Pc
m=1

p

p

λk,m (xi,k ,zi,m )∩λk,r (xi,k ,zi,r )> K1

K
Ph

h s=1

p

, (27)

λs,r (xi,s ,zi,r )}

where the second line approximates the probability of the event of interest using its
frequency of occurrence in the n posterior samples. As we demonstrate in the following
section, we thus obtain an detailed characterization of the dynamics of neural spiking
that accompany learning.
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.
21

45

C: Trial r at time k
Cue
A: Baseline during
conditioning
B: Average response
at time k
during habituation

41
36

Trials

31
26
C

21

A

16
11

B

6
1
-500

0

500

1000

1500

Time (ms)

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 16 through 45 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.

4
4.1

Applications
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 exper-

22

iment consists of 45 trials, each of which lasts 2 s. The conditioned stimulus becomes
active 1 s into a trial, while the aversive stimulus becomes active after trial 15. 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 16, along with the period from all
trials before the conditioned stimulus is presented. We assume that the rate of spiking
of the neuron is λA = 60 Hz. Region B consists of the period from trials following trial
15 after the conditioned stimulus is presented. The rate of spiking of the neuron in this
region is λB = 20 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 ≈ 1 to ≈ 2. Figure 4(b) provides a more
detailed characterization of the neural spiking dynamics. With probability close to 1,
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.

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.

23

45
41
36
Region A:

Trials

31

λA = 60 Hz

26
21
16

Region B:
λB = 20 Hz

11
6
1
-1000

-500

0

500

1000

Time (ms)

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 45 trials, each of which lasts 2 s. The conditioned stimulus
becomes active 1 s into a trial, while the aversive stimulus becomes active from the 16th
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 16, along
with the period from all trials before the conditioned stimulus is presented. We assume
that the rate of spiking of the neuron is λA = 60 Hz. Region B consists of the period
from trials following trial 15 after the conditioned stimulus is presented. The rate of
spiking of the neuron in this region is λB = 20 Hz.
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 et

24

Trials

45
41
36
31
26
21
16
11
6
1

45
41
36
31
26
21
16

1.0
0.8
0.6
0.4
0.2
0

-1000 -500

2

1

0

500 1000

Time (ms)
(b)

0

Rate
(Hz)

Trial Effect
60
40
20
0
-1000

-500

0

500

1000

Time (ms)
(a)

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 ≈ 1 to ≈ 2. (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 1,
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)
al., 2017). Briefly, an observer mouse observes a demonstrator receive conditioned
stimulus-shock pairings through a perforated transparent divider. The experiment consists of 45 to 50 trials, divided into two phases. During the first 15 trials of the exper25

iment, termed the habituation period, both the observer and the demonstrator simply
hear an auditory conditioned stimulus. From the 16th 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 ≈ 400 ms
following conditioned stimulus presentation. Accounting for this increase in withintrial 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 ≈< 1 (indicative of suppression, as can be seen through the sparseness of the raster during trials
1 through 5) to a peak average value of ≈ 4 at trial 23. 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.

26

Trials

45
41
36
31
26
21
16
11
6
1

45
41
36
31
26
21
16

1.0
0.8
0.6
0.4
0.2
0
-500

0

500 1000 1500

Time (ms)
(b)

8 6 4 2 0

Rate
(Hz)

Trial Effect
15
10
5
0

-500

0

500

1000

1500

Time (ms)
(a)

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 ≈ 400 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 ≈< 1 (indicative of suppression, as can be seen through the sparseness of the raster during trials 1 through 5) to a
peak average value of ≈ 4 at trial 23. This increase, however, does not persist as in the
case of the simulated data (Figure 4).

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
27

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 ≈ 400 ms following conditioned
stimulus presentation) and that the extent of the condition is highest first between trials
18 and 24 and then between trials 31 and 41. Panel (a), and panel (b) in particular,
suggest that conditioning is accompanied by intricate dynamics in neural modulation.
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 ≈ 400 ms following conditioned stimulus presentation) and that the extent of the condition is highest first between trials 18 and 24 and
then between trials 31 and 41.
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 ≈ 1. 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.

28

Trials

50
46
41
36
31
26
21
16
11
6
1

50
46
41
36
31
26
21
16

1.0
0.8
0.6
0.4
0.2
0
-500

6

4

2

0

500 1000 1500

Time (ms)
(b)

0

Rate
(Hz)

Trial Effect
5
0
-500

0

500

1000

1500

Time (ms)
(a)

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 ≈ 1. This indicates that conditioning does not result
in a significant increase in spiking rate.

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

29

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.
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 16
through 45 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 5, 6, 10 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 5, 6, 10 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
30

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 = 0 ms).

Learning

Neuron in

Neuron in

Neuron in

Neuron in

Figure 5

Figure 6

Figure 10

Figure 11

617

1316

202

20

16

34

18

18

time (ms)
Learning
trial
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 617 ms after conditioned stimulus onset and at
trial 16. The cortical unit from Figure 6 only exhibits a significant change in neural
spiking 1316 ms following the conditioned stimulus and at trial 34. 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 617 ms after conditioned stimulus
onset and at trial 16. The cortical unit from Figure 6 only exhibits a significant change
31

in neural spiking 1316 ms following the conditioned stimulus and at trial 34. 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 et al., 2008).
We briefly describe the experiment here and refer the reader to (Wirth et al., 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.

32

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 {λ̂pk,r ∆}K,R
k=1,r=1 (in Hz) of the neuron at time trial r
and time k 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 condition-

Trials

ing period. That’s why, it does not make sense to generate plots such as Figure 5(b).
55
51
46
41
36
31
26
21
16
11
6
1

160
140
120
100
80
60
40
20
3 2 1 0

-300 0

Rate
(Hz)

Trial Effect

12001500

Time (ms)
(b)

60
40
20
0

-300

500

0

500

1200 1500

Time (ms)
(a)

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.

33

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 {λ̂pk,r ∆}K,R
k=1,r=1 (in Hz) of the hippocampal neuron at time trial r and time k 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 et al., 2012; Yuan
et al., 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,

34

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 3 to 4 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 statespace 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 timedependent 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 et al., 2013)

35

(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 Directors New Innovator Award DP2-DK-102256-01 (NIDDK).

36

Appendix
The SMuRF model cannot be converted easily to a standard statespace model
We focus on the simple case when ρx = ρz = 1, and αx = αz = 0



 x =x
+  ,  ∼ N (0, σ 2 ), k = 1, · · · , K
k

k−1

k

k





 zr = zr−1 + δr , δr ∼ N (0, σδ2 ), r = 1, · · · , R
Let t = (r − 1) × K + k, r = 1, · · · , R, k = 1, · · · , K. The index t 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
st = Ast−1 + vt ,

(28)

where st ∈ R2 . Let st,1 and st,2 denote the first and second components of st respectively. We ask that st ∈ R2 because of the two dimensions present in the SMuRF
model. Allowing the dimensionality of st to increase up to K 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 K + R, while the dimensionality of the state sequence in Equation 28 is 2 × (K × R). For there to be an equivalence, the sequence st
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!
37





xt−(d Kt e−1)×K 
 
 ∈ R2 . The quantity t gives the trial index r correLet st = 
K


zd t e
K
sponding time index t. The within-trial index corresponding to index t is then obtained
by substracting (r − 1) × K from
 t. 

 

x1 
x1 
 and sK+1 =  . In general, st,1 = st0 ,1 = xk0
Note, for instance, that s1 = 
 
 
r1
r2
for some 1 ≤ k0 ≤ K if and only if t > t0 s.t. t − t0 = p × K for some integer p,
where we assume without loss of generality that t > t0 . That is, two different indices t
and t0 share the same within-trial component if and only if they are apart by an integer
multiple of K. Stated otherwise, the first component of st exhibit circular symmetry!
Therefore, for Equation 28 to hold, st−1,1 must equal st,1 , which is not possible because
t and t − 1 are not apart by an integer multiple of p!
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 st to dimension K +1, 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 x and z in the Gibbs sampler.
Since wk,r |xk , zr is drawn from a PG distribution, we can write the log pdf of
wk,r |xk , zr as,

38

!



∞
x2
(2i+1)2
X
k wk,r
(2i
+
1)
xk + zr
−
−
2
+ log
(−1)i p
log p(wk,r |xk , zr ) = log cosh
e 8wk,r
2
2
2πw
k,r
i=1
(29)
The complete data likelihood of the SMuRF model is,
p(∆N, x, w; θ) = p(∆N|x, w)p(w|x)p(x)
=

K Y
R
Y
k=1 r=1

(30)

{p(∆Nkr |xk , zr )p(wk,r |xk , zr )}

K
Y
k=1

p(xk |xk−1 ; σ2 )

R
Y
r=1

p(zr |zr−1 ; σδ2 )
(31)

The log of the complete data likelihood is therefore,
log p(∆N, x, w; θ)
=
=

R
K X
X

{log p(∆Nk,r |xk , zr ) + log p(wk,r |xk , zr )} + π(x, z)

k=1 r=1
R 
K X
X


∆Nk,r log

k=1 r=1

exk +zr
1 + exk +zr




+ (1 − ∆Nk,r ) log

1
1 + exk +zr

(32)






+ log cosh

∞
X

xk + zr
2

!#
(xk +zr )2 wk,r
(2i+1)2
(2i
+
1)
−
−
2
+ log
e 8wk,r
+ π(x, z)
(33)
(−1)i p
2
2πw
k,r
i=1


K X
R 
X
1 + exk +zr
xk +zr
=
∆Nk,r (xk + zr ) − log(1 + e
) + log
xk +zr
2e 2
r=1
k=1
!#
∞
(2i+1)2
(xk +zr )2 wk,r X
(2i
+
1)
−
2
+ log e−
(−1)i p
e 8wk,r
+ π(x, z)
(34)
2
2πw
k,r
i=1
K X
R 
X
xk + zr (xk + zr )2 wk,r
=
∆Nk,r (xk + zr ) − log(2) −
−
2
2
k=1 r=1
!#

∞
K 
(2i+1)2
X
X
1
(xk − xk−1 )2
2
i (2i + 1) − 8wk,r
+ log
e
+
log(2πσ ) −
+ π(x, z)
(−1) p
2
2
2
2σ
2πw

k,r
i=1
k=1
(35)
= K log(2) + π(x, z) +

K X
R 
X
k=1 r=1

∆Nkr (xk

xk + zr (xk + zr )2 wk,r
+ zr ) −
−
2
2

39





+

K X
R
X

(
log

k=1 r=1

∞
X

2

(2i + 1) − (2i+1)
(−1)i p
e 8wk,r
2
2πwk,r
i=1

!)
(36)

where
K 
X
1

π(x, z) =

2

log(2πσ2 )

k=1
R 
X

+

r=1

(xk − ρx xk−1 − αx ux,k )2
−
2σ2



1
(zr − ρz zr−1 − αz uz,k )2
log(2πσδ2 ) −
2
2σδ2

(37)


From the complete data log likelihood, we see that

K X
R 
X
xk + zr (xk + zr )2 wk,r
r
log p(x|∆N, z, w; θ) ∝
∆Nk (xk + zr ) −
−
+ π(x, z)
2
2
k=1 r=1
(38)
∝−
where ∆Ñk,r =

∆Nk,r − 12
wk,r

R
K X
X
1 (∆Ñk,r − (xk + xr ))2
k=1 r=1

2

1/wk,r

+ π(x, z),

(39)

.

Therefore, we can rewrite the augmented model as a linear Gaussian state space model
as stated in Theorem 2.




xk = ρx xk−1 + αx ux,k + k , k ∼ N (0, σ2 )



R
−1 !


P
ỹk = xk + v˜k , v˜k ∼ N 0,
wk,r

r=1




R

 ỹ = ∆Ñ = x − K − P x w

k

k

k

2

r

k,r

r=1

r
Let Hk,r = ∆Ñ1r , . . . , ∆Ñk−1
denote the history of the observed process up-to and

including k − 1. We can now write
p(x|∆N, z, w; θ) ∝

K
Y
k=1

p(∆N |xk )p(xk |Hk )

(40)

The forward filtering equations for this linear Gaussian state space model are as follows.
xk|k−1 = ρx xk−1|k−1 + αx ux,k

(41)

2
2
σk|k−1
= ρ2x σk−1|k−1
+ σ2

(42)

40

xk|k = ρx xk|k−1
R
P
+
1+

(43)

wk,r

r=1
R
P
r=1

2
σk|k−1


2
wk,r σk|k−1

R
P



αx ux,k


+
 P
R
2
wk,r σk|k−1

r=1

R
2

∆Nkr

r=1

R
P



− −
zr wk,r

r=1
R

− ρx xk|k−1 

P
wk,r
r=1

(44)

2
=
σk|k


1+

2
σk|k−1

R
P
2
wk,r σk|k−1

(45)

r=1

2
After running the forward filtering algorithm, we obtain xK|K and σK|K
from the final
2
). Now we can treat xK
iteration of the filter. We can then draw xK ∼ N (xK|K , σK|K

as the new observations and use the Kalman filter again to draw samples for xK−1 , and
repeat this process iteratively for xK−1 , ..., x1 . The new observation equation reads,



 xk = ρx xk−1 + k , k ∼ N (0, σ2 )


 xk+1 = xk + k
From Bayes Rule we have,
p(xk |xk+1 , Hk ) =

p(xk+1 |xk )p(xk |Hk )
p(xk+1 |Hk )

(46)

Denote the densities of xk |xk+1 , Hk as
∗

2
xk |xk+1 , Hk ∼ N (xk|k ∗ , σk|k
)

Then the update equations are,
log p(xk |Hk ) ∝ log p(xk+1 |xk ) + log p(xk |Hk−1 )
xk|k ∗ = ρx xk|k−1 +

41

2
σk|k

σ2

(xk+1 − ρx xk|k−1 )

(47)
(48)

2 ∗
σk|k

=

2
σ2 σk|k−1

(49)

2
σ2 + σk|k−1
∗

2
With this backward-sampling algorithm, we can draw xk ∼ N (xk|k ∗ , σk|k
), where

i = K −1, ..., 1. The forward-filtering and backward-sampling algorithm are symmetric
for xk and zr .

Initialization of the EM algorithm and the Gibbs sampler
We initialize the Monte-Carlo EM algorithm with values for σ2 and σδ2 obtained by applying the one-dimensional state space model from (Smith & Brown, 2003) to the raster
data aggregated across either trials or time. We initialize the Gibbs sampler using trajectories drawn from posterior distribution of the state in the one-dimensional state-space
models (Smith & Brown, 2003) used to initialize σ2 and σδ2 . The Gibbs sampler draws
5000 samples for x, z, and w at every iteration. The algorithm reaches convergence
when the absolute change in σ2 and σδ2 is less than a certain threshold (10−5 ).

Ability of SMuRF model to identify learning time and trial in simulated data
The results of our analysis of the cortical data in Section 4 demonstrate that learning
of a contingency by a neuron is a dynamic process that cannot be easily quantified in
terms of a static time and trial of learning. We also demonstrated (Table 1) how to use
inferences from the SMuRF model to identify a learning time and a trial.
Here, we use simulated data to determine the ability of the SMuRF model to identify
learning time and trial when learning is accompanied by sustained changes in neural
42

spiking following conditioned stimulus onset and during conditioning. In particular, we
assess the sensitivity of our method to the extent of the change in neural spiking rate
following conditioned stimulus onset and during conditioning.
We simulated neural spike raster data in the same manner as described in the Simulation Studies component of our Applications section (Section 4). As in said section,
the raster is divided into two regions (Figure 3). We assume that the rate of spiking of
the neuron in Region B is fixed and equal λB = 20 Hz. We vary the rate of spiking
λA of the neuron in Region A from 20 to 45 Hz in 5 Hz increments. For each value of
λA , we simulated 10 independent rasters and determine the learning time and trial as in
Table 1. We use the average over the 10 rasters as the learning time/trial pair. When our
method detects no change, we declare the learning time and trial as the last time and
trail pair in the simulated data, i.e. 1000 ms and trial. Figures 8(a) and 8(b) show the
averages of the identified learning times and trials as a function of the ratio

λA
.
λB

The true

learning time is at 0 ms with respect to conditioned stimulus onset, and true learning
trial is trial 16. The figures demonstrate that the inference performed from the SMuRF
model is able to detect the true learning time and trial when the rate in Region A is 1.8
and 2 times larger than that in Region B. Moreover, the lower the ratio

λA
,
λB

the larger

the delay. The intuitive reason why it is easier to determine the learning trial is that, for
a given trial, there are many more observations, compared to the number of trials for a
give time instant.

43

Average Learning Trial

Average Learning Time (ms)

1000
800
600
400
200
0
1.0

1.2

1.4

1.6

1.8

2.0

2.2

45
40
35
30
25
20
15
1.0

1.2

1.4

1.6

λA
λB

λA
λB

(a)

(b)

1.8

2.0

2.2

Figure 8: Plot of average (a) learning time and (b) learning trial identified by SMuRF
inference as a function of ratio of neural spiking rate following and preceding learning.
The learning time and trial are defined as the first time and trial pair when then empirical probability (Equation 27) 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) is larger than 95%. The average is
taken over 10 independently simulated rasters for each of the values for the ration

λA
.
λB

The figure demonstrate the ability of inference performed using the SMuRF model to
reliably identify learning.

Robustness of SMuRF model to the presence of error trials
We simulated neural spike raster data in the same manner as described in the Simulation Studies component of our Applications section (Section 4). We picked three
consecutive trials, starting from trial 21, to be a error trials in which all of the observations were 0. Note that the data were simulated in the same manner as in Figure 4,
44

except for the presence of the error trials. Figure 9 shows the result of applying the
SMuRF model to these simulated raster data. The presence of the error trials does not

Trials

affect our remarks for Figure 4.
45
41
36
31
26
21
16
11
6
1

45
41
36
31
26
21
16

1.0
0.8
0.6
0.4
0.2
0

-1000 -500

3

2

1

0

500 1000

Time (ms)
(b)

0

Rate
(Hz)

Trial Effect
80
60
40
20
0

-1000

-500

0

500

1000

Time (ms)
(a)

Figure 9: (a) Simulated neural spike raster with three consecutive error trials, 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. Despite the presence of the error trials, the left panel
of the figure shows the cross-trial effect which, following conditioning, increases above
its average initial value of ≈ 1 to ≈ 2.

45

Figure 9: (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. Despite the presence of the error trials, we can see that with probability
close to 1, 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

Trials

(Figure 2 A) and the average rate at the same time (Figure 2 B).
45
41
36
31
26
21
16
11
6
1

45
41
36
31
26
21
16

1.0
0.8
0.6
0.4
0.2
0
-500

0

500 1000 1500

Time (ms)
(b)

6 4 2 0

Rate
(Hz)

Trial Effect
8
6
4
2
0

-500

0

500

1000

1500

Time (ms)
(a)

Figure 10: (a) Raster, along with SMuRF within and cross-trial components, from a
cortical neuron that exhibits a transient 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. The bottom panel suggests that this
neuron exhibits a delayed response to the conditioned stimulus, beginning at ≈< 100
ms following conditioned stimulus presentation. The response then decreases and is
followed by a slight increase at ≈ 700 ms.
46

Figure 10: 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 ≈ 1 to a peak average value of ≈ 2 between trial ≈ 18 at trial 30. (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 provides a more detailed account of the the intricate dynamics in neural modulation that accompanies conditioning

Trials

to the stimulus for this neuron.
45
41
36
31
26
21
16
11
6
1

45
41
36
31
26
21
16

1.0
0.8
0.6
0.4
0.2
0
-500

4

2

0

500 1000 1500

Time (ms)
(b)

0

Rate
(Hz)

Trial Effect
30
20
10
0
-500

0

500

1000

1500

Time (ms)
(a)

Figure 11: (a) Raster, along with SMuRF within and cross-trial components, from a
cortical neuron that exhibits a transient 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.

47

Figure 11: The bottom panel suggests that this neuron exhibits a delayed response to
the conditioned stimulus, beginning at ≈ 20 ms following conditioned stimulus presentation. The response is sustained until ≈ 800 ms, and then decreases. 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 ≈ 1 to a peak average value of ≈< 2 between trial ≈ 18 at trial 30. (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 provides a more detailed account of the the intricate dynamics
in neural modulation that accompanies conditioning to the stimulus for this neuron.

References
Allsop, S. A., Vander Weele, C. M., Wichmann, R., & Tye, K. M. (2014). Optogenetic
insights on the relationship between anxiety-related behaviors and social deficits.
Frontiers in behavioral neuroscience, 8, 241.
Allsop, S. A., Wichmann, R., Mills, F., Felix-Ortiz, A. C., Viennel, A., Beyeler, A.,
. . . Tye, K. M. (2017). Observational learning requires cortico-amygdala transfer of
socially-derived information. Submitted to Cell.
Andrews, D. F., & Mallows, C. L. (1974). Scale mixtures of normal distributions.
Journal of the Royal Statistical Society. Series B (Methodological), 99–102.
Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems.
Journal of the Royal Statistical Society. Series B (Methodological), 192–236.

48

Casella, G. (2001). Empirical bayes gibbs sampling. Biostatistics, 2(4), 485–500.
Casella, G., & George, E. I. (1992). Explaining the gibbs sampler. The American
Statistician, 46(3), 167–174.
Choi, H. M., & Hobert, J. P. (2013). The polya-gamma gibbs sampler for bayesian
logistic regression is uniformly ergodic. Electronic Journal of Statistics, 7, 2054–
2064.
Chopin, N., Jacob, P. E., & Papaspiliopoulos, O. (2013). Smc2: an efficient algorithm
for sequential analysis of state space models. Journal of the Royal Statistical Society:
Series B (Statistical Methodology), 75(3), 397–426.
Czanner, G., Eden, U. T., Wirth, S., Yanike, M., Suzuki, W. A., & Brown, E. N. (2008).
An analysis of between-trial and within-trial neural spiking dynamics. J. Neurophysiol., 99(5), 2672–2693.
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from
incomplete data via the em algorithm. Journal of the royal statistical society. Series
B (methodological), 1–38.
Frühwirth-Schnatter, S. (1994). Data augmentation and dynamic linear models. Journal
of time series analysis, 15(2), 183–202.
Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models
using pólya–gamma latent variables. Journal of the American statistical Association,
108(504), 1339–1349.

49

Rad, K. R., & Paninski, L. (2010). Efficient, adaptive estimation of two-dimensional
firing rate surfaces via gaussian process methods. Network: Computation in Neural
Systems, 21(3-4), 142–168.
Scott, J., & Pillow, J. W. (2012). Fully bayesian inference for neural models with
negative-binomial spiking. In Advances in neural information processing systems
(pp. 1898–1906).
Smith, A. C., & Brown, E. N. (2003). Estimating a state-space model from point
process observations. Neural Computation, 15(5), 965–991.
Truccolo, W., Eden, U. T., Fellows, M. R., Donoghue, J. P., & Brown, E. N. (2005). A
point-process framework for relating neural spiking activity to spiking history, neural
ensemble, and extrensic covariate effects. J. Neurophysiol., 93, 1074–1089.
Wirth, S., Yanike, M., Frank, L. M., Smith, A. C., Brown, E. N., & Suzuki, W. A.
(2003). Single neurons in the monkey hippocampus and learning of new associations.
Science, 300(5625), 1578–1581.
Yuan, K., Girolami, M., & Niranjan, M. (2012). Markov chain monte carlo methods
for state-space models with point process observations. Neural Computation, 24(6),
1462–1486.
Zammit-Mangion, A., Sanguinetti, G., & Kadirkamanathan, V. (2012). Variational estimation in spatiotemporal systems from continuous and point-process observations.
IEEE Transactions on Signal Processing, 7(60), 3449–3459.

50

