II-Detection of non-stationary higher-order spike correlation11

2

Hideaki Shimazaki , Shun-ichi Amari , Emery N. Brown , and Sonja Grün62

1

1Theoretical Neuroscience Group, RIKEN Brain Science Institute, Wako-shi Saitama 351-0198 Japan
2
Massachusetts General Hospital, Depart. of Brain and Cognitive Sciences MIT, Harvard-MIT Division of Health Sciences and Technology, Cambridge MA USA
Contact: shimazaki@brain.riken.jp

The log-linear model provides a well-defined measure of higher-order correlation basedon information geometry [5].

Synchrony between two neurons were appropriately attributed toeither a chance coincidence or a pairwise correlation.
Constant rates + Time-dependent correlation

The Log-linear model of two binary variables

Neuron 1x1
0
10x2001

P-coordinate
p00

11

p11

θ12ψ ( θ) : cumulant generating function

log p 00 = −ψ

θ1( )
2-4

η2
η14

log p (x | θ ) = −ψ + ∑ θi xi + ∑ θij xi x j + θ123 x1 x2 x3

3 Neuronsi< jiN-th order correlationlog p (x | θ ) = −ψ + ∑ θi xi + ∑ θij xi x j + ⋯ + θ1⋯N x1 x2 x3 ⋯ xN

N Neurons

0.4

000.4

Time [s]00.4Time [s]03 Neurons
Time-dependence of triplewise correlation was estimatedwith an optimized covariance matrix and autoregressive parameters.

log p (x | θ) = −ψ + ∑ θi xi + ∑ θij xi x j + ∑ θijk xi x j xk
i< j

θ3(3)

diag(Q)

5

yt+1

yt

Observation

ξ t ∼ N(0, Q )θt = F θt −1 + ξ tPriorp (y1:T | θ 1:T ) p (θ1:T | w )p (y1:T | w ) EvidenceGoal:θ1 ∼ N(μ, Σ)

wheretθt = θ1t , …,θ Nt ,θ 12t ,θ 13t , …,θ 123, …,θ 1t⋯N  ′ty123

M-step: Hyper-parameter optimization, given the posterior.
Optimize the hyper-parameters by maximizes
the lower bound of evidence (Q-function).

{# of joint spikes of Neuron 1,2,3 occured at time t} / nn: the number of trials.

Filtering/smoothing and hyper-parameter optimization
E-step: Filtering/Smoothing
One-step predictioin

θt|t −1 = Fθt −1|t −1

Covariance

Wt|t −1 = FWt −1|t −1F′ + Q

observation
Prediction

One-step prediction distribution

Recursive Nonlinear Filtering O n e -step p red iction

What is the Higher Order Correlation?

t-1

θt|t = θt|t −1 + nWt|t −1 y t − ηt|t Expectation−1t |tW

Covariance1−1t |t −1=Wt

Recursive Filter distribution

O b se rva tion (e rro r)+ nI FO n e -step p red iction

Likelihood
FilterF ishe r In fo rm a tionp (θt | y1:t ) =Priorp (y t | θt , y1:t −1 ) p (θt | y1:t −1 )p (y t | y1:t −1 )(Bayes’ rule)

Fixed Interval Smoothing
Expectation

θt|T = θt|t + At θ t +1|T − θ t +1|t Covariance

Wt|T = Wt|t + At Wt +1|T − Wt +1|t  At′T1∝ exp  n (y Tt θt −ψ t )− (θt −θ t|t −1 ) Wt|−t −11 (θt − θ t|t −1 )
2At = Wt|t F T Wt +−1|1t

Neurons are independent.

O n e -step p red iction(Normal distrib ution )O b se rva tion at time t(Log-linear distribution)

Synchrony may appear by chance.
32
B

Neuron 1 and 2 are positively correlated.
Triplets may appear by chance.

Neuron 1

C

Neuron 2

Three neurons are connected with pairwise correlations.
Triplets may appear by chance.

Neuron 3
D
triplet by chanceexcess triplet

F1

A triplewise correlation is added.
Ecess triplets are generated.

Positive higher order correlation (HOC) indicates excess synchrony that can not be
explained by the lower order correlations.

M-step: Hyper-Parameter Optimizationl +1lll
Q-funtion q (w | w )= E log p (θ1:T , y1: T | w )| y1: T , w 

Optimization of Initial ValuesΣ = W1|T + (θ1|T − μ)(θ1|T − μ)′μ = θ 1|T

Optimization of covariance matrix
T(T − 1)Q = ∑ [Wt|T − Wt −1,t|T F ′ − FWt −′ 1,t|T + FWt −1|T F ′t =20.5(3)θ123

0+ (θt|K − F θt −1|T )(θ t|K − F θt −1|T )′ ]Optimization of AR(1) parameters

T T′F =  ∑ (Wt −1,t|T + θt|T θt −1|T )  ∑ (Wt −1|T + θ t −1|T θt′−1|T )
 t =2  t =2
−1q (wl +1|w )lq (wItterative Maximization of Q-funtion

θt = Fθ t −1 + G1ut1 + G 2ut2 + ξ t

State Equation

0.02
0−0.02−0.04−0.06

Stim 1

Stim 1:

ut1 1 at the onset of a signal; 0 otherwise.

12401

Triplewise Correlation

0Time [s]20.5Time [s]3[First column] 100 trials of simulated parallel spike train. Each pattern of spikes was counted with 1[ms] resolution. Red dots indicate synchronous spikes (triplets).
[Second column] Log-linear parameters. Shaded area indicates 99% confidence interval. [Third column] Optimized hyper-parameters.

Model Selection

401

Stim 2Ĝ1 Ĝ 22

Stim 2Estimated efect of stimuli to the interaction parameters
Individual

PairTriple

0.4

0

To validate that the inclusion of the triplewise correlation improves
the goodness-of-fit, we computed ABICs for the hierarchical state-space models.
The Akaike Bayesian information criterion (ABIC) is computed as
ABIC = −2 log p(y1:T |w) + 2 dim w,

Filter distribution

0Adiag(F)3Filter distribution

Expectation

0.5

Incorporation of the cue signals into the model may be useful for detecting
assemblies related to behavior.

0.8

via Recursive filtering/smoothing algorithm

0.4

−0.01

20.5

0.3

Modeling external signals’ effects

E-step: Estimate of the posterior distribution, given hyper-parameters.

ty t =  y1t , …, y Nt , y12t , y13t , …, y123
, …, y1t⋯N  ′Sufficient statistics(Synchrony rate)
Steinmetz P.N., Roy A., Fitzgerald P.J., Hsiao, S.S., Johnson K.O. & Niebur E. Nature, 2000, 404, 187-190(in Salinas E. & Sejnowski T.J. Nat Rev Neurosci, 2001, 2, 539-550)

Dθ23400.2

0.1

Solution:w = [μ,Σ ,Q, F ]Hyper-parameterw∗ = arg max ∫ p ( y1:T , θ1:T | w )d θ1:T .
w

t =1Parameterp (θ1:T | y1:T , w )∗Obtain the Posteriorp (y1:T | θ1:T )= ∏ exp  n (yt′ θt −ψ t )

Observation Equation

0AR(1) Matrix(3)

3−2Time

It has been shown that spike correlation with ms precision i) is modulated in time andii) occurs at behaviorally relevant instances.

