A simple indicator of nonstationarity of firing rate in spike trains

Journal of Neuroscience Methods 163 (2007) 181–187

A simple indicator of nonstationarity of firing rate in spike trains Boris Gour´evitch ∗ , Jos J. Eggermont Department of Physiology and Biophysics, Department of Psychology, University of Calgary, Calgary, Alta., Canada Received 4 January 2007; received in revised form 21 February 2007; accepted 21 February 2007

Abstract The assessment of stationarity of firing rate in neural spike trains is important but is often performed only visually. Facing the growing amount of neural data generated by multi-electrode recording, there is a need for an automatic method to identify and disqualify spike trains with highly nonstationary firing rates. In this report, we propose a simple test of nonstationarity, associated with an indicator quantifying the degree of nonstationary in a spike train. This method is compared to the Mann–Kendall test of trend detection and the Runs test on simulated and real spike trains. © 2007 Elsevier B.V. All rights reserved. Keywords: Nonstationarity; Spike trains; Anderson–Darling; Interspike intervals

1. Introduction Spike trains provide for information transfer between neurons. Consequently, spike trains emitted by individual neurons are used as the basis for estimating mutual information or other measures of dependence. Several types of nonstationarity or natural variability have been identified in spike trains, which can lead to misleading results: firing rate (FR) throughout the duration of experiment (Grun et al., 2001; Zhu et al., 2003), firing rate across trials (Grun et al., 2003), response latency variability across trials (Nawrot et al., 2003), spike waveforms during experiments (Penev et al., 2001). This paper deals with the nonstationarity of FR throughout the duration of experiment. These nonstationarities generally have two causes: either preparation-related ones, such as the emergence or loss of neural units due to small brain movements and changes in anesthesia level (Sokal et al., 2000) or stimulus-related ones causing fluctuations in FR (Averbeck and Lee, 2003). This paper focuses on spontaneous activity for long stationary stimuli that should result in long-term stationary spike trains. Most signal analysis methods require that the data are at least weakly stationary for the FR (Brown et al., 2004). The concept of weak stationarity in neural spike trains is defined as a stability of average firing rate and variance of FR over time. Breaking these rules ∗ Corresponding author at: Department of Psychology, University of Calgary, 2500 University Drive N.W., Calgary, Alta., Canada. Tel.: +1 403 220 7747 (office) ; fax: +1 403 282 8249. E-mail address: [email protected] (B. Gour´evitch).

0165-0270/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jneumeth.2007.02.021

