ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

Single-trial estimation of stimulus and spike-history
eﬀects on time-varying ensemble spiking activity ofmultiple neurons: a simulation study
Hideaki Shimazaki
RIKEN Brain Science Institute, Wako, Saitama, Japan
E-mail: shimazaki@brain.riken.jp
Abstract. Neurons in cortical circuits exhibit coordinated spiking activity, and can produce
correlated synchronous spikes during behavior and cognition. We recently developed amethod for estimating the dynamics of correlated ensemble activity by combining a modelof simultaneous neuronal interactions (e.g., a spin-glass model) with a state-space method(Shimazaki et al. 2012 PLoS Comput Biol 8 e1002385). This method allows us to estimate
stimulus-evoked dynamics of neuronal interactions which is reproducible in repeated trials underidentical experimental conditions. However, the method may not be suitable for detecting
stimulus responses if the neuronal dynamics exhibits significant variability across trials. In
addition, the previous model does not include eﬀects of past spiking activity of the neurons
on the current state of ensemble activity. In this study, we develop a parametric method for
simultaneously estimating the stimulus and spike-history eﬀects on the ensemble activity fromsingle-trial data even if the neurons exhibit dynamics that is largely unrelated to these eﬀects.
For this goal, we model ensemble neuronal activity as a latent process and include the stimulus
and spike-history eﬀects as exogenous inputs to the latent process. We develop an expectationmaximization algorithm that simultaneously achieves estimation of the latent process, stimulus
responses, and spike-history eﬀects. The proposed method is useful to analyze an interaction ofinternal cortical states and sensory evoked activity.

1. Introduction
Neurons in the brain make synaptic contacts to each other and form specific signaling networks.
A typical cortical neuron receives synaptic inputs from 3000 − 10000 other neurons, and makes
synaptic contacts to several thousands of other neurons. They send and receive signals using
pulsed electrical discharges known as action potentials, or spikes. Therefore, individual neurons
in a circuit can be activated in a coordinated manner when relevant information is processed. In
particular, nearly simultaneous spiking activity of multiple neurons (synchronous spikes) occurs
dynamically in relation to a stimulus presented to an animal, the animal’s behavior, and theinternal state of the brain (attention and expectation) [1–5].
Recently, it was reported that a model of synchronous spiking activity that accounts for spike
rates of individual neurons and interactions between pairs of neurons can explain ∼ 90% of thesynchronous spiking activity of a small subset (∼ 10) of retinal ganglion cells [6, 7] and corticalneurons [8] in vitro. This model is known as a maximum entropy model or an Ising/spin-glassmodel in statistical physics. However, since the model assumes stationary data, it is not directly
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distributionof this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Published under licence by IOP Publishing Ltd
1

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

applicable to non-stationary data recorded from awake behaving animals. In these data sets,
spike-rates of individual neurons and even interactions among them may vary across time.
In order to analyze time-dependent synchronous activity of neurons, we recently developed
a method for estimating the dynamics of correlations between neurons by combining the modelof neuronal interactions (e.g., the Ising/spin-glass model) with a state-space method [9, 10]. Inclassical neurophysiological experiments, neuronal activity is repeatedly recorded under identicalexperimental conditions in order to obtain reproducible features in the spiking activity acrossthe ‘trials’. Typically, neurophysiologists estimate average time-varying firing rates of individual
neurons in response to a stimulus from the repeated trials [11,12]. In the same fashion, the statespace method in [9, 10] aims to estimate the dynamics of the neuronal interactions, including
higher-order interactions, that occurs repeatedly upon the onset of externally triggered events.
When this method is applied to three neurons recorded simultaneously from the primary motor
cortex of a monkey engaged in a delayed motor task (data from [2]), it was revealed that theseneurons dynamically organized into a group characterized by the presence of a higher-order(triple-wise) interaction, depending on the behavioral demands to the monkey [12].
However, neurophysiological studies in the past decades revealed that spiking activity of
individual neurons is subject to large variability across trials due to structured ongoing activityof the networks that arises internally to the brain [13–15]. In these conditions, the method
developed in [12] would not eﬃciently detect the stimulus responses because a signal-to-noise
ratio may be small even in the trial-averaged activity. Although statistical methods for detecting
responses of individual neurons from single-trial data have been investigated [16–19], no methods
are available for estimating synchronous responses of multiple neurons to a stimulus in a single
trial when these neurons are subject to the activity that is largely unrelated to the stimulus.
In the analysis of single-trial data, it is critical to consider dependency of the current activityof neurons on the past history of their activity. A neuron undergoes an inactivation period
known as a refractory period after it generates an action potential. Therefore, a model neuron
significantly improves its goodness-of-fit to data if it captures this biophysical property [20, 21].
In addition, estimating the dependency of the current activity level of a neuron on past spiking
history of another neuron allows us to construct eﬀective connectivity of the network within an
observed set of neurons [22,23]. Including the spike-history eﬀects in the models of synchronous
ensemble activity is thus an important topic, and investigated also in [24] in the framework of
a continuous-time point process theory.
In this study we construct a method for simultaneously estimating the stimulus and spikehistory eﬀects on ensemble spiking activity when the activity of these neurons is dominatedby ongoing activity. For this goal, we extend the previously developed state-space modelof neuronal interactions: We model the ongoing activity, i.e., time-varying spike rates andinteractions, of neurons as a latent process, and include the stimulus and spike-history eﬀects on
the activity as exogenous inputs to the latent process. We develop an expectation-maximization(EM) algorithm for this model, which eﬃciently combines construction of a posterior densityof the latent process and estimation of the parameters for stimulus and spike-history eﬀects.
The method is tested using simulated spiking activity of 3 neurons with known underlyingarchitecture. We provide an approximation method for determining inclusion of these exogenous
inputs in the model and a surrogate method to test significance of the estimated parameters.
2. Methods
In this study, we analyze spike sequences simultaneously obtained from N neurons. Fromthese spike sequences, we construct binary spike patterns at discrete time steps by dividing
the sequences into disjoint time bins with an equal width of ∆ ms (in total, T bins). The width∆ determines a permissible range of synchronous activity of neurons in this analysis. We let Xit
be a binary variable of the i-th neuron (i = 1, 2, . . . , N ) in the t-th time bin (t = 1, 2, . . . , T ).