0.01

0θ126p (θ1:T | y1:T , w ) =TC(3)θ13( )0Q

0−5

Likelihood

Estimated Rates and Pairwise Correlations from 200 trials

Covariance Matrix−4i< j <kBayes’ Theorem

46Time [s]θ1(3)Neuron 2

Hidden State

3

Negatively Correlated−2

Triali< jTrial
θt+1

θtyt-1

Excess synchrony appeared with attention.

166 [Hz]

Log-linear parameters
0State-space Analysis

State Equation

Excess synchrony appeared with expectation.

28Positively Correlated

0.025Riehle A., Grun S., Diesmann M. & Aertsen A. Science, 1997, 278, 1950-1953(in Engel et. al. Nature Reviews Neuroscience, 2001, 2, 704-716)

True Rates and Pairwise Correlations

Neuron 1

Firing Rate

5log p (x | θ) = −ψ + ∑ θi xi + ∑ θij xi x j

4

Evidence of Dynamical Correlation

0.5

Time [s]Neuron 2

Precise spike coordination may appear due to the coordinated network activity.
Composition of the spike coordination may dynamically change.

00.4

Time [s]0Neuron 1

Posterior

Synchornous Spikes

0.5

7log p (x | θ) = −ψ + ∑ θi xi

3Time-dependent Network Activities

7-8-2

Pairwise correlation

To N-tuple parallel spike sequences, one can fit the sub-model with up to r-th order correlations:2556-7−4i< ji

θt-1

5-7

0

0

Tripplewise correlation

13

Neuron 8

Neuron 4

θ12( )

99% confidence interval−2Trial2

4-7

4Pairwise Correlation

22

State-Space Model of Parallel Spike Trains

Parallel Spike Sequences

43

excess doublet triplet by chance

Neuron 2

θ2(3)Timedoublet by chanceθ2−4

Neuron 7

The last parameter of the log-linear model provides the higher-order correlation.

i2BNeuron 3(2)

θ12A state-space framework allows regression of a time-dependent system to spike data.

A-4(2)

4Amari, IEEE Trans. Inf. Theory, 2001

Tripplewise-correlation model

13-720i

1

2-70[First column] 50 trials of simulated parallel spike trains. Each pattern of spikes was counted with 1[ms] resolution. Red dots indicate synchronous spikes(doublets). [Second column] Log-linear parameters. Shaded area indicates 99% confidence interval.

Introduction

1-72Pairwise Correlation

The mixed coordinate (η1 , η2 , θ12 ) is an orthogonal coordinate.

Pairwise-correlation modelExperiment

