Abstract Journal of Neuroscience Methods 161 (2007) 155–165 Detecting correlation changes

Journal of Neuroscience Methods

161 (2007) 155–165

Detecting correlation changes in electrophysiological data

Jianhua Wu a,Keith Kendrick b,Jianfeng Feng a,c,?

a Department of Computer Science,University of Warwick,Coventry CV47AL,UK

b Laboratory of Cognitive and Behavioural Neuroscience,The Babraham Institute,Cambridge CB24AT,UK

c Department of Mathematics,Hunan Normal University,Changsha410081,P.R.China

Received14March2006;received in revised form16October2006;accepted22October2006

Abstract

A correlation multi-variate analysis of variance(MANOV A)test to statistically analyze changing patterns of multi-electrode array(MEA) electrophysiology data is developed.The approach enables us not only to detect signi?cant mean changes,but also signi?cant correlation changes in response to external stimuli.Furthermore,a method to single out hot-spot variables in the MEA data both for the mean and correlation is provided.Our methods have been validated using both simulated spike data and recordings from sheep inferotemporal cortex.

? 2006 Elsevier B.V. All rights reserved.

Keywords:Multielectrode data;Hierarchical multi-variate analysis of variance;Mean;Variance

1.Introduction

The neural activities induced by a stimulus,such as an odor,an image or a sound,are routinely recorded by multi-electrode arrays(MEA,see Victor et al.,1994;Gross et al., 1995;Rutten et al.,2001;Feng,2004).Usually hundreds or even thousands of local?eld potentials and neuronal activi-ties are simultaneously recorded at a sample frequency between 1000and2000Hz(Mehring et al.,2003;Todd,2005).Such an experiment could last a few hours and produce data sets of hundreds of GBytes.Several analytical software packages(e.g., Neuroexplorer,MATLAB)are available for editing and visual-izing data,and involving in time–frequency analysis(Masimore et al.,2004).However,there is very little that has been proposed for detailed statistical analysis and subsequent visual represen-tations of changing spatial and temporal patterns of information recorded across an electrode array.Although many methods have been proposed to analyze single cell activity,such techniques are not directly applicable to multi-neurons activities(Shaw et al., 1999;Horton et al.,2005).

The?rst step of analyzing MEA data,in particular,the local ?eld potential data and spike trains,is to single out which vari-ables(neurons or electrodes)are sensitive or responsive to the stimulus(odor,image or sound).For such a task,multivariate

?Corresponding author.Tel.:+442476573788.

E-mail address:jianfeng.feng@https://www.wendangku.net/doc/6e15039807.html,(J.Feng).analysis of variance(MANOV A)(Johnson and Wichern,2002) is probably the?rst choice to test the difference between two populations with spatial and temporal patterns.When we apply MANOV A to MEA data,we have demonstrated its pros and cons in Horton et al.(2005).A modi?ed version of MANOV A,which we call adaptive MANOV A(Wu et al.,2006),has been proposed to deal with MEA data sets.We have shown that the adap-tive MANOV A cannot only perform the statistical testing task, but also accurately detect the hot-spot variables,i.e.,variables responding to an external stimulus.

However,adaptive MANOV A,as in the traditional MANOV A,is designed to test the changes in mean activities. As we all know,the mean activity of a neuronal network could remain unchanged when a stimulus is presented,but it could somehow reorganize itself to show more or less coherent activ-ities.A typical example is that two neurons?re independently, but synchronize themselves whenever a stimulus is applied.The mean?ring rates before-and during-stimulus are identical.The question we ask here is whether adaptive MANOV A could detect the coherence(correlation)changes in the recorded physiolog-ical data.Here coherence changes imply that the covariance matrix between variables before-and during-stimulus is differ-ent.Pattern changes include both mean activity(mean vector) and coherent activity(covariance matrix)changes.

A prerequisite for the applicability of adaptive MANOV A is that the covariance matrices among variables are equal.The fact is that for local?eld potentials and spike trains data,the correla-tions between variables are likely to change due to the stimulus,

0165-0270/$–see front matter? 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jneumeth.2006.10.017

156J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165

as mentioned above.Hence the adaptive MANOV A is not appli-cable when the covariance matrices vary during the stimulus.

In this paper,we propose methods to test pattern changes in physiological data,especially MEA data including both local ?eld potentials and spike trains.We extend adaptive MANOV A to detect the correlation matrix changes in response to an external stimulus.In doing this,the Fisher’s z-transformation and Edge-worth expansion(Bhattachary and Ghosh,1978;Blinnikov and Moessner,1998)are introduced.The transformation helps us to transform the distribution of covariance,which is not Gaussian in general,into Gaussian distribution,which is required in a MANOV A analysis.The expansion provides us with a criterion to assess how close is the joint distribution of transformed covari-ance and a multivariate normal distribution.After overcoming these dif?culties,the adaptive MANOV A can be directly applied to test the changes of the covariance matrix.The detected hot-spot correlations reveal the changed pair-variables in response to the stimulus.

The method is?rst applied to two types of simulated spike data,which demonstrates the necessity of correlation MANOV A test.We then apply the method to the local?eld potentials recorded while animals performed an operant discrimination task between different pairs of sheep faces in order to obtain a food reward.The paper?rst reviews the mean MANOV A test, and then develops correlation MANOV A test based on it.An algorithm for correlation MANOV A test is proposed in Section 2.A full version of MANOV A is introduced in Section2as well. Our approach is then tested by using simulated spike data and recordings from sheep’s vision experiments.