2

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

Here a time bin containing ‘1’ indicates that one or more spikes exist in the time bin whereas‘0’ indicates that no spike exists in the time bin. The binary pattern of N neurons at timet ]′ . The prime indicates the transposition operation to
bin t is denoted as Xt = [X1t , X2t , . . . , XN
the vector. The entire observation of the discretized ensemble spiking activity is represented as
X1:T = [X1 , X2 , . . . , XT ].
2.1. The model of time-varying simultaneous interactions of neurons
We analyze the ensemble spike patterns using time-dependent formulation of a joint probability
mass function for binary random variables. Let xi be a binary variable, namely xi = {0, 1}. The
joint probability mass function of N -tuple binary variables, x = [x1 , x2 , . . . , xN ], at time bin t(t = 1, 2, . . . , T ) can be written in an exponential form as
⎡
⎤
#
#
t
t
p(x|θ t ) = exp ⎣θit xi +θij
xi xj + · · · + θ1···N
x1 · · · xN − ψ(θ t )⎦ .
(1)
i

i<j

t , θt , . . . , θt′
Here θ t = [θ1t , θ2t , . . . , θ12131···N ] summarizes the time-dependent canonical parametersof the exponential family distribution. The canonical parameters for the interaction terms,t (i, j = 1, . . . , N ), represent time-dependent instantaneous interactions at time bin te.g., θijamong
the neurons denoted in its subscript. ψ(θ t ) is a log normalization parameter to satisfy&p(x|θ t ) = 1.
Using a feature vector that captures simultaneous spiking activities of subsets of the neurons,f = [f1 , f2 , . . . , f12 , f13 , . . . , f1···N ]′ , wherefi (x) = xi ,fij (x) = xi xj ,..
.

i = 1, · · · , N
i<jf1···N (x) = x1 · · · xN ,'(the probability mass function (Eq. 1) is compactly written as p(x|θ t ) = exp θ ′t f (x) − ψ(θ t ) .
The expected occurrence rates of simultaneous spikes of multiple neurons is given by a vector
η t = E [f (x) |θ t ], where expectation is performed using p(x|θ t ).
Eq. 1 specifies the probabilities of all 2N spike patterns by using 2N − 1 parameters.
One reasonable approach to reduce the number of parameters is to select and fix interesting
features in the spiking activity, and construct a probability model that maximizes entropy.
For example, maximization of entropy of the spike patterns given the low-order features,f = [f1 , f2 , . . . , fN , f12 , f13 , . . . , fN −1,N ]′ , yields a spin-glass model that is similar to Eq. 1, butdoes not include interactions higher than the second order. While it is important to explore acharacteristic feature vector to neuronal ensembles, here we note that the method developed inthis study does not depend on the choice of the vector, f . Below, we denote d as the number of
elements in the vector, f .
Given the observed ensemble spiking activity X1:T , the likelihood function of θ 1:T =[θ 1 , θ 2 , . . . , θ T ] is given asp (X1:T |θ 1:T ) =

T
)
t=1

exp[θ ′t f (Xt ) − ψ(θ t )],assuming conditional independence across the time bins.
equation of our state-space model.

3(2)

Eq.2 constitutes an observation

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

2.2. Inclusion of stimulus and spike-history eﬀects in the state model
The main focus of attention in this study is modeling of a process for the time-dependent
canonical parameters, θ t , in Eq. 1. We model their evolution as a first-order auto-regressive(AR) model. The eﬀects of the stimulus and spike history are included as exogenous inputs to
the AR model (an ARX model). In its full expression, the state model is written asθ t = Fθ t−1 + GSt +p#

Hi Xt−i + ξ t ,

(3)

i=1

