Mathematical details, supplemental to the paper: "Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves."Wellcome Trust (SJE; grant number 083205) and EPSRC (CSC

Cutts, C. S.

Eglen, S. J.

Apollo (Cambridge)

Detecting pairwise correlations in spike trains: anobjective comparison of methods and application tothe study of retinal wavesMathematical detailsCatherine S Cutts, Stephen J EglenAugust 13, 20141 Correlation index is confounded by firing rateIn this section we derive Equation 2 from the paper which is used to show that thecorrelation index is confounded by firing rate. To make this document self-contained wereproduce the formula for the correlation index iA,B between spike trains A and B. Here,the vectors a and b represent the spike times of neurons A and B; ai is the ith spike intrain A and bj is the jth spike in train B.iA,B =NA,BTNANB2∆t(1)whereNA = |a| (total number of spikes of A in recording),NB = |b|,T = length of recording,∆t = synchronicity window,and NA,B is the number of spike pairs where a spike from train A falls within ±∆t of aspike from train B:NA,B =NA∑i=1NB∑j=11∆t(|ai − bj |)1where1∆t(x) =1 if x ≤ ∆t0 otherwiseNo one statistical model is able to recreate the full range of observed spiking patterns,so a general proof that the correlation index is dependent on firing rate is not possible.Instead we assume a simple model for firing times and use this as a counter example tothe statement “the correlation index is in general independent of firing rate”.The model is as follows: the spikes of neurons A and B form a Poisson process withrates λA, λB respectively. A set proportion of these spikes occur in both trains (i.e. eachtrain fires a spike synchronously with the other), these occur with rate λS ≤ λA, λB .We derive an expression for the correlation index under this model and show that it isdependent on firing rate.The recording time T and synchronicity window ∆t are independent parameters (i.e.they cannot be evaluated in terms of firing rates). The expected values of NA and NBare as follows:NA =λATNB =λBTThe calculation of the remaining term, NA,B , is non-trivial. We use the following resultsfor Poisson processes (Kingman, 1992):• Superposition: If X and Y are two independent Poisson processes with respectiverates λX and λY then X+Y is a Poisson process with rate λX + λY .• Decomposition: To decompose this composite process into two processes withidentical statistics to X and Y (i.e. two Poisson processes with rates λX and λY ),it suffices to generate a series of Poisson events with rate λX +λY and assign eachevent as type X with probability λxλx+λy and as type Y with probabilityλyλx+λy.It follows from the above that the synchronous spikes in our model form a Poisson pro-cess with rate λS . The expected number of synchronous spikes is therefore λST . This2contributes to NA,B since each of these is one spike pair where a spike from A falls within±∆t of a spike from B. There is also a second contribution to NA,B from cases wherespikes from different processes (i.e. any combination except A and A, and B and B) fallwithin ±∆t of each other. We denote the expected number of such spike pairs C andthen:NA,B = λST + C (2)To evaluate C we consider the system as a decomposition of a master Poisson processwith rate λ = λA + λB − λS . Events are assigned as A only (A) with probability λA−λSλ ,B only (B) and synchronous (S) with probabilities given similarly. By considering themaster process we can work out the probability of spikes from each combination of trainsfalling within ±∆t of each other by calculating the probability that any spikes fall within±∆t of each other and then multiplying this by the probability of them being from thecorrect combinations.We note that multiple spikes can occur within ±∆t of each other, not just pairs andso our calculations need to be recursive, that is we have to consider the case when nsuccessive inter-spike-intervals (ISIs) sum to less than ∆t and sum the contributions fromall n. We also have to be careful not to double-count events: if we consider all the spikepairs formed when n inter-spike-intervals sum to less than ∆t then most of those spikepairs will have been covered in the previous cases n = 1 . . . n − 1. The only pair notalready counted is the pair formed by the first and last spikes. Therefore, for arbitrary n,the expected number of new spike pairs (Cn) is:Cn = In(∆t)PnSn (3)whereIn(∆t) =probability that n successive ISIs sum to ≤ ∆t,Pn =probability that outer two spikes are not both A or both B (can be synchronous),Sn =number of n successive ISIs in the recording.3andC =∞∑n=1CnWe find each quantity in turn:In(∆t): Let the ith ISI in the master train be Zi. Since the master train is a Poissonprocess with rate λ:Zi ∼ Exponential(λ) (4)andk∑i=1Zi ∼ Erlang(k, λ) (5)(Kingman, 1992) where the cumulative distribution function of the Erlang distribution is:ErlangC(z; k, λ) = 1−k−1∑n=01n!e−λz(λz)n (6)The probability of n ISIs summing to less than ∆t is then given by:In(∆t) = ErlangC(∆t;n, λ) (7)Pn: The combinations of spikes which will contribute a spike pair to NA,B are as follows:AB(BA),AS(SA),BS(SB). Combinations in brackets occur with the same probability astheir counterpart outside the brackets. The combination SS will contribute two spikepairs to NA,B so to account for this we multiply its probability by two (meaning that Pnis not strictly a probability). To simplify the notation we introduce the following:λA′ =λA − λSλB′ =λB − λSand Pn is then given by:Pn =2λ2S + 2λSλA′ + 2λSλB′ + 2λA′λB′λ2(8)4This is independent of n so we drop the subscript, Pn = P . This simplifies toP =2λAλB(λA + λB − λC)2 (9)Sn: The number of n successive ISIs in the recording is λT −n (λT is the total numberof spikes in the recording).We substitute these three terms into Equation 3 to find Cn and sum over all n to findC:C =λT∑n=1(λT − n)ErlangC(∆t;n, λ)P (10)for simplicity, we leave this in terms of P . Then,NA,B = λST +λT∑n=1(λT − n)ErlangC(∆t;n, λ)P (11)Having found NA,B , the correlation index is:ia,b =(λST +λT∑n=1(λT − n)ErlangC (∆t;n, λ)P)T(λAT ) (λBT ) 2∆t(12)This expression is clearly dependent on rate which is what we set out to show. Wenow simplify it to gain an understanding about the strength of the dependency and howit varies with rate.We first simplify C:C =λT∑n=1(λT − n)ErlangC(∆t;nλ) =λT∑n=1(λT − n)(1−n−1∑i=01i!e−λ∆t(λ∆t)i) (13)After rearrangement involving rewriting the double sums into single sums:C =λT (λT − 1)2(1−e−λ∆tλT∑n=01n!(λ∆t)n)+(λ2T∆t−λ∆t+λ2∆t22)e−λ∆tλT∑n=01n!(λ∆t)n(14)5writinge−λ∆tλT∑n=01n!(λ∆t)n = e−λ∆t(eλ∆t −∞∑n=λT+11n!(λ∆t)n) (15)= 1− e−λ∆t∞∑n=λT+11n!(λ∆t)n (16)we assume that1− e−λ∆t∞∑n=λT+11n!(λ∆t)n ≈ 1 (17)which is true if the following hold:• λ∆t ≤ 1. This holds within experimentally observed ranges since λ is O(1) Hz andtypically ∆t = 0.05− 0.1 s.• λT is sufficiently large. This also holds within experimentally observed ranges,for instance if λ = 1 Hz, ∆t = 0.1 s and T = 4 s then Equation 16 is equal to0.9999999 (7 d.p). Since T is typically O(1000) s, in practice this approximationwill be even more accurate.If the above holds, thenC ≈ λ2T∆t− λ∆t− λ2∆t22(18)Using this approximation:iAB ≈ (λST + (λ2T∆t− λ∆t− λ2∆t22)P )1λAλBT2∆t(19)substituting in our expression for P and simplifying gives us the form of Equation 2 in thepaper:iA,B =λS2∆tλAλB+ 1− 1(λA + λB − λS)T −∆t2T(20)In the paper we use the sub-case of auto-correlation λ = λA = λB = λS to illustratethe dependence on firing rate. In this case:iA,B =1λ(12∆t− 1T)+(1− ∆t2T)(21)6Both the term in 1λ and the second term are positive since ∆t ≤ 1 (within exper-imentally observed ranges) and T = O(1000) s. Therefore increasing λ decreases thecorrelation index.72 Normalisation of the Spike Time Tiling CoefficientIn this section we justify the form of the normalisation of the Spike Time Tiling Coefficient.To make this self contained, we reproduce the formula for the Spike Time Tiling CoefficientC from Figure 1 of the paper:C =12(PA − TB1− PATB)+12(PB − TA1− PBTA)wherePA =The proportion of spikes in train A which fall within ±∆t of any spike in train B,PB =the proportion of spikes in train B which fall within ±∆t of any spike in train A,TA =the proportion of the recording time which falls within ±∆t of any spike in train A,TB =the proportion of the recording time which falls within ±∆t of any spike in train B.The rationale behind this is that we would expect the proportion of spikes in A whichfall within ±∆t of any spike in B by chance to be equal to the proportion of the recordingtime which falls within ±∆t of any spike from B. The quantity PA−TB will therefore bepositive/negative if this happens more/less than expected by chance which is indicativeof positive/negative correlation. It does not, however, possess the following necessaryqualities:• N1: The measure should be symmetric i.e. C(A,B) = C(B,A).• N4: The measure should be bounded taking a value of +1 when the spike trains areidentical, with a value of zero corresponding to no correlation and -1 correspondingto anti-correlation.We normalise this quantity to ensure it has property N4 and then take the mean of thetwo directed variants to ensure it satisfies property N1 (hence the factor 12 ).For simplicity, we consider one of the directed variants and define:8f(PA, TB) =PA − TB1− PATBIn order to fulfil N4, the following is needed (where 0 ≤ PA ≤ 1, 0 ≤ TB ≤ 1):1. f(1, TB) = 1.This ensures that auto-correlations always take the value +1.2. f(0, 1) = −1. Ensures that the most anti-correlated patterns take value -1.3. f(PA, TB) = 0, if PA = TB (PA 6= 1, TB 6= 1). Ensures that uncorrelated spikestake value zero.4. f(PA, TB) > 0 if PA > TB and f(PA, TB) < 0 if PA < TB . Ensures that themeasure is positive/negative for positive/negative correlation.The following two properties must also hold since they are implicitly implied by thenotion of a correlation measure:5. f(PA, TB) must be monotonically increasing in PA for PA 6= 1, TB 6= 1.6. f(PA, TB) must be monotonically decreasing in TB for PA 6= 1, TB 6= 1.i.e. if the correlation increases, so does f , if it decreases, so does f .Properties One,Two and Three are trivially satisfied. We note that for PA 6= 1, TB 6= 1,1− PATB > 0 so property Four is also satisfied. Note also:∂f∂PA=1− T 2B(1− PATB)2which is positive within the ranges of PA and TB , so f increases monotonically with PA(property Five) and∂f∂TB=P 2A − 1(1− PATB)2which is negative within the ranges of PA and TB , so f decreases monotonically with TB(property Six).Therefore f satisfies all required properties and the normalisation is appropriate.9Alternative normalisationsWe note that this is not the only possible normalisation of the Spike Time Tiling Coeffi-cient. Other normalisations exist which have all the properties listed above, for instance2− PA − TB . Using another normalisation factor which has these properties will notchange whether the Spike Time Tiling Coefficient possess all the necessary and desiredproperties. However, the sensitivity of the measure to changes in values of PA and TBdoes change for different normalisations. This is important as it will affect how big achange in firing patterns is needed to affect a noticeable change in the Spike Time TilingCoefficient which should be as sensitive as possible to changes within the experimentallyobserved ranges of PA and TB . The version in the paper was chosen as we found it tobe sensitive to small changes in these values.To illustrate this, we briefly compare the sensitivity of f with that of the followingdifferently normalised Spike Time Tiling Coefficient, g:g(PA, TB) =PA − TB2− PA − TBIt is important that the measure is sensitive to changes within the experimentallyobserved ranges of values. The values of PA typically take all values in the range 0− 1,but the distribution is highly skewed towards small values and values of TB are usuallysmall (Figure 1). We therefore want our coefficient to be as sensitive as possible tochanges within this range and since within it f(PA, TB) is more sensitive than g(PA, TB)(Figure 2), f is preferred.ReferencesKingman C (1992) Poisson Processes, Vol. 3 Clarendon Press.Kirkby LA, Feller MB (2013) Intrinsically photosensitive ganglion cells contribute toplasticity in retinal wave circuits. Proc Natl Acad Sci USA 110:12090–12095.10PA (equivalently PB) Probability0.0 0.2 0.4 0.6 0.8 1.00.00.10.20.30.40.5TA (equivalently TB) Probability0.00 0.02 0.04 0.06 0.080.000.050.100.150.200.25A BFigure 1: Experimentally observed distributions of PA and PB are skewed towardslow values and values of TA and TB are low. A: Histogram of the values of PA (orequivalently PB) from a recording of spontaneous retinal activity in P5 wild type mousefrom Kirkby and Feller (2013). Note the local maximum at PA or PB = 1, which are thevalues of the auto-correlations. B: Histogram of the values of TA (or equivalently TB)from the same recording. This recording has duration 998.1 s, 44 electrodes and meanfiring rate 0.28 Hz. In both panels ∆t = 0.1 s as in the original publication.11PAT B0.00.20.40.60.81.00.0 0.2 0.4 0.6 0.8 1.0PA − TB1 − PATB0.0 0.2 0.4 0.6 0.8 1.0PA − TB2 − PA − TB−1.0−0.50.00.51.0Figure 2: The form of normalisation affects the sensitivity of the Spike Time TilingCoefficient to changes in PA and TB. A: Heatmap of f (one half of the Spike TimeTiling Coefficient with denominator 1−PATB) for 0 ≤ PA ≤ 1, 0 ≤ TB ≤ 1. Changes incolour show changes in contours of f . The sensitivity of f (and therefore the sensitivityof the Spike Time Tiling Coefficient) to changes in the values of PA and TB can be seenby the widths of the contours- narrower contours show higher sensitivity. Since typicallyTB is small and the distribution of PA is skewed towards small values, the sensitivityof the Spike Time Tiling Coefficient in this region (marked with a black dotted-line) ismost crucial. B: Same as A but with f replaced by g (one half of the Spike Time TilingCoefficient with denominator 2− PA − TB).12

Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves. Mathematical Details

http://www.repository.cam.ac.uk/handle/1810/245681

Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves. Mathematical Details

Abstract

Similar works

Full text

Available Versions

Apollo (Cambridge)