Copyright 2009 IEEE. Published in the IEEE 2009 International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2009),
scheduled for April 19 - 24, 2009 in Taipei, Taiwan Personal use of this material is permitted. However, permission to reprint/republish
this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists,or to reuse any copyrighted component of this work in other works, must be obtained from the IEEE. Contact: Manager, Copyrights andPermissions / IEEE Service Center / 445 Hoes Lane / P.O. Box 1331 / Piscataway, NJ 08855-1331, USA. Telephone: + Intl. 908-562-3966.

STATE-SPACE ANALYSIS ON TIME-VARYING CORRELATIONSIN PARALLEL SPIKE SEQUENCES
Hideaki Shimazaki1 , Shun-ichi Amari1 , Emery N. Brown2,3 , and Sonja Grün1
1

Theoretical Neuroscience Group, RIKEN Brain Science Institute, Wako-shi, Saitama, Japan
2
Anesthesia and Critical Care, Massachusetts General Hospital, Boston, MA, USA
3
Harvard-MIT Division of Health Sciences and Technology, Cambridge, MA, USA
ABSTRACTA state-space method for simultaneously estimating timedependent rate and higher-order correlation underlying parallel spike sequences is proposed. Discretized parallel spikesequences are modeled by a conditionally independent multivariate Bernoulli process using a log-linear link function,
which contains a state of higher-order interaction factors. A
nonlinear recursive filtering formula is derived from a logquadratic approximation to the posterior distribution of thestate. Together with a fixed-interval smoothing algorithm,
time-dependent log-linear parameters are estimated. The
smoothed estimates are optimized via EM-algorithm suchthat their prior covariance matrix maximizes the expectedcomplete data log-likelihood. In addition, we perform model
selection on the hierarchical log-linear state-space models
to avoid over-fitting. Application of the method to simultaneously recorded neuronal spike sequences is expected tocontribute to uncover dynamic cooperative activities of neurons in relation to behavior.
Index Terms— State space methods, Point processes,
Generalized linear model, Correlation, Information geometry
1. INTRODUCTIONClassical studies in neurophysiology are based on the ideathat stimulus information is encoded in the firing rates of single neurons. Alternatively, precise spike coordination is discussed as an indication of coordinated network activity inform of a cell assembly [1] relevant for information processing. Activation of a cell assembly predicts higher-order correlation (HOC) between the spiking activities of its member
neurons [2]. Supportive evidence for this concept was provided by existence of excess synchrony among neuronal spiking activities occurring dynamically in relation to behavioral
context [3, 4, 5]. However, available approaches for correlation analysis do not allow for identifying time-dependent
HOCs to trace active assemblies.
To characterize higher-order interaction, the log-linear
model is an useful tool because it provides a well-definedmeasure of correlation based on information geometry [6].

978-1-4244-2354-5/09/$25.00 ©2009 IEEE

3501

Although its natural (canonical) parameters are not orthogonal, HOCs can be extracted in a quasi-orthogonal mannerfrom a mixture of dually affine Riemannian coordinates [7, 8].
Former studies performed a regression analysis on parallel
spike trains using either a log-linear model [7, 8, 9], or the
log-linear model considering up to pairwise interaction only
(maximum entropy model) [10] to characterize neuronal data.
The existing approaches, however, assume stationarity, a condition that is typically not fulfilled in neuronal spike datafrom awake behaving animals.
The state-space method [11] was suggested as a framework to model a time-dependent system by representing its
parameters (states) as a Markov process. It allows to estimate a filtered/smoothed posterior distribution of the timedependent state conditional on observed data. The approximation method for non-Gaussian point process observation
in a recursive filtering algorithm was successfully applied toneuronal spike data [12, 13, 14]. Existing state-space models
incorporated ensembles of spike histories into the univariate
spike response model [14, 15], however, without consideringcorrelations between the spike sequences.
In this contribution, we provide a method for estimating
the dynamics of HOCs by combining the log-linear modelwith a state-space analysis. For that estimation, we combined
a fixed-interval smoother [11] with a nonlinear recursive filtering algorithm which we derived from a log-quadratic approximation of the posterior distribution. To obtain the most
predictive model, we compared goodness-of-fit of hierarchical log-linear state-space models with different order of interactions by using the Akaike’s Bayesian information criterion(ABIC) [16]. While inclusion of increasingly higher-order interaction terms improves the model accuracy, the estimationof higher-order parameters may suffer large variance due tothe paucity of synchronous spikes in the data. This trade-off is
optimally resolved with the model that minimizes the ABIC.
The model complexity is thus selected based on the sample
size as well as the prominence of higher-order structure.
An earlier version of this paper was presented in abstract
form [Shimazaki, Brown, Grün, Statistical Analysis of Neuronal Data 4, Pittsburgh, May 2008].