for t = 2, . . . , T . Here the matrix F (d × d matrix) is the first order auto-regressive parameter.
ξ t (d × 1 matrix) is a random vector independently drawn from a zero-mean multivariate normal
distribution with covariance matrix Q (d × d matrix) at each time bin. The state process starts
with an initial value θ 1 that follows a normal distribution with mean µ (d × 1 matrix) andcovariance matrix Σ (d × d matrix), namely θ 1 ∼ N (µ, Σ). Below, we describe details of theexogenous terms.
The second term represents responses to external signals, or stimuli, St , which are observed
concurrently with the spike sequences. The vector St is a column vector of ns external signals
at time bin t. The each element is the value of an external signal at time bin t. If an externalsignal is represented as a sequence of discrete events, we denote the corresponding element of
St by ‘1’ if an event occurred within time bin t and ‘0’ otherwise. Multiplying St by the matrix
G (d × ns matrix) produces weighted linear combinations of the external signals at time bin t.
The third term represents the eﬀects of spiking activity during the previous p time bins, Xt−i(i = 1, . . . , p), on the current activity. The matrix Hi (d × N matrix) represents the spikehistory eﬀects of spiking activity in the previous time bin t − i on the state at time bin t. The
spike-history eﬀects are collectively denoted as H ≡ [H1 , H2 , . . . , Hp ] (d × N p matrix).
Eq. 3 constitutes a prior density of the latent process in our state-space model. We denote theset of parameters in the prior distribution, called hyper-parameters, as w ≡ [F, G, H, Q, µ, Σ].
In this study, we refer to w as a parameter. In addition, we simplify Eq. 3 as θ t =Fθ t−1 + Uut + ξ t , where ut is a single column vector constructed by stacking the stimulus vector
and spike-history vectors at time bin t in a row, i.e., ut = [St ; Xt−1 ; Xt−2 ; . . . ; Xt−p ]. Similarly,
we define a matrix U as U = [G, H]. With this simplification, the prior density defined in Eq. 3*
is written as p(θ 1:T |w) = p(θ 1 |µ, Σ) Tt=2 p(θ t |θ t−1 , F, U, Q), where the transition probability,p(θ t |θ t−1 , F, U, Q), is given as a normal distribution with mean Fθ t−1 + Uut and covariancematrix Q.
3. Estimation of stimulus responses and spike-history eﬀects
We estimate the parameter, w, based on the principle of maximizing a (log) marginal likelihood
function. Namely, we select the parameter that maximizes+
l (w) = log p (X1:T , θ 1:T |w) dθ 1:T .
(4)
For this goal, we use the expectation-maximization (EM) algorithm [25–27]. In this method,
we iteratively obtain the optimal parameter w∗ that maximizes the lower bound of the above
log marginal likelihood. This alternative function, known as the expected complete data log-

4

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

likelihood (a.k.a., q-function), is computed asq (w∗ |w) ≡ E [log p (X1:T , θ 1:T |w∗ ) |X1:T , w]
=

T
#
,
t=1

- d
1
Eθ ′t f (Xt ) − Eψ (θ t ) − log 2π − log det Σ∗221(T − 1) d(T − 1)− E[(θ 1 − µ∗ )′ Σ∗−1 (θ 1 − µ∗ )] −log 2π −log det Q∗
2
2
2
T
1#
−
E[(θ t − F∗ θ t−1 − U∗ ut )′ Q∗−1 (θ t − F∗ θ t−1 − U∗ ut )].
2

(5)

t=2

The expectation, E[!|X1:T , w], in Eq. 5 is performed using the smoother posterior density of thestate obtained by a nominal parameter, w, namelyp (θ 1:T |X1:T , w) =p (X1:T |θ 1:T ) p (θ 1:T |w).
p (X1:T |w)(6)

In particular, Eq. 5 can be computed using the following expected values by the posterior
density: The smoother mean θ t|T = E [θ t |X1:T , w], the smoother covariance matrix Wt|T =E[(θ t − θ t|T )(θ t − θ t|T )′ |X1:T , w], and the lag-one covariance matrix, Wt,t−1|T = E[(θ t −θ t|T )(θ t−1 − θ t−1|T )′ |X1:T , w]. These values are obtained using the approximate recursive
Bayesian filtering/smoothing algorithm developed in [9, 10] (See Appendix A and Eqs. A.5,
A.6, and A.7 therein).
In the EM-algorithm, we obtain the parameter that maximizes the q-function by alternating
the expectation (E) and maximization (M) steps. In the E-step, we obtain the above expectedvalues in Eq. 5 by the approximate recursive Bayesian filtering/smoothing algorithm using afixed w (Appendix A). In the M-step, we obtain the parameter, w∗ , that maximizes Eq. 5. The
resulting w∗ is then used in the next E-step. Below, we derive an algorithm for optimizing theparameter at the M-step.
For the state model that includes the auto-regressive parameter and stimulus and/or spike∂∗history eﬀects, these parameters are estimated simultaneously. From ∂F∗ q (w |w) = 0, we
obtain

F∗

T .
#
t=2

T
T .
/
/
#
#
Wt−1,t|T + θ t−1|T θ ′t−1|T + U∗
ut θ ′t−1|T =Wt−1,t|T + θ t|T θ ′t−1|T .
t=2