leads to biased results in numerous techniques widely used in spike train analysis, such as firing rate related techniques, peristimulus time histograms (PSTH/JPSTH (Gerstein, 1988) or cross-correlation histograms (Perkel et al., 1967). Nonstationarities may sometimes even be the only cause of correlation between spike trains (Brody, 1999). Consequently, it is crucial to a priori check that neural data do not show too strong nonstationarities, even if in practice few papers mention it in their methods. When mentioned, nonstationarities in spike trains have been identified by: (1) visual inspection of interspike intervals (Benoit and Chataignier, 1973); (2) estimation of the mean, variance or coefficient of variation (CV) of the FR per defined time bin or trial (Poggio and Viernstein, 1964; Levine and Troy, 1986; Habets et al., 1987; Ramakers et al., 1990); (3) using the Side test of Griffith and Horn (Griffith and Horn, 1966; Nakahama et al., 1971, 1972); (4) using the Wald–Wolfowitz Runs test (Bradley, 1968) on interspike intervals (ISI) (Correia and Landolt, 1977; Starr et al., 2006) or FR on small bins (Horwitz and Newsome, 2001); (5) using change-point detectors, such as the Cusum algorithm (Ellaway, 1978; Belisle et al., 1998). However, visual inspection becomes a very long and largely irrelevant task with the increasing number of recording electrodes resulting in a growing amount of neural data. Moreover, methods based on variance or coefficient of variation implicitly assume Gaussianity of the data and often arbitrarily define analysis intervals or bins. The Runs test and the Side test, which are in agreement (Rodden and Landolt, 1974), act as autocorrelation detectors, which implies that they should

182

B. Gour´evitch, J.J. Eggermont / Journal of Neuroscience Methods 163 (2007) 181–187

be able to detect nonstationarities. They are not very powerful though as we will show. Finally, the change-point detectors do not take into account the various types of fluctuations in FR that may occur in a recorded neural spike train. Zhu and colleagues proposed other ways to compensate for the effects of nonstationarity, such as “dividing the data into a set of short but relatively stationary segments or removing the time-varying average from the data” (Zhu et al., 2003). These methods, of course, only are applicable when studies are restricted to FR analysis. As a consequence, there is a need for a simple and automatic preprocessing method for the detection of significant nonstationarities in spike trains, i.e. of trends in the fluctuations of interspike intervals. In addition to the Runs test, another reference test for trend detection is the Mann–Kendall (MK) test (Mann, 1945; Kendall, 1975) whose power is similar to the Spearman’s Rho test (Yue et al., 2002). However, it is restricted to monotonic trends. In this paper, we propose a simple test for nonstationarity in the FR of spike trains based on the comparison of the cumulative spike arrival times with an ideal uniform cumulative distribution function (CDF). A quantification of nonstationarity is extracted from this test. We compare our method’s features to the MK test and the Runs test on simulations and real data. A Matlab code is provided as supplementary material.

In our case where FU (x) = x and because FZ (x) is a step function, some manipulations lead to the following expression

2. Methods

 AN =

2.1. Stationarity test Let tn be the occurrence time of the nth spike in a spike train of N spikes recorded between 0 and T seconds. No hypothesis is made about the spike train. Let Zn =tn /T be the normalized occurrence time of the spikes. Under the assumption of weak stationarity, the mean FR and the variance of FR should remain constant over time. Consequently, (Zn )n=1,. . .,N should be a set of uniform values on [0 1]. Intuitively, it means that the empirical cumulative distribution function (CDF) of this set, FZ (Zn )=n/N, should oscillate closely around the CDF of the uniform law, FU (x) = x, x ∈ [0 1], as illustrated in Fig. 1. In the absence of any hypothesis about the spike train, we use the results from statistical theory to derive a maximal acceptable variation of FZ (x) around FU (x): in the general case of any empirical CDF FZ (x) built from a sample (Zn )n=1,. . .,N following any given distribution U, the asymptotic distribution (N → ∞) of the following statistic  AN =

0

1

(FZ (x) − FU (x))2 dx x(1 − x)

Fig. 1. Empirical cumulative distribution function (CDF) FZ (x) of the spike arrival times Zn normalized between 0 and 1 for a spike train of N spikes oscillates around the CDF FU (x) of the uniform distribution between 0 and 1. The area between the two curves is shaded.

0

1

N

(FZ (x) − x)2 1 (2n − 1) dx = −N − x(1 − x) N n=1

ln (Zn (1 − ZN−n+1 )) .

(2)

Though the distribution of AN has a very complicated closed form, and so only a few tabulated values have been used for 50 years, Marsaglia and Marsaglia recently proposed a very accurate approximation (5 digits) of the asymptotic CDF (Marsaglia and Marsaglia, 2004), which leads to P(AN < z) = z−1/2 e−1.2337141/z (2.00012 + (0.247105 −(0.0649821 − (0.0347962 −(0.0116720 − 0.00168691z)z)z)z)z)

(3)

for 0 < z < 2, with |error| < 0.000002, and P(AN < z) = exp(−exp(1.0776 − (2.30695 − (0.43424 −(0.082433 − (0.008056−0.0003146z)z)z)z)z)) (4)

(1)

has been established in Anderson and Darling (1952) as a goodness-of-fit test of a sample to a probability distribution. It is usually known as the Anderson–Darling (AD) test. Eq. (1) shows that AN compares the area between the two CDFs (Fig. 1) and thus quantifies the variation of FZ (x) around FU (x).

for 2 < z < ∞, with |error| < 0.0000008. The p-value associated with the proposed stationarity (S) test (null hypothesis: the spike train is stationary) is p = 1 − P(AN < z). Given the very small values of p that will be manipulated in the following, we found it more convenient to use a nonstationarity indicator NS = −ln(p) as a “surprise” value for AN (Palm, 1981). A Matlab 6 code for this test is provided in the supplementary material.

B. Gour´evitch, J.J. Eggermont / Journal of Neuroscience Methods 163 (2007) 181–187

183

Fig. 2. The two types of spike trains (SP1, SP2) used in simulations. Each dot represents a spike. Experiment time is represented along the abscissa; dots are plotted randomly on the vertical scale. For SP1 (a, b), a Gamma process of 90 s at a rate of 10 spikes per second is followed by another Gamma process of 90 s with a firing rate of 3 spikes per second. For SP2 (c, d), two Gamma processes of 90 s at 10 spikes per second surround another Gamma process of 90 s with a firing rate of 3 spikes per second. Different shape parameters for the Gamma distribution are used, from (a, c) A = 0.2 with a very bursty behavior to (b, d) A = 4 with more regular interspike intervals.

2.2. Mann–Kendall test Let ISIn = tn+1 − tn be the nth ISI. The statistics associated to the MK test is TMK =

N  i−1 

sign(ISIi − ISIj )

(5)

i=2 j=1

where sign(x) is 1, 0, −1 if x > 0, x = 0, x < 0, respectively. Under the null hypothesis of a monotonic trend in ISI values, we have asymptotically    (N − 1)(N − 2)(2N + 3) TMK → N 0, . (6) 18 2.3. Runs test We divide the recording of T seconds in M intervals of q seconds. Let Fm be the number of spikes in the mth interval. Let Xm be 1 if Fm is greater than the median of (Fm )m=1,. . .,M , 0 otherwise. A “run” is a sequence of adjacent equal symbols in (Xm )m=1,. . .,M . For example, the sequence 0011100011110 is composed of 5 runs, 3 of them are made of “0” and the others are made of “1”. Let M0 be the number of “0” in (Xm )m=1,. . .,M , and M1 the number of “1”. If R is the number of runs in (Xm )m=1,. . .,M , then under the assumption of the randomness of (Xm )m=1,. . .,M , we have asymptotically    M + 2M0 M1 2M0 M1 (2M0 M1 − M) R→N , (7) M M 2 (M − 1) 2.4. Simulations In order to evaluate the properties of our proposed test, we simulate two spike trains composed of renewal processes following a Gamma law γ(A, 1/rA) of density fγ (x) =

(rA)A xA−1 e−rxA Γ (A)

value in the interval [5 15] during 90 s. For the second spike train (SP2), a 90 s train with r ∈ [5 15] is inserted between two trains of 90 s each with r = 10 sp/s (Fig. 2). For each train, A is chosen between 0.2 and 4. A is the inverse square of the coefficient of variation. The Poisson process corresponds to A = 1. The greater the value of A, the more regular the spike train is. A Gamma process is generally a better representation of a neural spike train than a Poisson process, since it allows for the presence of a relatively refractory period (Kuffler et al., 1957; Stein, 1965; Troy and Robson, 1992; Mukherjee and Kaplan, 1998; Baker and Lemon, 2000). For each value of A and r, 1000 datasets of spike trains are generated and the p-values for S, MK and Runs tests are computed. The results are averaged over these 1000 datasets. For each value of A and r, we therefore estimate the power of the tests by the proportion of datasets for which the p-value < 0.05. The power of a test is the probability of rejecting the null hypothesis of stationarity when the alternative hypothesis of nonstationarity is true. This quantifies the sensitivity of the tests to the variations of the parameters A and r. The specificity of the tests – do not reject the null hypothesis of stationarity when it is true – can be quantified by 1-power for the case r = 10 sp/s, i.e. when no nonstationarity is introduced. For the Runs test, we chose q = 3 s as the value which gave the best power to the test. 2.5. Real data The S, MK and Runs tests are also applied to spontaneous activity spike trains recorded from the primary auditory cortex of the anesthetized cat. The length of the recording is 900 s for each dataset. This dataset was recorded for a previous study (Tomita and Eggermont, 2005) where the protocol is detailed. The spike trains represent multiunit recordings. The lowest FR for our examples being 0.7 sp/s, we used the value q = 10 s for the Runs test to get at least 7 spikes per bin in average for every recording. 3. Results

(8)

where r is the average FR, A a shape parameter for the ISI distribution and Γ is the Gamma function. For the first spike train (SP1), r = 10 sp/s during 90 s then r is suddenly modified to a

3.1. Simulations For simulated spike train 1 (SP1), the MK test performs slightly better than the S test for regular data (A > 1) but its

184

B. Gour´evitch, J.J. Eggermont / Journal of Neuroscience Methods 163 (2007) 181–187

Fig. 3. Simulated spike train SP1: detection of a permanent change in FR in the second half of the spike train. For the S test, (a) power of the test at 5% risk, i.e. probability to reject the null hypothesis of stationarity and (d) average p-value of the test as a function of the r and A values defining alternative hypotheses; contour lines for several NS values are shown. Results for the MK test as shown in (b and e) and for the Runs test in (c and f) are in the same format. Specificity of the tests can be quantified by 1-power for the case r = 10 sp/s.

power (probability to reject the null hypothesis of stationarity) is dramatically decreasing for A < 1 (Fig. 3b and e). In contrast, whatever the value of A, the power of the S test (Fig. 3a) is always close to 1 when r is above 12 sp/s or below 8 sp/s, and the average NS value remains almost constant (Fig. 3d). This suggests that the ability of the test to detect large nonstationarities remains high even for a very bursty spike train. This nice behavior has a downside: the risk of the first kind (probability to wrongly reject the null hypothesis of stationarity, i.e. the results in Fig. 3a–c for r = 10) slightly increases when A < 1 for the S test whereas it remains constant for the MK test. In other words, the specificity of the S test, which is better than that of the MK test when A > 1, becomes worse than for the MK test when A < 1. This is explained by the fact that the lower A is, the more bursty

the spike train is and the more questionable the definition of stationarity is. The Runs test has a similar behavior as the MK test but always performs worse (Fig. 3c and f). For simulated spike train 2 (SP2), the FR fluctuation is nonmonotonic since the discontinuity in FR occurs in the middle of the train. Consequently, the MK test has a constant p-value whatever the values of A and r and is not able to detect the FR fluctuation (Fig. 4b and e). In contrast, both S and Runs tests show similar results for this train and the SP1 train (Fig. 4a, c, d and f). The S test is slightly less selective than for the SP1 train but remains more powerful than the Runs test for A < 2, i.e. the most common case in real spike trains. The same decrease of specificity of the S test than for the SP1 train is observed for A < 1.

Fig. 4. Simulated spike train SP2: detection of a temporary change in FR in the middle of the spike train. For the S test, (a) power of the test at 5% risk, i.e. probability to reject the null hypothesis of stationarity and (d) average p-value of the test as a function of the r and A values defining alternative hypotheses; contour lines for several NS values are shown. Same format for the MK test in (b and e) and for the Runs test in (c and f). Specificity of the tests can be quantified by 1-power for the case r = 10 sp/s.

B. Gour´evitch, J.J. Eggermont / Journal of Neuroscience Methods 163 (2007) 181–187

185

Fig. 5. Examples of four nonstationary spike trains (a–d) and a stationary one (e) recorded in auditory cortex. The p-value for S, MK and Runs tests and NS statistics are shown on the left. The function FZ (x) − FU (x), on which the Anderson–Darling test is based (see Eq. (1)), is represented on the right, above the dot-raster representation of the spike train (time in abscissa, random value in ordinate).

3.2. Real data Fig. 5 presents four examples of nonstationary cortical spike trains recorded in a cat’s primary auditory cortex, along with one spike train considered as stationary on a visual criterion, using the provided dot-raster representation of the spike train. Nonstationarities in spike trains of Fig. 5a and c likely ensue from the too close position of the electrode to a unit, irritating the neuron and provoking a sudden bursting and a subsequent temporary spike loss. For all these spike trains, the ISI distribution is well fitted with a Gamma law, and the shape parameter A is under 0.8. The p-values for the S test are quite low, except for the apparent stationary case (Fig. 5e). At risk 5%, the MK and Runs tests fail to detect the nonstationarities for cases c and

d. However, the results of the Runs test heavily depend on the choice for q (Fig. 6). 4. Discussion To our knowledge, the most powerful methods to detect nonstationarities in spike trains derive either from the Runs and Mann–Kendall tests or from the maximum likelihood method leading to the Cusum algorithm (Ellaway, 1978). We showed that, compared to the S test, the Runs test performs poorly for our real neuronal spike train data and generally for simulations with A < 2. Moreover, the choice of the interval length q can dramatically affect the results of the Runs test. In contrast to the S test, the Mann–Kendall test does not detect any

186

B. Gour´evitch, J.J. Eggermont / Journal of Neuroscience Methods 163 (2007) 181–187

Fig. 6. Case (c) of Fig. 5: p-value for the Runs test as a function of the value for q.

non-monotonic trend in the FR fluctuations. Although efficient, the Cusum method acts as a change-point detector (Belisle et al., 1998), and it is not obvious how to apply it as a general FR-fluctuation detector. However, in a certain way, our method could be seen as a global nonparametric Cusum method: the function FZ (x) − FU (x) is indeed proportional to the current FR minus the average FR, on which the Cusum is based. The AD statistic acts as a transform of this function to obtain a level of significance. Our method uses the only implicit assumption that the spike times arrivals are somewhat regular, the ideal case being a constant ISI. An apparent weakness of the method is that very irregular ISI may lead to false alarms even if the FR is globally constant, as it locally deviates from the implicit assumption. For example, an ideal neural spike train, such as a realization of a Poisson process, has a strong intrinsic stationarity since the ISI distribution remains the same over time. However, the choice of a small analysis window may give rise to strong variations of local average FR and variance. Consequently, the p-value of our test slightly decreases for stationary irregular spike trains (risk of the first kind) as shown in Figs. 3a and 4a, i.e. its specificity decreases. However, the most important behaviour to expect from our test is that global nonstationarities lead to lower pvalues than local nonstationarities. This is one quality of the proposed test that maintains some very low p-values for highly nonstationary process regardless the shape of the ISI distribution (Figs. 3 and 4), showing its robustness to local nonstationarities. In addition, the influence of local variations of interspike intervals decreases with the length of the spike train. In contrast, for both MK and Runs tests, nonstationarities become difficult to detect when the variance of the ISI increases. Because of this behavior, the following limitations apply: (1) we recommend a strict significance level for the S test, for example NS ≥ 5, which is equivalent to a p-value of 0.007. Notice that due to the accuracy of the approximation of P(AN < z) to 10−6 , NS statistic values can not be very precise above 14. The AD test is generally considered as more powerful than the similar Kolmogorov–Smirnov (KS) or Cramer–von Mises tests (Marhuenda et al., 2005). Moreover, the convergence of P(AN < z) towards the asymptotic distribution (N →√ ∞) is about 0.044/N, much faster than for the KS test (0.278/ N), which would had prevented us to safely use the asymptotic approxima-

tion for small values of N; (2) henceforth, N should be bigger than 1000; (3) in addition, we recommend not to use this test if the ratio between the variance of ISI and the length of spike trains is too high. For instance, for A = 1 and a Poisson process, our test should probably not be used for spike trains shorter than 100 s. Finally, the motivation for the present study came from the difficulty to manage the growing amount of spike data due to the increasing number of recording electrodes. It is common to use spontaneous activity to get reference values or to investigate some functional relations between neurons (Tomita and Eggermont, 2005). Our method is also applicable to responses to stationary stimuli. However, visual inspection of 16, 32 or more spike trains of several minutes duration for every experiment is time consuming. Moreover, the nonstationarity of the example of Fig. 5d may be quite difficult to detect visually. The proposed nonstationary test associated with the NS statistics may represent a fast and automatic way to detect any trend in the FR fluctuations: fixing at first one or several thresholds for NS, then remove or further scrutinize suspicious recording channels. Acknowledgements This research was supported by a CIHR-New Emerging Team grant, the Natural Sciences and Engineering Council of Canada (NSERC), the Alberta Heritage Foundation of Medical research, and the Campbell McLaurin Chair for Hearing Deficiencies. The authors would like to thank Greg Shaw for his valuable comments on the basic ideas of this paper. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.jneumeth.2007.02.021. References Anderson TW, Darling DA. Asymptotic theory of certain ‘goodness of fit’ criteria based on stochastic processes. Ann Math Stat 1952;23:192–212. Averbeck BB, Lee D. Neural noise and movement-related codes in the macaque supplementary motor area. J Neurosci 2003;23(20):7630–41. Baker SN, Lemon RN. Precise spatiotemporal repeating patterns in monkey primary and supplementary motor areas occur at chance levels. J Neurophysiol 2000;84(4):1770–80. Belisle P, Joseph L, MacGibbon B, Wolfson DB, du Berger R. Change-point analysis of neuron spike train data. Biometrics 1998;54(1):113–23. Benoit O, Chataignier C. Patterns of spontaneous unitary discharge in thalamic ventro-basal complex during wakefulness and sleep. Exp Brain Res 1973;17(4):348–63. Bradley J. Distribution-free statistical tests. NJ: Prentice-Hall: Englewood Cliffs; 1968. Brody CD. Correlations without synchrony. Neural Comput 1999;11(7):1537–51. Brown EN, Kass RE, Mitra PP. Multiple neural spike train data analysis: stateof-the-art and future challenges. Nat Neurosci 2004;7(5):456–61. Correia MJ, Landolt JP. A point process analysis of the spontaneous activity of anterior semicicular canal units in the anesthetized pigeon. Biol Cybern 1977;27(4):199–213. Ellaway PH. Cumulative sum technique and its application to the analysis of peristimulus time histograms. Electroencephalogr Clin Neurophysiol 1978;45(2):302–4.

B. Gour´evitch, J.J. Eggermont / Journal of Neuroscience Methods 163 (2007) 181–187 Gerstein GL. Information flow and state in cortical neural networks: interpreting multineuron experiments. In: von Seelen W, Leinhos U, Shaw G, editors. Organization of neural networks: structures and models. VCH: Weinheim; 1988. p. 53–75. Griffith JS, Horn G. An analysis of spontaneous impulse activity of units in the striate cortex of unrestrained cats. J Physiol 1966;186(3):616–34. Grun S, Diesmann M, Aertsen A. Unitary events in multiple single-neuron spiking activity. II. Nonstationary data. Neural Comput 2001;14:81–119. Grun S, Riehle A, Diesmann M. Effect of cross-trial nonstationarity on jointspike events. Biol Cybern 2003;88(5):335–51. Habets AM, Van Dongen AM, Van Huizen F, Corner MA. Spontaneous neuronal firing patterns in fetal rat cortical networks during development in vitro: a quantitative analysis. Exp Brain Res 1987;69(1):43–52. Horwitz GD, Newsome WT. Target selection for saccadic eye movements: prelude activity in the superior colliculus during a direction-discrimination task. J Neurophysiol 2001;86(5):2543–58. Kendall MG. Rank correlation methods. London: Charles Griffin; 1975. Kuffler SW, Fitzhugh R, Barlow HB. Maintained activity in the cat’s retina in light and darkness. J Gen Physiol 1957;40(5):683–702. Levine MW, Troy JB. The variability of the maintained discharge of cat dorsal lateral geniculate cells. J Physiol 1986;375:339–59. Mann HB. Non-parametric tests against trend. Econometrica 1945;13:245–59. Marhuenda Y, Morales D, Pardo MC. A comparison of uniformity tests. Statistics 2005;39(4):315–28. Marsaglia G, Marsaglia J. Evaluating the Anderson–Darling distribution. J Stat Soft 2004;9(2):1–5. Mukherjee P, Kaplan E. The maintained discharge of neurons in the cat lateral geniculate nucleus: spectral analysis and computational modeling. Vis Neurosci 1998;15(3):529–39. Nakahama H, Ishii N, Yamamoto M. Markov process of maintained impulse activity in central single neurons. Kybernetik 1972;11(2):61–72. Nakahama H, Ishii N, Yamamoto M, Saito H. Stochastic properties of spontaneous impulse activity in central single neurons. Tohoku J Exp Med 1971;104(4):373–409.

187

Nawrot MP, Aertsen A, Rotter S. Elimination of response latency variability in neuronal spike trains. Biol Cybern 2003;88(5):321–34. Palm G. Evidence, information, and surprise. Biol Cybern 1981;42(1):57–68. Penev PS, Dimitrov AG, Miller JP. Characterization of, and compensation for the nonstationarity of spike shapes during physiological recordings. Neurocomputing 2001;38–40:1695–701. Perkel DH, Gerstein GL, Moore GP. Neuronal spike trains and stochastic point processes. II. Simultaneous spike trains. Biophys J 1967;7(4):419–40. Poggio GF, Viernstein LJ. Time series analysis of impulse sequences of thalamic somatic sensory neurons. J Neurophysiol 1964;27:517–45. Ramakers GJ, Corner MA, Habets AM. Development in the absence of spontaneous bioelectric activity results in increased stereotyped burst firing in cultures of dissociated cerebral cortex. Exp Brain Res 1990;79(1):157–66. Rodden BE, Landolt JP. Note on the side test of Griffith and Horn. Kybernetik 1974;14(4):245–7. Sokal DM, Mason R, Parker TL. Multi-neuronal recordings reveal a differential effect of thapsigargin on bicuculline- or gabazine-induced epileptiform excitability in rat hippocampal neuronal networks. Neuropharmacology 2000;39(12):2408–17. Starr PA, Turner RS, Rau G, Lindsey N, Heath S, Volz M, et al. Microelectrodeguided implantation of deep brain stimulators into the globus pallidus internus for dystonia: techniques, electrode locations, and outcomes. J Neurosurg 2006;104(4):488–501. Stein RB. A theoretical analysis of neuronal variability. Biophys J 1965;91:173–94. Tomita M, Eggermont JJ. Cross-correlation and joint spectro-temporal receptive field properties in auditory cortex. J Neurophysiol 2005;93(1):378–92. Troy JB, Robson JG. Steady discharges of X and Y retinal ganglion cells of cat under photopic illuminance. Vis Neurosci 1992;9(6):535–53. Yue S, Pilon P, Cavadias G. Power of the Mann–Kendall and Spearman’s rho tests for detecting monotonic trends in hydrological series. J Hydrol 2002;259(1–4):254–71. Zhu L, Lai YC, Hoppensteadt FC, He J. Probing changes in neural interaction during adaptation. Neural Comput 2003;15(10):2359–77.