State-Space Analysis of Time-Varying Higher-Order SpikeCorrelation for Multiple Neural Spike Train Data
Hideaki Shimazaki1*, Shun-ichi Amari1, Emery N. Brown2,3,4, Sonja Grün1,5,6
1 RIKEN Brain Science Institute, Wako-shi, Saitama, Japan, 2 Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge,
Massachusetts, United States of America, 3 Department of Anesthesia and Critical Care, Massachusetts General Hospital, Boston, Massachusetts, United States of America,
4 Division of Health Sciences and Technology, Harvard Medical School/Massachusetts Institute of Technology, Cambridge, Massachusetts, United States of America,
5 Institute of Neuroscience and Medicine (INM-6), Computational and Systems Neuroscience, Research Center Jülich, Jülich, Germany, 6 Theoretical Systems
Neurobiology, RWTH Aachen University, Aachen, Germany

Abstract
Precise spike coordination between the spiking activities of multiple neurons is suggested as an indication of coordinatednetwork activity in active cell assemblies. Spike correlation analysis aims to identify such cooperative network activity bydetecting excess spike synchrony in simultaneously recorded multiple neural spike sequences. Cooperative activity is
expected to organize dynamically during behavior and cognition; therefore currently available analysis techniques must be
extended to enable the estimation of multiple time-varying spike interactions between neurons simultaneously. In
particular, new methods must take advantage of the simultaneous observations of multiple neurons by addressing theirhigher-order dependencies, which cannot be revealed by pairwise analyses alone. In this paper, we develop a method for
estimating time-varying spike interactions by means of a state-space analysis. Discretized parallel spike sequences aremodeled as multi-variate binary processes using a log-linear model that provides a well-defined measure of higher-order
spike correlation in an information geometry framework. We construct a recursive Bayesian filter/smoother for the
extraction of spike interaction parameters. This method can simultaneously estimate the dynamic pairwise spikeinteractions of multiple single neurons, thereby extending the Ising/spin-glass model analysis of multiple neural spike train
data to a nonstationary analysis. Furthermore, the method can estimate dynamic higher-order spike interactions. To validate
the inclusion of the higher-order terms in the model, we construct an approximation method to assess the goodness-of-fit
to spike data. In addition, we formulate a test method for the presence of higher-order spike correlation even in
nonstationary spike data, e.g., data from awake behaving animals. The utility of the proposed methods is tested using
simulated spike data with known underlying correlation dynamics. Finally, we apply the methods to neural spike data
simultaneously recorded from the motor cortex of an awake monkey and demonstrate that the higher-order spikecorrelation organizes dynamically in relation to a behavioral demand.
Citation: Shimazaki H, Amari S-i, Brown EN, Grün S (2012) State-Space Analysis of Time-Varying Higher-Order Spike Correlation for Multiple Neural Spike Train
Data. PLoS Comput Biol 8(3): e1002385. doi:10.1371/journal.pcbi.1002385
Editor: Olaf Sporns, Indiana University, United States of America
Received May 24, 2011; Accepted December 28, 2011; Published March 8, 2012
Copyright: ß 2012 Shimazaki et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported in part by JSPS Research Fellowships for Young Scientists (HS), R01 MH59733 (ENB), DP1 OD003646 (ENB), R01 MH071847(ENB), and RIKEN Strategic Programs for R&D (SG). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of themanuscript.
Competing Interests: The authors have declared that no competing interests exist.
* E-mail: shimazaki@brain.riken.jp

Over the years, various statistical tools have been developed toanalyze the dependency between neurons, with continuous
improvement in their applicability to neuronal experimental data
(see [24–26] for recent reviews). The cross-correlogram [27] was
the first analysis method for detecting the correlation between
pairs of neurons and focused on the detection of stationary
correlation. The joint-peri stimulus time histogram (JPSTH)
introduced by [11,28] is an extension of the cross-correlogramthat allows a time resolved analysis of the correlation dynamics
between a pair of neurons. This method relates the joint spiking
activity of two neurons to a trigger event, as was done in the peristimulus time histogram (PSTH) [29–31] for estimating the timedependent firing rate of a single neuron. The Unitary Event
analysis method [25,32,33] further extended the correlation
analysis to enable it to test the statistical dependencies between
multiple, nonstationary spike sequences against a null hypothesisof full independence among neurons. Staude et al. developed a test

Introduction
Precise spike coordination within the spiking activities of multiple
single neurons is discussed as an indication of coordinated network
activity in the form of cell assemblies [1] comprising neuronal
information processing. Possible theoretical mechanisms and conditions for generating and maintaining such precise spike coordination
have been proposed on the basis of neuronal network models [2–4].
The effect of synchronous spiking activities on downstream neurons
has been theoretically investigated and it was demonstrated that theseare more effective in generating output spikes [5]. Assembly activity
was hypothesized to organize dynamically as a result of sensory input
and/or in relation to behavioral context [6–10]. Supportive
experimental evidence was provided by findings of the presence ofexcess spike synchrony occurring dynamically in relation to stimuli[11–14], behavior [14–19], or internal states such as memory
retention, expectation, and attention [8,20–23].
PLoS Computational Biology | www.ploscompbiol.org

1

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

time bins, and are modeled as a multivariate Bernoulli distribution
using a multinomial logit link function. The dependencies among
the binary units are modeled in the generalized linear framework
by introducing undirected pairwise and higher-order interaction
terms for instantaneous, synchronous spike events. This statisticalmodel is referred to as the ‘log-linear model’ [44,45], or Ising/spin-glass model if the model contains only lower-order interactions. The latter is also referred to as the maximum entropy model
when these parameters are estimated under the maximum
likelihood principle. In contrast to the former, biologically-inspired
network-style modeling, the latter approach using a log-linear
model was motivated by the computational theory of artificial
neural networks originating from the stationary distribution of a
Boltzmann machine [46,47], which is in turn the stochastic
analogue of the Hopfield network model [48,49] for associative
memory.
The merit of the log-linear model is its ability to provide a welldefined measure of spike correlation. While the cross-correlogramand JPSTH provide a measure of the marginal correlation of two
neurons, these methods cannot distinguish direct pairwise
correlations from correlations that are indirectly induced throughother neurons. In contrast, a simultaneous pairwise analysis basedon the log-linear model (an analogue of the Ising/spin-glass
analysis in statistical mechanics) can sort out all of the pairdependencies of the observed neurons. A further merit of the loglinear model is that it can provide a measure of the ‘pure’ higherorder spike correlation, i.e., a state that can not be explained bylower-order interactions. Using the viewpoint of an informationgeometry framework [50], Amari et al. [44,45,51] demonstrated
that the higher-order spike correlations can be extracted from thehigher-order parameters of the log-linear model (a.k.a. the natural
or canonical parameters). The strengths of these parameters areinterpreted in relation to the lower-order parameters of the dual
orthogonal coordinates (a.k.a. the expectation parameters). The
information contained in the higher-order spike interactions of aparticular log-linear model can be extracted by measuring the
distance (e.g., the Kullback-Leibler divergence) between thehigher-order model and its projection to a lower-order model
space, i.e., a manifold spanned by the natural parameters whose
higher-order interaction terms are fixed at zero [44,52–54].
Recently, a log-linear model that considered only up to pairwiseinteractions (i.e., an Ising/spin-glass model) was proposed as amodel for parallel spike sequences. Its adequateness was shown by
the fact that the firing rates and pairwise interactions explained
more than *90% of the data [55–57]. However, Roudi et al. [54]demonstrated that the small contribution of higher-order correlations found from their measure based on the Kullback-Leibler
divergence could be an artifact caused by the small number ofneurons analyzed. Other studies have reported that higher-ordercorrelations are required to account for the dependencies between
parallel spike sequences [58,59], or for stimulus encoding [53,60].
In [60], they reported the existence of triple-wise spike correlationsin the spiking activity of the neurons in the visual cortex andshowed their stimulus dependent changes. It should be noted,
though, that these analyses assumed stationarity, both of thefiring rates of individual neurons and of their spike correlations.
This was possible because the authors restricted themselves todata recorded either from in vitro slices or from anesthetizedanimals. However, in order to assess the behavioral relevance of
pairwise and higher-order spike correlations in awake behaving
animals, it is necessary to appropriately correct for time-varying
firing rates within an experimental trial and provide an algorithm
that reliably estimates the time-varying spike correlations within
multiple neurons.

Author Summary
Nearly half a century ago, the Canadian psychologist D. O.
Hebb postulated the formation of assemblies of tightly
connected cells in cortical recurrent networks because ofchanges in synaptic weight (Hebb’s learning rule) byrepetitive sensory stimulation of the network. Consequently, the activation of such an assembly for processing
sensory or behavioral information is likely to be expressedby precisely coordinated spiking activities of the participating neurons. However, the available analysis techniques
for multiple parallel neural spike data do not allow us toreveal the detailed structure of transiently active assemblies as indicated by their dynamical pairwise and higherorder spike correlations. Here, we construct a state-spacemodel of dynamic spike interactions, and present a
recursive Bayesian method that makes it possible to trace
multiple neurons exhibiting such precisely coordinatedspiking activities in a time-varying manner. We also
formulate a hypothesis test of the underlying dynamic
spike correlation, which enables us to detect theassemblies activated in association with behavioral events.
Therefore, the proposed method can serve as a useful tool
to test Hebb’s cell assembly hypothesis.

method (CuBIC) that enables the detection of higher-order spikecorrelation by computing the cumulants of the bin-wise population
spike counts [34,35].
In the last decade, other model-based methods have beendeveloped that make it possible to capture the dependency among
spike sequences by direct statistical modeling of the parallel spikesequences. Two related approaches based on a generalized linear
framework are being extensively investigated. One models the
spiking activities of single neurons as a continuous-time point
process or as a discrete-time Bernoulli process. The point processintensities (instantaneous spike rates) or Bernoulli success probabilities of individual neurons are modeled in a generalized linear
manner using a log link function or a logit link function,
respectively [36–38]. The dependency among neurons is modeledby introducing coupling terms that incorporate the spike history ofother observed neurons into the instantaneous spike rate [38–40].
Recent development in causality analysis for point process data[41] makes it possible to perform formal statistical significance testsof the causal interactions in these models. Typically, the models
additionally include the covariate stimulus signals in order toinvestigate receptive field properties of neurons, i.e., the relations
between neural spiking activities and the known covariate signals.
However, they are not suitable for capturing instantaneous,synchronous spiking activities, which are likely to be induced by anunobserved external stimulus or a common input from anunobserved set of neurons. Recently, a model was proposed todissociate instantaneous synchrony from the spike-history dependencies; it additionally includes a common, non-spike driven latentsignal [42,43]. These models provide a concise description of themultiple neural spike train data by assuming independent spiking
activities across neurons conditional on these explanatory
variables. As a result, however, they do not aim to directly model
the joint distribution of instantaneous spiking activities of multiple
neurons.
In contrast, an alternative approach, which we will follow andextend in this paper, directly models the instantaneous, joint
spiking activities by treating the neuronal system as an ensemble
binary pattern generator. In this approach, parallel spikesequences are represented as binary events occurring in discretized
PLoS Computational Biology | www.ploscompbiol.org

2

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

We consider the presence of excess spike synchrony, inparticular the excess synchrony explained by higher-ordercorrelation, as an indicator of an active cell assembly. If some ofthe observed neurons are a subset of the neurons that comprise
an assembly, they are likely to exhibit nearly completely
synchronous spikes every time the assembly is activated. It may
be that such spike patterns are not explained by mere pairwise
correlations, but require higher-order correlations for explanationof their occurrence. One of the potential physiological mechanisms for higher-order correlated activity is a common input froma set of unobserved neurons to the assembly that includes theneurons under observation [25,61–63]. Such higher-orderactivity is transient in nature and expresses a momentarysnapshot of the neuronal dynamics. Thus, methods that arecapable of evaluating time-varying, higher-order spike correlations are crucial to test the hypothesis that biological neuronal
networks organize in dynamic cell assemblies for information
processing. However, many of the current approaches based on
the log-linear model [44,45,53,55,56,61,62,64,65] are notdesigned to capture their dynamics. Very recently two approaches were proposed for testing the presence of non-zero pairwise[66] and higher-order [67] correlations using a time-dependent
formulation of a log-linear model. In contrast to these methods,
the present paper aims to directly provide optimized estimates of
the individual time-varying interactions with confidence intervals.
This enables to identify short lasting, time-varying higher-ordercorrelation and thus to relate them to behaviorally relevant time
periods.
In this paper, we propose an approach to estimate the dynamic
assembly activities from multiple neural spike train data using a‘state-space log-linear’ framework. A state-space model offers ageneral framework for modeling time-dependent systems byrepresenting its parameters (states) as a Markov process. Brown
et al. [37] developed a recursive filtering algorithm for a point
process observation model that is applicable to neural spike train
data. Further, Smith and Brown [68] developed a paradigm forjoint state-space and parameter estimation for point process
observations using an expectation-maximization (EM) algorithm.
Since then, the algorithm has been continuously improved andwas successfully applied to experimental neuronal spike data fromvarious systems [38,69–71] (see [72] for a review). Here, we
extend this framework, and construct a multivariate state-spacemodel of multiple neural spike sequences by using the log-linear
model to follow the dynamics of the higher-order spikeinteractions. Note that we assume for this analysis typicalelectrophysiological experiments in which multiple neural spiketrain data are repeatedly collected under identical experimental
conditions (‘trials’). Thus, with the proposed method, we deal
with the within-trial nonstationarity of the spike data that is
expected in the recordings from awake behaving animals. We
assume, however, that dynamics of the spiking statistics withintrials, such as time-varying spike rates and higher-orderinteractions, are identical across the multiple experimental trials(across-trial stationarity).
To validate the necessity of including higher-order interactionsin the model, we provide a method for evaluating the goodness-offit of the state-space model to the observed parallel spike sequencesusing the Akaike information criterion [73]. We then formulate ahypothesis test for the presence of the latent, time-varying spikeinteraction parameters by combining the Bayesian modelcomparison method [74–76] with a surrogate method. The lattertest method provides us with a tool to detect assemblies that are
momentarily activated, e.g., in association with behavioral events.
We test the utility of these methods by applying them to simulated
PLoS Computational Biology | www.ploscompbiol.org

parallel spike sequences with known dependencies. Finally, weapply the methods to spike data of three neurons simultaneously
recorded from motor cortex of an awake monkey and demonstratethat a triple-wise spike correlation dynamically organizes inrelation to a behavioral demand.
The preliminary results were presented in the proceedings of theIEEE ICASSP meeting in 2009 [77], as well as in conference
abstracts (Shimazaki et al., Neuro08, SAND4, NIPS08WS,
Cosyne09, and CNS09).

Results
General formulation
Log-linear model of multiple neural spike sequences.