2.Methods

Traditional multi-variate analysis of variance(MANOV A)is used to investigate whether two mean vectors are signi?cantly different or not(Johnson and Wichern,2002).In dealing with physiological data MANOV A could be extremely useful when the dimension of the population vectors is small.But as we have mentioned above,for MEA data we usually face a huge amount of data with a large dimensionality(i.e.,a dimension>1000). MANOV A in its textbook version is proved not to be the ideal choice to handle high dimensional population dataset,especially for complex biological data such as local?eld potentials(LFPs) recorded from MEA,Gene Microarrays,multi Spike trains and fMRI data(functional magnetic resonance imaging)(Horton et al.,2005).Observing from electrophysiological recordings,it is possible that information seems to be carried by a small number of signi?cant channels or neurons rather than the whole set of channels or neurons.The dif?culty we are facing is how to pick up the hot-spots(variables with signi?cantly changed activities) from the whole group of variables.An adaptive MANOV A and a corresponding HOTTOR algorithm to resolve the issue has been proposed in Wu et al.(2006),and the method is proved to be very useful in detecting the hot-spots.But a restriction in the paper is that the correlation matrix of the variables of pre-and during-stimulus must be identical.

Fig.1demonstrates the situation with mean changes(left) and variance changes(right).When the covariance is identical for pre-and during-stimulus,the mean MANOV A test could be successfully applied to test the difference.When the covariance is different,mean MANOV A test is no long applicable no mat-ter whether the mean values are the same or not.Fig.1(right) shows that the covariance of variables is signi?cantly changed but the mean values are identical.As we all know,synchroniza-tion between neurons could play a signi?cant role in information processing and a proper statistical approach to test and single out changed pair-neurons is of vital importance.

In this section,we will?rst introduce the adaptive MANOV A, which we call mean MANOV A test.Then,we will discuss the method to test the changes of the correlation matrix of recorded variables,which is termed correlation MANOV A test.An

algo-Fig.1.Raster plot of100neuron activity generated from the integrate-and-?re models.(Left)The mean?ring rate of neurons in the left and right hand of the vertical thick line is different,i.e.,at around400ms,a stimulus is turned on(a constant current is added in the model).(Right)The mean?ring rate of neurons in the left and right hand of the vertical thick line is identical,but their correlation changes,i.e.,at around400ms,a stimulus is turned on(inputs are correlated).

J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165157 rithm(HOTTOR)to pick up changed variables is proposed.

Finally,a full version of correlation MANOV A is introduced.

2.1.Models

Let X1=(x11,...,x1m)(say,pre-stimulus LFPs),X2=(x21,

...,x2m)(say,during-stimulus LFPs)be two random matrices.

The spatial and temporal pattern of the recordings can be repre-

sented by the rows and columns of the matrices.For example,

for LFP data,x1i is the one pre-stimulus LFP variable recorded

during time period[0,Nh],where h is the sampling resolution

and N is the total sample size.Correspondingly,we have

1={φ1ij?(x1i?Ex1i)(x1j?Ex1j)}i,j=1,...,m

and

2={φ2ij?(x2i?Ex2i)(x2j?Ex2j)}i,j=1,...,m

the random covariance matrix of X1and X2.

The mean MANOV A test is designed to?nd all variables

between X1and X2,which are signi?cantly different.In fact,a

pre-requirement to perform the mean MANOV A test is that

E 1=E 2(2.1)

In other words,before carrying out a mean MANOV A test,

we have to test whether Eq.(2.1)is true or not,which is the main

concern of the current paper.

In Fig.1,raster plot of100neuron activities is presented.

Spikes are generated from the integrate-and-?re model with

independent or correlated Poisson inputs.For spikes data and

for a given time window[kT w,(k+1)T w]?[0,T],where T w is

the recording time window,de?ne

x ik=N[kT

w,(k+1)T w]

T w

,i=1,2

where N[kT

w,(k+1)T w]is the number of spikes recorded in the time

window,k ranges from0to T/T w?1,and T is the total length of recordings.It is known that for i=1,2the joint distribution of X i is normal(Johnson and Wichern,2002).For example,in both?gures in Fig.1,the spikes in the left hand side of the thick vertical line could be spontaneous activity(pre-stimulus), which yields X1.The right hand side of the thick vertical line gives us X2.The data in Fig.1is generated according to the fact that EX1=EX2,but E 1=E 2(left)and EX1=EX2,but E 1=E 2(right).Mean MANOV A can handle the situation in Fig.1(left)but not in Fig.1(right).

2.2.Mean MANOVA test

For the completeness of our current paper,we?rst present a brief review of the mean MANOV A.For population1,let us assume that the?rst n variables have mean1,and the rest have mean0.For population2,the?rst n variables have means c1, c2,...,c n,respectively,and the rest have means0.In other words,n random variables take different values due to,say,dif-ferent treatments or stimuli.The hot-spot is{1,2,...,n}.In the literature,the most common way to compare the difference between two populations is MANOV A.In MANOV A,we?rst construct a statistical quantity called Wilks lambda.The validity

of the hypothesis that μ1:=EX1(the mean activity of popula-

tion1)= μ2:=EX2(the mean activity of population2)or not