θ1( )

-2

θ2( )

Neuron 2i

Hypothetical Dynamic Cooperative Activity of Neurons

Neuron 60Rate model (Independent model)

To study cooperative neural network activity, we develop state-space method toestimate the time-dependent correlation structures embedded in parallel spike trains.

True Pairwise correlation
Estimated from n=200 trials
Neuron 22

η12 = E [x1 x2 ] = p11log p11 = −ψ + θ1+ θ 2+ θ12Pairwise Correlations between Neuron 7 and others

0-2

η2 = E [x2 ] = p01 + p11

log p 01 = −ψ + θ2p01

Neuron 54

η1 = E [x1 ] = p10 + p11

log p 10= −ψ + θ1p10

Neuron 1

Log-linear parameters0η -coordinateθ -coordinate

Neuron 1

Log-linear parameters

50 Trials

log p (x1, x 2 )= − ψ + θ1x1+ θ2x 2+ θ12x1x 2State-space method applied to a sub-model is useful for large-scale data analysis.

Time-dependent rates + Constant correlation

50 Trials

ate Bernoulli process using a log-linear
link function [4-8]. A nonlinear recursive filtering formula is derived from a logquadratic approximation to the posterior
distribution of the state. Together with afixed-interval smoothing algorithm, timedependent log-linear parameters are estimated. The smoothed estimates are optimized via EM-algorithm such that theirprior covariance matrix maximizes the expected complete data log-likelihood. In
addition, we perform model selection on
the hierarchical log-linear models to avoid
over-fitting. Application of the method
to simultaneously recorded neuronal spikesequences is expected to contribute to
uncover dynamic cooperative activities ofneurons in relation to behavior.

8 Neurons

50 Trials

Classical studies in neurophysiology arebased on the idea that stimulus information is encoded in the firing rates of single
neurons. Alternatively, precise spike coordination is discussed as an indication ofcoordinated network activity and may beexpressed by higher-order correlation between the simultaneous spiking activitiesof neurons. As shown in earlier studies [1],correlation between neurons may modulating in time in relation to the behavioral
demand.
Here, we simultaneously estimate thetime-dependent rate and correlation underlying multiple-neuron spiking activities
by means of state-space analysis [2,3].
We model discretized parallel spike trains
by a conditionally independent multivari-

2 Neurons

50 Trials

Log-linear model

50 Trials

Abstract

where log p(y1:T |w) =PT

t=1 log p(yt|y1:t−1, w) ∼ nPT′θ
(yt=1 t t|t−1 − ψt|t−1).

Note that the marginal likelihood takes account of the parameter uncertainty only
(its approximation is known as BIC).
Rate model Pair-corr. model Triple-corr. model
r=1
r=2
r=3
n=5
2085∗
2097
2144
7913
7698∗
7728
n = 20
19263
18542
18540∗
n = 5036811
36781∗
n = 100 38478
73330
73238∗
n = 200 76683
The complexity of the model needs to be selected based on the sample size and theprominence of higher-order structure.

Correspondingparameters

θ1

θ2

θ3

θ4

θ12 θ13 θ23 θ14 θ24 θ24 θ123 θ124 θ134 θ234

Summary and Outlook
We developed a method for estimating time-varying rate and higher-order correlation
structure in parallel spike sequences. To our knowledge, it is the first method that canresolve the time-dependent spike rates and ’well-defined’ spike correlation measuressimultaneously. The method is thus applicable to simultaneously recorded neuronal
spike sequences recorded from an awake behaving animal. Such an application isexpected to provide us with new insights into dynamic assembly activities, their compositions, and behavioral relevance.
The presented method will be summarized in
Shimazaki H., Amari S., Brown E. N., and Grün S. “State-space Analysis on Timevarying Correlations in Parallel Spike Sequences”, Proc. IEEE ICASSP2009 in press

Jensenl +1| w )≤ log p (y1:T | wl +1 )lComplete data log-likelihood

References[1] Riehle et al. (1997) Spike synchronization and rate modulation differentially involved in motor cortical function, Science, 278, 1950-1953.
[2] Smith and Brown (2003) Estimating a state-space model from point process observations. Neural Comput., 15, 965-991.
[3] Eden et al. (2004)Dynamic analysis of neural encoding by point process adaptive filtering., Neural Comput, 16, 971–998.
[4] Martignon et al. (2000) Neural coding: higher-order temporal patterns in the neurostatistics of cell assemblies, Neural Comput, 12 2621

[5] Amari (2001) Information geometry on hierarchy of probability distributions, IEEE Trans. Inf. Theory 47(5) 1701-1711[6] Nakahara and Amari (2002) Information-geometric measure for neural spikes. Neural Comput. 18: 1259-1267.
[7] Schneidman et al. (2006) Weak pairwise correlations imply strongly correlated network states in a neural population, Nature, 440, 1007.
[8] Shlens et al. (2006) The structure of multi-neuron firing patterns in primate retina, J Neurosci, 26, 8254–8266.