We consider an ensemble spike pattern of N neurons. The state ofeach neuron is represented by a binary random variable, Xi(i~1, . . . ,N) where ‘1’ denotes a spike occurrence and ‘0’ denotes
silence. An ensemble spike pattern of N neurons is represented bya vector, X~ðX1 ,X2 , . . . ,XN Þ, with the total number of possible
spike patterns equal to 2N . Let p(x), where x~ðx1 ,x2 , . . . ,xN Þ
and xi ~1 or 0 (i~1, . . . ,N), represent the joint probability massfunction of thePN-tuple binary random variables, X . Because of
the constraint p(x)~1, the probabilities of all the spike patterns
are specified using 2N {1 parameters. In information geometry[44,50], these 2N {1 parameters are viewed as ‘coordinates’ thatspan a manifold composed of the set of all probability
distributions, fp(x)g. In the following, we consider twocoordinate systems.
The logarithm of the probability mass function can be expandedaslog p(x)~Xihi xi zX

hij xi xj z    zh1N x1    xN {y(h), ð1Þ

ivj

where h~(h1 ,h2 , . . . ,h12 ,h13 , . . . ,h1N )’. In this study, a prime
indicates the transposition operation to a vector or matrix.
y(h)~{logP p(f0, . . . ,0g) is a log normalization parameter tosatisfyp(x)~1. The log normalization parameter, y(h), is afunction of h. Eq. 1 is referred to as the log-linear model. The
parameters hi , hij ,…, h1N are the natural parameters of an
exponential family distribution and form the ‘h-coordinates’. The
h-coordinates play a central role in this paper.
The other coordinates, called g-coordinates, can be constructedusing the following expectation parameters:
gi ~E ½Xi ,gij ~E Xi Xj ,..
.

i~1,    ,Nivj
ð2Þg1N ~E ½X1    XN :
These parameters represent the expected rates of joint spike
occurrences among the neurons indexed by the subscripts. We
denote the set of expectation parameters as a vector
g~(g1 ,g2 , . . . ,g12 ,g13 , . . . ,g1N )’.
Eqs. 1 and 2 can be compactly written using the following‘multi-indices’ notation. Let Vk be collections of all the k-element
subsets of a set with N elements (i.e., a k-subset of N elements):
V1 ~f1,2, . . . ,N g, V2 ~f12,13, . . .g, V3 ~f123,124, . . .g, etc.
We use the ‘multi-indices’ I[fV1 , . . . ,VN g to represent thenatural parameters and expectation parameters as hI and gI ,respectively. Similarly, let us denote the interaction terms in Eq. 1
as3

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlationfi ðxÞ~xi ,fij ðxÞ~xi xj ,..
.
f1N ðxÞ~x1    xN :do so, we parameterize
the natural parameters as h t ~0½ht1 , . . . ,ht12 , . . . ,ht1r , . . . in discrete time steps t(~1,2, . . . ,T).
We denote the dimension of h t for an rth-order model as
P
N. The time-dependent log-linear model up to thed~ rk~1kr-th order interactions (Er ) at time t is then defined asi~1,    ,Nivj
ð3Þpðxjht Þ~exp½

Using fI ðxÞ (I[fV1 , . P. . ,VN g), Eqs. 1 and 2 can now be compactly
written as log p(x)~ I hI fI ðxÞ{y(h) and gI ~E½fI (X), respectively.
The higher-order natural parameters in the log-linear model
represent the strength of the higher-order spike interactions.
Amari et al. [44,45,50,51] proved that the h- and g-coordinates
are dually ‘orthogonal’ coordinates and demonstrated that thenatural parameters that are greater than or equal to the rth-order,
fhI g (I[fVr , . . . ,VN g, rƒN), represent an excess or paucity ofhigher-order synchronous spikes in the fgI g (I[fVr , . . . ,VN g)coordinates. To understand the potential influence of the higherorder natural parameters on the higher-order joint spike event
rates, let us consider a log-linear model in which the parameters
higher than or equal to the rth-order vanish: hI ~0(I[fVr , . . . ,VN g). In the dual representation, the higher-orderjoint spike event rates, fgI g (I[fVr , . . . ,VN g), are chancecoincidences expected from the lower-order joint spike event
rates, fgI g (i[fV1 , . . . ,Vr{1 g). From this, it follows that the nonzero higher-order natural parameters that are greater than orequal to the rth-order of a full log-linear model represent theexcess or paucity of the higher-order joint spike event rates, fgI g(I[fVr , . . . ,VN g), in comparison to their chance rates expectedfrom the lower-order joint spike event rates, fgI g(I[fV1 , . . . ,Vr{1 g).
However, in this framework, the excess or scarce synchronous
spike events of r subset neurons reflected in gI (I[Vr ) are caused not
only by the non-zero rth-order interactions fhI g (I[Vr ), but alsoby all non-zero higher-order interactions fhI g (I[fVr , . . . ,VN g).
Therefore, one cannot extract the influence of pure rth-order spikeinteractions (rvN) on the higher-order synchrony rates unlessthe parameters higher than the rth-orders vanish, hI ~0(I[fVrz1 , . . . ,VN g). To formally extract the rth-order spikeinteractions, let Er denote a log-linear model whose parameters
higher than the rth-order are fixed at zero. By successively adding
higher-order terms into the model, we can consider a set of
log-linear models forming a hierarchical structure as
E1 5E2 5    5EN . Here, Er 5Erz1 denotes that Er is a submanifold of Erz1 in the h-coordinates. In Er , the set of non-zero
rth-order natural parameters explains the excess or paucity of rthorder synchronous spike events, gI (I[Vr ), in comparison to theirchance occurrence rates expected from the lower-order marginal
rates, gI (I[fV1 , . . . ,Vr{1 g). In a full model EN , the lastparameter in the log-linear model, h1N , represents the pure
Nth-order spike correlation among the observed N neurons. The
information geometry theory developed by Amari and others[44,45,50,51] provides a framework for illustrating the duality andorthogonal relations between the h- and g-coordinates and singlingout the higher-order correlations from these coordinates using
hierarchical models. In the Methods subsection ‘Mathematicalproperties of log-linear model’ we describe the known properties ofthe log-linear model utilized in this study, and give references[44,45,50,51] for further details on the information geometry
approach to spike correlations.

I[fV1 ,...,Vr g

htI fI ðxÞ{y(h t ):ð4Þ

The corresponding time-varying expectation parameters at time t
0

are written as gt ~½gt1 , . . . ,gt12 , . . . ,gt1r , . . . , where each element(I[fV1 , . . . ,Vr g) is given as
gtI ~E ½fI ðxÞjh t :

ð5Þ

In neurophysiological experiments, neuronal responses arerepeatedly obtained under identical experimental conditions(‘trials’). Thus, we here consider the parallel spike sequences
repeatedly recorded simultaneously from N neurons over n trials
and align them to a trigger event. We model these simultaneous
spike sequences as time-varying, multivariate binary processes byassuming that these parallel spike sequences are discretized in time
into T bins of bin-width D. Let Xt,l ~(X1t,l ,X2t,l , . . . ,XNt,l ) be theobserved N-tuple binary variables in the t-th bin of the l-th trial,
whereby a bin containing ‘1’ values indicates that one or more
spikes occurred in the bin whereas ‘0’ values indicate that no spike
occurred in the bin. We regard the observed spike pattern Xt,l attime t as a sample (from a total of n trials) from a joint probability
mass function. An efficient estimator of gtI is the observed spikesynchrony rate defined for the t-th bin as
ytI ~

n
 
1X
fI Xt,ln l~1

ð6Þ

for I[fV1 , . . . ,Vr g. The synchrony rates up to the r-th order,0

yt ~½yt1 , . . . ,yt12 , . . . ,yt1r , . . . , constitute a sufficient statistic forthe log-linear model up to the r-th order interaction.
Using Eqs. 4 and 6, the likelihood of the observed parallel spikesequences is given asnpðy1:T jq1:T Þ~

T

P P exp½

X
I[fV1 ,...,Vr gl~1 t~1
T~

P exp½n(

X
I[fV1 ,...,Vr g

t~1
T

~

0
t t

 
htI fI Xt,l {y(h t )htI ytI {y(h t ))ð7Þ

P exp½n(y h {y(h )),tt~1

with y1:T ~fy1 ,y2 , . . . ,yT g and h 1:T ~fh 1 ,h 2 , . . . ,h T g. Here, we
assumed that the observed spike patterns are conditionallyindependent across bins given the time-dependent naturalparameters, and that the samples across trials are independent.
Our prior assumption about the time-dependent naturalparameters is expressed by the following state equationh t ~Fh t{1 zjt ,ð8Þ

for t~2, . . . ,T. Matrix F (d|d matrix) contains the first order
autoregressive parameters. jt (d|1 matrix) is a random vector
drawn from a zero-mean multivariate normal distribution withcovariance matrix Q (d|d matrix). The initial values obey

State-space log-linear model of multiple neural spikesequences. To model the dynamics of the spike interactions, we

extend the log-linear model to a time-dependent formulation. To
PLoS Computational Biology | www.ploscompbiol.org

X

4

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlationh 1 *N ðm,SÞ, with m (d|1 matrix) being the mean and S (d|dmatrix) the covariance matrix. The parameters F, Q, m, and S are
called hyper-parameters. In the following, we denote the set of
hyper-parameters by w: w~½F,Q,m,S.
Given the likelihood (Eq. 7) and prior distribution (Eq. 8), we
aim at obtaining the posterior density
pðh 1:T jy1:T ,wÞ~pðy1:T jh 1:T Þpðh 1:T jwÞ,pðy1:T jwÞ

ð9Þusing Bayes’ theorem. The posterior density provides us with themost likely paths of the log-liner parameters given the spike data
(maximum a posteriori, MAP, estimates) as well as the uncertainty
in its estimation. The posterior density depends on the choice ofthe hyper-parameters, w. Here, the hyper-parameters, w, except
for the covariance matrix S for the initial parameters, areoptimized using the principle of maximizing the logarithm of themarginal likelihood (the denominator of Eq. 9, also referred to as
the evidence),l ðwÞ~logðpðy1:T ,h 1:T jwÞdh 1:T :ð10Þ

For non-Gaussian observation models such as Eq. 7, the exact
calculation of Eq. 10 is difficult (but see the approximate formulain the Methods section). Instead, we use the EM algorithm [68,78]
to efficiently combine the construction of the posterior density andthe optimization of the hyper-parameters under the maximum
likelihood principle. Using this algorithm, we iteratively construct
a posterior density with the given hyper-parameters (E-step), andthen use it to optimize the hyper-parameters (M-step). To obtain
the posterior density (Eq. 9) in the E-step, we develop a nonlinear
recursive Bayesian filter/smoother. The filter distribution issequentially constructed by combining the prediction distribution
for time t based on the state equation (Eq. 8) and the likelihood
function for the spike data at time t, Eq. 4. Figure 1 illustrates therecursive filtering process in model subspace Er . In combination
with a fixed-interval smoothing algorithm, we derive the smooth
posterior density (Eq. 9). The time-dependent log-linear parameters (natural parameters) are estimated as MAP estimates of the
optimized smooth posterior density. In the Methods section, pleaserefer to the subsection on ‘Bayesian estimation of dynamic spikeinteractions’ for the derivation of the optimization method for
hyper-parameters, along with the filtering/smoothing methods.
We summarize a method for estimating the dynamic spikeinteractions in Table 1.

Figure 1. Geometric view of recursive filtering in subspace Er .
Each point in this figure represents a probability distribution, p(x), of an
N-tuple binary variable, x~(x1 ,x2 , . . . ,xN ). The underlying timedependent model is represented by white circles in the space of EN .
The dashed lines indicate projections of the underlying models to themodel subspace, Er . The maximum a posteriori (MAP) estimates of theunderlying models projected on subspace Er were obtained recursively:
Starting from the MAP estimate at time t{1 (filter estimate, red circle),
the model at time t is predicted based on the prior knowledge of thestate transition, Eq. 8 (blue arrow, prediction; black cross, a predicteddistribution). The maximum likelihood estimate (MLE, black circle) for
the spike data at time t derived by Eq. 4 is expected to appear near the
projection point of the underlying model at time t in Er . The filter
distribution at time t is obtained by correcting the prediction by theobservation of data at time t (black arrow). The filter estimation at time t
is used for predicting the model at time tz1 and so on. This recursive
procedure allows us to retain past information while tracking theunderlying time-dependent model based on the current observation.
doi:10.1371/journal.pcbi.1002385.g001

section, we obtain the joint posterior density of the log-linear
parameters (Eq. 9). The posterior density provides, not only themost likely path of the log-liner parameters (MAP estimates), butalso the uncertainty in its estimation. The credible interval allows
us to examine whether the pairwise spike correlation is statisticallysignificant (but see the later section on ‘‘Testing spike correlationin nonstationary spike data’’ for the formal use of the joint
posterior density for testing the existence of the spike correlation in
behaviorally relevant time periods). Second, an EM algorithm
developed in the proposed method optimizes the smoothness ofthe estimated dynamics of the pairwise correlation (i.e.,
optimization of the hyper-parameters, w, in the state equation,
Eq. 8). By the automatic selection of the smoothness parameter, we
can avoid the problem of spurious modulation in the estimateddynamic spike correlation caused by local noise, or excessivesmoothing of the underlying modulation.
Figure 2A displays an application of our state-space method to 2
parallel spike sequences, Xit,l (i~1,2, t~1, . . . ,T, and l~1, . . . ,n),which are correlated in a time-varying fashion. The data are
generated as realizations from a time-dependent formulation of a

Application of state-space log-linear model to simulated
spike data
Estimation of time-varying pairwise spike interaction.

To demonstrate the utility of the developed methods for the
analysis of dynamic spike correlations, we first consider anonstationary pairwise spike correlation analysis. For this goal,
we apply the state-space method to two examples of simulated
spike data, with N~2 neurons. The dynamic spike correlation
between two neurons can be analyzed by conceptually simpler
histogram-based methods, e.g., a joint peri-stimulus time
histogram (JPSTH). However, even for the pair-analysis, theproposed method can be advantageous in the following two
aspects. First, the proposed method provides a credible interval (aBayesian analogue of a confidence interval). Using the recursive
Bayesian filtering/smoothing algorithm developed in the Methods
PLoS Computational Biology | www.ploscompbiol.org

5

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Table 1. Method for estimating dynamic spike interactions.

I

II

Preprocessing parallel spike data
(1)

Align parallel spike sequences from N neurons at the onset of external clock that repeated n times.

(2)

Construct binary sequences Xt,l (t~1, . . . ,T and l~1, . . . ,n) using T bins of width D from the spike timing data.

(3)

Select rð~1,2, . . . ,N Þ, the order of interactions included in the model. At each bin, compute the joint spike event rates up to the rth order yt , using Eq. 6.

Optimized estimation of time-varying log-linear parameters(1)

Initialize the hyper parameters: w~½F,Q,m,S.{(2)E-step: Apply the recursive Bayesian filter/smoother to obtain posterior densities.
(i)

Filtering: For t~1,2, . . . ,T, recursively obtain
the one-step prediction density, pðht jy1:t{1 ,wÞ&N h tjt{1 ,Wtjt{1 , using Eqs. 25 and 26,

the filter density, pðht jy1:t ,wÞ&N htjt ,Wtjt , using Eqs. 31 and 32.

(ii)

Smoothing: For t~T{1,T{2, . . . ,1, recursively obtain
the smooth density, pðh t jy1:T ,wÞ&N h tjT ,WtjT , using Eqs. 34 and 35.

(3)M-step: Optimize the hyper-parameters.

(4)

Repeat (2) and (3) until the iterations satisfy a predetermined convergence criterion.{

Update the hyper-parameters, F and Q, using Eqs. 38 and 39, and m~h 1jT .

{In this study, we initialized the hyper-parameters using F~I, Q~0:05I, m~0 and used a fixed diagonal covariance matrix for an initial density, as S~0:1I, unlessspecified otherwise in the main text or figure captions.
{In our algorithm, we computed the approximate log marginal likelihood, l(w), of the model using Eq. 45. We stopped the EM algorithm if the increment of the log
marginal likelihood was smaller than 0:1.
doi:10.1371/journal.pcbi.1002385.t001

despite the presence of time-varying rates of synchronous spike
events (Figure 2B right, bottom).

full log-linear model of 2 neurons (Figure 2A left, repeated trials:
n~50; duration: T~400 bins of width D). Here, the underlying
model parameters, ht1 , ht2 , and ht12 (Figure 2A right, dashed lines),
are designed so that the individual spike rates are constant(gt1 ~0:0384 ½spike=D and gt2 ~0:0194 ½spike=D), while the spikecorrelation between the two neurons, ht12 , varies in time, i.e.,
across bins (synchronous spike events caused by the timedependent correlation are marked as blue circles in Figure 2Aleft). While the bin-width, D, is an arbitrary value in this simulation
analysis, the bin-width typically selected in spike correlation
analyses is on the milli-second order. If the bin-width is D~1 ms,
the individual spike rates of simulated neurons 1 and 2 are 38.4 Hz
and 19.4 Hz, respectively. The correlation coefficient calculated
from these parallel spike sequences is 0.0763. These values arewithin the realistic range of values obtained from experimentally
recorded neuronal spike sequences. By applying the state-spacemethod to the parallel spike sequences, we obtain the smooth
posterior density of the log-linear parameters. The right panels in
Figure 2A display the MAP estimates of the log-linear parameters(solid lines). The analysis reveals the time-varying pairwiseinteraction between the two neurons (Figure 2A right, bottom).
The gray bands indicate 99% credible intervals from the posterior
density, Eq. 9. We used marginal posterior densities, p(htI jy1:T ,w)(I[fV1 ,V2 g), to display the credible intervals for the individual loglinear parameters. The variances of the individual marginaldensities were obtained from the diagonal of a covariance matrixof the smooth joint posterior density (Eq. 35).
Figure 2B shows an application of the method to parallel spikesequences of time-varying spike rates. Here, the underlying model
parameters are constructed so that the two parallel spike sequences
are independent (ht12 ~0 for t~1, . . . ,T), while the individual
spike rates vary in time (i.e., across the bins). The observation ofsynchronous spike events (Figure 2B left, blue circles) confirms thatchance spike coincidences frequently and trivially occur at higher
spike rates. The analysis based on our state-space method reveals
that virtually no spike correlation exists between the two neurons,

PLoS Computational Biology | www.ploscompbiol.org

Simultaneous estimation of time-varying pairwise spikeinteractions. In this subsection, we extend the pairwise

correlation analysis of 2 neurons to the simultaneous analysis ofmultiple pairwise interactions in the parallel spike sequences
obtained from more than 2 neurons.
Figure 3 demonstrates an application of our method tosimulated spike sequences of N~8 neurons. As an underlying
model, we construct a time-dependent log-linear model of 8
neurons with time-varying rates and pairwise interactions (htI ,I[V1 ,V2 , duration: T~500 bins). The higher-order log-linear
parameters are set to zero, i.e., no higher-order interactions areincluded in the model. Figure 3A displays snapshots of the
dynamics of the parameters of the individual spike rates, gtI (I[V1 ),and pairwise interactions, htI (I[V2 ), at t~100,200,300,
400 and 500. Figure 3B shows the parallel spike sequences (50out of 200 trials are displayed) simulated on the basis of this model.
The spikes involved in the pairwise, synchronous spike events
between any two of the neurons (in total: 28 pairs) aresuperimposed and marked with blue circles. Figure 3C displays
snapshots of the simultaneous MAP estimates of the pairwiseinteractions, htI (I[V2 ), of a log-linear model of 8 neurons applied
to the parallel spike train data. In addition, the spike rates were
estimated from the dual coordinates, i.e., gtI (I[V1 ). The results
demonstrate that the simultaneous estimation of time-varying,multiple pairwise interactions can be carried out by using a statespace log-linear model with up to pairwise interaction terms.
Estimation of time-varying triple-wise spike interaction.

Another important aspect of the proposed method is its ability toestimate time-varying higher-order spike interactions that cannot be
revealed by a pairwise analysis. To demonstrate this, we apply thestate-space log-linear model to N~3 parallel spike sequences byconsidering up to a triple-wise interaction (i.e., the full log-linear
model). Spike data (Figure 4A) are generated by a time-dependentlog-linear model (Figure 4C, dashed lines) repeatedly in n~100

6

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Figure 2. Estimation of pairwise interactions in two simulated parallel spike sequences. (A) Application of the state-space log-linear
model to parallel spike sequences with time-varying spike interaction. (Left) Based on a time-dependent formulation of the log-linear model (dashedlines in the right panels represent the model parameters), N~2 parallel spike sequences, Xt,l , are simulated repeatedly for n~50 trials (duration:
T~400 bins). The two panels show dot displays of the spike events of the variables, X1t,l or X2t,l (t~1, . . . ,T and l~1, . . . ,n). The observed
synchronous spike events across the two spike sequences within the same trials are marked by blue circles. (Right) Smoothed estimates of the loglinear parameters, h t (solid lines, red: pairwise interaction; blue and green: the first order), estimated from the data shown in the left panels. The graybands indicate the 99% credible interval from the posterior density of the log-linear parameters. The dashed lines are the underlying time-dependent
model parameters used for the generation of the spike sequences in the left panels. (B) Application of the state-space log-linear model to
independent parallel spike sequences with time-varying spike rates. Each panel retains the same presentation format as in A.
doi:10.1371/journal.pcbi.1002385.g002

order joint spike events, the automatic smoothness optimization
method selects hyper-parameters that make the trajectories of thehigher-order log-linear parameters stiff in order to avoid statisticalfluctuation caused by a local noise structure. Given the limited
number of trials available for data analyses, these observations
show the necessity of a method to validate inclusion of the higherorder interaction terms in the model.
Additionally, we observe that in the later period of spike data
(300–500 bins), the dynamics of the estimated triple-wise spikeinteraction do not follow the underlying trajectory faithfully: The
underlying trajectory falls on outside the 99% credible interval.
Similar results are sometimes observed when an autoregressive
parameter, F, in a state model is optimized (Eq. 8). In contrast,
when we replace the autoregressive parameter with an identity
matrix (i.e., F~I, where I is the identity matrix), the credibleintervals become larger. Therefore, such observations do not
typically occur. Thus, we also need a method for validating theinclusion of the autoregressive parameter in the state model using

trials. Figure 4C displays the MAP estimates (solid lines) of the loglinear parameters from the data shown in Figure 4A. Here, nonzero parameter ht123 represents a triple-wise spike correlation, i.e.,
excess synchronous spikes across the three neurons or absence ofsuch synchrony compared to the expectation if assuming pairwise
correlations. The gray band is the 99% credible interval from themarginal posterior density, p htI jy1:T ,w , for I[fV1 ,V2 ,V3 g.
The credible interval of the higher-order (triple-wise) log-linear
parameter in the bottom panel of Figure 4C appears to be larger
than those of the lower-order parameters. In general, the observed
frequency of simultaneous spike occurrences decreases as thenumber of neurons that join the synchronous spiking activities
increases (note that the marginal joint spike occurrence rate, Eq. 2,
is a non-increasing function with respect to the order ofinteraction, i.e., gI §gJ if the elements of I are included in J).
Thus, the estimation variance typically increases for the higherorder parameters, as seen in the bottom panel of Figure 4C.
Related to the above, because of the paucity of samples for higher-

PLoS Computational Biology | www.ploscompbiol.org

7

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Figure 3. Simultaneous estimation of pairwise interactions of 8 simulated neurons. (A) Snapshots of the underlying model parameters of atime-dependent log-linear model of N~8 neurons containing up to pairwise interactions (duration: T~500 bins) at t~100,200,300,400 and 500th
bins. No higher-order interactions are included in the model. Each node represents a single neuron. The strength of a pairwise interaction between
the ith and jth neurons, htij (i,j~1, . . . ,8), is expressed by the color as well as the thickness of the link between the neurons (see legend at the right ofpanel B). A red solid line indicates a positive pairwise interaction, whereas a blue dashed line represents a negative pairwise interaction. The
underlying spike rates of the individual neurons, gti (i~1, . . . ,8), are coded by the color of the nodes (see color bar to the right of panel A). (B) Dot
displays of the simulated parallel spike sequences of 8 neurons, Xt,l , sampled repeatedly for n~200 trials from the time-dependent log-linear model
shown in A. For better visibility, only the first 50 trials are displayed (l~1, . . . ,50). Synchronous spike events between any two neurons (28 pairs intotal) are marked by blue circles. (C) Pairwise analysis of the data illustrated in B (using all n~200 trials) assuming a pairwise model (r~2) of 8
neurons. The snapshots at the t~100,200,300,400 and 500th bins show smoothed estimates of the time-varying pairwise interactions, htI (I[V2 ), andthe spike rates, gtI (I[V1 ). For this estimation, we use S~10{4 I for the prior density of initial parameters. The scales are identical to the one in panelA.
doi:10.1371/journal.pcbi.1002385.g003(r~1, . . . ,N). While the inclusion of increasingly higher-orderinteraction terms in the model improves its accuracy whendescribing the probabilities of 2N spike patterns, the estimation of
the higher-order log-linear parameters of the model may sufferfrom large statistical fluctuations caused by the paucity ofsynchronous spikes in the data, leading to an erroneous estimationof such parameters. This problem is known as ‘over-fitting’ themodel to the data. An over-fitted model explains the observed
data, but loses its predictive ability for unseen data (e.g., spikesequences in a new trial under the same experimental conditions).
In this case, the exclusion of higher-order parameters from themodel may better explain the unseen data even if an underlying
spike generation process contains higher-order interactions. The

an objective criterion. Detailed analyses of these topics will begiven in the next section using the example of 3 simulated neurons
displayed in Figure 4. In the above example, the state-space loglinear model with an optimized F provides a better overall fits tothe spike data than the model using F~I, despite an inaccurate
representation in part of its estimation. However, for the purposeof testing the spike correlation in a particular period of spike data,
we recommend using an identity matrix as an autoregressive
parameter, i.e. F~I, in the state model.

Selection of state-space log-linear model
For a given number of neurons, N, we can construct state-spacelog-linear models that contain up to the rth-order interactions
PLoS Computational Biology | www.ploscompbiol.org

8

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Figure 4. Estimation of triple-wise interaction from simulated parallel spike sequences of 3 neurons. (A) Dot displays of the simulated
spike sequences, Xt,l , which are sampled repeatedly for n~100 trials from a time-dependent log-linear model containing time varying pairwise andtriple-wise interactions (duration: T~500 bins; see the dashed lines in C for the model parameters). Each of the 3 panels shows the spike events foreach of the 3 variables, X1t,l ,X2t,l ,X3t,l (t~1, . . . ,T and l~1, . . . ,n), as black dots. Synchronous spike events across the 3 neurons as detected inindividual trials are marked by blue circles. (B) Observed rates of joint spike events, ytI (I[fV2 ,V3 g). (Top) Observed rates of the synchronous spike
events between all possible pair constellations as specified by index I (I[V2 ). (Bottom) Observed rate of the synchronous spikes across all 3 neurons,
yt123 . (C) Smoothed estimates of the time-varying log-linear parameters, h t . The three panels depict the smoothed estimates (solid lines) of the loglinear parameters, htI , of the different orders (I[fV1 ,V2 ,V3 g), as obtained from the data shown in A and B (top and middle: the first and second orderlog-linear parameters; bottom: triple-wise spike interaction, ht123 ). The gray bands indicate the 99% credible interval of the marginal posterior
densities of the log-linear parameters. The dashed lines indicate the underlying time-dependent parameters used for the generation of the spikesequences.
doi:10.1371/journal.pcbi.1002385.g004

by using the AIC. We apply the model selection method to the
spike train data of N~3 simulated neurons. The data are
generated by the time-varying, full log-linear model that contains a
non-zero triple-wise interaction terms shown in Figure 4C (dashedlines). The AICs are computed for hierarchical state-space loglinear models, i.e., for models of interaction orders up to r~1,2,3.
To test the influence of the data sample size on the model
selection, we vary the number of trials, n, used to fit thehierarchical log-linear models. The results are shown in Table 2
for n~2,5,20,100,200. For a small number of trials (n~5), amodel without any interaction structure (r~1) is selected. For
larger numbers of trials, models with larger interaction orders are
selected. For n~100 and 200, the full log-linear model (r~3) isselected.
Below, we examine whether the AIC selected a model thatminimizes the KL risk function by directly computing its
approximation using the known underlying model parameters.
First, Table 3 shows how often a specific order is selected by the
AIC by repeatedly applying the method to different samples
generated from the same underlying log-linear parameters(Figure 4C, dashed lines). We examine two examples: One inwhich a sample is composed of n~5 trials (left) and the other ofn~100 trials (right). We repeatedly compute the AICs of statespace models of different orders (r~1,2,3) applied to 100 data
realizations (of the respective number of trials). We then count how

model that has this predictive ability by optimally resolving thebalance between goodness-of-fit to the observed data and themodel simplicity is obtained by maximizing the cross-validated
likelihood or minimizing the so-called information criterion. In
this section, we select a state-space model that minimizes the
Akaike information criterion (AIC) [73], which is given asð
AIC~{2log

pðy1:T ,h1:T jwÞdh 1:T z2k:

ð11ÞThe first term is the log marginal likelihood, as in Eq. 10. The
second term that includes k is a penalization term. The AIC uses
the number of free parameters in the marginal model (i.e., thenumber of free parameters in w) for k. Please see in the Methods
subsection ‘Selection of state-space model by information criteria’
for an approximation method to compute the marginal likelihood.
Selecting a model that minimizes the AIC is expected to beequivalent to selecting a model that minimizes the expected (oraverage) distance between the estimated model and unknown
underlying distribution that generated the data, where the‘distance’ measure used is the Kullback-Leibler (KL) divergence.
The expectation of the KL divergence is called the KL riskfunction.
Selection from hierarchical models. Here, we examine thevalidity of including higher-order interaction terms in the model
PLoS Computational Biology | www.ploscompbiol.org

9March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Table 2. AICs for different numbers of trials.

The number of trials

Table 4. Model orders selected by different information
criteria for different numbers of trials.

Model order
r~1

r~2

r~3

n~2

775.15*807.146841.96n~5

1957.4

1890.6*

1922

n~20

7857.27565.1*

The number of trials
n~2

n~5

n~20

n~100

n~200

AIC

1

2

2

3

3

7585.8

BIC

2

2

2

2

3

1

1

1

2

2

n~100

37791

36264

36231*

PDIO

n~200

75443

72366

72283*

KL-risk

1

2

3

3

3

MSE

1

1

3

3

3

This table displays the AICs of a state-space log-linear model with increasing
interaction orders: an independent model (r~1), pairwise model (r~2), and fullmodel (r~3), applied to simulated N~3 spike sequences with an increasing
number of trials, n, considered in the analysis. The spike data are identical to thatshown in Figure 4A. The asterisk indicates the model that minimizes the AIC.
doi:10.1371/journal.pcbi.1002385.t002

The state-space log-linear models of different orders (r~1,2,3) are applied tosamples of the N~3 spike sequences of n repeated trials generated from atime-dependent full log-linear model (indicated by the dashed lines in
Figure 4C). Three data-driven information criteria, AIC, BIC, and PDIO, arecomputed for the fitted state-space models of the different orders, r~1,2,3.
The count for the model order r that minimizes the respective criteria isdetermined by repeating the process for 100 repetitions as in Table 3. The mostfrequently selected model order, r, is displayed for each of the information
criteria and for the different numbers of trials, n. For comparison, we also show
the order of interactions that minimizes the KL risk function (KL-risk) and meansquared error (MSE). We approximated the KL-risk and MSE as follows. At eachbin, we compute the KL-divergence (Eq. 21), between a full underlying loglinear model of N~3 neurons and the estimated log-linear model whoseparameters are given by the MAP estimates of the r th-order model. The totalsum of the all KL-divergences from T~500 bins is used
Pas the distance
between the two (time-dependent) models: i.e., KL~ Tt~1 D½q(ht ),p(htjT ),where the function D½q,p is given in Eq. 21. ht represents the underlying loglinear parameters used to generate the data. htjT is its estimate from one
sample composed of n trials. The parameters higher than the r th-order that are
not included in the model are set to zero. The KL-risk function is estimated as
the average of the KL-divergences of 100 realizations of the spike data, eachcomposed of n trials. To obtainthe MSE, we first computed the sum of thePsquared errors (SE) : SE~ Tt~1 jjh t {h tjT jj2 , using one sample composed of ntrials. The MSE is then estimated as the average of the SEs over 100 samples.
doi:10.1371/journal.pcbi.1002385.t004

often a model of order r is selected by minimizing the AIC. For
comparison, the table includes the outcomes from other criteriasuch as the Bayesian information criterion (BIC) [79,80] and thepredictive divergence for indirect observation models (PDIO) [81],which are suggested for models containing latent variables. Please
see the Methods section for the details of these criteria. Next,
Table 4 displays the most frequently selected model (fromr = 1,2,3) by the various information criteria when they areapplied to 100 data realizations as a function of the number oftrials in each data set, n~2,5,20,100,200 (see Table 3 for theoutcomes of n~5 and n~100).
We compare the outcomes of the AIC and two otherinformation criteria with the model order that minimizes the KL
risk function (KL-risk). For this goal, we include in Table 4 themodel order that minimizes the KL-risk between the underlying
log-linear model (Figure 4C, dashed lines) and estimates of the rth-order model. In addition to the KL-risk, we calculate the meansquared error (MSE) between the underlying model parametersand corresponding estimates. Please see in the caption of Table 4
how we compute the KL-risk and MSE for this analysis. We findthat the KL-risk and MSE select the same model, except for thecase of n~5 trials. In comparison to the KL-risk, the BIC tends toselect models with an excessively higher-order of interaction (overfitting) for a small number of trials (n~2) and tends to choose

lower-order models for n~20 and 100. The PDIO mostly selects
lower-order models. In contrast, the AIC follows the selection of
the KL-risk minimization principle, except for n~20, where itshows a conservative choice compared to the KL-risk.
We repeat the same analysis for spike data generated from anunderlying model that contains up to pairwise interactions, but does
not contain the triple-wise interaction term. The purpose of thisanalysis is to show that the methods do not select models withexcessively higher orders of interaction than those actually containedin the data. To construct such an underlying model, we project thefull model (shown in Figure 4C) onto the subspace of a pairwise loglinear model, E2 . The projection model does not contain any triplewise correlation, while the 1st and 2nd order expectation parameters(gI for I[fV1 ,V2 g) are the same as those of the full model that wasused to generate the data in the analysis of Tables 2–4. Table 5
displays the most frequently selected model orders by the AIC, BIC,
PDIO, along with the selections by the KL-risk and MSE. We findthat the pairwise model is the most frequently selected model under
all of the criteria for the samples with a large number of trials(n~100,200). Under this condition, only the AIC among the otherdata-driven methods follows the KL-risk selection. These results lead
us to the conclusion that the AIC is a reliable measure to assess the
goodness of fit of the state-space log-linear model.
Selection of state transition model. In addition tovalidating the inclusion of the order of spike interactions, we
examine state models (Eq. 8) with different conditions for thehyper-parameters by the AIC, using as an example the same spiketrain data of 3 neurons with n~100 trials (displayed in Figure 4A).
The tested state models are (i) a time-independent model in which

Table 3. Models selected using different information criteria.

n~5 trials

n~100 trials

r~1

r~2

r~3

r~1

r~2

r~3

AIC

29

71*

0

0

3

97*

BIC

15

85*

0

0

99*

1

PDIO

46*

25

29

3

60*

37

The state-space log-linear models containing interactions up to the rth-order(r~1,2,3) are applied to the data from N~3 simultaneous spike sequences.
The spike data is generated from a time-dependent full log-linear model (seedashed lines in Figure 4C) repeatedly for n trials; either n~5 (left) or n~100(right) trials. For this data set, we compute three information criteria (AIC, BIC,
and PDIO) and find the order of the model that minimizes these information
criteria. We repeated the selection of the model order 100 times, using each of
the criteria and using the spike data that contains respective number of trials(n~5 or n~100). The count of the order of spike interactions that minimizes
the applied information criteria is increased by 1, and finally expressed as afrequency. The asterisk marks the most frequently selected model.
doi:10.1371/journal.pcbi.1002385.t003

PLoS Computational Biology | www.ploscompbiol.org

10

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlationnecessity of dealing with the within-trial nonstationarity that is
typically present in the data from awake behaving animals.
However, we know from other correlation analysis approaches
that, if not well corrected for, nonstationary spike data bear the
potential danger of generating false outcomes [25,83]. Here wedeal with the within-trial nonstationarity of the data by using thestate-space log-linear model while assuming identical dynamic
spiking statistics across trials (across-trial stationarity). In order tocorrectly detect the time-varying correlation structure within trials,
we apply to the state-space log-linear model a Bayesian modelcomparison method based on the Bayes factor (BF) [74–76], andcombine it with a surrogate approach. The BF is a likelihood ratio
for two different hypothetical models of latent signals, e.g., in ourapplication, different underlying spike correlation structures.
Using the BF, we determine which of the two spike correlation
models the spike data supports. By computing the BF for aparticular task period in a behavioral experiment, we can test
whether the assumed correlation structure appears in association
with the timing of the animal’s behavior. In the following, we
denote a specific task period of interest by the time period ½a,b.
In this study, the BF, B12 (ya:b ), is defined as the ratio of themarginal likelihoods of the observed spike patterns, ya:b , in thetime period ½a,b under different models, M1 or M2 , assumed forthe hidden state parameters,

Table 5. Model orders selected using different information
criteria for model without triple-wise spike interaction.

The number of trials
n~2

n~5

n~20

n~100

n~200

AIC

1

2

2

2

2

BIC

2

2

2

2

2

PDIO

1

1

1

2

2

KL-risk

1

2

2

2

2

MSE

1

2

2

2

2

Similar to Table 4, the state-space log-linear models (r~1,2,3) are repeatedly
applied to 100 samples of N~3 spike sequences with n trials. In contrast toTable 4, here a sample of spike sequences is generated from a pairwise (r~2)
time-dependent log-liner model of N~3 neurons. This underlying model is a
projection of the full log-liner model examined in Table 4 (dashed lines in
Figure 4C) to the pairwise model space. The frequencies of the r th order modelthat minimizes AIC, BIC, and PDIO are counted by repeatedly applying models
with different orders (r~1,2,3) to 100 samples (each with n trials). The mostfrequently selected models are displayed for the different criteria and for
different numbers of trials, n, as in Table 4. The rows for the KL-risk and MSE
display the models that minimize the estimates of the KL-risk and MSE.
doi:10.1371/journal.pcbi.1002385.t005

the hyper-parameters are fixed as Q~0 and F~I, where I is an
identity matrix; (ii) a random walk model (F~I), in which only thecovariance matrix, Q, is optimized via the EM algorithm; and (iii)
a 1st-order autoregressive model, in which Q and F are both
optimized. With these settings, the state-space log-linear model ofcase (i) becomes a stationary log-linear model, which has beenfrequently employed in spike data analyses, often in the form of a
maximum entropy model [55–57,82]. In contrast, the state-spacemodels of cases (ii) and (iii) are nonstationary because the joint
distribution of the spike pattern changes in time as a result of thetime-varying log-linear parameters. The interaction parameters
themselves are assumed to be nonstationary in the state model incase (ii), and can be either stationary or nonstationary in the statemodel in case (iii). To see this, let li be the eigenvalues of F, i.e.,the solutions for det(F{lI)~0. For a stationary process, theeigenvalues have to satisfy jli jv1, otherwise the process isnonstationary. This excludes the case F~I.
The AICs of the full model with the various state-equations
become progressively smaller for increasingly complex state-spacemodels (38525 for case (i), 36263 for case (ii), and 36231 for case(iii)). The fact that conditions (ii) and (iii) are better fitted to the data
than condition (i) confirms that the proposed method of timevarying spike correlation analysis performs better than an analysis
based on a stationary log-linear model for this data. In addition, thefact that condition (iii) better fits the data than condition (ii) supports
the use of autoregressive models (note: The estimation in Figure 4B
is done with state model (iii)). We find that the fitted F (to the data in
Figure 4A) yields an eigenvalue larger than 1, indicating that theunderlying state process (T~500 bins) is modeled as a nonstationary process. We consider that the nonstationary state model isselected because of the relatively short observation period, during
which a nonstationary trend for the parameters can appear even if
the state process is stationary in the long run.

B12 (ya:b )~ð12Þ

By successively conditioning the past, the BF is computed by themultiplication of the bin-by-bin one-step BF given at time t asPb

B12 (ya:b )~B (y jy). Here, the bin-by-bin BF at timet~a 12 t a:t{1t, B12 (yt jya:t{1 ), can be calculated as (see the Methods subsection,‘Bayesian model comparison method for detecting spike correlation’),Ð

S1 pðh t jy1:t ,wÞdh t

B12 (yt jy1:t{1 )~ Ð

S2 pðh t jy1:t ,wÞdh t

Ð
S

=Ð 1

pðh t jy1:t{1 ,wÞdh t

S2 pðh t jy1:t{1 ,wÞdh t

,

ð13Þ

where Si is the space of the interaction parameters, h t , for themodel, Mi (i~1,2). In Eq. 13, pðh t jy1:t ,wÞ is the filter density andpðh t jy1:t{1 ,wÞ is called the one-step prediction density, both ofwhich are obtained in the Bayesian recursive filtering algorithm
developed in the Methods section (cf. Eqs. 25, 26 and Eqs. 31, 32).
Therefore, the bin-by-bin BF at time t, B12 (yt jy1:t{1 ), is the ratioof the odds (of opposing models) found by observing the spike train
data up to time t (filter odds, the numerator in Eq. 13) to the odds
predicted from y1:t{1 without observing the data at time t(prediction odds, the denominator in Eq. 13). Thus, an unexpectedsynchronous spike pattern that significantly updates the filter odds
for the interaction parameters from their predicted odds gives rise
to a large absolute value for the BF. Because the posterior densities
are approximated as a multivariate normal distribution in ourfiltering algorithm, the BF at time t can be easily computed byusing normal distribution functions. Please see the subsection,‘Bayesian model comparison method for detecting spike correlation’, in the Methods section for the derivation of Eq. 13 anddetailed analysis of the BF.
The BF becomes larger than 1 if the data, ya:b , support modelM1 as opposed to M2 as an underlying spike correlation structure
and becomes smaller than 1 if the data support model M2 as
opposed to M1 . Alternatively, it is possible to use the logarithm of
the BF, known as the ‘weight of evidence’ [75] which becomes
positive if the data support model M1 as opposed to model M2

Test for presence of spike correlation in nonstationary
spike sequences
One of the goals of a time resolved analysis of spike correlation
is to discover dynamical changes in the correlated activities ofneurons that reflect the behavior of an animal. This implies the
PLoS Computational Biology | www.ploscompbiol.org

pðya:b jM1 Þ:pðya:b jM2 Þ

11March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

assumptions about the first and second order log-linear parameters;
thus, they are allowed to be real numbers. These parameters areintegrated out in Eq. 13. In this simulation study, we independently
applied a state-space log-linear model to each of the four periods
and computed the respective BFs. Figure 5B displays the groundtruth of the time-dependent triple-wise interaction parameter, ht123 ,of the full log-linear model. The gray areas indicate periods where
the triple-wise spike interaction, ht123 , is positive: Model M1 is true inthis period. In the remaining periods (white), the model of a negative
triple-wise spike interaction, M2 , is true.
We compare the observed BF for the tested correlation models
to the BFs computed from surrogate data sets resampled under anull hypothesis of no tested order of spike interaction. Specifically,
to test the existence of a triple-wise spike correlation, the surrogate
BFs are computed from resampled spikes generated from a model
containing no triple-wise spike interaction. For this, we first apply
a pairwise state-space log-linear model (r~2) to the observed spikedata to derive estimates of the time-varying spike rates anddynamic pair-interactions. Then we generate 1000 surrogate
samples of N~3 parallel spike sequences, where each sample iscomposed of n~100 trials, using the fitted pairwise state-spacelog-linear model. Thus, the resampled spike sequences exhibit thesame time-varying spike rates and pairwise correlations as theoriginal data, but the triplet coincidences occur on a chance level,
as expected by the individual spike rates and pairwise correlations
among the neurons. The surrogate BFs of a triple-wise spikecorrelation (using S1 : ht123 w0, S2 : ht123 ƒ0) are computed foreach surrogate data set by applying the full model (r~3).
Figure 5C shows the observed BF (red vertical line and red
triangle) and the cumulative distribution function (CDF) of thesurrogate BFs (solid black lines) for each of the four time periods,
I–IV. In each of the graphs, the gray area indicates the 95%
confidence interval derived from the distribution of 1000
respective surrogate BFs. If the observed BF falls outside theconfidence interval, the null hypothesis that no tested order of
spike interaction exists is rejected. Further, by the two-tailed test, if
the observed BF falls outside the confidence interval into the red
area of the distribution on the right (i.e., significantly positive BF
values), the data supports model M1 as an underlying correlation
structure, whereas if it falls into the blue area on the left, weconclude that it is a significantly negative BF, i.e., M2 is supported.
If the observed BF falls into the 95% confidence interval, the null
hypothesis that no tested order of spike correlation exists is not
rejected. We first look at the results for periods II and IV in
Figure 5C. In these periods, the null hypothesis of no triple-wise
spike correlation is rejected, and the observed BF correctlysuggests a positive triple-wise correlation model (M1 ), reflecting thefact that the underlying triple-wise interaction term is positive
throughout the two periods (see periods II and IV in Figure 5B). In
contrast, in periods I and III, where a negative triple-wisecorrelation model (M2 ) is true, the null hypothesis of no triple-wisecorrelation cannot be rejected.
It is possible that the synchronous spike events observed among
multiple neurons are detected as simultaneous increases in pairinteraction terms using the log-linear model of the second order(r~2), without having to take higher-order spike interactions into
consideration. Paradoxically, this issue becomes relevant when atriple-wise spike correlation is a dominant factor in the generationof synchronous spike events because the projection of a full model(r~3) that contains a positive triple-wise interaction term to apairwise model induces simultaneous increases in the pairinteraction terms in that projected model (r~2). In such a case,an analysis of the higher-order spike correlations might just add
redundant information about an animal’s behavior to lower-orderand negative in the opposite situation. Below, we display the
results for the BF in bit units (logarithm of the BF to base 2), i.e.,
the weight of evidence. By sequentially computing the bin-by-bin
BF, we can obtain the weight of evidence in a period ½a,b as
the summation of the local weight of evidence: log2 B12 (y1:T )~Pbt~a log2 B12 (yt jy1:t{1 ).
An intuitive interpretation of the BF values is provided in theliterature [74,76]. For example, in [76], a BF (weight of evidence)from 1.6 to 4.3 bit was interpreted as ‘positive’ evidence in favor ofM1 against M2 . Similarly, a BF from 4.3 to 7.2 bit was interpreted
as ‘strong’ evidence, and a BF larger than 7.2 bit was found to be‘very strong’ evidence in favor of M1 against M2 . While theclassical guidelines are useful in practical situations, they aredefined subjectively. Thus, in this study, in order to objectively
analyze the observed value of the BF, we combine the Bayesianmodel comparison method with a surrogate approach. In thissurrogate method, we test the significance of the observed BF forthe tested spike interactions by comparing it with the surrogate
BFs computed from the null-data generated by destroying only thetarget spike interactions while the other structures such as thetime-varying spike-rates and lower-order spike interactions are
kept intact.
The BF in a behaviorally relevant sub-interval ½a,b can becomputed from the optimized state-space log-linear model fitted tothe entire spike train data in ½0,T. Here, for the purpose of testing
spike correlation in the sub-interval, we recommend to use F~I inthe state model because the autoregressive parameters are
optimized for entire spike train data, which are not necessarily
optimal for the sub-interval. Similarly, a typical trial-based
experiment is characterized by discrete behavioral or behaviorally
relevant events, e.g. movement onset after a go signal or a cue
signal for trial start, etc. Thus, on top of the expected smooth timevarying change in the spike-rate and spike-correlation, suddentransitions may be expected in their temporal trajectories. Because
we use time-independent smoothing parameters (i.e., hyper
parameter Q in Eq. 8) that were optimized to entire data (seethe EM algorithm in the Methods section), such abrupt changes
may not be captured very well. This may cause a false detection or
failure in the detection of the spike correlation at the edge of a task
period. For such data, we suggest applying the Bayesian modelcomparison method to state-space models which are independently fitted to each of the task periods (or relatively smooth subintervals within each task period).
Detecting triple-wise spike correlation in simulated
nonstationary spike sequences. In this subsection, we

examine a method to test for the presence of a higher-order(triple-wise) spike correlation using simulated spike train data withknown spike interaction dynamics. In this context, we again use
the simulated spike data of N~3 neurons (of length 500 bins withn~100 trials) shown in Figure 4A, which was generated by thetime-varying model, as shown in Figure 4C (dashed lines). In thefollowing, we assume for this data an underlying experimental
protocol that can be segmented into behaviorally-relevant time
periods, I–IV, as indicated in Figure 5A, as for example a period
before the trial starts, a preparatory period, a movement period,
etc. (e.g., see [8,19] for a typical experimental protocol).
In each of the time periods, we compare opposing models for thehidden log-linear parameters of N~3 neurons, using a full model(r~3). In one hypothetical model, M1 , we assume that the triplewise interaction term is positive, S1 : ht123 w0. The time (bin) index t
expresses that the BF is computed under the same models for all of
the time steps, t, in the respective time period. In the other model
M2 , we assume that the triple-wise interaction term is smaller thanor equal to zero, S2 : ht123 ƒ0. Neither model makes specific
PLoS Computational Biology | www.ploscompbiol.org

12

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Figure 5. Detection of spike correlations and their relation to pseudo experimental protocol (simulation study). (A) Sketch of anassumed experimental time course composed of four epochs (I–IV), e.g., of different behavioral task conditions. Epochs I–III have a duration of
100 bins, period IV has a duration of 200 bins. (B) Time-varying triple-wise spike interaction parameter of the underlying model (cf. Figure 4C, ht123 )used for the simulation of spike data (N~3 neurons, n~100 trials) during the time course outlined in A. The gray areas indicate the time intervals inwhich the triple-wise interaction is positive: ht123 w0. (C) Hypothesis testing for a triple-wise spike correlation based on surrogates. In each time period(I–IV), we perform a test on the Bayes factor (BF) resulting from the original data. The observed BF (Eq. 12), marked by a red line and triangle, iscomputed as evidence of a positive triple-wise interaction, M1 : ht123 w0, as opposed to a zero or negative triple-wise interaction, M2 : ht123 ƒ0. We thencompare the ‘observed BF’ with cumulative distribution functions (CDFs, solid lines) for the BFs derived from surrogate data sets generated from amodel containing only up to pairwise interactions (r~2). In all of the CDFs, the gray area indicates the 95% confidence interval of the distribution. If
the observed BF falls into the lower tail of the distribution (blue area), M2 is supported; if it falls into the upper tail (red marked area), M1 is
supported. (D) Time courses of the underlying pair-interaction parameters of a pairwise log-linear model (r~2), htI (I[V2 ). These underlying
parameters are obtained by projecting the full underlying model (r~3) in Figure 4 to the pairwise model space, E2 . The gray areas indicate the time
periods in which all of the pairwise interactions are positive. (E) Similar tests as in C, but for the BF computed as evidence for the presence of
simultaneous positive pairwise interactions, M1 : VI[V2 ,htI w0, as opposed to the absence of such an assembly, M2 : AI[V2 ,htI ƒ0. (F) Compact
visualization of the test results from C and E. The colored bars show which hypotheses are supported in the different time segments (red and blue)
and where the null hypothesis cannot be rejected.
doi:10.1371/journal.pcbi.1002385.g005(Figure 4C dashed lines) to the probability space of the secondorder model (E2 ). Note that in period II, the pair-interaction terms
simultaneously increase because the effect of the triple-wiseinteraction in E3 is projected to this sub-space model (E2 ) (seeFigure 4C, middle panel, for comparison, and observe that there isno such increase in the pair-interaction terms of the full log-linear
model). The gray areas indicate the period where model M1(simultaneously positive pairwise interactions) is true. Figure 5E
displays the results for testing the simultaneous increases in thepair interactions. We again look at periods II and IV. In period II,
the null hypothesis of independent spike sequences is rejected, andmodel M1 is supported, reflecting the fact that the projected
pairwise interactions are positive for most of this time period. In
contrast, the same null hypothesis cannot be rejected in period IV
because neither model M1 nor M2 alone can support the data inthis period: Models M1 and M2 are true in the first and second
halves of period IV, respectively (see Figure 5D).

analyses. Thus, we repeat an analysis similar to the one above
under this alternative hypothesis, using a model that contains up to
pairwise interactions (r~2). In this test, the BF (Eq. 12) iscomputed using the following assumed models on the hidden loglinear parameters: M1 , all three neurons simultaneously exhibit
positive pairwise interactions (S1 : htI w0 for all I[V2 ); and M2 , atleast one of the pairwise interactions is not positive (S2 : htI ƒ0 for
at least one I[V2 ). The first order log-linear parameters are again
assumed to be real values. Surrogate data sets are obtained byresampling spike sequences from the first order state-space loglinear model (r~1) fitted to the data. The surrogate BFs for themodel of simultaneously positive pairwise interactions, as opposed
to its complementary model, are computed for each of thesurrogate data sets by applying a pairwise model (r~2). Figure 5D
displays the ground truth of the dynamics of the pairwiseinteraction terms, htI (I[V2 ), of the pairwise log-linear model.
These are obtained by projecting the full underlying model

PLoS Computational Biology | www.ploscompbiol.org

13

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

constructed binary sequences, Xt,l , from the spike times of theneurons by binning the data using a bin-width of D~3 ms.
Figure 6C displays the occurrence rates for the individual spiking
activities (upper panel) and the pair and triple synchronous spike
events (middle and bottom panels) detected in the bins with awidth of D~3 ms.
We apply a full log-linear model of N~3 neurons with up to a
triple-wise interaction term to the binary spike data shown in
Figure 6C. Figure 6D displays the estimated dynamics of theinteraction parameters in the log-linear model by using themethod summarized in Table 1. The thick bold lines indicate the
MAP estimates, i.e., the most probable paths, of the interaction
parameters. The gray bands and thin solid lines mark the 95%
credible interval (the Bayesian analogue of the confidence interval)
computed from a marginal posterior density. The last parameterof the log-linear model, ht123 , indicates the triple-wise interaction
among the three neurons in the MI. Strong positive (or negative)
value of this term indicates that the three neurons are triple-wisecorrelated, i.e. dependent in a manner that cannot be explained bypair correlations of the neurons. Specifically, a positive triple-wiseinteraction value means that the occurrence of the synchronous
events among the three neurons is more frequent than the chancelevel expected from the observed individual rates and pairwise
spike correlations among them.
The analysis of the MI neurons using the state-space log-linear
model reveals that the triple-wise interaction, ht123 , during the PP
gradually increases, with additional local peaks at the ESs that also
increase in heights towards the end of the PP (Figure 6D, bottom
panel). This result is consistent with the results found by Riehle
and collaborators for a subset of the neurons of the same data set(Figure 2 in [8]). In this previous study neuron 1 and 2 (neuron id2, 3) were analyzed for excess spike synchrony using the Unitary
Events analysis [32,33]. The analysis revealed that the two
neurons exhibit a modulation of significant excess spike synchrony,with peaks at the ESs (despite the first) and at the RS.[8,32,84].
That observation was interpreted as evidence that the neurons
cooperate to prepare for motor action and facilitate the efficiency
for movement execution [8,84]. Similarly, occurrences ofsynchronous spike events of the three neurons roughly coincidewith the expected events (ESs) (Figure 6C, bottom panel);accordingly, our result shows that the triple-wise interaction,
ht123 , is also locked to the ESs, however decreases at RS. Note thatthe triple-wise interaction is not only determined by the frequency
of synchronous spike events of all three neurons but is also
determined by frequencies of other observed spike patterns acrosstrials. However, for typical neural spike train data with low spike
rates and low rates of synchronous spiking of pairs of neurons, theoccurrence of synchronous spike events of three neurons cansignificantly increase the triple-wise interaction. In the subsection‘Bayesian model comparison method for detecting spike correlation’ in the Methods section we explore the contribution ofdifferent spike patterns in different spiking scenarios to theevidence of a triple-wise spike correlation in simulated data.
Figure 6E displays modulation of the interaction term, ht123 , using
the binary data constructed with bins of smaller (D~2 ms, upper
panel) and larger (D~5 ms, lower panel) widths. The top of eachpanel shows the timing of the synchronous spike events of all threeneurons. The sample sizes for synchronous spike events across all
three neurons for the 2 ms bin-width are greatly reduced fromthose observed for a bin-width of 3 ms. Thus estimated dynamics
is much less structured: The data prevents us from detecting
existence of a triple-wise spike correlation using the proposed
statistical test (see below for the application of the test method).
With larger bins of a width of 5 ms, the precise locking of the

Figure 5F summarizes the results for periods I–IV using bars incorresponding colors. The results of our tests reflect the dynamics
of the underlying parameters in Figures 5B and D, and we findthat in period IV, only the test for a triple-wise spike correlation
detected the presence of interactions among all three neurons in
the data. This is because, in this time period, the dynamics of theunderlying triple-wise interaction correlate well with the assumed
behavioral time table, whereas the dynamics of the simultaneous
pairwise interactions do not. In summary, it is possible that thehigher-order analysis allows us to discover the correlated activities
of multiple neurons associated with behavioral events, which may
not be revealed by pairwise analyses.
Detecting triple-wise spike correlation in neural spikedata from behaving monkey. Finally, we demonstrate theapplication of our method on simultaneous spike recordings fromthe primary motor cortex (MI) of an awake, behaving monkey.
The data were recorded by Alexa Riehle and her colleagues [8] totest the hypothesis that neuronal cooperativity is involved in theplanning of motor actions. Therefore, Riehle et al. designed abehavioral task in which different durations of preparationintervals were provided to the monkeys before they had toperform an arm movement and touch a target on a screen. The
detailed time table of the task is as follows (cf. Figure 6A). A trial
started by a signal indicating to the monkey that he may initiate
the task by pressing a button. After initiating the trial, the monkey
had to wait for 1000 ms until a preparatory signal (PS) was
delivered, indicating to the monkey that now the preparationinterval started. After appearance of a response signal (RS) at 600,
900, 1200, or 1500 ms (randomly selected with equal probability),
the monkey had to move his arm and touch the target position.
The reaction times of the monkey showed a dependence on the
duration of the preparatory period (PP): The longer the PPs theshorter the reaction times [8,19,84]. The conditional probability
for the occurrence of the RS increases during the longer PP, inparticular at the times when the RS may occur. The hypothesis forthe underlying network activity realizing reduced reaction times
with longer PPs was that the system is better prepared for therequested movement by increasingly enhancing the synchrony inthe relevant network to facilitate the response. Indeed, Riehle et al.
showed that cortical neurons modulate the degree of excess
synchronous spiking activities independently from their individual
firing rates in conjunction with the occurrence of the expected RSs[8,19,84]. An open question that could not be conclusively
answered in that study was if larger groups of neurons than pairs
coordinate their activity in relation to the expected events. Here,we approach this question by applying our newly developed
method and will indeed demonstrate the existence of higher-order(triple-wise) spike interaction in relation to the behavioral task.
To that end, we analyze data which were in part analyzed by
the Unitary Events analysis method (neuron id 2 and 3; shown in
Figure 2 in [8]). Figure 6B shows this particular data set thatconsists of three simultaneously recorded neurons (neuron id 2, 3,
and 5; Here we denote them as neuron 1, 2, and 3) observed
during trials (n~36) of the longest preparatory period of 1500 ms.
As in Riehle et al. (Figure 2 in [8]) we align the trials at the PS, andanalyze the data for the interaction parameters as a function oftime from 200 ms before the PS and 1800 ms after the PS. This
time segment is composed of three behavioral epochs: an interval
before the PS (200 ms), the PP with three expected signals (ESs, at
600, 900 and 1200 ms), and a period after the RS at 1500 ms
including the reaction time and partly the movement (reaction andmovement time, RT-MT) (see Figure 6A). The average firing ratesof the neurons during the PP are 29.4, 12.9, and 41.9 Hz forneurons 1, 2, and 3 (neuron id 2, 3, and 5) respectively. We
PLoS Computational Biology | www.ploscompbiol.org

14

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

Figure 6. Analysis of experimental spike data using the state-space log-linear model. (A) Experimental time table of delayed-response
hand movement task. The experiments were designed and conducted by Riehle and her colleagues [8,84]. During the preparatory period (PP,1500 ms) that starts with the preparatory signal (PS), the presentation of the response signal (RS) was expected at three distinct moments at 600, 900,
and 1200 ms (expected signals, ESs). Here the RS occurs finally after the longest possible delay of 1500 ms. After the RS, the requested movement
was executed (reaction and movement time, RT-MT). See [8,84] and the text for the detailed experimental protocol. (B) Dot displays of the spikesequences (duration: 2 s, sampling resolution 1 ms) of three neurons simultaneously recorded from the primary motor cortex (MI). The spike

PLoS Computational Biology | www.ploscompbiol.org

15

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

sequences are aligned at the onset of the PS (n~36 trials). The synchronous spike events across the 3 neurons detected in individual trials (detection
in bins of D~3 ms width) are marked by blue circles. (C) Observed rates of joint spike events, ytI (I[fV1 ,V2 ,V3 g). All of the events are detected in bins
with a width of D~3 ms. (Top) Observed rates of the spike occurrence of the individual neurons (I[V1 ). (Middle) Observed rates of the synchronous
spike events between two neurons specified by index I (I[V2 ). (Bottom) Observed rate of the synchronous spike events of all 3 neurons, yt123 . (D)Estimation of the time-varying log-linear parameters of 3 neurons, htI (I[fV1 ,V2 ,V3 g), from the binary data shown in C, according to the method
summarized in Table 1. The panels depict the MAP estimates of the log-linear interaction parameters (solid lines; from top to bottom, the first andsecond order log-linear parameters and a triple-wise spike interaction, ht123 ). The gray bands indicate the 95% credible interval. In this analysis, weused an identity matrix as an autoregressive parameter, i.e. F~I, in the state model. The covariance matrix of the initial parameter was fixed toS~0:05I. (E) The smoothed estimate of a triple-wise spike interaction, ht123 , is computed from binary data constructed using a bin-width of 2 ms
(Top) and a bin-width of 5 ms (Bottom). The top of each panel shows the timing of the synchronous spike events of all three neurons.
doi:10.1371/journal.pcbi.1002385.g006

spike correlation should be stronger in the late stage of the PP than
before. To test this idea, we divided the PP into earlier and laterstages of the PP (each with a 750 ms duration, see the bottom of
Figure 7A), and investigate whether the three MI neurons exhibit atriple-wise spike correlation in these periods. In addition, we selecta duration of 300 ms after the onset of the RS, during which theanimal starts to initiate the motor action (reaction and movement
time, RT-MT). The upper panels of Figure 7C display the results
obtained with bins of width D~3 ms. The weights of evidence forthe triple-wise spike correlation models are 7.16, 11.2, and 0.98 bitin the early PP, late PP, and RT-MT periods respectively: i.e.,‘very strong’ evidence is found at the late stage of the PP. The null
hypothesis is rejected in this late period of the PP whereas the samenull hypothesis is not rejected in the earlier period of PP or theperiod after the RS. In the late PP, the significantly large positive
BF indicates the existence of a positive triple-wise spikecorrelation. The bottom panels of Figure 7C display the same
analyses, but with a larger (D~5 ms) bin-width. We obtain thesame results as in the analysis with a bin-width of 3 ms. When weanalyze the data using a smaller bin-width (D~2 ms), the BFs are
not significant in any of the periods (not shown here) because of
the small samples of synchronous spikes. In addition, as in thesimulation study in the previous subsection, we tested whether theobserved synchronous spiking activities are explained merely bysimultaneous increases in the pairwise interaction terms of thesecond order log-linear model. We did not find any evidence forsuch simultaneous increases in the pairwise interactions from thedata for all three periods (bin-width D~3 ms) (not shown here).
Thus only by the application of the higher-order analysis we wereable to detect the task-dependent changes in the joint interactionsof all three neurons.
The critical assumptions made in this study are independence
and identical sampling across the trials (across-trial stationarity).
However, in the spike sequences of neuron 2 (neuron id 3) shownin Figure 6B, we notice an increase in the firing rates across trials.
Higher firing rates yield a higher chance rate for synchronous
spiking events. Thus, underestimation of the spike rate caused byaveraging across the entire trials might induce spurious estimationof higher-order spike interaction. Indeed, synchronous spike eventsof the three neurons are mostly observed in later trials (e.g., trials
19–36, see blue dots in Figure 6B) when the firing rates werehigher than in the earlier trials. Hence, we repeat the same
analysis by using only the latter half of the trials (trials 19–36).
Figure 7D displays the results. We still observe a significantly
positive BF in the late PP. We note that the estimated dynamics ofthe log-linear parameters using only the latter half of the trials are
not changed much from those using all trials, i.e., Figure 6D. In
the analysis using the first half of the trials (1–19, not shown), the
BF is not significant during the late PP. Kilavik et al. [19] reported
that ‘‘during practice, the temporal structure of synchrony wasshaped, with synchrony becoming stronger and more localized intime during late experimental sessions’’. We observe similar effects
in our triple-wise spike correlation analysis.

synchronous spike events across the three neurons to ESs is nolonger apparent. However, the state-space log-linear model reveals
a gradual increase in the strength of the triple-wise interaction
until the end of the PP.
To strengthen the findings that a triple-wise interaction termincreases during the PP, we test for the presence of a triple-wise
spike correlation in the PP using the Bayes factor (marginal
likelihood ratio, Eq. 12). We compute the BF for the opposing
models using a full model (r~3) of N~3 neurons for the hidden
log-linear parameters. In one hypothetical model, M1 , we assume
that the triple-wise interaction term is positive, S1 : ht123 w0, in theother model, M2 we assume that the triple-wise interaction term issmaller than or equal to zero, S2 : ht123 ƒ0. A large positive valueof the BF indicates that the spike data in that period supports
model M1 as opposed to model M2 , i.e., the existence of excess
synchronous spike events among the three neurons in comparison
with the chance rate given by spike rates and pairwise correlations.
A large negative value of the BF shows support of model M2 as
opposed to model M1 , i.e., a paucity of synchronous spikes for thethree neurons. Figure 7A shows the BF computed bin-by-bin (Eq.
13) in bit units. In the bottom panel of Figure 7A, we indicate thebehavioral periods for which we test for the presence of triple-wise
spike correlation (the PP and RT to MT; in later analysis we divide
the PP into early and late stages). The evidence for model M1 , as
opposed to M2 , in the behaviorally relevant time periods isobtained by summing the log of the bin-by-bin BF in that period(cf. Eqs. 12 and 13). The BF in the PP is found to be 18.08 bit,which is interpreted as ‘very strong’ evidence for presence of apositive triple-wise spike correlation by the classical guideline [76].
Next, we test the observed value of the BF in the PP (18.08 bit) bycomparing it with surrogate BFs. These are computed from null datathat contain no triple-wise spike interactions, while keeping the timevarying structure of the pairwise correlations and the individual
spike rates the same as those observed in the original spike train data.
The construction of the null data follows the same procedure as in
the simulation study shown before: We apply a pairwise state-spacelog-linear model (r~2) to the spike data and then generate 1000
surrogate samples (each with n~36 trials) of N~3 parallel spikesequences using the fitted pairwise state-space log-linear model.
Figure 7B (left panel) displays the surrogate distribution along withthe observed BF in the PP. The observed positive BF falls out of the
95% confidence interval, suggesting the presence of a positive triplewise interaction as an underlying model for the spike train data inthis period. In the right panels of Figure 7B, we display the results ofthe same analysis, but using spike train data binned by a larger(D~5 ms, top panel) and a smaller (D~2 ms, bottom panel) binwidth. We obtained almost the same results as in the analysis with abin-width of 3 ms as for the analysis with the larger bin-width (5 ms).
However, we could not reach the same results for the smaller bin
width (2 ms) because the small count of synchronous spike events
made it impossible to reject the null hypothesis.
If the higher-order dependency among the three neurons isrelated to motion preparation, the evidence for the triple-wise
PLoS Computational Biology | www.ploscompbiol.org

16

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

triple-wise spike interaction, M1 : ht123 w0, as opposed to a zero or
negative triple-wise interaction, M2 : ht123 ƒ0. The evidence for modelM1 as opposed to M2 in a behaviorally relevant time period is obtainedby summing the log of the bin-by-bin BF in that period (cf. Eq. 12).
(Bottom) Two behavioral periods (preparatory period, PP, and reaction
and movement time, RT-MT) tested for presence of a triple-wise spikecorrelation. In addition, to examine the evolution of the triple-wise
spike correlation in the PP, the PP is divided into early and late stages at
the middle of the PP. (B) (Left) The observed BF for an entire PP (marked
by a red line and triangle) computed using Eq. 12 (bin-width: 3 ms). We
then test the ‘observed BF’ using a distribution function (solid line) ofnull BFs derived from surrogate data sets generated from a fitted model
containing only up to pairwise interaction terms (r~2). The gray area
indicates the 95% confidence interval of the distribution. Verticaldashed lines indicate the 90% confidence interval. If the observed BF
falls into the lower tail of the distribution (blue area), M2 is supported, if
it falls into the upper tail (red marked area), M1 is supported. (Right)
The same analysis in the left panel, but with binary data constructedusing larger (5 ms, upper panel) and smaller (2 ms, lower panel) binwidths. (C) (Top) Test of observed BFs (bin-width 3 ms) computed in
distinct periods: early PP (0–750 ms), late PP (750–1500 ms), and RT-MT(1500–1800 ms). (Bottom) The same test as in the top panel, except thatbinary data using a bin-width of 5 ms were used. (D) Test of observed
BFs (bin-width: 3 ms) in each period using spike data from only the last
half of the 36 trials (trials 19–36).
doi:10.1371/journal.pcbi.1002385.g007

Discussion
In this study, we introduced a novel method for estimating
dynamic spike interactions in multiple parallel spike sequences bymeans of a state-space analysis (see Methods for details). By
applying this method to nonstationary spike train data using the
pairwise log-linear model, we can extend the stationary analysis of
the spike train data by the Ising/spin-glass model to within-trial
nonstationary analysis (Figure 3). In addition, our approach is not
limited to a pairwise analysis, but can perform analyses of timevarying higher-order spike interactions (Figure 4). It has been
discussed whether higher-order spike correlations are important tocharacterize neuronal population spiking activities, assuming
stationarity in the spike data [54–59]. Based on the state-spacemodel optimized by our algorithm, we developed two methods tovalidate and test its latent spike interaction parameters, inparticular the higher-order interaction parameters, which may
dynamically change within an experimental trial. In the first
method, we selected the proper order for the spike interactions
incorporated in the model under the model selection frameworkusing the approximate formula of the AIC for this state-spacemodel (In Methods, ‘Selection of state-space model by information
criteria’). This method selects the model that best fits the data
overall across the entire observation period. The selected model
can then be used to visualize the dynamic spike interactions or fora performance comparison with other statistical models ofneuronal spike data. However, more importantly, the detailed
structure of the transient higher-order spike interactions needs tobe tested locally in time, particularly in conjunction with the
behavioral paradigm. To meet this goal, we combined theBayesian model comparison method (the Bayes factor) with asurrogate method (In Methods, ‘Bayesian model comparison
method for detecting spike correlation’). The method allows us totest for the presence of higher-order spike correlations andexamine its relations to experimentally relevant events. We
demonstrated the utility of the method using neural spike train
data simultaneously recorded from primary motor cortex of an
awake monkey. The result is consistent with, and further extended
the findings in the previous report [8]: We detected an increase intriple-wise spike interaction among three neurons in the motor

Figure 7. Detection of triple-wise spike correlation of MI
neurons. (A) (Top) The bin-by-bin Bayes factor (BF) for a model of
triple-wise spike interaction computed locally in time in bit units, Eq. 13.
The bin-by-bin BF is computed from Eq. 13, using the state-space loglinear model fitted to the spike data (a total of 36 trials, 2 s binned using
3 ms bin-width; cf. Figure 6C). The BF computes evidence of a positive

PLoS Computational Biology | www.ploscompbiol.org

17

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

compound Poisson point process (CPP) model for the pooled spiketrain [34,35]. This method can detect an assembly from a few
occurrences of synchronous spike events to which many neurons
belong to, typically by using the lower-order cumulants of theobserved spike counts. In contrast, the information geometry
measure for the higher-order spike correlation used in this study
aims to represent the correlated state that cannot be explained bylower-order interactions. Consequently, the information geometry
measure extracts the relative strength of the higher-orderdependency to the lower-order correlated state. Therefore, thepresence of positive higher-order spike correlations does notnecessarily indicate that many neurons spike synchronously
whenever they spike because such activities can be induced bypositive pairwise spike correlations alone [65,85,86] (see also
Figure 8A and B in the Methods section). In contrast, thecumulant-based correlation method by Staude et al. [34] infers thepresence of ‘higher-order correlation’ for such data by determining
the presence of high amplitudes in the jump size distribution of thecortex during the preparatory period in a delayed motor task,
which was also tightly locked to the expected signals. Although theanalysis was done for a limited number of neurons, smaller thanthe expected size of an assembly, it demonstrates that the
nonstationary analysis of the higher-order activities is useful toreveal cooperative activities of the neurons that are organized inrelation to behavioral demand. Of course, further analysis isrequired to strengthen the findings made above including a metaanalysis of many different sets of multiple neurons recorded under
the same conditions.
In this study, we adopted the log-linear model to describe thehigher-order correlations among the spiking activities of neurons.
There are, however, other definitions for ‘higher-order spikecorrelation’. An important alternative concept is the definition
based on cumulants. Using the cumulants of an observed countdistribution from a spike train pooled across neurons, Staude et al.
developed an iterative test method that can detect the existence ofa high amplitude in the jump size distribution of the assumed

Figure 8. Analysis of stationary spike correlations of N~3 simulated neurons using Bayes factor. (A) Sketch of different time periods andthe underlying models used for the generation of parallel spike sequences: (I) Model of independent spiking (hti ~{2:2, htij ~0, ht123 ~0 for
t~1, . . . ,250 and 1ƒivjƒ3); (II) Model of simultaneously positive pairwise interactions, without a triple-wise interaction (hti ~{2:77, htij ~1:57,
ht123 ~0 for t~251, . . . ,500); (III) Model of triple-wise interaction, with negative pair interactions (hti ~{2:09, htij ~{2:69, ht123 ~10 for t~501, . . . ,750).
(B) Raster display of three parallel spike sequences, Xt,1 , within one example trial. Each spike is colored according to the spike pattern in which itappears: Spikes occurring in triplets are shown in red, spikes within doublets (all types) are marked in blue, and spikes not involved in any synchrony
pattern are shown in black. (C) The bar plots demonstrate the Bayes factors (BFs), Eq. 48, for each of the time periods I–III in bit units. The upper panel
shows the average BF when testing simultaneously positive pairwise interactions (Test 1), averaged across 200 realizations (+1SD). A positive value
for the log of the BF supports the model for the presence of simultaneously positive pairwise interactions, M1 : VI[V2 ,htI w0, while a negative value
supports the absence of such an assembly, M2 : AI[V2 ,htI ƒ0. The BFs per sample and time period are computed by applying a pairwise state-spacelog-linear model (r~2) independently to the three periods. We use a state model with F~I. The lower panel shows the BF for the positive triple-wise
spike interaction (Test 2), M1 : ht123 w0, as opposed to a zero or negative triple-wise spike interaction, M2 : ht123 ƒ0, in each of the three periods. The
BFs are computed from a full state-space log-linear model (r~3). (D) Bar plot of the bin-by-bin BFs (Eq. 52) sorted by the different spike patterns in
the three periods (from top to bottom). The contributions to the BFs (shown in C) from the different spike patterns are sorted and displayed using theindicated spike patterns (000, 100, 110 and 111) as representative examples. The gray bars indicate the average BFs of simultaneously positive
pairwise interactions, with the average computed for the spike patterns observed in 200 realizations in the respective periods, while the dark graybars indicate the average BFs for the triple-wise spike correlation.
doi:10.1371/journal.pcbi.1002385.g008

PLoS Computational Biology | www.ploscompbiol.org

18

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

massively parallel spike sequences that can be analyzed by othermethods such as those based on the statistics pooled acrossneurons. Thus, we consider it to be important to combine thedetailed analysis method proposed in this study with other state-ofthe-art analysis techniques in practical applications. For example,
test methods based on population measures such as the Unitary
Event method and cumulant-based inference method [34,35]
allow us to detect the existence of statistically dependent neurons
in massively parallel spike sequences. If the null-hypothesis ofindependence among those neurons is not rejected in thesemethods, we can exclude those neurons from any further detailed
analysis of their dynamics using the methods proposed in thisstudy.
Several critical assumptions made in the current frameworkneed to be addressed. First, it was assumed in constructing the
likelihood (Eq. 7) that no spike history effect exists in thegeneration of a population spike pattern. Second, we assumed
the use of identically and independently distributed samples acrosstrials when constructing the likelihood (Eq. 7). The first
assumption may appear to be strong constraint given the factthat individual neurons exhibit non-Poisson spiking activities [94].
However, as in the case of the estimation of the firing rate of asingle neuron, the pooled spike train across the (independent) trials
is assumed to obey a Poisson point process because of the general
limit theorem for point processes [30,31,95,96]. This is because
most of the spikes in the pooled data come from independent
different trials. They are thus nearly statistically independent fromeach other, even if the individual processes are non-Poisson.
Similarly, in our analysis, we used statistics from a pooled binary
spike train, assuming independence across trials: The occurrencesof joint spikes in the binary data pooled across trials are mostly
independent of each other across bins. Because these joint spike
occurrences are sparse (i.e., they rarely happen closely to eachother in the same trial), it is even more feasible to assume theirstatistical independence across bins. Third, however, while poolingindependent and identical trials (the second assumption) may
validate the first assumption of the independence of the samples
across bins, that assumption of independently and identically
distributed samples across trials has itself been challenged [97,98]
and is known to be violated in some cases, e.g., by drifting
attention, ongoing brain activity, adaptation, etc. It is possible thatthe trial-by-trial jitter/variation in the spike data causes spurious
higher-order spike correlation. Thus, as discussed in the section on
the application of our methods to real neuronal data, it isimportant to examine the stationarity of the spike train data acrosstrials. Note that, not only the firing rates, but also the spikesynchrony can be shaped on a longer time scale by repeatedly
practicing a task [19]. In fact, the current analysis method can beused to examine the long-term evolution of pairwise and higherorder spike interactions across trial sessions by replacing the role ofa bin with a trial, assuming within-trial stationarity. It will be achallenge to construct a state-space log-linear model that
additionally applies a smoothing method across trials (see [98]for such a method for a point process model).
The present method is left with one free parameter, namely thebin-width D. The bin-width determines the permissible temporal
precision of synchronous spike events. Very large bin-widths resultin binary data that are highly synchronized across sequences, whilevery small bin-widths result in asynchronous multiple spikesequences. In the latter case, we might overlook the existing
dependency between multiple neural spike sequences due todisjunct binning [99] (but see [57,100] that aim to overcome sucha problem by modeling the spike interactions across different
consecutive time bins). Within our proposed modeling framework,

assumed CPP model. Yet another important tool for analyzing
higher-order dependency among multiple neurons is the copulafunction, a standardized cumulative distribution function used tomodel the dependence structure of multiple random variables (see[87–89] for an analysis of neurophysiological data using the
copula, including an analysis for higher-order dependency [89]).
In summary, it should be remembered that the analysis method
used for the higher-order dependency among neuronal spikes
inherits its goal from the assumed model for spike generation as
well as a parametric measure defined for the ‘higher-order’ spikecorrelation [34,62].
Although we face a high-dimensional optimization problem inour settings, we are able to successfully obtain MAP estimates of
the underlying parameters because of the simplicity of theformulation of the state-space model: The use of the log-concave
exponential family distributions [50,90] in both the state andobservation models guarantees that the MAP estimates can beobtained using a convex optimization program. At each bin, themethod numerically solves a nonlinear filter equation to obtain themode of the posterior state density (the MAP estimates, see Eqs.
28, 29, and 30 in Methods). With only a few (3–8) NewtonRaphson iterations, the solution reaches a plateau (the incrementsof all the elements of the updated state space vector are smaller
than 10{5 ). The entire optimization procedure can be performed
in a reasonable amount of time: On a contemporary standard
laptop computer it takes no longer than 30 seconds to obtainsmooth estimates of a full log-linear model for N~3 neurons(T~500 bins, Figure 4), which includes 100 EM iterations. The
method is even faster when approximating the posterior modeusing the update formula Eq. 30 without any iterations, using the
one-step prediction mean as an initial value. This fast approximation method suggested in [69] could even be utilized in a realtime, on-line application of our filter (the filtering method applied
to a single trial, n~1, using predetermined hyper-parameters) at
the cost of estimation accuracy.
The pairwise analysis can be applied up to about N~12
neurons simultaneously to derive time-dependent pair interactions.
However, the current version of the algorithm does not scale to alarger number of neurons because the number of spike patterns
that need to be considered suffers from a combinatorial explosion.
The major difficulty arises from the coordinate transformationfrom the h-coordinates to the g-coordinates that appear in the
non-linear filter equation (Eq. 31 in Methods). The coordinatetransformation is required in this equation to calculate theinnovation signal, i.e., the difference between the observed
synchrony rates, yt , and the expected synchrony rates (gcoordinates) based on the model. We numerically derived theexact g-coordinates by marginalizing the 2N dimensional joint
probability mass function computed from the h-coordinates.
Thereby, a full knowledge of the probability mass function isrequired even if the model considers only the lower-orderinteractions. Because this is a common problem in the learningof artificial neural networks [46,47,91], sampling algorithms suchas the Markov chain Monte Carlo method have been developed toapproximate the expectation parameters, h, without having tocompute the partition function [92]. The inclusion of suchmethods allows us to analyze the time-varying low order spikeinteractions from a larger number of parallel spike sequences.
Recent progress [93], e.g., in the mean field approach and/or theminimum probability flow learning algorithm for an Ising model,
may allow us to further increase the number of neurons that canbe treated in this nonstationary pairwise analysis. Nonetheless, themethod presented here, which aims at a detailed analysis of the
dynamics in higher-order spike interactions, may not easily scale to
PLoS Computational Biology | www.ploscompbiol.org

19

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

which focuses on instantaneous higher-order spike correlations, it
is important to catch the innate temporal precision of the neuronal
population under investigation using the appropriate bin-width.
Thus, the choice may be guided by the biophysical properties of
the neurons. However, it may be of advantage to derive the binwidth in a data-driven manner. For example, in the context of an
encoding problem, the proper bin-width can be chosen based on
the goodness-of-fit test for single neuron spiking activities [101],
conditional on the spiking activities of the other neurons [40]. For
questions about the relation of coincident spiking to stimulus/
behavior, the bin-width may be selected based, for example, on thepredictive ability of an external signal. For this goal, it is importantto search the optimal bin-width using elaborate methods such as
those developed in the context of the Unitary Event analysis
method [99] (see [25] for a review of related methods).
A substantial number of studies have demonstrated that
stimulus and behavioral signals can be decoded simply based on
the firing rates of individual neurons. At the same time, it has been
discussed whether spike correlations, particularly higher-order
spike correlations, are necessary to characterize neuronal population spiking activities [54–58] or to encode or decode
information related to stimuli [53,60,102]. At this point in time,a smaller number of dedicated experiments have supported theconceptual framework of information processing using neuronalassemblies formed by neurons momentarily engaged in coordinated activities, as expressed by temporally precise spikecorrelations (see [6,7,9,10] for reviews of these experiments).
Nevertheless, it is possible that the current perspective on thissubject has been partly formed by a lack of proper analysis
approaches for simultaneously tracing time-varying individual
pairwise spike interactions, and/or their higher-order interactions.
Indeed, we demonstrated by the time-resolved higher-orderanalysis that three cortical neurons coordinated their spiking
activities in accordance with behaviorally relevant points in time.
Thus our suggested analysis methods are expected to be useful toreveal the dynamics of assembly activities and their neuronalcomposition, as well as for testing their behavioral relevance. We
hope that these methods help shed more light on the cooperative
mechanisms of neurons underlying information processing.

transformation from one to the other is given by the Legendre
transformation [40,50]. From Eq. 14, the log normalization
function, y(h), is written asP"#
hI fI ðxÞXy(h)~logeI:

The first derivative of the log normalization function, y(h), withrespect to hI (I[fV1 , . . . ,VN g), provides the expectation parameter, gI ~E½fI (X):
Ly(h)
~gI :
LhI

1Q(g)~LI LJ y(h)~gIJ

ð19Þ

for I,J[fV1 , . . . ,VN g. gIJ is known as the Fisher metric withrespect to the natural parameters. Eqs. 17 and 19 are important
relations used in this study to construct a non-linear filtering
equation for a dynamic estimate of the natural parameters because
we approximate the log-linear model (Eq. 14) with a precision of
up to a (log) quadratic function (cf. Eqs. 28 and 29). Similarly, thehigher-order derivatives yield higher-order multivariate cumulants. For example, the third-order derivative yields the third order
cumulant, LI LJ LK y(h)~C?IJ,K , where C?IJ,K :gI|J|K {gI gJ|K
{gJ gI|K {gK gI|J z2gI gJ gK .
The pseudo distance between two different distributions, q(x)and p(x) is defined using the Kullback-Leibler (KL) divergence

ð14Þ

and
ð15Þ

where fI ðxÞ is a feature function, here representing an interaction
among the neurons indicated by the multi-index, I (Eq. 3).
The h- and g-coordinates are dually flat coordinates in the
exponential family probability space [44,50], and the coordinate
PLoS Computational Biology | www.ploscompbiol.org

ð18Þ

Eqs. 17 and 18 complete the Legendre transformation from hcoordinates to g-coordinates [44,50]. The Legendre transformation transfers the functional relationship of h and y(h) to theequivalent relation in the dual coordinates, g and Q(g). The
LQ(g)inverse transformation is given by Eq. 18 and
~hI .
LgI
Using the log normalization function, we can obtain themultivariate cumulants of p(x) with respect to the randomvariables, fI (X). The cumulant generating function of the
exponential family distribution is given as K ðuÞ~y(hzu){y(h). Let us compactly write the partial derivative with respectL
to hI (i.e.,) as LI . Then, the first order cumulant is given as
 LhI

L
K ðuÞ ~LI y(h)~gI , as shown in Eq. 17. In general, theLu
u~0cumulants of the exponential family distribution are given by
the derivatives of the log normalization function. Thus, thesecond derivative of y(h) yields the second-order cumulant,gIJ :E½(fI (X){gI )(fJ (X){gJ )~gI|J {gI gJ (by the cup, |, wemean the multi-index representation of an union of the elements ofthe two multi-indices, e.g., if I~12 and J~234, then
I|J~1234):

I

gI ~E½fI (X),p(x) log p(x)~h’g{y(h):In this subsection, we review the known mathematicalproperties of a log-linear model for binary random variables.
These properties are used in constructing recursive filtering/smoothing formulas in the next section. Using the multi-index,
I[fV1 , . . . ,VN g (see the Results subsection ‘Log-linear model ofmultiple neural spike sequences’), the probability mass function(Eq. 1), p(x), where x~ðx1 ,x2 , . . . ,xN Þ and xi ~1 or 0
(i~1, . . . ,N), and the expectation parameters (Eq. 2) arecompactly written as
hI fI ðxÞ{y(h)XxMathematical properties of log-linear model

X

ð17Þ

Let Q(h) be the negative entropy of the distribution:

Methodslog p(x)~ð16ÞxD½q,p:X

q(x) log q(x)=p(x)

ð20Þ

x

We represent distribution q(x) by using g-coordinates as gðqÞ, and20

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlationp(x) by using h-coordinates as h ð pÞ. Here, we used h ð pÞ for the hparameters of p(x) (and gðqÞ for g-parameters of q(x)) in order todifferentiate it from the representation of distribution q(x) in the hcoordinates (and the representation of p(x) in the g-coordinates).
Then, the KL-divergence between the two distributions, q(x) andp(x), is computed as [44,50]
D½q,p~y(hð pÞ)zQ(gðqÞ){h ð pÞgðqÞ:at time t given the spike data up to and including time t, whereas
the backward algorithm constructs a posterior density at time t
given the entire data. The posterior density allows us to computethe maximum a posteriori (MAP) estimate or Bayes estimator andprovides uncertainty for the estimate. In the following, h tjs andmean, E ½h t jy1:s ,w, and covariance,
Wtjs denote the conditional0E½(h t {h tjs )(ht {h tjs ) jy1:s ,w. The filter mean and covariance aredenoted as h tjt and Wtjt , respectively. The mean and covariance ofa smooth posterior density are denoted as h tjT and WtjT ,respectively.
We first compute the one-step prediction density,
pðh t jy1:t{1 ,wÞ. This is the conditional density of the state at time
t given the observation of parallel spike sequences up to time t{1.
The one-step prediction density is written using the ChapmanKolmogorov equation as [37,69,105]

ð21Þ

Optimized estimation of dynamic spike interactions
We develop a non-linear recursive Bayesian filtering/smoothingalgorithm and its optimization method in order to trace
dynamically changing spike interactions from parallel spikesequences. To reach this goal, we use the expectation-maximization (EM) algorithm [68,73,103,104], which is known to efficiently
combine the construction of the posterior density of a state (thenatural parameters) and the optimization of the hyper-parameters.
This method maximizes the lower bound of the log marginal
likelihood, Eq. 10. Using Jensen’s inequality and nominal hyperparameters, w~½F,Q,m,S, the lower bound of the log marginal
likelihood with hyper-parameters w ~½F ,Q ,m ,S  is given byl(w )~logðpðh 1:T jy1:T ,wÞðpðh t jy1:t{1 ,wÞ~ pðh t jh t{1 ,wÞpðh t{1 jy1:t{1 ,wÞdht{1 :Here the transition probability, pðh t jht{1 ,wÞ, is a multivariate
normal distribution with mean Fh t{1 and covariance Q, as
defined in the state equation, Eq. 8. For an initial prior, pðh 1 Þ, we
use a normal distribution with mean h 1j0 ~m and covariance
W1j0 ~S. The other distribution, pðh t{1 jy1:t{1 ,wÞ, in Eq. 24 is thefilter density at time t{1. In the next paragraph, the filter density
will be obtained by approximating it with a normal distribution
whose mean and covariance are denoted as h t{1jt{1 andWt{1jt{1 . Under this condition, the one-step prediction density(Eq. 24) again becomes a normal distribution whose mean, h tjt{1 ,and covariance, Wtjt{1 , are given as [37,68]

pðy1:T ,h 1:T jw Þ
dh 1:Tpðh 1:T jy1:T ,wÞðpðy1:T ,h 1:T jw Þ
dh 1:T§ pðh 1:T jy1:T ,wÞlog
pðh 1:T jy1:T ,wÞð22Þ~Qðw jwÞzH½pðh 1:T jy1:T ,wÞ:
Here, H represents a negative entropy. The maximization of thelower bound with respect to w is equivalent to maximizing the
expected complete data log-likelihood in Eq. 22, known as the Qfunction:h tjt{1 ~Fh t{1jt{1 ,0

Wtjt{1 ~FWt{1jt{1 F zQ:


ð24Þ

ð25Þ

ð26Þ



Qðw jwÞ~E ½log pðy1:T ,h 1:T jw Þjy1:T ,w~nTXThe filter density, pðh t jy1:t ,wÞ, is the conditional distribution of
the state given the observation of parallel spike sequences up to
time t. Using the likelihood function and one-step prediction
density, the filter density is given by Bayes’ theorem as0

yt Eht {Eyðh t Þ

t~1

d
1
{ log 2p{ log det S2201{ E½ðh 1 {m Þ S{1 ðh 1 {m Þ
2
ðT{1ÞdðT{1Þ{log 2p{log det Q
2
2
T
0
1X
{
E½ðh t {F h t{1 Þ Q{1 ðh t {F h t{1 Þ:
2 t~2

ð23Þ
pðht jy1:t ,wÞ~!exp½n0yt ht {yðh t Þð27Þ
0 {1 1{ h t {h tjt{1 Wtjt{1 h t {h tjt{1 :
2This posterior density is a complicated function with respect to thenatural parameter, h t . Here, we apply a Gaussian approximation
to the posterior density using the Laplace method [37,104,106,107]: The filter mean, h tjt , is identified with a mode of
the posterior density as h tjt ~argmaxq log pðhjy1:t ,wÞ, and thefilter covariance, Wtjt , is determined from the Hessian of the log

The expectation in the above equation is read as E½ jy1:T ,w. We
maximize the Q-function by alternating the expectation (E) andmaximization (M) steps. In the E-step, we obtain the expectedvalues with respect to h t in Eq. 23 using a fixed w. In the M-step,
we obtain the hyper-parameter, w , that maximizes Eq. 23. The
resulting w is then used in the next E-step. The details of eachstep are now given as follows.
E-step: Bayesian recursive filter/smoother. The E-step iscomposed of filtering and smoothing algorithms conducted byforward and backward recursions, respectively. The forward
algorithm sequentially constructs a posterior density of the state
PLoS Computational Biology | www.ploscompbiol.org

pðyt jh t ,y1:t{1 Þpðh t jy1:t{1 ,wÞpðyt jy1:t{1 ,wÞposterior at the mode as Wtjt ~{½++log pðhjy1:t ,wÞjh~h {1 .
tjt

The approximate posterior mode is obtained using the iterative
procedure of a gradient ascent method or the Newton-Raphson
method using the gradient and Hessian matrix. The gradient andHessian of the log of the posterior density at h is calculated as{1+f:+log pðhjy1:t ,wÞ~nðyt {+yðh ÞÞ{Wtjt{1h{h tjt{1 , ð28Þ21

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

H:++f~{n++yðh Þ{W{1
tjt{1 :

Because the density functions in the recursive formula wereapproximated as a normal distribution, we follow the fixed-intervalsmoothing algorithm [37,68,105] established for a Gaussian stateand observation equation. Starting from h TjT and WTjT , which
are obtained from the filtering algorithm, we obtain the smoothedmean and covariance,

ð29Þ

Note that, as in Eqs. 17 and 19, the first derivative of the log
normalization function, yðh Þ, with respect to h provides the dual
coordinates, g: +yðh Þ~g. Furthermore, the second derivative
yields the Fisher metric: ++yðh Þ~G. In this study, we adopt the
Newton-Raphson method:
h new ~h old {EH{1 +f:

ð30Þ

The gradient and Hessian, Eqs.28 and 29, are evaluated using an
old value, h old . Here, E is a learning coefficient that was introducedin the context of ‘natural’ gradient search algorithm [91]. For thesmall population size analyzed in this study, we used E~1.
However, it is recommended that a smaller positive value be
selected for the analysis of a larger system to avoid numericalinstability. Because both the likelihood function and prior density
are logarithmically concave functions, the posterior is also logconcave. Thus, this optimization problem is convex, whichguarantees a unique solution for the filter estimate [90]. The
optimized natural parameter is selected as the filter mean, h tjt . The
filter covariance is approximated as Wtjt ~{H{1 by using h tjt . It
is also possible to use a simple gradient ascent method to obtainthe mode, and then compute the Hessian matrix at the mode.
The above method is equivalent to solving the following
nonlinear recursive filter formulas:
h tjt ~h tjt{1 znWtjt{1 (yt {gtjt ),

ð31Þ

{1
W{1
tjt ~Wtjt{1 znGtjt :

ð32Þ

ð


 0
WtjT ~Wtjt zAt Wtz1jT {Wtz1jt At :

ð35Þ

0

At ~Wtjt F W{1tz1jt :

ð36Þfor t~T{1,T{2, . . . 1. The lag-one covariance smoother,
Wt{1,tjT , which appears in the Q-function, is obtained using themethod of De Jong and Mackinnon [108]:0
Wt{1,tjT ~E½(h t{1 {h t{1jT )(h t {h tjT ) y1:T 

ð37Þ

~At{1 WtjT :M-step: Optimization of hyper-parameters. At the Mstep, we optimize hyper-parameter w given the posterior density
under the principle of maximizing the Q-function. FromLQðw jwÞ~0, the update rule of the covariance matrix, Q,
LQ
is obtained as
T
1 X
0
0
0
½WtjT {Wt{1,tjT F {FWt{1,tjT zFWt{1jT F 
T{1 t~2
ð38Þ
T 

0
1 X
z
h tjT {Fh t{1jT h tjT {Fh t{1jT :T{1 t~2Q ~Here, htjT , WtjT , and Wt{1,tjT are the smoother mean andcovariance, and the lag-one covariance matrix given by Eqs. 34,
LQðw jwÞ~0, the35, and 37, respectively. Similarly, fromLFauto-regressive parameter, F, is updated according to
"
#"
#{1
T 
T 
X
 X

:
Wt{1,tjTzh tjT h’t{1jT
Wt{1jTzh t{1jT h’t{1jT
:ð39Þ
F~


t~2

t~2

The mean of the initial distribution is updated with m ~h 1jT fromLQðw jwÞ~0. The covariance matrix, S, is not updated;Lm
instead, we use a fixed matrix, S. Alternatively, the covariancematrix of the initial distribution can be updated according to
LS ~W1jT z(h 1jT {m)(h 1jT {m)’ from  Qðw jwÞ~0, while theLSmean, m, is fixed. Both methods work well with appropriate
choices for the fixed parameters. In this study, we updated themean vector, m, of the initial normal distribution, and used a fixeddiagonal matrix for its covariance matrix, S. It was also suggestedpðhtz1 jy1:T ,wÞpðh tz1 jh t ,wÞdh tz1 : ð33Þ
pðh tz1 jy1:t ,wÞPLoS Computational Biology | www.ploscompbiol.org

ð34Þ

with

Eq. 31 was obtained from +log pðh t jy1:t ,wÞ~0. Eqs. 31 and 32
are recursively computed in combination with the one-step
prediction equations, Eqs. 25 and 26, for t~1, . . . ,T. Because
h tjt and gtjt are dual representations of the same probability
distribution, Eq. 31 is a nonlinear equation and needs to be solved
as suggested above (Eqs. 28, 29, 30). As pointed out for a point
process adaptive filter [37,69], the residual, yt {gtjt , in Eq. 31 acts
similarly to an innovation vector of a standard Kalman filter. The
same error signal between the observed synchrony rates andexpected synchrony rates is utilized in training the Boltzmann
machine [46,47,91]. The innovation term corrects the one-step
prediction mean, h tjt{1 , i.e., an expected state by the priordistribution. The degree of the correction is determined by the
number of repeated trials, n, and uncertainty of the predictedstate, Wtjt{1 . The hyper-parameters of the prior density
significantly affect the latter gain (see Eq. 26), and thereby thesmoothness of the estimated processes. The filter covarianceequation, Eq. 32, describes the reduction of the prediction
uncertainty, Wtjt{1 , by observing the parallel spike sequences attime t, with the amount determined by the number of repeatedtrials, n, and the Fisher information, Gtjt . Please see Figure 1 for ageometric view of the recursive Bayesian filter.
Finally, we compute the smooth posterior density using the
backward recursive algorithm [68,105,107],pðh t jy1:T ,wÞ~pðht jy1:t ,wÞh tjT ~h tjt zAt h tz1jT {h tz1jt ,22March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

to use a stationary mean and covariance of an unconstrained stateprocess (Eq. 8) as parameters of the initial prior distribution[68,109].
The equilibrium condition from Eq. 8 yields m~Fm andX~F
F0 zQ. The solutions are obtained as m ~0 and

l(w)&

S

T
X
t~1

vec(g )~½I{F6F{1 vec(Q) (p. 121 and p. 426 in [109], p. 112
in [110]). In the latter, 6 denotes the Kronecker product (tensor
product) and the vec operator creates a single column vector froma matrix by stacking its column vectors. However, the stationary
condition is not always satisfied. As we found that the fitted F
sometimes violates the stationary assumption, we did not adopt
this approach in the current study.

~

TX

The method developed in the previous subsection is applicable
to a full log-linear model, as well as a model that considers anarbitrary order of interactions. In order to select the most
predictive model among the hierarchical models in accordancewith the observed spike data, we select the state-space model thatminimizes the Akaike information criterion (AIC) for a model withlatent variables [73,105]:
ð40Þ

~

T
X

ð
log

2 PDIO~{2 l ðwÞz2 tr

p(yt jh t )p(h t jy1:t{1 ,w)dh t :We make a log quadratic approximation to evaluate the integral.
To accomplish this, we denoteð0

exp½q(h t )dht , ð42Þ

0


1
h t {h tjt{1 W{1tjt{1 h t {h tjt{1 : ð43Þ2

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
  
exp½qðh t Þdh t & (2p)d jWtjt j exp q h tjt :B12 (ya:b )~ð44Þ

pðya:b jM1 ,wÞ:pðya:b jM2 ,wÞð48Þ

The BF becomes larger than 1 if the data, ya:b in a time period½a,b, supports model M1 as opposed to model M2 , and becomes
smaller than 1 if the data supports model M2 as opposed to model
M1 . The BF can be computed from the one-step prediction andBy applying Eqs. 42, 43 and 44 to Eq. 41, the log marginal
likelihood is approximated as
PLoS Computational Biology | www.ploscompbiol.org

ð47Þ

In this subsection, we formulate a method for detecting thehidden structure of spike interaction by means of a Bayesian modelcomparison based on the Bayes factor (BF) [74–76]. The BF is aratio of likelihoods for the observed data, based on two different
assumptions about the hidden states (model M1 and M2 ). Here we
reiterate the definition of the BF for model M1 as opposed tomodel M2 used in this paper (cf. Eq. 12):

The Laplace approximation of the integral in Eq. 42 is given as
[80]
ð

!
LM(w) {1,Lw0

Bayesian model comparison method for detecting spikecorrelation

withq(h t )~n yt h t {yðh t Þ {I{Here, we redefine w as a one-dimensional vector of free hyperparameters, while M(w) denotes the one-step operator of EM
iteration. To obtain the Jacobian matrix for the EM operator, wefollow the algorithm described in Meng and Rubin [111]. In thismethod, the Jacobian matrix was approximated using a numericaldifferentiation of the EM operator. By perturbing one hyperparameter and then computing a one-step EM iteration,
numerical differentiations of the hyper-parameters with respect
to the perturbed hyper-parameter were obtained. An entire
Jacobian matrix was approximated by repeating the process whilechanging the hyper-parameter to be perturbed.

ð41Þ

1p(yt jht)p(ht jy1:t{1 ,w)dh t ~ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃd(2p) jWtjt{1 jð46Þ

Shimodaira’s predictive divergence for indirect observation models(PDIO) [81] is given as

t~1

ð

T 

1X
log det Wtjt {log det Wtjt{1 :2 t~1BIC~{2 l(w)zk log n:log p(yt jy1:t{1 ,w)

t~1

ð45Þ

We confirmed that the log-quadratic approximation provided abetter estimate of the marginal likelihood than the first orderapproximation used in [77] by comparing them with a Monte
Carlo approximation of the integral in Eq. 42. We select the statespace log-linear model that minimizes the AIC (Eq. 40), where thelog marginal likelihood is approximated using Eq. 45.
For comparison with the AIC, we compute two otherinformation criteria that employ different forms of the penalization
term. The Bayesian information criterion [79,80] (also known as
Schwartz’s criterion) are obtained by replacing the penalization
term of Eq. 40, k, with k log n:Here, l ðwÞ is the log marginal likelihood (Eq. 10) and k is thenumber of free hyper-parameters in the prior distribution. For the
rth order model,
the number of natural parameters is given by
Xr  N d~. The number of free parameters in the prior
k~1
kdistribution is computed as k~d 2 zd(dz1)=2zd, where eachterm corresponds to the number of free parameters in F, Q, and m.
Note that the AIC applied to the state-space model is sometimes
referred to as the Akaike Bayesian information criterion (ABIC)[73]. In the following, we derive the approximation method toevaluate the AIC for the state-space log-linear model.
The log marginal likelihood, Eq. 10, can be written as
T
X

 
0
n yt h tjt {y h tjt

T
h


0 i
1X
{
tr W{1tjt{1 h tjt {h tjt{1 h tjt {h tjt{1
2 t~1

Selection of state-space model by information criteria

l(w)~

T 

1X
log det Wtjt {log det Wtjt{1
2 t~1

t~1

z

AIC~{2l ðwÞz2kq(htjt )z23

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

analysis where the BF is applied to stationary spike sequences ofN~3 simulated neurons. To demonstrate the kind of correlation
scenarios that is distinguished by the BF, we generate stationary
spike sequences using log-linear models of different spikecorrelation structures in three distinct time periods (I–III,
Figure 8A). These are: (I) an independent model (htI ~{2:2 forI[V1 and htI ~0 for I[fV2 ,V3 g at t~1, . . . ,250), (II) a model
containing simultaneously positive pairwise interaction terms,
without a triple-wise interaction (htI ~{2:77 for I[V1 , htI ~1:57
for I[V2 , and ht123 ~0 at t~251, . . . ,500), and (III) a model
containing a positive triple-wise interaction, with negative pairwiseinteractions (htI ~{2:09 for I[V1 , htI ~{2:69 for I[V2 , ht123 ~10
at t~501, . . . ,750). A sample of the spike sequences from a single
trial (n~1) is shown in Figure 8B. In all of the models, the firstorder log-linear parameters, hI (i[V1 ), are adjusted such that the
spike rates of the individual neurons are gI ~0:1½spikes=D (I[V1 ).
Therefore, the time segments cannot be distinguished on the basisof the firing rates. In addition, the model in period III is designed
such that the projection to E2 yields zero pairwise interactions, i.e.,
the projection yields an independent model. In this model, thesecond-order joint synchrony rates, gI ~0:01 (I[V2 ), are at the
chance level expected from the individual spike rates, gI ~0:1
(I[V1 ). However, the model expresses excess joint spike synchrony
between all three neurons, g123 ~0:0092, which is larger than the
expected rate (0.001). The presence of a positive triple-wisecorrelation, combined with the absence of marginal pairwise
correlations, leads to a very sparse number of excess synchronous
spike triplets, as illustrated by the dot displays in period III in
Figure 8B. In contrast, a period containing only positive pairwise
correlations without a triple-wise interaction (see Figure 8B, periodII) contains a relatively large number of triplets, which are,however, the mere expression of excess pairwise correlations.
We calculate the BFs to test for the presence of the two different
correlation structures, i.e., the same tests discussed in the Results
section: Test 1 checks for the presence of simultaneously positive
pairwise interactions and Test 2 checks for the presence of a triplewise spike correlation. In one test (Test 1) we ask whether the threeneurons simultaneously exhibit positive pairwise interactions byapplying a model that contains up to pairwise interactions (r~2) toeach of the three time periods. In that case, the BF (Eq. 48) is
computed using a model which assumes that all of the pairwise
terms are positive: M1 , specified as S1 : htI w0 for all I[V2 , as
opposed to a model where at least one of the pairwise interactions
is not positive: M2 , specified as S2 : htI ƒ0 for at least one I[V2 .
Time (bin) index, t, expresses that the BF is computed under thesame models for all time steps t in the respective time period. In
both models, we do not make specific assumptions for the firstorder log-linear parameters; thus, they are allowed to be real
numbers. In the other test (Test 2), we ask whether the neurons
exhibit excess triple-wise synchrony by applying a full model (r~3)
to each of the periods. Here, the BF is computed with the opposing
models: M1 , specified as S1 : ht123 w0, and M2 , specified as
S2 : ht123 ƒ0. The lower-order log-linear parameters have again
real values.
Figure 8C displays the average BFs (in units of bits) computedfrom 200 samples (each sample contains n~1 trial) within each of
the periods I–III. In addition, we show the standard deviations of
the BFs (shown as +1 SD) derived from the BFs of the 200
realizations. The upper panel shows the result of Test 1, expressing
the weight of evidence for the presence of excess synchrony
realized by simultaneously positive pairwise interactions. We findsubstantial evidence for the simultaneously positive pairwiseinteractions in period II, but not in the other two periods. The
lower panel displays the evidence for the presence of a positive

filter density obtained in the method developed in the preceding
subsection. From Eq. 48, the BF can be rewritten asb

B12 (ya:b )~pðyt jy1:t{1 ,M1 ,wÞ:t 1:t{1 ,M2 ,wÞt~a

P pðy jy

ð49ÞLet us define the bin-by-bin BF for the spike data at time t as
B12 (yt jy1:t{1 )~pðyt jy1:t{1 ,M1 ,wÞ:pðyt jy1:t{1 ,M2 ,wÞð50Þ

Using Bayes’ theorem, we obtain [74–76]
pðMi jyt ,y1:t{1 ,wÞ~p(yt jy1:t{1 ,Mi ,w)p(Mi jy1:t{1 ,w)p(yt jy1:t{1 ,w)

ð51Þ

for i~1,2. Using Eq. 51, we can rewrite the bin-by-bin BF aspðM1 jy1:t ,wÞ pðM1 jy1:t{1 ,wÞ=pðM2 jy1:t ,wÞ pðM2 jy1:t{1 ,wÞÐÐ
S1 pðh t jy1:t ,wÞdh t
S pðh t jy1:t{1 ,wÞdh t
=Ð 1
,
~Ð
S pðh t jy1:t ,wÞdh t
S pðh t jy1:t{1 ,wÞdh t

B12 (yt jy1:t{1 )~

2

ð52Þ

2

where S1 and S2 denote spaces that the natural parameters
occupy, supported by models M1 and M2 , respectively. Here,pðh t jy1:t ,wÞ and pðh t jy1:t{1 ,wÞ are the filter density and one-step
prediction density, respectively. By sequentially computing Eq. 52,
we can obtain the BF with respect to the sub-interval ½a,b asbP

B12 (ya:b )~B (y jy).
t~a 12 t 1:t{1A test with the following models in the sub-interval ½a,b allows
us to detect a momentarily active cell assembly of more than two
neurons by the presence of simultaneously positive rth-order spikeinteractions. In the rth-order log-linear model of N neurons, thenatural parameters, htI (I[Vr ), represent the r th-order spikeinteractions among neurons denoted in index I. We examine
whether a subset of N’ neurons among the total N neurons
simultaneously exhibit positive r th-order interactions. Let V’1 bethe subset of N’ neurons from N neurons, e.g., V’1 ~f2,3,5,8g
from V1 ~f1,2, . . . ,10g if N’~4 and N~10. Let V’r be an rsubset from V’1 , e.g. V’2 ~f23,25,28,35,38,58g if r~2. Then, themodel where the subset neurons simultaneously exhibit positive
rth-order interactions (M1 ) and its complementary hypothesis (M2 )are represented as
S1 : htI w0 for all I[V’rð53Þ

S2 : htI ƒ0 for at least one of I[V’rð54Þ

for t~a,az1, . . . ,b. The remaining parameters are real: htI [R(I[fV1 , . . . ,Vr{1 g and I[Vr . excluding V’r ). These parameters
are integrated out in Eq. 52. The above definition of the assemblyis a clique [12], a subset in which each neuron is connected toevery other neuron through the positive rth-order interactions.
Depending on the assembly structure one wishes to uncover, othermodels can be tested such as one in which neurons bounded in
non-exclusive manner.
Analysis of Bayesian model comparison method using
simulated stationary spike sequences. In this subsection, we

analyze the Bayesian model comparison method by using
simulated stationary spike sequences. Figure 8 illustrates an
PLoS Computational Biology | www.ploscompbiol.org

24

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

triple-wise correlation as evaluated by Test 2. The existing positive
triple-wise interaction present in period III is well detected by thelarge BF, and the very low BFs in period I and II correctly detect
the absence of triple-wise correlation.
In Figure 8D, we examine how much each individual spikepattern contributes to the calculation of the BFs discussed
above. The weight of evidence (log of the BF) in an observation
period ½a,b is obtained by summing the pieces of evidence
computed locally at time t using Eq. 52, as log2 B12 (ya:b )~P
bt~a log2 B12 (yt jya:t{1 ). Here we elucidate the contributions ofeach individual spike pattern to the BF by sorting the bin-by-bin
evidence, log2 B12 (yt jya:t{1 ), with each spike pattern. Figure 8D
displays the average BFs for each spike pattern in the three periods
in bit unit. The degree of the contribution by a specific spikepattern to the evidence for the tested spike interaction model varies
depending on the context in which the spike pattern is observed.
For example, the magnitudes of the BFs for a triple-wisecorrelation (Test 2) found by observing a spike triplet (111) vary
greatly across the three periods (see dark gray bars at 111 inperiods I, II, and III). In period III, where pair-synchronous spikes
are rarely observed because of the existence of negative pairwiseinteraction terms, the observation of triplet spikes provides
substantial evidence of a triple-wise correlation. In contrast, inperiod II, where pair-synchronous spikes frequently occur becauseof the existence of positive pairwise interactions, the chance triplet
spikes do not provide substantial evidence of a triple-wisecorrelation. These results come from the fact that an unexpectedsynchronous spike pattern that significantly alters and updates thefilter density from its one-step prediction density gives rise to alarge absolute BF value (Eq. 52). By the same logic, the BFs for the
000 and 100 patterns (and 010 and 001, not shown) are small
because they do not substantially change the posterior densities.
However, they should not be neglected because these patterns areabundant as a result of the low firing rates. For example, in period
I in Figure 8C, the accumulation of the small weight of evidence
from the abundant spike patterns (such as 000 and 100) offsets thelarge weight of evidence induced by a few chance coincidencessuch as triplet spikes (111), thus producing virtually zero weight ofevidence for this period.

Acknowledgments
We thank B. Staude, S. Koyama, Z. Chen, D. Ba, K. Sadeghi, and T.
Toyoizumi for the valuable discussions, and the Diesmann Laboratory forproviding the computing environment. We also thank Prof. Alexa Riehle
for kindly allowing the use of her data in this study.

Author Contributions
Conceived and designed the experiments: HS SA ENB SG. Performed theexperiments: HS SG. Analyzed the data: HS SG. Contributed reagents/
materials/analysis tools: HS SG. Wrote the paper: HS SA ENB SG.

References
1.
2.
3.
4.

5.
6.
7.

8.

9.
10.
11.

12.
13.

14.

15.

16.

17.

18.

19.

Hebb DO (1949) The Organization of Behavior: A Neuropsychological
Theory. New York: Wiley.
Abeles M (1991) Corticonics: Neural Circuits of the Cerebral Cortex.
Cambridge; New York: Cambridge University Press.
Diesmann M, Gewaltig M, Aertsen A (1999) Stable propagation ofsynchronous spiking in cortical neural networks. Nature 402: 529–533.
Kumar A, Rotter S, Aertsen A (2010) Spiking activity propagation in neuronal
networks: reconciling different perspectives on neural coding. Nat Rev
Neurosci 11: 615–627.
Kuhn A, Aertsen A, Rotter S (2003) Higher-order statistics of input ensembles
and the response of simple model neurons. Neural Comput 15: 67–101.
Gerstein GL, Bedenbaugh P, Aertsen MH (1989) Neuronal assemblies. IEEE
Trans Biomed Eng 36: 4–14.
Fujii H, Ito H, Aihara K, Ichinose N, Tsukada M (1996) Dynamical cell
assembly hypothesis - theoretical possibility of spatio-temporal coding in thecortex. Neural Netw 9: 1303–1350.
Riehle A, Grün S, Diesmann M, Aertsen A (1997) Spike synchronization andrate modulation differentially involved in motor cortical function. Science 278:
1950–1953.
Sakurai Y (1996) Population coding by cell assemblies–what it really is in thebrain. Neurosci Res 26: 1–16.
Harris KD (2005) Neural signatures of cell assembly organization. Nat Rev
Neurosci 6: 399–407.
Aertsen AM, Gerstein GL, Habib MK, Palm G (1989) Dynamics of neuronalfiring correlation: modulation of ‘‘effective connectivity’’. J Neurophysiol 61:
900–917.
Berger D, Warren D, Normann R, Arieli A, Grün S (2007) Spatially organized
spike correlation in cat visual cortex. Neurocomputing 70: 2112–2116.
Maldonado P, Babul C, Singer W, Rodriguez E, Berger D, et al. (2008)
Synchronization of neuronal responses in primary visual cortex of monkeys
viewing natural images. J Neurophysiol 100: 1523–1532.
Ito H, Maldonado PE, Gray CM (2010) Dynamics of stimulus-evoked spiketiming correlations in the cat lateral geniculate nucleus. J Neurophysiol 104:
3276–3292.
Ahissar E, Vaadia E, Ahissar M, Bergman H, Arieli A, et al. (1992)
Dependence of cortical plasticity on correlated activity of single neurons and on
behavioral context. Science 257: 1412–1415.
Vaadia E, Haalman I, Abeles M, Bergman H, Prut Y, et al. (1995) Dynamics ofneuronal interactions in monkey cortex in relation to behavioral events. Nature
373: 515–518.
Ishikane H, Gangi M, Honda S, Tachibana M (2005) Synchronized retinal
oscillations encode essential information for escape behavior in frogs. Nat
Neurosci 8: 1087–1095.
Fujisawa S, Amarasingham A, Harrison MT, Buzsáki G (2008) Behaviordependent short-term assembly dynamics in the medial prefrontal cortex. Nat
Neurosci 11: 823–833.

PLoS Computational Biology | www.ploscompbiol.org

20.
21.

22.

23.

24.
25.
26.
27.
28.
29.
30.
31.
32.

33.
34.
35.

36.
37.

25

Kilavik BE, Roux S, Ponce-Alvarez A, Confais J, Grün S, et al. (2009) Longterm modi–cations in motor cortical dynamics induced by intensive practice.
J Neurosci 29: 12653–12663.
Sakurai Y (1996) Hippocampal and neocortical cell assemblies encode memory
processes for different types of stimuli in the rat. J Neurosci 16: 2809–2819.
Steinmetz PN, Roy A, Fitzgerald PJ, Hsiao SS, Johnson KO, et al. (2000)
Attention modulates synchronized neuronal firing in primate somatosensory
cortex. Nature 404: 187–190.
Sakurai Y, Takahashi S (2006) Dynamic synchrony of firing in the monkey
prefrontal cortex during working-memory tasks. J Neurosci 26:
10141–10153.
Takahashi S, Sakurai Y (2009) Sub-millisecond firing synchrony of closely
neighboring pyramidal neurons in hippocampal CA1 of rats during delayed
non-matching to sample task. Front Neural Circuits 3: 9.
Brown EN, Kass RE, Mitra PP (2004) Multiple neural spike train data analysis:state-of-the-art and future challenges. Nat Neurosci 7: 456–461.
Grün S (2009) Data-driven significance estimation for precise spike correlation.
J Neurophysiol 101: 1126–1140.
Grün S, Rotter S, eds. (2010) Analysis of Parallel Spike Trains. New York:
Springer.
Perkel DH, Gerstein GL, Moore GP (1967) Neuronal spike trains andstochastic point processes. II. simultaneous spike trains. Biophys J 7: 419–440.
Gerstein GL, Perkel DH (1969) Simultaneously recorded trains of action
potentials - analysis and functional interpretation. Science 164: 828–830.
Gerstein GL, Kiang NYS (1960) An approach to the quantitative analysis ofelectrophysiological data from single neurons. Biophys J 1: 15–28.
Shimazaki H, Shinomoto S (2007) A method for selecting the bin size of a time
histogram. Neural Comput 19: 1503–1527.
Shimazaki H, Shinomoto S (2010) Kernel bandwidth optimization in spike rateestimation. J Comput Neurosci 29: 171–182.
Grün S, Diesmann M, Aertsen A (2002) Unitary events in multiple singleneuron spiking activity: I. detection and significance. Neural Comput 14:
43–80.
Grün S, Diesmann M, Aertsen A (2002) Unitary events in multiple singleneuron spiking activity: II. nonstationary data. Neural Comput 14: 81–119.
Staude B, Rotter S, Grün S (2010) Cubic: cumulant based inference of higher-ordercorrelations in massively parallel spike trains. J Comput Neurosci 29: 327–350.
Staude B, Grün S, Rotter S (2010) Higher-order correlations in non-stationary
parallel spike trains: statistical modeling and inference. Front Comput Neurosci
4: 16.
Chornoboy ES, Schramm LP, Karr AF (1988) Maximum likelihood
identification of neural point process systems. Biol Cybern 59: 265–275.
Brown EN, Frank LM, Tang D, Quirk MC, Wilson MA (1998) A statisticalparadigm for neural spike train decoding applied to position prediction fromensemble firing patterns of rat hippocampal place cells. J Neurosci 18:
7411–7425.

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

38.

39.

40.

41.

42.
43.

44.
45.
46.
47.
48.
49.
50.
51.
52.
53.

54.

55.

56.

57.

58.

59.
60.

61.

62.

63.

64.
65.
66.
67.

68.
69.

70.

Stevenson IH, Rebesco JM, Hatsopoulos NG, Haga Z, Miller LE, et al. (2009)
Bayesian inference of functional connectivity and network structure fromspikes. IEEE Trans Neural Syst Rehabil Eng 17: 203–213.
71. Truccolo W, Hochberg LR, Donoghue JP (2010) Collective dynamics inhuman and monkey sensorimotor cortex: predicting single neuron spikes. Nat
Neurosci 13: 105–111.
72. Chen Z, Barbieri R, Brown EN (2010) State-space modeling of neural spiketrain and behavioral data. In: Oweiss KG, ed. Statistical Signal Processing for
Neuroscience and Neurotechnology. Burlington: Elsevier Science. pp 175–218.
73. Akaike H (1980) Likelihood and the bayes procedure. Trabajos deEstadı́stica y de Investigación Operativa 31: 143–166.
74. Jeffreys H (1961) Theory of Probability. Oxford: Clarendon Press.
75. Good IJ (1985) Weight of Evidence: A brief survey. In: Bernardo JM, De
Groot MH, Lindley DV, Smith AFM, eds. Bayesian statistics 2. New York:North Holland. pp 249–269.
76. Kass RE, Raftery AE (1995) Bayes factors. J Am Stat Assoc 90: 773–795.
77. Shimazaki H, Amari SI, Brown EN, Grün S (2009) State-space analysis on
time-varying correlations in parallel spike sequences. In: Proceedings of IEEE
International Conference on Acoustics, Speech, and Signal Processing, 2009,
ICASSP-2009. pp 3501–3504.
78. Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood fromincomplete data via the EM algorithm. J Roy Stat Soc B Met 39: 1–38.
79. Schwarz G (1978) Estimating the dimension of a model. Ann Stat 6: 461–464.
80. Rissanen J (2009) Information and Complexity in Statistical Modeling. New
York: Springer.
81. Shimodaira H (1994) A new criterion for selecting models from partiallyobserved data. In: Cheeseman P, Oldford RW, eds. Selecting models from data
- artificial intelligence and statistics IV. New York: Springer-Verlag. pp 21–30.
82. Shlens J, Field GD, Gauthier JL, Greschner M, Sher A, et al. (2009) The
structure of large-scale synchronized –ring in primate retina. J Neurosci 29:
5022–5031.
83. Louis S, Gerstein G, Grün S, Diesmann M (2010) Surrogate spike train
generation through dithering in operational time. Front Comput Neurosci 4:
127.
84. Riehle A, Grammont F, Diesmann M, Grün S (2000) Dynamical changes andtemporal precision of synchronized spiking activity in monkey motor cortex
during movement preparation. J Physiol Paris 94: 569–582.
85. Bohte SM, Spekreijse H, Roelfsema PR (2000) The effects of pair-wise andhigher order correlations on the firing rate of a post-synaptic neuron. NeuralComput 12: 153–179.
86. Staude B, Grün S, Rotter S (2010) Higher order correlations and cumulants.
In: Grün, Rotter, eds. Analysis of parallel spike trains Springer New York. pp
253–280.
87. Jenison RL, Reale RA (2004) The shape of neural dependence. NeuralComput 16: 665–672.
88. Berkes P, Wood F, Pillow J (2009) Characterizing neural dependencies withcopula models. In: Koller D, Schuurmans D, Bengio Y, Bottou L, eds.
Advances in NIPS 21. pp 129–136.
89. Onken A, Grünewälder S, Munk MHJ, Obermayer K (2009) Analyzing shortterm noise dependencies of spike-counts in macaque prefrontal cortex using
copulas and the flashlight transformation. PLoS Comput Biol 5: e1000577.
90. Paninski L (2005) Log-concavity results on gaussian process methods forsupervised and unsupervised learning. In: Saul LK, Weiss Y, Bottou L, eds.
Advances in NIPS 17. pp 1025–1032.
91. Amari S (1998) Natural gradient works efficiently in learning. Neural Comput
10: 251–276.
92. Murray I, Ghahramani Z (2004) Bayesian learning in undirected graphicalmodels: approximate MCMC algorithms. In: Proceedings of the 20th Annual
Conference on Uncertainty in Artificial Intelligence. UAI-2004. Arlington:
AUAI Press. pp 392–399.
93. Schaub MT, Schultz SR (2011) The Ising decoder: reading out the activity oflarge neural ensembles. J Comput Neurosci 32: 101–118. 101007/s10827-0110342-z.
94. Shinomoto S, Kim H, Shimokawa T, Matsuno N, Funahashi S, et al. (2009)
Relating neuronal firing patterns to functional differentiation of cerebral
cortex. PLoS Comput Biol 5: e1000433.
95. Daley D, Vere-Jones D (1988) An Introduction to the Theory of Point
Processes. New York: Springer-Verlag.
96. Kass RE, Ventura V, Brown EN (2005) Statistical issues in the analysis ofneuronal data. J Neurophysiol 94: 8–25.
97. Brody CD (1999) Correlations without synchrony. Neural Comput 11:
1537–1551.
98. Czanner G, Eden UT, Wirth S, Yanike M, Suzuki WA, et al. (2008) Analysis of
between-trial and within-trial neural spiking dynamics. J Neurophysiol 99:
2672–2693.
99. Grün S, Diesmann M, Grammont F, Riehle A, Aertsen A (1999) Detecting
unitary events without discretization of time. J Neurosci Methods 94: 67–79.
100. Marre O, Boustani SE, Frégnac Y, Destexhe A (2009) Prediction of
spatiotemporal patterns of neural activity from pairwise correlations. Phys
Rev Lett 102: 138101.
101. Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM (2001) The timerescaling theorem and its application to neural spike train data analysis. NeuralComput 14: 325–346.

Truccolo W, Eden UT, Fellows MR, Donoghue JP, Brown EN (2005) A point
process framework for relating neural spiking activity to spiking history, neuralensemble, and extrinsic covariate effects. J Neurophysiol 93: 1074–1089.
Okatan M, Wilson MA, Brown EN (2005) Analyzing functional connectivity
using a network likelihood model of ensemble neural spiking activity. NeuralComput 17: 1927–1961.
Pillow JW, Shlens J, Paninski L, Sher A, Litke AM, et al. (2008) Spatiotemporal correlations and visual signalling in a complete neuronal population.
Nature 454: 995–999.
Kim S, Putrino D, Ghosh S, Brown EN (2011) A granger causality measure forpoint process models of ensemble neural spiking activity. PLoS Comput Biol 7:
e1001110.
Kulkarni JE, Paninski L (2007) Common-input models for multiple neuralspike-train data. Network: Comp Neural Sys 18: 375–407.
Vidne M, Ahmadian Y, Shlens J, Pillow JW, Kulkarni J, et al. (2011) Modelingthe impact of common noise inputs on the network activity of retinal ganglion
cells. J Comput NeurosciIn press. DOI: 101007/s10827-011-0376-2.
Amari SI (2001) Information geometry on hierarchy of probability distributions. IEEE T Inform Theory 47: 1701–1711.
Nakahara H, Amari SI (2002) Information-geometric measure for neuralspikes. Neural Comput 14: 2269–2316.
Ackley D, Hinton G, Sejnowski T (1985) A learning algorithm for boltzmann
machines. Cognitive Sci 9: 147–169.
Amari S, Kurata K, Nagaoka H (1992) Information geometry of boltzmann
machines. IEEE T Neural Networ 3: 260–271.
Amari S (1972) Learning patterns and pattern sequences by self-organizing netsof threshold elements. IEEE Trans Comput 21: 1197–1206.
Hopfield J (1982) Neural networks and physical systems with emergent
collective computational abilities. Proc Natl Acad Sci U S A 79: 2554–2558.
Amari S, Nagaoka H (2000) Methods of Information Geometry. Providence:
American Mathematical Society.
Amari S (2009) Measure of correlation orthogonal to change in firing rate.
Neural Comput 21: 960–972.
Schneidman E, Still S, Berry MJ, Bialek W (2003) Network information andconnected correlations. Phys Rev Lett 91: 238701.
Montani F, Ince RAA, Senatore R, Arabzadeh E, Diamond ME, et al. (2009)
The impact of high-order interactions on the rate of synchronous discharge
and information transmission in somatosensory cortex. Philos Transact A Math
Phys Eng Sci 367: 3297–3310.
Roudi Y, Nirenberg S, Latham PE (2009) Pairwise maximum entropy models
for studying large biological systems: when they can work and when they can’t.
PLoS Comput Biol 5: e1000380.
Schneidman E, Berry MJ, Segev R, Bialek W (2006) Weak pairwise
correlations imply strongly correlated network states in a neural population.
Nature 440: 1007–1012.
Shlens J, Field GD, Gauthier JL, Grivich MI, Petrusca D, et al. (2006) The
structure of multineuron firing patterns in primate retina. J Neurosci 26:
8254–8266.
Tang A, Jackson D, Hobbs J, Chen W, Smith JL, et al. (2008) A maximum
entropy model applied to spatial and temporal correlations from cortical
networks in vitro. J Neurosci 28: 505–518.
Santos GS, Gireesh ED, Plenz D, Nakahara H (2010) Hierarchical interaction
structure of neural activities in cortical slice cultures. J Neurosci 30:
8720–8733.
Yu S, Yang H, Nakahara H, Santos GS, Nikolic D, et al. (2011) Higher-orderinteractions characterized in cortical activity. J Neurosci 31: 17514–17526.
Ohiorhenuan IE, Mechler F, Purpura KP, Schmid AM, Hu Q, et al. (2010)
Sparse coding and high-order correlations in fine-scale cortical networks.
Nature 466: 617–621.
Martignon L, Hasseln HV, Grün S, Aertsen A, Palm G (1995) Detecting
higher-order interactions among the spiking events in a group of neurons. Biol
Cybern 73: 69–81.
Martignon L, Deco G, Laskey K, Diamond M, Freiwald W, et al. (2000)
Neural coding: higherorder temporal patterns in the neurostatistics of cell
assemblies. Neural Comput 12: 2621–2653.
Macke JH, Opper M, Bethge M (2011) Common input explains higher-ordercorrelations and entropy in a simple model of neural population activity. Phys
Rev Lett 106: 208102.
Gütig R, Aertsen A, Rotter S (2003) Analysis of higher-order neuronalinteractions based on conditional inference. Biol Cybern 88: 352–359.
Amari SI, Nakahara H, Wu S, Sakai Y (2003) Synchronous firing and higherorder interactions in neuron pool. Neural Comput 15: 127–142.
Long II J, Carmena J (2011) A statistical description of neural ensemble
dynamics. Front Comput Neurosci 5: 52.
Kass RE, Kelly RC, Loh WL (2011) Assessment of synchrony in multiple
neural spike trains using loglinear point process models. Ann Appl Stat 5:
1262–1292.
Smith AC, Brown EN (2003) Estimating a state-space model from point
process observations. Neural Comput 15: 965–991.
Eden UT, Frank LM, Barbieri R, Solo V, Brown EN (2004) Dynamic analysisof neural encoding by point process adaptive filtering. Neural Comput 16:
971–998.

PLoS Computational Biology | www.ploscompbiol.org

26

March 2012 | Volume 8 | Issue 3 | e1002385

State-Space Analysis of Spike Correlation

107. Koyama S, Paninski L (2010) Efficient computation of the maximum aposteriori path and parameter estimation in integrate-and-fire and moregeneral state-space models. J Comput Neurosci 89–105: 1–2.
108. De Jong P, Mackinnon MJ (1988) Covariances for smoothed estimates in state
space models. Biometrika 75: 601–602.
109. Harvey AC (1989) Forecasting Structural Time Series Models and the Kalman
Filter. Cambridge; New York: Cambridge University Press.
110. Durbin J, Koopman SJ (2001) Time Series Analysis by State Space Methods.
Oxford; New York: Oxford University Press.
111. Meng X, Rubin D (1991) Using EM to obtain asymptotic variance-covariancematrices: the SEM algorithm. J Am Stat Assoc 86: 899–909.

102. Oizumi M, Ishii T, Ishibashi K, Hosoya T, Okada M (2010) Mismatched
decoding in the brain. J Neurosci 30: 4815–4826.
103. Shumway R, Stoffer D (1982) An approach to time series smoothing andforecasting using the em algorithm. J Time Ser Anal 3: 253–264.
104. Smith AC, Frank LM, Wirth S, Yanike M, Hu D, et al. (2004) Dynamic
analysis of learning in behavioral experiments. J Neurosci 24: 447–461.
105. Kitagawa G (1987) Non-Gaussian state-space modeling of nonstationary timeseries. J Am Stat Assoc 82: 1032–1041.
106. Fahrmeir L (1992) Posterior mode estimation by extended kalman filtering for
multivariate dynamic generalized linear models. J Am Stat Assoc 87: 501–509.

PLoS Computational Biology | www.ploscompbiol.org

27

March 2012 | Volume 8 | Issue 3 | e1002385