ICASSP 2009

2. METHODS
2.1. Log-linear model of higher-order interactions
We model N simultaneous spike sequences as a multivariate Bernoulli process. The processes, repeated over n trials, are discretized into T bins of bin size Δ. Let Xt,l =t,l) be N -tuple binary variables at t-th bin(X1t,l , X2t,l , . . . , XN
of l-th trial, where 1 denotes a bin filled with one or more
spikes, and 0 denotes no spike. Let x = (x1 , x2 , . . . , xN )represent 2N spike patterns, and Ωk be a k-subset of N elements: Ω1 = {1, 2, . . . , N } , Ω2 = {12, 13, . . .}, Ω3 =
{123, 124, . . .}, etc. We define the interaction terms acrossthe variables as Fi (x) = xi , Fij (x) = xi xj , Fijk (x) =xi xj xk (1 ≤ i < j < k ≤ N ), etc. The log-linear model up
to the r-th order interaction (r ≤ N ) is defined as [6](r)θit Fi (x)(1)
log pt (x) = −ψt +is the first order autoregressive parameters. The initial value
obeys θ1 ∼ N (μ, Σ); μ (d × 1 matrix) and Σ (d × d matrix).
We denote the set of hyper-parameters w = [F, Q, μ, Σ].
The joint probability distribution is obtained by combining the observation equation (Eq. 4) and the state equation(Eq. 5). The optimization principle for w is chosen as the
maximization of the marginal log-likelihood,(6)l (w) = log p (y1:T , θ 1:T |w) dθ 1:T .
We make use of the expectation-maximization (EM) algorithm known to efficiently combining the estimation of theposterior density and the optimization of the hyper-parameter.
Instead of Eq. 6, the method maximizes its lower bound,
Q (w|w∗ ) = E [log p (y1:T , θ 1:T |w) |y1:T , w∗ ]i∈{Ω1 ,...,Ωr }with ψt being a normalization parameter. θit is a natural pa-

rameter referring to an interaction among variables indicatedttin theindex i. Let θ t = [θ1t , . . . , θ12, . . . , θ1···r, . . .] be arN
d = k=1 k dimensional vector. Our goal is to obtain theestimate of θ t from the observed parallel spike sequences.
The log-linear model belongs to the exponential family
which defines a dually flat space [6]. Its e-affine coordinate isgiven as the natural parameter θ t , its m-affine coordinates astt, . . . , η1···r, . . .] ,the expectation parameter η t = [η1t , . . . , η12
where
(2)
ηit = E [Fi (x) |θ t ]for i ∈ {Ω1 , . . . , Ωr }. An efficient estimator of ηit is the synchrony rate defined for the t-th bin asyit =n
1   t,l 
Fi X
n

(3)

l=1

for i ∈ {Ω1 , . . . , Ωr }. Synchrony rates up to the r-th order,tt, . . . , y1···r, . . .] , constitute sufficient statisyt = [y1t , . . . , y12
tic of the log-linear model up to the r-th order interaction.
2.2. State-space analysis on the log-linear model
Assuming conditionally independent observations, the likelihood of observed spike sequences is given as
p (y1:T |θ 1:T ) =T

exp [n (yt θ t − ψt )] ,

(4)

=n

√
√
−d log 2π − log det Σ1− E (θ 1 − μ) Σ−1 (θ 1 − μ)2√
− (T − 1) d log 2π − (T − 1) log det Q
T−for t = 2, . . . , T . A random vector ξ t (d × 1 matrix) is drawn
from a zero-mean multivariate normal distribution with a covariance matrix Q (d × d matrix). A matrix F (d × d matrix)

3502

1

E (θ t − Fθ t−1 ) Q−1 (θ t − Fθ t−1 ) , (7)2 t=2by alternating expectation (E) and maximization (M) steps. In
the E-step, we obtain the expected values with respect to θ t
in Eq. 7 using w∗ . In the M-step, the hyper-parameter w thatmaximizes Eq. 7 is obtained. The obtained w is then used forw∗ in the next E-step. The details of each step are given as
follows.
E-step Mean and covariance of the one-step prediction
density p (θ t |y1:t−1 ) for the Gaussian state transition is given
as [11, 12, 13]
θ t|t−1 = Fθ t−1|t−1 ,Wt|t−1 = FWt−1|t−1 F + Q.

(8)
(9)

Here θ t|s and Wt|s are the conditional mean E [θ t |y1:s ] andcovariance E[(θ t − θ t|s )(θ t − θ t|s )|y1:s ]. The filter distribution is given by the Bayes’ theorem asp (yt |θ t , y1:t−1 ) p (θ t |y1:t−1 )p (yt |y1:t−1 )∝ exp[n (yt θ t − ψt ) −1 1
θ t − θ t|t−1 ].
− θ t − θ t|t−1 Wt|t−1
2

p (θ t |y1:t ) =with y1:T = {y1 , y2 , . . . , yT } and θ 1:T = {θ 1 , θ 2 , . . . , θ T }.
Our prior assumption is the following state equation(5)(yt Eθ t − Eψt )

t=1

t=1

θ t = Fθ t−1 + ξ t ,

T


(10)

Here we make a Gaussian approximation to the posterior by Laplace’s method: the posterior mean, θ t|t , is ob∂
log p (θ t |y1:t ) = 0; the covariance is given
tained from ∂θt2

as −[ ∂θ∂t ∂θ log p (θ t |y1:t ) |θt =θt|t ]−1 . Note that the firsttTrial, lA

Estimated log-linear parameters

0

θ3t

−2

t

θ2

−4

Variable X t,l
1

−6
6

5
20
50

t

θ23

0
t
θ12

t,l

Variable X 2

100

Trial, l

t

θ1

50

100

Trial, l

B

Spike Sequences
5
20

−6

5
20

t

θ13
t

θ123

6

50
0

t,l

Variable X 3

100
1

Bin t [Δ]

−4

500

1

Bin t [Δ]500

Fig. 1. (A) Simulated parallel spike sequences, Xt,l . According to a time-dependent formulation of the log-linear model (seedashed lines in B for the model parameters), N = 3 parallel spike sequences were simulated repeatedly for n = 100 trials(duration: T = 500 bins of width Δ = 1[ms]). Each of the 3 panels show the spike events for the 3 variables X1t,l , X2t,l , X3t,l(t = 1, . . . , T and l = 1, . . . , n). Synchronous spikes across all 3 simultaneous sequences (detection in bins of 1ms width) aremarked with circles. Horizontal lines indicate trial numbers l = 5, 20, 50, 100. (B) Smoothed estimates of the log-linear parameters, θ t|T . The different panels show the smoothed estimates (solid lines) of the log-linear parameters θit of different orders(i ∈ {Ω1 , Ω2 , Ω3 }) derived from the data shown in A. The gray band indicates the 99% credible interval of the marginalized
posterior density of the state (see text). The dashed lines show the true time-dependent parameters used for the generation of
the spike sequences.
derivative of the moment generating function provides the ext
pectation parameter ∂ψ∂θ t = η t . The second derivative is themethod of De Jong and Mackinnon [17]:Wt−1,t|T = E[ (θ t−1 − θ t−1|T )(θ t − θ t|T ) y1:T ](16)= At−1 Wt|T .

2

ψt
t
t
t
Fisher metric ∂θ∂ t ∂θ
 = Gt . Here (Gt )ij = ηωi ∪ωj −ηωi ηωj ,twhere ωi is an index of the ith element of η t . Hence we obtain the equations for the filtered mean and covariance asθ t|t = θ t|t−1 + nWt|t−1 (yt − η t|t ),
−1
−1
Wt|t
= Wt|t−1+ nGt .

(11)
(12)

Eq. 11 is solved by the Newton-Raphson method. Together
with the one-step prediction formulae Eqs. 8 and 9, we evaluate the nonlinear recursion formulae Eqs. 11 and 12 for t =2, . . . , T , by using the initial values θ 1|0 = μ and W1|0 = Σ.
Starting from θ T |T and WT |T obtained by the filtering
algorithm, the following fixed-interval smoothing algorithm[11, 12, 13] provides the smoothed mean and covariance. For
t = T − 1, T − 2, . . . 2, we compute(13)
θ t|T = θ t|t + At θ t+1|T − θ t+1|t , 
Wt|T = Wt|t + At Wt+1|T − Wt+1|t At .
(14)

M-step Given θ t|T , Wt|T , Wt−1,t|T , the update of w is
given by the following equations. The covariance matrix Q
and the auto-regressive parameter F are updated according to
Q=

T
1 
[Wt|T − Wt−1,t|T F
T − 1 t=2

− FWt−1,t|T + FWt−1|T F ]

+

T


1 
θ t|T − Fθ t−1|T θ t|T − Fθ t−1|T (17)T − 1 t=2

and
T


F=

Wt−1,t|T + θ t|T θ t−1|T

t=2

·

T


−1

Wt−1|T + θ t−1|T θ t−1|T

.

(18)

t=2

with
−1
At = Wt|t F Wt+1|t.

(15)

The lag-one covariance smoother Wt−1,t|T is obtained by the3503

The mean of the initial distribution is updated with μ = θ 1|T .
The covariance matrix Σ is not updated. Instead, we use anominal diagonal matrix as Σ.

Table 1. The ABIC of the r-th order log-linear model applied
to data of different number of trials n (data shown in Fig.1-A).
The asterisk indicates the model that minimizes the ABIC.
n=5
n = 20n = 50n = 100n = 200

r=1
2085∗
7913
19263
38478
76683

r=2
2097
7698∗
18542
36811
73330

r=3
2144
7728
18540∗
36781∗
73238∗

4. CONLUSION
We developed a method for identifying the time-varying
higher-order correlation structure in parallel spike sequences.
Its application to simultaneous recordings of neuronal activity
is expected to provide us with new insights into the dynamics of assembly activities, their composition and behavioral
relevance.
Acknowledgements This work was supported in part by JSPS Research Fellowships for Young Scientists (HS), R01 MH59733 / DP1
OD 003646 (EB), and RIKEN Strategic Programs for R&D (SG).
5. REFERENCES

2.3. Model Selection
Let us define a log-linear state-space model with up to ther-th order interaction by Er . Since Er is a submanifold of
Er+1 , the log-linear models naturally form hierarchical structure, E1 ⊂ E2 ⊂ · · · ⊂ EN . Comparison of the hierarchicalmodels is performed by computing the ABIC [16]:ABIC = −2l (w) + 2 dim w.

(19)

The marginal log-likelihood is computed by using the onestep prediction formula for the mean Eq. 8 as l(w) =TT
t=1 log p(yt |y1:t−1 , w) ≈ n
t=1 (yt θ t|t−1 − ψt|t−1 ).
3. RESULTS
We applied the method to N = 3 parallel spike sequenceswith n = 100 repetitions (Fig.1-A), generated by a timedependent log-linear model (Fig.1-B, dashed lines). Figure
1-B (solid lines) shows the estimates of the log-linear parameters from the data in Figure 1-A. The gray band is the
99% credible interval of the marginalized posterior distribution p (θit | y1:T ) for i ∈ {Ω1 , Ω2 , Ω3 }. Note that parametert
indicates triplewise correlation, i.e. excess synchrony
θ123
compared to expectation given pairwise correlations.
While the full model (E3 ) is an unbiased model for the
triple binary sequences, the fitted parameter for the triplewiset
may suffer large variance (over-fitting) due
correlation θ123
to the paucity of triplet synchrony. Thus a submodel withr < 3 can be close to the generative model in terms of the
Kullback-Leibler risk. To validate that the inclusion of thetriplewise correlation improves the goodness-of-fit, we computed the ABICs for the hierarchical models (r = 1, 2, 3) asshown in Table 1. To test the influence of the sample sizeof the data upon the model selection, we varied the numberof trials n used to fit the hierarchical log-linear models. For
small number of trials (n = 5), the model without correlation structure was selected. With increasing trial numbers,
models with increasing correlation orders were selected. For
n = 50, 100, and 200, the full model predicts best.

3504[1] D. O. Hebb, The Organization of Behavior: A Neuropsychological Theory, Wiley, New York, 1949.
[2] S. Grün et al. “Impact of higher-order correlations on coincidence distributions of massively parallel data.,” in Lect Notes
Comput Sc, vol. 5286, pp. 96-114, Nov 2008[3] E. Vaadia et al. “Dynamics of neuronal interactions in monkey
cortex in relation to behavioral events,” Nature, vol. 373, no.
6514, pp. 515–518, Feb 1995.
[4] A. Riehle et al. “Spike synchronization and rate modulationdifferentially involved in motor cortical function,” Science, vol.
278, no. 5345, pp. 1950–1953, Dec 1997.
[5] S. Grün et al.“Unitary events in multiple single-neuron spiking
activity: II. nonstationary data.,” Neural Comput, vol. 14, no.
1, pp. 81–119, Jan 2002.
[6] S.-I. Amari and H. Nagaoka, Methods of Information Geometry , American Mathematical Society, 2000.
[7] S.-I. Amari, “Information geometry on hierarchy of probability distributions,” IEEE Trans. Inf. Theory, vol. 47, no. 5, pp.
1701–1711, Jul 2001.
[8] H. Nakahara and S.-I. Amari, “Information-geometric measurefor neural spikes.,” Neural Comput, vol. 14, no. 10, pp. 2269–
2316, Oct 2002.
[9] L. Martignon et al. “Detecting higher-order interactions among
the spiking events in a group of neurons.,” Biol Cybern, vol.
73, no. 1, pp. 69–81, Jun 1995.
[10] E Schneidman et al. “Weak pairwise correlations imply
strongly correlated network states in a neural population.,” Nature, vol. 440, no. 7087, pp. 1007–1012, Apr 2006.
[11] G. Kitagawa and W. Gersch, Smoothness Priors Analysis ofTime Series, Springer-Verlag, New York, 1996.
[12] E. N. Brown et al. “A statistical paradigm for neural spike train
decoding applied to position prediction from ensemble firing
patterns of rat hippocampal place cells.,” J Neurosci, vol. 18,
no. 18, pp. 7411–7425, Sep 1998.
[13] A. C. Smith and E. N. Brown, “Estimating a state-space modelfrom point process observations.,” Neural Comput, vol. 15, no.
5, pp. 965–991, May 2003.
[14] W. Truccolo et al. “A point process framework for relatingneural spiking activity to spiking history, neural ensemble, andextrinsic covariate effects.,” J Neurophysiol, vol. 93, no. 2, pp.
1074–1089, Feb 2005.
[15] M. Okatan et al. “Analyzing functional connectivity using a
network likelihood model of ensemble neural spiking activity.,” Neural Comput, vol. 17, no. 9, pp. 1927–1961, Sep 2005.
[16] H. Akaike, “Likelihood and the Bayes procedure,” in Bayesian
Statistics, University Press, pp. 143–166, 1980.
[17] P. De Jong and M. J. Mackinnon, “Covariances for smoothed
estimates in state space models,” Biometrika, vol. 75, pp. 601–
602, Sep 1988.