is assessed based on Wilks lambda.

In order to evaluate the Wilks lambda,we need the covari-

ance matrices of both populations.In this section,the covariance

matrix =E 1=E 2is assumed to be equal and generated randomly(each upper triangle element is i.i.d.N(0,σ2),since

the covariance matrix must be symmetric and positive de?nite).

MANOV A is de?ned based on the Wilks lambda ,as in

Horton et al.(2005).The signi?cance score of A={1,...,m}is

given by

S(A)=?

N?

m+2

2

log( (A))/χ2m(α)(2.2)

whereχ2m(α)is the upper(100α)th percentile of a Chi-square distribution with m degrees.The element of set A is the index of the variable,i.e.,the i th element of A stands for the i th variable for both populations. (A)is computed according to Appendix A.For any subset B∈{1,...,m},the signi?cance score could be de?ned accordingly.Obviously,a subset B is signi?cantly changed if its score S(B)is larger than unity.

By de?ning the signi?cance score,a few interesting prop-erties are found for the signi?cance score.First of all,the signi?cance score is a decreasing function of m when m>n.Sec-ondly the signi?cance score is an increasing function of n when m is?xed.Thirdly,the signi?cance score is a decreasing func-tion with respect toρ(ρis the correlation coef?cient between variables)with all other parameters?xed.Fourthly,the deriva-tive of signi?cance score with respect to a subset B={1, (i)

is not continuous at i=n.It is the crucial property to develop the following HOTTOR algorithm to detect the hot-spot.

Fig.2presents a simple simulated example to demonstrate the properties we discussed above.In this example,we

assume

Fig.2.Upper panel:3-D signi?cance score is plotted against the subset A and c∈[0,1].There is an in?exion point for all three?gures in the upper panel,and this point could be identi?ed clearly at the lower panel at point30.

158J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165 all cis equal to c for simplicity and the covariance matrix is

=?

??

??

??

??

??

??

?

??

?

1ρ

..

.

ρ1

?

??

?

n×n

?

??

?

ρ1ρ1

..

.

ρ1ρ1

?

??

?

n×(m?n)?

??

?

ρ1ρ1

..

.

ρ1ρ1

?

??

?

(m?n)×n

?

??

?

1ρ

..

.

ρ1

?

??

?

(m?n)×(m?n)

?

??

??

??

??

??

??

(2.3)

where we assume that|ρ1| |ρ|,re?ecting the structure in means.There are total100variables for each population,and only?rst30have signi?cant change from population1to pop-ulation2.The sample size is200.

Three dimensional signi?cance score is plotted against the subset A and values of c∈[0,1]for three different cases,ρ<0,ρ=0,ρ>0.Obviously,for?xed subset A and c,the interesting point here is that the signi?cance score forρ<0is always larger than that forρ>0.This agrees with one of our long-term working hypotheses:it is much easier to detect the changes in a system with negatively correlated units than with positively correlated or independent units.

The two dimensional signi?cance score is also presented with c?xed.The signi?cance score reaches maximum when only part

of the signi?cantly changing variables are included.In?exion occurs when exactly all the signi?cantly changing variables are involved in the test.This means in?exion rather than the maxi-mum signi?cance score can be regarded as the criterion to detect the whole signi?cance variables.

After the signi?cance score is de?ned and the proper algo-rithm is employed,the detection of the hot-spots is quite routine. But as we mentioned at the beginning of the section,it is a pre-requirement to satisfy Eq.(2.1)in mean MANOV A.In the next subsection we turn our attention to the assumption Eq.

(2.1),which is of physiologically interests as well,as mentioned before.

2.3.Correlation MANOVA test

For two populations,which both follow multi-normal dis-tributions,there only exist the following relations listed below since the distributions are fully determined by their means and variances.

(1)Equal means and variances?same population.

(2)Different means,equal variances?mean MANOV A test

applicable.(3)Different means,equal variances?mean MANOV A test

not applicable.

(4)Equal means,different variances?mean MANOV A test

not applicable.

If the variances are different,mean MANOV A test is not applicable.To determine whether the two populations are sig-ni?cantly different or not,we have to compare the covariance matrices 1and 2.From Kenny and Keeping(1962)and Eckhardt(1984),we have the exact form of the distribution of the correlation(Eqs.(2.5)and(2.6)).Actually,we will go even a step further to de?ne

φkij(N)=

N

ξ=1

(x ki(ξ)?Ex ki(ξ))(x kj(ξ)?Ex kj(ξ))

N

ξ=1

(x ki(ξ)?Ex ki(ξ))

N

ξ=1

(x kj(ξ)?Ex kj(ξ))

(2.4) where x ki(ξ)is theξth element of the vector x ki,k=1,2,i

If Eφkij=0,the probability density function ofφkij(N)is

p(r)=r

ν

1?r2

～t?distribution,ν=N?2.(2.5)

If Eφkij=0,let I=((ν?2)/2)and J=((ν?3)/2),the proba-bility density function ofρis

p(r)=?

???

???

???

???

1?

2

√

π

Γ((ν+1)/2)

Γ(ν/2)

I

k=0

(?1)k

I!

(I?k)!k!

|r|2k+1

2k+1

,forνeven;

1?

2

√

π

sin?1|r|+|r|

J

k=0

(2k)!!

(2k+1)!!

(1?r2)k+1/2

,forνodd.

(2.6)

If we de?ne

Z kij(N)=tanh?1φkij(N),(2.7)

the variable Z kij(N)is approximately normally distributed with a

meanμZ(k,i,j,N)=tanh?1Eφkij and varianceσ2Z(k,i,j,N)=

1/(N?3).Actually,the transformation de?ned by Eq.(2.7)is

usually called Fisher’s z transformation.

Now it is clear that every element of the correlation matrices

follows a normal distribution with the same variance if N is kept

as a constant.The issue now is whether the joint distribution of

ˉZ k(N)={Z kij(N)}

i

k=1,2(2.8)

still follows a normal distribution.However,it is possible that

the joint distribution does not follow a multivariate normal dis-

tribution even if every element is normally distributed.

In order to determine whether the joint distribution of

the elements follows a normal distribution or not,the Edge-

worth expansion of the density p(ˉρ),ˉρ=[ρ1,...,ρd]with

d=m(m?1)/2,is introduced.Edgeworth expansion approxi-

mates a probability distribution in terms of its cumulants and

Hermite polynomials.Let p n be the probability density func-

tion,and it is chosen to be a normal distribution with mean and

J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165159

variance equal to the mean and variance of {ˉZ

k (N )},k =1,2.Then we can expand the probability density function p (ˉρ)to

p (ˉρ)≈p n (ˉρ)??1+13! ξ,η,ζ

k ξ,η,ζh ξηζ(ˉρ)+14!

ξ,η,ζ,τ

k ξ,η,ζ,τ

h ξηζτ(ˉρ)+O (ˉρ)?

?(2.9)

with h ξηζbeing the ξηζth Hermite polynomial,ξ,η,ζ∈{1,...,

d },and k ξ,η,ζth

e corresponding standardized cumulant de?ned by k ξ,η,ζ=k ξηζ/δξδηδζ,with k ξηζthe cumulants ξ,η,ζ,τ∈{1,...,d }and δζis the standard deviation o

f the ζth random vari-able.Similarly,h ξ,η,ζ,τis the ξηζτth Hermite polynomial and the correspondin

g standardized cumulant k ξ,η,ζ,τis de?ned by k ξ,η,ζ,τ=k ξηζτ/δξδηδζδτwit

h k ????the cumulant.

In Eq.(2.9),the probability density function is approximated by a normal distribution plus a residual term,which could be fully determined by the cumulants of the given variables.p n (ˉρ)is the main term for the approximation,but the magnitude of the second term is important in determining the approxima-tion accuracy.To estimate the difference between p (ˉρ)and its corresponding normal distribution p n (ˉρ),let us de?ne Q =1?13!

ξ,η,ζ

k ξ,η,ζmin k =1,2h ξηζ({μZ (k,i,j,N )}i

(2.10)as the measurement of the approximation accuracy.In the sequel,if Q is greater than a signi?cant similarity level (by running the simulation,85%can be regarded as the signi?cant similar-ity level),then we can conclude that the joint distribution p (ˉρ)approximately follows a multivariate normal distribution.

To determine the magnitude of 13! ξ,η,ζk ξ,η,ζh ξηζ(ˉρ),we

need some more basic knowledge of cumulant (Kenny and Keeping,1962;Abramowitz and Stegun,1972),which is explained in Appendix B .

The joint distribution of the transformed correlations could be calculated through Eq.(2.9)when the cumulants are estimated from the empirical data.If the joint distribution approximates a multivariate normal distribution well,then Mean MANOV A test is applicable to the correlations.

The procedure for correlation MANOV A model can be as (1)Calculate the correlation matrix 1and 2for population

1and 2,respectively.

(2)Apply inverse hyperbolic tangent transformation to each

φkij (N ).

(3)Calculate the joint distribution of the transformed cor-relations,compare the distribution with corresponding multivariate normal distribution (Eqs.(2.9)and (2.10)).(4)Apply mean MANOV A test to the transformed correlation

variables.The key point for correlation MANOV A test is whether the joint distribution of the transformed correlations follows a mul-tivariate normal distribution.As we discussed,every single correlation transformed by Fisher’s z transformation approxi-mately follows a Gaussian distribution,but the joint distribution may not follow multivariate normal distribution.The Edgeworth

expansion is applied to test whether the joint distribution can be approximated by a Gaussian distribution.As mentioned,if the approximation is good enough (more than 85%similarity),the joint distribution could be regarded as a multivariate normal dis-tribution and then the mean MANOV A test is directly applicable to the transformed correlation variables.The critical similarity value 58%is determined by applying the method to simulated data:the ?nal results do not change if the similarity of the dis-tribution with Gaussian distribution is greater or equal to 85%.The application procedure is fully described in the following section.2.4.Algorithm

The whole algorithm is divided into three parts.The ?rst part is unique for correlation MANOV A test,the other two parts are the same as the mean MANOV A test.

2.4.1.Part I:transformation part

First we introduce the notation involved in the algorithm.Suppose we have two populations X 1=(x 11,x 12,···,x 1m )and X 2=(x 21,x 22,···,x 2m )with same number of variables m and sample size N .The elements ρkij (h )of correlation matrix k are calculated by Eq.(2.11),the dimension of 1is m ×m ×(N ?M +1).M is the number of samples employed to calculate the correlation.h increases from 1to N ?M +1.

ρkij (h )= h +M ?1

q =h

x ki (q )x kj (q ) h +M ?1q =h x ki (q ) h +M ?1

q =h x kj (q )(2.11)Let z kij (h )=tanh ?1ρkij (h )(Fisher’s z transform).After the

transformation,the Edgeworth expansion is employed to mea-sure the similarity of the joint distribution with Gaussian distribution.If the distribution is not multivariate normal dis-tribution,correlation MANOV A test is not applicable (see Section 4).If the distribution is approximately multivariate normal distribution,we de?ne Y k be the collection of all the upper diagonal elements of {z kij }.The dimension of Y k is (m (m ?1)/2)×(N ?N +1),and can be regarded as two new populations for further analysis.

2.4.2.Part II:monotonic part

First we compute the single unit pairs,which produce maxi-mum or minimum signi?cance score between population Y 1and Y 2,and let these two pairs be (i min ,j min )and (i max ,j max ).Let A be the whole set containing {(i ,j )},where i ranges from 1to m ?1and j from 2to m ,A q be the set containing signi?cantly changing variables.The initial value for A q is (i max ,j max ).

Then selecting a unit pair,say (i 1,j 1),from set A /{(i min ,j min ),(i max ,j max )},we calculate S (A q ∪{(i 1,j 1)})(refer to Eq.(2.2)).If S ({(i 1,j 1)})

160J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165

signi?cance score decreases.The monotonic step stops when the signi?cance score reaches its maximum.This can be achieved by comparing all the signi?cance scores.

2.4.

3.Part III:reference part

After?nding the unit pairs,which positively contribute to the signi?cance score,now we have to deal with those units,which do not positively contribute to the signi?cance score,but they change their mean values(refer to the units between solid line and dash line in Fig.2,lower panel).At this step,we assume that at least one pair,which does not change its activity is picked up,say(i min,j min).We use the score of this unit as a reference to assess whether a remaining pair,say(i,j),changes its activity or not.Denote the obtained set as A q+?A q.If it changes,i.e., S(A q+∪{(i min,j min)})

If it does not change,S(A q+∪{(i min,j min)})≥S(A q+∪{(i,j)}), we then simply choose a new pair to check.Continuous the procedure above until we check all units.We then have all the signi?cantly changing pairs A q+.

The procedure of the algorithm could be expressed as a?ow chart,seen in Fig.3.The?ow chart is divided into three parts, which represent the three steps.Each part should be?nished ?rst,then next part starts.The red arrow is the only relation between part I and part II,and part II and part III.The meaning of the symbols are explained below:

?A:the whole set{(i,j)},i ranges from1to m?1and j from 2to m;

?(i min,j min):the unit pair with minimum signi?cance score;?(i max,j max):the unit pair with maximum signi?cance score;?A q:the set containing indexes of the signi?cantly changing variables;

?A p:the set containing indexes of the signi?cantly and insignif-icantly changing variables.

2.5.Recording from inferotemporal cortex

Electrophysiological data were acquired from the inferotem-poral cortex of awake behaving sheep(Nicol et al.,2005; Tate et al.,2005).All experimental procedures involving ani-mals were conducted in strict accordance with guidance on the operation(The Animals Act,1986).Sheep,under halothane (?uothane)anaesthesia,were implanted with two chronic64-channel MEAs in the right and left inferotemporal cortices, respectively.Recordings were made from the unanaesthetised sheep while they performed an operant discrimination task in which different pairs of sheep faces were presented and a correct panel-press response elicited a food reward.Local event-related LFPs are recorded simultaneously from128chronically implanted electrodes.LFP time series data are sampled at2000 samples/s from around3s prior to stimulus onset to3.5s after the stimulus onset in each trial of a session,and stored as16-bit numbers.

3.Results

3.1.Simulated spike data

Before applying correlation MANOV A test to real data,we ?rst test the model by using spike data generated using

MATLAB Fig.3.Algorithm?ow chart.The transformation part will produce two new populations Y1and Y2.Part II and part III are processed based on these two populations. The initial values for A q and A p are{(i max,j max)}and{(i min,j min)},respectively.The initial values for A q+and A p+are the?nal values of A q and A p,respectively. (i1,j1)and(i,j)are randomly selected from set A/{(i min,j min),(i max,j max)}and A p,respectively.

J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165161

7.01(The MathWorks,Natick,MA,USA).Two types of spike data are generated.The?rst type is that the covariance matrix is?xed for pre-and during-stimulus,and the mean?ring rate for a small number of neurons varies between pre-and during-stimulus.Another type is that the mean?ring rate is?xed,but the covariance matrix of a small number of neurons varies between pre-and during-stimulus.

The raster plot of the data sets is presented in Fig.4.The data is generated103s before stimulus,and103s during-stimulus.In Fig.4,only2s data for pre-and during-stimulus are included. In total there are300neurons,but only10neurons respond to the stimulus.We randomize the300variables such that we cannot tell,which variable is responsive to the stimulus.The randomized data set is presented in Fig.4(bottom panel).For the randomized data,the signi?cantly changed neuronal activi-ties are obscured by the large number of unresponsive neuronal activities.

In our analysis,sampling interval t=100ms and total record-ing time T=1000s are used.There are N=T/t=104samples for each neuron.The correlation MANOV A test is?rst applied to the data set on the left hand side of Fig.4(left panels).A total d=300(300?1)/2=44850correlation variables are cal-culated for both pre-and during-stimulus.When facing such a huge number of variables,we need to divide the variables into groups arti?cially,and apply mean MANOV A to each divided group.After Fisher’s z transformation,10sets Y1i and Y2i

are Fig.4.Raster plot of generated data from the integrate-and-?re model.In the upper left panel,a case with different means is illustrated.The covariance between population1and population2is identical for pre-and during-stimulus(stimulus onset is indicated by the vertical line).The mean?ring rate is changed for10neurons between pre-and during-stimulus.In the upper right panel,a case with different covariance is illustrated.The mean?ring rates for all variables are the same for pre-and during-stimulus.The covariance matrix between10neurons is changed from pre-stimulus to during-stimulus.Bottom panels:identical to the upper panels but neuron labels(from1to300)are randomized.The neurons with signi?cantly changing activities are mixed with other unresponsive neurons.

162J.Wu et al./Journal of Neuroscience Methods

161 (2007) 155–165

Fig.5.Results of mean MANOV A and correlation MANOV A.Upper panels:left hand panel depicts average ?ring rate for pre-stimulus;right hand panel presents average ?ring rate for during-stimulus.The red regions indicate neurons with signi?cantly changing activities.Bottom panels:left hand panel is the average correlation for pre-stimulus;right hand panel shows average correlation for during-stimulus.The red regions indicate neurons with signi?cantly changing activities.(For interpretation of the references to color in this ?gure legend,the reader is referred to the web version of the article.)

constructed,the dimension for Y 1i and Y 2i is 4485×9000(i =1,...,10),which means that each transformed correlation variable has a sample size of 9000.By applying the mean MANOV A test to Y 1i and Y 2i ,after 4485selection steps,no signi?cantly chang-ing variables are detected,which con?rms that the correlation matrix does not have any signi?cant change when the stimulus is applied.The mean MANOV A is then directly applied to 300mean activities.The average ?ring rates for pre-and during-stimulus are presented in Fig.5(upper panel),and the regions marked red show the signi?cantly changed neurons.All the sig-ni?cantly changing neurons (10neurons)are detected by mean MANOV A test.

For the data set of Fig.5(right upper panel),the average ?ring rate is a constant.It is not surprising that there is no detected changed activity when mean MANOV A test is directly applied (the maximal signi?cance score is 0.72).The correla-tion MANOV A test is applied to the correlation matrix.There are a total of d =300(300?1)/2=44850correlation variables.These variables are again divided into 10groups.After perform-ing Fisher’s z transformation,we have 10sets of Y 3i and Y 4i ,all with a dimension of 4485×9000.Then mean MANOV A test is employed to detect the signi?cantly changed variables in each groups.Only 2500correlations among the total of 44850vari-ables are depicted in Fig.5(bottom panel,left and right).The red regions in right bottom panel are the variables that have signi?-cantly changed activities.All the changed correlations are fully detected by our correlation MANOV A.3.2.Application to recordings

Translating the experiment results into our setup,we have collected data recorded from the left inferotemporal cortex:X 1=(x 1,1,x 1,2,...,x 1,64)pre-stimulus,and X 2=(x 2,1,x 2,2,...,x 2,64)during-stimulus.For x i ,j ,i is referred to pre-or during-stimulus,j is the channel number.The similar data have been recorded from the right inferotemporal cortex.As discussed at the beginning of the paper,the issue we face now is to pick up changed variables from the large experimental data set.

Before applying the algorithm,pre-processing the data is nec-essary.The data is standardized:each variable is divided by its square root of variance.As before, 1and 2are the sam-ple correlation matrix with dimension 64×64×6000.Fisher’s z transformation is applied to the sample correlation matrix.Then by employing Edgeworth expansion,the similarity of transformed joint distribution of the correlation matrix to a multi-normal distribution is calculated.The similarity factor for pre-and during-stimulus are 88.46and 93.28%respectively,which imply that the joint distribution can be regarded as a multivariate normal distribution and mean MANOV A test is applicable to the covariance matrix.By applying the algorithm to the transformed correlation matrices,the signi?cantly changed correlations are detected and presented in Fig.6.The circle is the template of the multi-electrode array,and each number represents one electrode.The red numbers represent the variables,which are responsive to the stimulus.Fig.6presents the results of the correlations for both left and right cortex.The number of signi?cantly changed variables detected by correlation MANOV A test is 10and 14out of 64channels for the left and right cortex,respectively.

The channels with a signi?cantly changed correlations are connected by lines in Fig.6.The green lines indicate that the correlation increases during the stimulus,but the blue lines imply the reduction of the correlation during the stimulus.The values of correlations for pre-and during-stimulus between these variables are further listed in Tables 1and 2.Among

J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165

163

Fig.6.Results of correlation MANOV A test.The green connection lines between variables show that the correlation between channels increase during the stimulus,but the blue lines indicate the decrease of the correlation during the stimulus.The values of correlations among signi?cantly changed channels are listed in Tables1and2. (For interpretation of the references to color in this?gure legend,the reader is referred to the web version of the article.)

Table1

Correlation list for left visual cortex

Variable pair Correlation(pre-stimulus)Correlation(during-stimulus) 20,270.340.58

27,340.240.49

22,290.130.42

20,36?0.040.51

38,450.31?0.10

52,550.110.45

Table2

Correlation list for right visual cortex

Variable pair Correlation(pre-stimulus)Correlation(during-stimulus) 12,190.210.52

12,21?0.120.31

21,290.630.18

23,390.020.47

27,360.180.59

27,34?0.07?0.42

38,390.240.56

46,530.41?0.00

53,590.080.38

changed channels,correlations increase more than decrease.The increased correlations tell us that more regions are coherently more active or less active during stimulus.

4.Discussion

We have presented a novel approach:correlation MANOV A test based upon our previous work(Horton et al.,2005;Wu et al.,2006),to analyze the electrophysiological data recorded by multi-electrode arrays.By applying the method to both simu-lated and biological data(spikes and LFPs),we conclude that the approach is successful in detecting variables,which are sta-tistically signi?cantly responsive to an external stimulus.As mentioned before,we face increasingly huge amounts of data in electro-physiological experiments and our method could play a crucial role in tackling such data sets.4.1.Correlation detection versus synchronization

The proposed methods can successfully detect the pattern changes in response to the stimulus,and they also pick up the variables that are statistically signi?cantly responsive to the stimulus.We have demonstrated that the correlation MANOV A test can detect synchronization or de-synchronization between variables.The test detects pair-variables(neurons)that have signi?cantly changed during stimulus.If positive correlation between two variables increases,we could conclude that the variables during stimulus are more synchronized than pre-stimulus.

Of course correlation is a more general notation than syn-chronization.For example,if we assume that two spike trains are(perfectly)synchronized at spike to spike level,the corre-lation between the two spike trains is unity.If,in statistical sense,two out of10spikes are synchronized,then the corre-lation between the two spike trains is around2/10=0.2.We could say that the two spike trains are partially synchronized.If the second spike train is delayed by a fraction of the sampling interval,then the two spike trains are no longer(perfectly)syn-chronized.However,the correlation between two spike trains is unchanged.Hence the detection of correlation changes could be more meaningful than the detection of synchronization,which is more appealing to a visual inspection.

During the detection of correlation and synchronization, timescale is a very important factor.For example,hippocampal place cells with similar place?elds will?re in a highly corre-lated fashion at macroscopic(>200ms)timescales.However,?ne-timescale(<120ms)dynamics are also tightly regulated such that there is little synchronous?ring(O’Keefe and Recce, 1993;Skaggs et al.,1996).These observations suggest that correlation and synchronization occur at multiple temporal lev-els.So the interaction of observation timescales needs to be addressed during the detection of correlation and synchroniza-tion.

Furthermore,spike trains could exhibit negative correlations, which is beyond desynchronization.Roughly speaking,nega-tive correlation between two variables implies that if one?res

164J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165

faster,the other should be slower and vice visa.For two peri-odic signals,a negative correlation between them means there is a phase shift between them.When two periodic signals are negatively correlated with a correlation coef?cient of?1,we know that they are in opposite phases.In terms of correlation MANOV A we are able to single out both negative(anti-phase)and positive(on-phase)changes in the recorded signals (neurons).

4.2.Instantaneous?ring rates versus interspike intervals

When dealing with spike trains,the results presented in the current paper are based upon the information of instantaneous ?ring rates.One might argue that it is not the optimal way to tackle the spike train data since interspike intervals might pro-vide more information than instantaneous?ring rates.Hence a natural question is whether we could apply our approach to interspike intervals.The answer is con?rmative,but it requires some more developments.Essentially,as we have emphasized in the paper,it is a prerequisite for applying MANOV A that the variables are normally distributed.With interspike intervals,we know that they are not normally distributed and,even worse,the exact distribution is unknown.Fortunately,in the past few years, we have successfully developed a non-parametric approach to accurately?t interspike interval distributions(Rossoni and Feng, 2005).Therefore,a combination of the results in the current paper(mean MANOV A and correlation MANOV A)and the results in Rossoni and Feng(2005)could lead to a successful application of the hierarchical MANOV A to interspike interval data.

4.3.Correlation MANOVA and causality

Due to the correlations(interactions)between variables, multi-variable analysis is the only appropriate way to deal with multi-variable recordings such as multi-electrode array data. Nevertheless,with multi-variables,usually we face the interac-tions between two or more groups of variables.An appropriate tool is required to statistically assess the changes between groups of variables.Correlation MANOV A as we brie?y summarized in the current paper serves this role.

In our approach,we consider the interaction between two groups of variables.The interaction(correction)is no-directional.An observation that the activity of neuron B is correlated with that of neuron C could mean that the input of neu-ron B comes from the output of neuron C,or that both of them receive a common input.To distinguish between these cases, we have to take into account the causality analysis(Chen et al., 2004;Albo et al.,2004).Again,causality analysis is usually per-formed on single to single variable bases and a multi-variable causality analysis is still lacking.

Acknowledgment

J.F.was partially supported by grants from UK EPSRC (EP/C51338X),(EP/D051916),and(GR/S30443).Appendix A

De?neˉx l=

N

j=1

x lj/N for l=1,2andˉx=

2

l=1

N

j=1

x lj/ (2N),where N is the number of samplings.Furthermore,let us de?ne sampling covariances

SSW=

2

l=1

N

j=1

(x lj?ˉx l)(x lj?ˉx l)

SSB=

2

l=1

N(ˉx l?ˉx)(ˉx l?ˉx)

Then Wilks’lambda(Johnson and Wichern,2002)is intro-duced as:

=

SSW

SSB+SSW

=

2(N?1)

SSB+2(N?1)

=

11 12

21 22

11+B 12

21 22

where =

11 12

21 22

and B=SSB/2(N?1).

Appendix B

In probability theory and statistics,the cumulant k n of the probability distribution of a random variable X are given by (Kenny and Keeping,1962)

E(e tX)=exp

∞

n=1

k n t n

n!

The correlation between cumulants and moments gives us an easy way to calculate cumulants from known moments.Cumu-lants are related to the moments by the following recursion formula:

k n=μ n?

n?1

k=1

n?1

k?1

k kμ n?k

The joint cumulant of several random variables X1,...,X n is k(X1,...,X n)=

π

B∈π

(|B|?1)!(?1)|B|?1E

i∈B

X i

whereπruns through the list of all partitions of{1,2,...,n}, and B runs through the list of all blocks of the partitionπ,|B|is the size of the set B.

References

Abramowitz M,Stegun IA,editors.Handbook of mathematical functions with formulas,graphs,and mathematical tables.9th ed.Dover;1972.

J.Wu et al./Journal of Neuroscience Methods 161 (2007) 155–165165

Albo Z,Di Prisco GV,Chen Y,Rangarajan G,Truccolo W,Feng JF,et al.Is partial coherence a viable technique for identifying generators of neural oscillations?Bio Cyber2004;90:318–26.

Bhattachary RN,Ghosh JK.On the validity of the formal Edgeworth expansion.

Ann Stat1978;6:434–51.

Blinnikov S,Moessner R.Expansions for nearly Gaussian distributions.Astron Astrophys Suppl Ser1998;130:193–205.

Chen Y,Rangarajan G,Feng JF,Ding M.Analyzing multiple nonlinear time series with extended Granger causality.Phys Lett A2004;324:26–35. Eckhardt DH.Correlations between global features of terrestrial?elds.Math Geol1984;16:155–71.

Feng JF,https://www.wendangku.net/doc/6e15039807.html,putational neuroscience:a comprehensive approach.Chap-man and Hall/CRC Press;2004.

Gross GW,Rhoades BK,Azzazy HM,Wu M-C.The use of neuronal networks on multielectrode arrays as biosensors.Biosens Bioelectron1995;10:553–67. Guidance on the operation of The Animals(scienti?c procedures) Act,1986.Online source:http://www.archive.of?https://www.wendangku.net/doc/6e15039807.html,/document/hoc/321/321.htm.

Horton P,Bonny L,Nicol A,Kendrick K,Feng JF.Applications of multi-variate analysis of variance(MANOV A)to multi-electrode array electrophysiology data.J Neurosci Methods2005;146:22–41.

Johnson RA,Wichern DW.Applied multivariate statistical analysis.5th ed.

Prentice Hall;2002.

Kenny JF,Keeping ES.Mathematics of statistics.3rd ed.New Jersey:Van Nostrand;1962.p.71–82.

Masimore B,Kakalios J,Redish AD.Measuring fundamental frequencies in local?eld potentials.J Neurosci Methods2004;138:97–105.

Mehring C,Rickert J,Vaadia E,Cardoso de Oliveira S,Aertsen A,Rotter S.

Inference of hand movements from local?eld potentials in monkey motor cortex.Nat Neurosci2003;6:1253–4.Nicol AU,Magnusson MS,Segonds-Pichon A,Tate A,Feng JF,Kendrick KM.

Local and global encoding of odor stimuli by olfactory bulb neural networks.

Program No.476.62005Abstract Viever/Itinerary Planner.Washington,DC: Society for Neuroscience;2005.

O’Keefe J,Recce ML.Phase relationship between hippocampal place units and the EEG theta rhythm.Hippocampus1993;3:317–30.

Rossoni E,Feng JF.Non-parametric approach to extract information from inter-spike intervals.J Neurosci Methods2005;150:30–40.

Rutten W,Mouveroux J,Buitenweg J,heida C,Ruardij T,Marani E,et al.

Neuroelectronic interfacing with cultured multielectrode arrays toward a cultured probe.Proc IEEE2001;89:1013–29.

Shaw F-Z,Chen R-F,Tsao H-W,Yen C-T.A multichannel system for recording and analysis of cortical?eld potentials in freely moving rats.J Neurosci Methods1999;88:33–43.

Skaggs W,McNaughton B,Wilson M,Barnes C.Theta phase precession in hip-pocampal neuronal populations and the compression of temporal sequences.

Hippocampus1996;6:149–72.

Tate A,Nicol A,Fischer H,Segonds-Pichon A,Feng JF,Magnusson M,et al.

Lateralised local and global encoding of face stimuli by neural networks in the temporal cortex.Program no.362.1.2005Abstract Viewer/Itinerary Planner.Washington,DC:Society for Neuroscience;2005.

Todd C,editor.Event-related potentials:a methods handbook.MIT Press;2005. Victor JD,Purpura K,Katz E,Mao B.Population encoding of spatial frequency, orientation,and color in macaque V1.J Neurophysiol1994;72:2151–

66.

Wu JH,Kendrick K,Feng JF.Detecting hot-spot in spatio-temporal patterns of complex biological data;2006.Online source:http://www.dcs.warwick.ac.

uk/～jian-hua/doc/hottor.html.