(7)

t=2

Here, θ t|T , Wt|T , and Wt−1,t|T are the smoother mean and covariance, and the lag-one covariance∂∗matrix given by Eqs. A.5, A.6, and A.7, respectively. Similarly, from ∂U∗ q (w |w) = 0, we obtain
F∗

T
#
t=2

θ t−1|T u′t + U∗

T
#
t=2

ut u′t =T#θ t|T u′t .

t=2

Hence, the simultaneous update rule for F∗ and U∗ is given as
/ &
1
' ∗
( 0 &T .
T′′
F U∗ =W+θθθut−1,t|T
t|T
t|T
t
t=2
t=2
t−1|T
/ &
2 & .
3−1
T
T
′
′
W
+
θ
θ
θ
u
t−1,t|T
t−1|T t−1|T
t=2
t=2 t−1|T t
.
&T
&T
′
′
t=2 ut θ t−1|T
t=2 ut ut
5

(8)

(9)

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

Here the inverse matrix on the r.h.s. is obtained by using the blockwise inversion formula:
,
-−1
,
-−1 3
4
5−1 2 −1
A B
A + A−1 B D − CA−1 B
CA−1 −A−1 B D − CA−1 B
,
-−1
,
-−1
=
.
C D
− D − CA−1 B
CA−1
D − CA−1 B∂∗
The covariance matrix, Q, can be optimized separately. From ∂Q∗ q (w |w) = 0, the update
rule of Q is obtained as
T

Q∗ =
+

1 #
[Wt|T − Wt−1,t|T F′ − FW′t−1,t|T + FWt−1|T F′ ]
T −1
1
T −1

t=2
T
#
t=2

,

θ t|T − Fθ t−1|T − Uut-,-′θ t|T − Fθ t−1|T − Uut .

(10)

∂
∗
Finally, the mean of the initial distribution is updated with µ∗ = θ 1|T from ∂µ∗ q (w |w) = 0.
The covariance matrix Σ for the initial parameters is fixed in this optimization.

4. Results
4.1. Simulation of a network of 3 neurons
In order to test the method, we simulate spiking activity of 3 neurons that possess
specific characteristics in spike generation and connectivity as follows (See Fig. 1A). (1) The
instantaneous firing rate of each neuron model depends on its own spike history in order toreproduce refractoriness in neuronal spiking activity. To achieve this, we adopt a renewal point
process model whose instantaneous inter-spike interval (ISI) distribution is given by an inverse
Gaussian distribution as a model of the stochastic spiking activity. (2) The firing rate of eachneuron model varies across time in order to reproduce the ongoing activity. For that purpose,
spike times of each neuron are generated from the renewal process by adding inhomogeneity
to the underlying rate using the time-rescaling method described in [20]. The underlying rateof the inhomogeneous renewal point process model is modulated using a sinusoidal function(frequency: 1 Hz; mean and amplitude: 30 spikes/s). This rate modulation is common to
the 3 neurons. (3) The neurons are activated by externally triggered stimulus inputs. To
realize the stimulus responses, we deterministically induce spikes at predetermined timings ofthe stimuli. We consider two stimuli, one (Stimulus 1) that induces a spike in Neuron 1, andthe other (Stimulus 2) that induces synchronous spikes in Neuron 2 and Neuron 3. The timings
of external stimuli are not related to the sinusoidal time-varying rate, but randomly selectedin the observation period (On average each stimulus happens once in 1 second). (4) Thereis feedforward circuitry in the network. We assume that Neuron 1 makes excitatory synaptic
contacts to Neuron 2 and Neuron 3. To realize this, 5ms after a spike occurs in Neuron 1, we
induce simultaneous spikes in Neuron 2 and Neuron 3 with a probability 0.5.
We simulate spike sequences with a length of 30 seconds using 1 ms resolution for numericaltime steps (An example of a short period (1 s) is shown in Fig. 1B). Figure 1C displays theinstantaneous spike-rates (conditional intensity functions of point processes) of Neuron 1 (Top,red line) and Neuron 2 & 3 (Bottom, green and blue lines) underlying the spiking activity in
Fig. 1B. The black lines indicate sinusoidal rate modulation common to all neurons. In addition,
spikes are induced in Neuron 1 at the onsets of Stimulus 1 (magenta triangles). Similarly,
simultaneous spikes of Neuron 2 and Neuron 3 are generated at the onsets of Stimulus 2 (cyantriangles). In the traces of instantaneous spike-rates in Fig. 1C, instantaneous increases causedby the stimuli and synaptic inputs are not displayed. The instantaneous spike-rate of a neuron
is reset to zero whenever a spike is induced in that neuron.

6

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

B

Stimulus 1

Neuron 1
Neuron 2A

Sinusoidal modulation

Neuron 3
10

Stimulus 2

10.5

11

Time [s]Stimulus 1

Neuron 2

5 ms delay

Neuron 3

Stimulus 2

Instantaneous rate (Neuron 1)
Sinusoidal modulation

100
50
Spike rate [event / s]

Neuron 1

C0
Stimulus 1
Instantaneous rate (Neuron 2)
Instantaneous rate (Neuron 3)

100
50
0
10

10.5
Stimulus 2

11Time [s]Figure 1. (A) Schematic diagram of a simulated network of 3 neurons. Neuron 1 makes
excitatory synaptic contacts to Neuron 2 and 3. Stimulus 1 excites Neuron 1 whereas Stimulus
2 excites Neuron 2 and 3 simultaneously. In addition, all neurons receive sinusoidal rate
modulation. (B) Simulated spiking activity of the network. A short period (1 s) of the total
30 s length is shown. The magenta and cyan triangles represent occurrence times of Stimulus
1 and 2, respectively. The gray bar highlights simultaneous spikes of Neuron 2 and 3 thatare causally induced 5 ms after a spike occurs in Neuron 1. (C) Instantaneous spike rates.
(Top) The red trace is the instantaneous spike rate of Neuron 1 simulated as an inhomogeneous
renewal point process whos instantaneousinter-spikeinterval is given by the inverse Gaussian
0
1
6 κ
2
κ
distribution (f (t; κ) = 2πt3 exp − 2t (t − 1) for x > 0, 0 for x < 0, κ = 1.8 for all neurons).
The inhomogeneous rate is modulated by the sinusoidal function (black solid line, frequency: 1
Hz; mean rate and amplitude: 30 spikes/s). (Bottom) Instantaneous spike rates of Neuron 2
and 3 (solid green line and dashed blue line, respectively).
4.2. Selection of a state model
We analyze the simulated ensemble activity by the proposed state-space model. For this goal,
we first construct binary spike patterns, X1:T , from the simulated spike sequences of 30 seconds(Note: spike times are recorded in 1 ms resolution) by discretizing them using disjoint time bins
with 2 ms width. We then apply state-space models to the binary data. The observation model
used here contains interactions up to the second order (a pairwise interaction model):
'
(
t
t
t
p(x|θ t ) = exp θ1t x1 + θ2t x2 + θ3t x3 + θ12x1 x2 + θ13
x1 x3 + θ23
x2 x3 − ψ(θ t ) .
(11)

For the state model, we consider 5 diﬀerent models that include a set of diﬀerent components in
Eq. 3. We select a model based on the framework of model selection in order to avoid over-fittingof a model to the data. Details of each state model are described as follows.
The first state model assumes F = I, where I is an identity matrix, and does not include any ofexogenous inputs. In this model, we optimize only the covariance matrix Q. The first model is7

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

denoted as [Q]. The second state model, denoted as [Q, F], is the first-order auto-regressivemodel. In this model, we optimize both the covariance matrix Q and the auto-regressive
parameter F. The third model, denoted as [Q, F, G], additionally includes the stimulus term
as exogenous inputs (Stimulus 1 and Stimulus 2). Both the matrix F and G are optimized
simultaneously in addition to Q. The fourth model includes both stimulus and spike-history
terms. In this model, the state model includes the history of spiking activity up to the last 6
time bins (p = 6). All parameters F, G, and H are optimized simultaneously in addition to
Q. This model is denoted as [Q, F, G, H6]. The structure of the fifth model is the same as the
fourth model, but contains the history of spiking activity up to the last 12 time bins (p = 12).
The last model is denoted as [Q, F, G, H12].
In order to select the most predictive model among them, we select the state-space modelthat minimizes the Akaike (Bayesian) information criterion (AIC) [28],AIC = −2l (w∗ ) + 2dim w∗ ,(12)

where w∗ is the optimized parameter in the Methods section. The (marginal) likelihood function
in Eq. 12 is obtained by a log-quadratic approximation, i.e, the Laplace method [10] (SeeAppendix B for the complete equation). Figure 2B displays decreases in AICs (∆AIC) of thelast four models from the AIC of the first model, [Q]. The larger the ∆AIC is, the better thestate-space model is expected to predict unseen data. For these data sets, inclusion of exogenous
inputs, in particular the spike history, significantly decreases the AIC. From this result, we selectthe state model that includes the stimulus response term and the spike-history terms up to theprevious 6 time bins.
4.3. Parameter estimation
We now look at the estimated parameters of the model selected by the AIC, namely[Q, F, G, H6]. Due to the limitation in the space, we do not display the estimated dynamics ofthe canonical parameters, θ t , by the recursive Bayesian method (See [9, 10] for the detailed
analysis on dynamics of θ t by this method). The estimated parameters, G and H, aresummarized in Fig. 3.
Here, in order to test the significance of the estimated parameters, we construct confidencebounds of the estimates, using a surrogate method. In this approach, we apply the same statespace model, [Q, F, G, H6], to the surrogate data set for the exogenous inputs. In the surrogate
data set, the onset times of external signals (Stimulus 1 and Stimulus 2) are randomized in the
1800

Figure 2. Comparison of state models by
the Akaike Information Criterion (AIC). The
state-space models with the following fivediﬀerent state models are comapared: [Q],
[Q, F],[Q,F,G], [Q,F, G, H6], and [Q,F,G, H12] (See details of the models for maintext). The reduction of the AIC of thelast four models from the AIC of the model[Q] (∆AIC) was repeatedly computed for 10
times. The height of the bar indicates theaverage ∆AIC. The error bar indicates ± 2
S.E. The numbers marked on each bar are
dimensions of the models (The number of freeparameters in the state model).

1600
1400

ΔAIC

1200
1000
800
600
400
200

105

117

225

333

0
[Q,F]

[Q,F,G] [Q,F,G,H6] [Q,F,G,H12]

Models8

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

Stimulus effects, G

A

Summed spike-history effects, H

B

Neuron 1 Neuron 2 Neuron 3

1.5
Stimulus 1
Stimulus 2

θ

0.5

(t)
1
(t)

1

θ2

(t)

θ3
0.5

0
(t)

θ 12
(t)

0

θ 13
(t)

θ 23

−0.5

−0.5
(t)

θ1

C

(t)

θ2

(t)

θ3

(t)

θ 12

(t)

θ 13

(t)

θ 23

D

(t)

Spike-history effects of Neuron i on θ i(Hp) ii 0.1(t)

Spike history effects on θ 23

(Hp) 6i
Neuron 1

0.4

0

Neuron 2

0.3−0.1

Neuron 3

0.2(t)

Neuron 1 to θ 1

−0.2

0.1

(t)

Neuron 2 to θ 20(t)

Neuron 3 to θ 3

−0.3

−0.1

−0.4

−0.2
2

4

6
8
lag [ms]

10

12

2

4

6
8
lag [ms]

10

12

Figure 3. Parameter estimation of the state-space model. (A) Eﬀects of Stimulus 1 and 2
on the canonical paramters (the first and second column of G). The vertical ticks on abscissa
indicate the 95% confidence bounds for each parameters obtained by the surrogate method. (B)Summed spike-history eﬀects. The matrices of spike-history eﬀects, Hp , are summed over thetime-lags and shown using color. (C) The eﬀect of a spike occurrence in Neuron i at p time(t)bins before the tth bin on θi (i = 1, 2, 3). (D) The eﬀect of a spike occurence in Neuron i(t)(i = 1, 2, 3) on the interaction parameter θ23 .
observation period. Similarly, we randomly select p = 6 bins from the past spiking activity toobtain surrogate spike history, instead of selecting the last consecutive 6 bins from time bint. Thus the estimated parameters, G and H, are not related to the structure specified in the
Section 4.1. We repeatedly applied the state-space model to the surrogate data (1000 times) toobtain the 95% confidence bound for the parameter estimation (vertical ticks in Fig. 3A, C andD).
Figure 3A displays eﬀects of the two stimuli, G, on the respective elements in θ t . First,(t)Stimulus 1 significantly increases θ1 whereas changes in the pairwise interactions by Stimulus
1 are relatively small, indicating that Stimulus 1 induces spikes in Neuron 1. On the contrary,
(t)
(t)
(t)
Stimulus 2 increases to θ2 , θ3 , and θ23 . In particular, the increase in the interaction parameter

9ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009(t)

θ23 by Stimulus 2 indicates that the presence of Stimulus 2 induces excess simultaneous spikes
in Neuron 2 and Neuron 3 more often than the chance coincidence expected
for the two neurons.
&
The spike-history eﬀects are summarized as a summed matrix, pi=1 Hi , shown in Fig. 3B.
(t)
Two major eﬀects are observed. First, the spike history of Neuron i significantly decreases θi(i = 1, 2, 3) (See diagonal of the first 3×3 matrix in Fig. 3B). Figure 3C displays the contribution(t)of a spike in Neuron i during the previous p time bins to the parameter θi . These components
primarily, albeit not exclusively, capture the renewal property of the simulated neuron models.
(t)
(t)
(t)
Second, the spike history of Neuron 1 increases θ2 , θ3 , and θ23 (See the first column in
Fig. 3B), indicating that spike interactions from Neuron 1 to Neuron 2 & 3. In particular, the(t)increase in θ23 due to the spikes in Neuron 1 during previous 1-3 time bins (Fig. 3D) indicatesthat the inputs from Neuron 1 induces excess synchronous spikes in the other two neurons withapproximately 2-6 ms delay.
5. Conclusion
We developed a parametric method for estimating stimulus responses and spike-history eﬀects
on the simultaneous spiking activity of multiple neurons when the ensemble themselves exhibit
ongoing activity. The method was tested by simulated multiple neuronal spiking activity withknown underlying architecture. We provided two methods to corroborate the fitted models.
First, based on the result in the preceding paper, we provided an approximate equation for thelog marginal likelihood (see Appendix B), which was used to select the most predictive statespace model. Second, we provided a method for obtaining confidence bounds of the estimatedparameters based on a surrogate approach.
Example spike sequences simulated in this study are overly simplified. Therefore, the methodneeds be tested using real neuronal spike data, e.g., from cultured neurons whose underlying
circuit is identified by electrophysiological studies. In practical applications, it is recommended
to utilize basis functions such as raised cosine bumps used in [23] in the exogenous terms inorder to capture the stimulus and spike-history eﬀects with a fewer parameters. In addition, an
appropriate bin size must be selected in order to obtain a meaningful result in the analysis of realdata. Since the bin size determines a permissible range of synchronous activity, a physiologicalinterpretation of the result depends on the choice of the bin size. It is thus recommended topresent results based on multiple diﬀerent bin sizes in order to confirm a specific hypothesis
in a study as shown in [10, 29]. Methods to overcome an artifact due to the disjoint binningare discussed in [30–34]. In future, inclusion of such advanced methods will allow us to detect
near-synchronous responses without sacrificing temporal resolution of the analysis.
Given that applicability of the method is confirmed in real data, the proposed method is useful
to investigate how ensemble activity of multiple neurons in a local circuit changes configurationsof their simultaneous responses (synchronous responses) to diﬀerent stimuli applied to an animal.
Further, it would be interesting to see diﬀerent eﬀects of the same stimulus on the ensemble
activity when an animal undergoes diﬀerent cortical states.
Acknowledgments
The present study is based on the modeling framework developed in [9, 10]. The author
acknowledges Prof. Shun-ichi Amari, Prof. Emery N. Brown, and Prof. Sonja Grün for theirsupport in construction of the original model. The author also thanks to Dr. Christopher L.
Buckley and Dr. Erin Munro for critical reading of the manuscript.

10ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

Appendix A. Construction of a posterior density by the recursive Bayesian
filtering/smoothing algorithm
A posterior density of the time-varying θ t , which specifies the joint probability mass function of
spike patterns at time bin t, are obtained by a non-linear recursive Bayesian estimation method
developed in [9, 10]. The method allows us to find a maximum a posteriori (MAP) estimate ofθ t and its uncertainty, namely the most probable paths of time-varying canonical parameters θ t
and their confidence bounds given the observed simultaneous activity of multiple neurons. The
estimation procedure completes by a forward recursion to construct a filter posterior density
and then by a backward recursion to construct a smoother posterior density. In this approach,the posterior densities are approximated as a multivariate normal probability density function.
In the forward filtering step, we first compute mean and covariance of one step prediction
density:θ t|t−1 = Fθ t−1|t−1 + Uut ,′Wt|t−1 = FWt−1|t−1 F + Q.

(A.1)
(A.2)

Then, a mean vector and covariance matrix of the filter posterior density, which is approximatedas a normal density, is given as
θ t|t = θ t|t−1 + nWt|t−1 (yt − η t|t ),

−1
−1
Wt|t
= Wt|t−1+ nJt|t ,

(A.3)
(A.4)

where η t|t = E [f (x) |X1:t , w] is the simultaneous spike rates at time bin t expected from thejoint probability mass function, Eq. 1, specified by θ t|t . Thus Eq. A.3 is a non-linear equation.
We solve Eq. A.3 by a Newton-Raphson method. It can be shown that the solution is unique.
The matrix Jt|t is a Fisher information matrix of Eq. 1 evaluated at θ t|t .
Finally, we compute mean and covariance of a smoother posterior density as,θ t|T = θ t|t + At θ t+1|T − θ t+1|t ,
(A.5)
,
- ′
Wt|T = Wt|t + At Wt+1|T − Wt+1|t At .
(A.6)

−1
with At = Wt|t F′ Wt+1|t
for t = T, T − 1, . . . , 2, 1. Namely, we start computing Eqs. A.5 andA.6 in a backward manner, using θ T |T and WT |T obtained in the filtering method at the initialstep. The lag-one covariance smoother, Wt−1,t|T , is obtained using the method of De Jong andMackinnon [35]:7
Wt−1,t|T ≡ E[ (θ t−1 − θ t−1|T )(θ t − θ t|T )′ 7 y1:T ] = At−1 Wt|T .
(A.7)

Appendix B. Approximate marginal likelihood function
The approximated formula of the log marginal likelihood (Eq. 4) was obtained in [10] as
T
T
#
, ′
, -- 1 #,l(w) ≈
n yt θ t|t − ψ θ t|t +log det Wt|t − log det Wt|t−1
2
t=1

−

1
2

t=1

T
#
t=1

0

,
-,
-′ 1
−1
tr Wt|t−1
θ t|t − θ t|t−1 θ t|t − θ t|t−1 .

(B.1)

Here we briefly provide the derivation (See [10] for details). The log marginal likelihood is
written as
+
T
T
#
#
l(w) =
log p(yt |y1:t−1 , w) =log p(yt |θ t )p(θ t |y1:t−1 , w)dθ t .
(B.2)
t=1

t=1

11

ELC International Meeting on Inference, Computation, and Spin Glasses (ICSG2013)
IOP Publishing
Journal of Physics: Conference Series 473 (2013) 012009
doi:10.1088/1742-6596/473/1/012009

The integral in the above equation is approximated as+1
p(yt |θ t )p(θ t |y1:t−1 , w)dθ t = 8(2π)d |Wt|t−1 |+8exp [q(θ t )] dθ t ≈ 8(2π)d |Wt|t |(2π)d |Wt|t−1 |' , -(
exp q θ t|t ,(B.3),-′ −1 ,
where q(θ t ) = n (yt′ θ t − ψ (θ t )) − 12 θ t − θ t|t−1 Wt|t−1
θ t − θ t|t−1 . To obtain the second
approximate equality,
approximation: the integral in Eq. B.3 is given as
8 we used the Laplace' , -(
9d
exp [q (θ t )] dθ t ≈ (2π) |Wt|t | exp q θ t|t . Here we note that a solution of q(θ t ) = 0 is
equivalent to the filter MAP estimate, θ t|t . By applying Eq. B.3 to Eq. B.2, we obtain Eq. B.1.
References[1] Vaadia E, Haalman I, Abeles M, Bergman H, Prut Y, Slovin H and Aertsen A 1995 Nature 373 515–518[2] Riehle A, Grün S, Diesmann M and Aertsen A 1997 Science 278 1950–1953[3] Steinmetz P N, Roy A, Fitzgerald P J, Hsiao S S, Johnson K O and Niebur E 2000 Nature 404 187–190[4] Fujisawa S, Amarasingham A, Harrison M T and Buzsáki G 2008 Nat Neurosci 11 823–833[5] Ito H, Maldonado P E and Gray C M 2010 J Neurophysiol 104 3276–3292[6] Schneidman E, Berry M J, Segev R and Bialek W 2006 Nature 440 1007–1012[7] Shlens J, Field G D, Gauthier J L, Grivich M I, Petrusca D, Sher A, Litke A M and Chichilnisky E J 2006
J Neurosci 26 8254–8266[8] Tang A, Jackson D, Hobbs J, Chen W, Smith J L, Patel H, Prieto A, Petrusca D, Grivich M I, Sher A,
Hottowy P, Dabrowski W, Litke A M and Beggs J M 2008 J Neurosci 28 505–518[9] Shimazaki H, Amari S I, Brown E N and Grün S 2009 Proc. IEEE ICASSP2009 pp 3501–3504[10] Shimazaki H, Amari S I, Brown E N and Grün S 2012 PLoS Comput Biol 8 e1002385[11] Shimazaki H and Shinomoto S 2007 Neural Comput 19 1503–1527[12] Shimazaki H and Shinomoto S 2010 J Comput Neurosci 29 171–182[13] Arieli A, Sterkin A, Grinvald A and Aertsen A 1996 Science 273 1868–1871[14] Tsodyks M, Kenet T, Grinvald A and Arieli A 1999 Science 286 1943–1946[15] Kenet T, Bibitchkov D, Tsodyks M, Grinvald A and Arieli A 2003 Nature 425 954–956[16] Nawrot M, Aertsen A and Rotter S 1999 Journal of Neuroscience Methods 94 81–92[17] Cunningham J, Yu B, Shenoy K, Sahani M, Platt J, Koller D, Singer Y and Roweis S 2008 Advances in
Neural Information Processing Systems 20 329–336[18] Czanner G, Eden U T, Wirth S, Yanike M, Suzuki W A and Brown E N 2008 J Neurophysiol 99 2672–2693[19] Yu B M, Cunningham J P, Santhanam G, Ryu S I, Shenoy K V and Sahani M 2009 J Neurophysiol 102
614–635[20] Brown E N, Barbieri R, Ventura V, Kass R E and Frank L M 2001 Neural Comput 14 325–346[21] Barbieri R, Quirk M C, Frank L M, Wilson M A and Brown E N 2001 J Neurosci Methods 105 25–37[22] Truccolo W, Eden U T, Fellows M R, Donoghue J P and Brown E N 2005 J Neurophysiol 93 1074–1089[23] Pillow J W, Shlens J, Paninski L, Sher A, Litke A M, Chichilnisky E J and Simoncelli E P 2008 Nature 454
995–999[24] Kass R E, Kelly R C and Loh W L 2011 Ann Appl Stat 5 1262–1292[25] Dempster A P, Laird N M and Rubin D B 1977 J Roy Stat Soc B Met 39 1–38 ISSN 00359246[26] Shumway R and Stoﬀer D 1982 J Time Ser Anal 3 253–264[27] Smith A C and Brown E N 2003 Neural Comput 15 965–991[28] Akaike H 1980 Bayesian Statistics ed Bernardo J M, Groot M H D, Lindley D V and Smith A F M (Valencia,
Spain: University Press) pp 143–166[29] Riehle A, Grammont F, Diesmann M and Grün S 2000 J Physiol Paris 94 569–582[30] Grün S, Diesmann M, Grammont F, Riehle A and Aertsen A 1999 J Neurosci Methods 94 67–79[31] Pipa G, Wheeler D W, Singer W and Nikolić D 2008 J Comput Neurosci 25 64–88[32] Grün S 2009 J Neurophysiol 101 1126–1140[33] Hayashi T and Yoshida N 2005 Bernoulli 11 359–379[34] Chakraborti A, Toke I M, Patriarca M and Abergel F 2011 Quant Finance 11 991–1012[35] De Jong P and Mackinnon M J 1988 Biometrika 75 601–602

12