Search tips
Search criteria 


Logo of springeropenLink to Publisher's site
Journal of Computational Neuroscience
J Comput Neurosci. 2010 August; 29(1-2): 327–350.
Published online 2009 October 28. doi:  10.1007/s10827-009-0195-x
PMCID: PMC2940040

CuBIC: cumulant based inference of higher-order correlations in massively parallel spike trains


Recent developments in electrophysiological and optical recording techniques enable the simultaneous observation of large numbers of neurons. A meaningful interpretation of the resulting multivariate data, however, presents a serious challenge. In particular, the estimation of higher-order correlations that characterize the cooperative dynamics of groups of neurons is impeded by the combinatorial explosion of the parameter space. The resulting requirements with respect to sample size and recording time has rendered the detection of coordinated neuronal groups exceedingly difficult. Here we describe a novel approach to infer higher-order correlations in massively parallel spike trains that is less susceptible to these problems. Based on the superimposed activity of all recorded neurons, the cumulant-based inference of higher-order correlations (CuBIC) presented here exploits the fact that the absence of higher-order correlations imposes also strong constraints on correlations of lower order. Thus, estimates of only few lower-order cumulants suffice to infer higher-order correlations in the population. As a consequence, CuBIC is much better compatible with the constraints of in vivo recordings than previous approaches, which is shown by a systematic analysis of its parameter dependence.

Keywords: Multichannel recordings, Spike train analysis, Higher-order correlations, Cell assembly hypothesis


More than 50 years after its first conception (Hebb 1949), the idea that entities of thought or perception are represented by the coordinated activity of (large) neuronal groups has lost nothing of its attraction (e.g., Wennekers et al. 2003; Harris 2005; Sakurai and Takahashi 2006). Now known as the “assembly hypothesis”, the concept of a single, unified principle underlying such diverse computational tasks as association and pattern completion (e.g., Palm 1982), binding (e.g., Singer and Gray 1995), language processing (e.g., Wennekers et al. 2006), or memory formation and retrieval (e.g., Pastalkova et al. 2008) has inspired numerous theoretical and experimental neuroscientists. However, whether or not the dynamic formation of cell assemblies constitutes a fundamental principle of cortical information processing remains a controversial issue of current research (e.g., Singer 1999; Shadlen and Movshon 1999; van Vreeswijk 2006). While initially mainly technical problems limited the experimental surge for support of the assembly hypothesis (state-of-the-art electrophysiological setups in the last century allowed to record only few neurons simultaneously), the recent advent of multi-electrode array and optical imaging techniques reveals fundamental shortcomings of available analysis tools (Brown et al. 2004).

Based on the efficacy of synchronized presynaptic spikes to reliably generate output spikes (Abeles 1982; König et al. 1995; Diesmann et al. 1999), the temporal coordination of spike timing is a commonly accepted signature of assembly activity (e.g., Gerstein et al. 1989; Abeles 1991; Singer et al. 1997; Harris 2005). Consequently, approaches to detect assembly activity have focused on the detection of correlated spiking. Statistically, coordinated spiking of large neuronal groups has been associated to higher-order correlations among the corresponding spike trains (Martignon et al. 1995, 2000), where “genuine higher-order correlations” are assumed, if coincident spikes of a neuronal group cannot be explained/predicted by the firing rate and pairwise correlations alone. Mathematical frameworks for the estimation of such correlations exist (Nakahara and Amari 2002; Martignon et al. 1995; Schneider and Grün 2003; Gütig et al. 2003), however, run into combinatorial difficulties, as they assign one correlation parameter for each group of neurons and thus require in the order of 2N parameters for a population of N simultaneously recorded neurons. Comparing the number of parameters to the available sample size of typical electrophysiological recordings (e.g. a population of N = 100 neurons implies ~1030 parameters while 100 s of data sampled at 10 kHz provide ~106 samples) illustrates the principal infeasibility of this approach. In fact, the estimation of such higher-order correlations runs into severe practical problems even for populations of N~10 neurons (Martignon et al. 1995, 2000; Del Prete et al. 2004).

Due to the limitations of multivariate approaches, most studies in favor of cortical cell-assemblies resort to pairwise interactions. And indeed, the existence and functional relevance of pairwise interactions has been demonstrated in various cortical systems and behavioral paradigms (e.g., Eggermont 1990; Vaadia et al. 1995; Kreiter and Singer 1996; Riehle et al. 1997; Kohn and Smith 2005; Sakurai and Takahashi 2006; Fujisawa et al. 2008). However, although pairwise analysis can indicate highly correlated groups of neurons (Berger et al. 2007; Fujisawa et al. 2008), knowledge of higher-order correlations is essential to conclude on large, synchronously firing assemblies (see examples in Fig. 4 for illustration). The importance to transcend pairwise descriptions of parallel spike trains is further underscored by the strong impact of higher-order correlations on the dynamics of cortical neurons (e.g., Bohte et al. 2000; Kuhn et al. 2003). Taken together, we believe that the increasing number of simultaneously observable neurons can only be exploited by analysis techniques that go beyond mere pairwise correlations but are nevertheless applicable given the limited sample size of electrophysiological data (Brown et al. 2004).

Fig. 4
CuBIC results for three example data sets with identical rates and pairwise correlations, but different higher-order correlations (each column shows one data set). Amplitude distributions (top, logarithmic y-scale) illustrate the correlation structure ...

In this study, we present a novel cumulant based inference procedure for higher-order correlations (CuBIC). The method is specifically designed to detect the presence of higher-order correlations under the constraints of the limited sample size typical of standard experimental paradigms. This is achieved by combining three basic ideas. First, CuBIC is based on the binned and summed spiking activity across all recorded neurons (population spike count). This transforms the multivariate problem to estimate all ~2N model parameters into a parsimoniously parametrized univariate problem, thereby dramatically reducing the required sample size. Furthermore, pooling the spiking activity avoids the need for sorting the multi-unit spike trains recorded on a single electrode into isolated single neuron spike trains. Correlations among the individual spike trains are measured by the cumulants of the population spike count. Importantly, the cumulant correlations of higher order used here do not conform to the higher-order parameters of the exponential family used by, e.g., Schneidman et al. (2006). Second, we use the compound Poisson process (e.g., Snyder and Miller 1991; Daley and Vere-Jones 2005) as a flexible, intuitive, and analytically tractable model of correlated spiking activity. Here, correlations among spike trains are modeled by the insertion of additional coincident events in continuous time (Holgate 1964; Ehm et al. 2007; Johnson and Goodman 2008; Staude et al. 2010). As a consequence, the conceptual difference between the discretely sampled data (population spike count) and the continuous time spike trains is an explicit feature of the present framework. And third, we formalize the observation that even small correlations of lower order can imply synchronized spiking of large neuronal groups (Amari et al. 2003; Benucci et al. 2004; Schneidman et al. 2006).

The first and second ideas are elaborated in Section 2. Based on the compound Poisson process, Section 3 thoroughly formalizes the third idea, where we analytically derive confidence intervals for a hierarchy of statistical hypothesis tests. The tests are then combined to compute a lower bound equation M1 for the order of correlations in a given data set . Importantly, the inferred lower bound equation M2 can considerably exceed the order of the estimated correlations. Thereby, CuBIC avoids the direct estimation of higher-order correlations, which is practically infeasible for orders of equation M3, but nevertheless reveals their presence. As shown by extensive Monte Carlo simulations (Section 4), using the third cumulant suffices to reliably detect existing correlations of order > 10 in large neuronal pools (here N~100 neurons with average count-correlation coefficients of c~0.01). Thus, CuBIC achieves an unprecedented sensitivity for higher-order correlations in scenarios with reasonable sample sizes, i.e. experiment durations of T < 100 s. The Discussion section relates the cumulant-based correlations used here to the higher-order interaction parameters of the exponential family (e.g., Martignon et al. 1995; Nakahara and Amari 2002; Shlens et al. 2006) and critically discusses the implications of the hypothesized compound Poisson process. Preliminary results have been presented previously in abstract form (Staude et al. 2007).

Measurement & model


Assume an observation of a large number of parallel spike trains. To measure correlation, we describe such a population as a succession of “patterns”, equation M4 (T denotes transpose), one pattern for every time bin of width h. The components of X are the binned, discretized spike trains, i.e. Xi(s) is the spike count of the ith neuron in the interval [s,s + h). Given these patterns, we define the population spike count Z(s) at s as the total spike count in the sth bin (Fig. 1)

equation M5
Fig. 1
Schema of the compound Poisson process and its measurement. Left half: spike event times (horizontal bars) of individual neurons x1(t),...,xN(t) and tick marks of the carrier process z(t) (top) with the associated amplitudes (numbers above the ticks), ...

The variable Z(s), from here on referred to as the “complexity” of the population at s, counts the number of spikes that fall into the interval [sh,(s + 1)h), irrespective of the neuron IDs that emitted these spikes. In the case where the Xi are binary (“1” for one or more spikes in the bin, “0” for no spike), Z(s) is simply the number of neurons that spike in the time slice. As opposed to most other frameworks for correlation analysis (e.g., Aertsen et al. 1989; Martignon et al. 1995; Grün et al. 2002a; Nakahara and Amari 2002; Shlens et al. 2006), however, the method presented in this study does not assume binary variables. As a consequence, our term “complexity” as the “number of spikes per bin in the population” does not comply with the “number of active neurons per bin” used in binary frameworks.

To regard the values of Z(s) in different bins as i.i.d. random variables, we here assume that Z(s) and Z(s + k) are independent for k  0 (zero memory), and that the distribution of Z(s) does not depend on the time bin s (stationarity). We name the resulting distribution of population spike counts

equation M6

the “complexity distribution” of the population. The validity of these assumptions with respect to real spike trains, and potential adaptations of CuBIC, are discussed in Section 5.3.

Despite ignoring the specific neuron IDs that contribute to the patterns X(s), the complexity distribution nevertheless contains information about the correlation structure of the population. For instance, if the counting variables {Xi}i [set membership] {1,...,N} are independent Poisson counts, the corresponding population spike count, being the sum of the independent Poisson variables, is again a Poisson variable, and thus the complexity distribution fZ is a Poisson distribution. Correlations change the relative probabilities for patterns of high and low complexities as compared to the independent case, as can be seen by comparing the complexity distribution of the correlated population and its independent Poisson fit in Fig. 1 (blue bars and dashed gray line, respectively; see also Grün et al. 2008a and Louis and Grün 2009). We quantify such deviations from independence by means of the cumulants κm[Z] of the population spike count Z.

Correlations and cumulants

Like the more familiar (raw) moments equation M7 of a random variable Z, the cumulants κm[Z] characterize the shape of its distribution (see Appendix A and, e.g., Stratonovich 1967; Gardiner 2003). Most common are the first two cumulants, the expectation and the variance. These admit simple expressions in terms of the moments: equation M8 and equation M9. Similar expressions for higher cumulants, however, are exceedingly complicated (see Stuart and Ord 1987 for explicit expressions for m  10).

The most important property of cumulants, and their advantage over the raw moments is that the cumulant of the sum of independent random variables is the sum of their cumulants. For m = 2 and equation M10, for instance, we have the well-known variance-covariance relationship

equation M11

Equation (1) shows that equation M12 if equation M13, which is the case if the Xi are (pairwise) uncorrelated. Cumulant correlations of higher order, the so-called “mixed” or “connected” cumulants, are a straightforward generalization of the covariance in exactly this sense: if the population has neither pairwise nor triple-wise cumulant correlations, then equation M14 (see Appendix A for a concise definition, and Gardiner 2003, for a general introduction). Higher-order cumulant correlations generalize the covariance also with respect to the fact that equation M15 implies equation M16. That is, if the multivariate random vector X decomposes into independent subgroups, then the cumulant correlations vanish (Streitberg 1990). For notational consistency, we will from now on stick to the cumulant notation, e.g., use “first/second cumulant” instead of the more familiar terms “mean/variance”.

A further consequence of Eq. (1) is that κ2[Z] is influenced only by the single process statistics (via equation M17) and pairwise correlations (via equation M18). No higher-order correlations contribute to the second cumulant. This holds also for higher cumulants, i.e. κm[Z] depends on correlations among the Xi of maximal order m. And finally, in the same way that κ2[Z] depends on pairwise correlations and the single process statistics, also κm[Z] does not measure pure mth order correlations, but depends on correlations of all orders up to m. While a correction of the second cumulant for the influence of the single process statistics would be straightforward (subtracting equation M19 in Eq. (1), see Appendix B), correcting higher cumulants for the influence of correlations of lower order is exceedingly complicated. We therefore employ a parametric model for Z, the compound Poisson process (see next section), the parameters of which can be interpreted straightforwardly in terms of higher-order correlations among the Xi.

We would like to stress that the higher-order correlations defined by cumulants differ strongly from the higher-order parameters of the exponential family used by, e.g. Martignon et al. (1995), Nakahara and Amari (2002), Shlens et al. (2006). The relationship between these two frameworks is discussed in more detail in Section 5.1 (see also Staude et al. 2010).


As opposed to the discretized, i.e. binned, population spike count Z(s) of the previous section, the proposed model operates in continuous, i.e. unbinned time. That is, we model the process equation M20, where equation M21 denotes the ith unbinned, continuous-time spike train (i = 1,...,N) with spike-event times equation M22. The model we propose for z(t) is that of a compound Poisson process (CPP)

equation M23

where the event times tj constitute a Poisson process, and the marks aj are i.i.d. integer-valued random variables, drawn independently for all tj. The marks aj determine the number of neurons that fire at time tj, and will be referred to as the “amplitude” of the event at time tj. The probability that an event has a specific amplitude is determined by the amplitude distribution fA, i.e. equation M24 (see Fig. 1). The Poisson process that generates the events tj is called the “carrier process” of the model and its rate ν is the “carrier rate”. Processes of this type are also referred to as generalized, or marked, Poisson processes (see e.g. Snyder and Miller 1991 for a general definition and Ehm et al. 2007 for an application to spike train analysis).

To interpret the CPP as the lumped process of a correlated population of N spike trains, and to utilize the CPP for the generation of artificial data, the events of z(t) are assigned to individual processes xi(t) (i = 1,...,N). The simplest model draws the aj process IDs that receive a spike at tj as a random subset from {1,...N}, independently for all event times tj (Fig. 1). This results in a homogeneous population of correlated Poisson processes, where all processes have the same rate, all pairs of processes have identical correlations, and so on. The SIP/MIP models presented in Kuhn et al. (2003) are special cases of this more general model (see also Appendix B for examples). As the inference procedure CuBIC presented in the remainder of this study is based solely on the summed activity, however, it is not affected by the details of the copying procedure. In particular, CuBIC does not presuppose that the population is homogeneous.

Relating measurement and model

To interpret the continuous-time model parameters ν and fA in terms of correlations among the counting variables Xi, we now relate the former to the cumulants of the population spike count Z. First of all, note that as the carriers process z(t) is a Poisson process, the population spike counts in different bins are i.i.d. For l = 1,...,N, define the processes yl(t) that determine all event times tj in Eq. (2) with given amplitude aj = l (compare Appendix A). Then the CPP z(t) admits the representation equation M25. As a consequence, the discretized population spike count satisfies equation M26, where the Yl are the counting variables obtained from the processes yl(t) using a bin size h. As the event times tj of z(t) follow a Poisson process and the subsequent amplitudes aj are independent, the yl(t) are independent Poisson processes. The rate νl of yl(t) is given by νl = fA(lν (l = 1,2,...). Hence, the Yl are Poisson variables with parameter νlh. As a consequence, all cumulants of Yl are identical and given by κm[Yl] = νlh (m = 1,2,...; Gardiner 2003). The scaling behavior of the cumulants equation M27 for all m (Mattner 1999) yields

equation M28

Using νl = fA(lν and equation M29, we finally have

equation M30

Equation (3) is the central equation of this study and requires a few more remarks. First, it implies that the first two cumulants of the population spike count, κ1[Z] and κ2[Z], are determined solely by the carrier rate ν and the first two moments of the amplitude distribution, μ1[A] and μ2[A]. As κ1[Z] and κ2[Z] determine firing rates and pairwise correlations in the population (Appendix B, Eqs. (24) and (26)), we conclude that all CPP models with identical carrier rate ν and identical μ1[A] and μ2[A] generate populations with identical rates and pairwise correlations. Differences in higher-order correlations can thus be modeled by choosing different higher moments for fA (see Section 2.4 for the precise relationship between the entries of the amplitude distribution and higher-order correlations in the population). This makes the CPP a very flexible and convenient model to generate artificial data sets with identical firing rates and pairwise correlations, yet different higher-order correlations (see Section 4 and Appendix B for examples, and Kuhn et al. 2003; Ehm et al. 2007; Staude et al. 2010 for alternative parametrizations).

Second, the moments of the strictly positive, integer-valued variable A increase with the order, i.e. satisfy μm[A]  μm + 1[A] for m  1. By Eq. (3), this implies that also the cumulants of Z increase with the order, i.e. κm[Z]  κm + 1[Z]. This is a reformulation of the fact that insertion of joint spikes as done in the CPP generates only positive correlations (e.g., Brette 2009; Johnson and Goodman 2008, Johnson and Goodman, unpublished manuscript), and it shows that certain combinations of cumulants cannot be realized by the CPP.


The conceptual difference between the measured data (the counting variables Xi and the complexity distribution fZ) and the parametric model that we assume to underlie this measurement (the CPP) is crucial for the remainder of this study. To illustrate this difference, note that fZ depends on the bin size h used in the discretization. Larger bin widths increase the probabilities for high complexities, while smaller bin widths increase the probability for empty bins. The parameters of the CPP, on the other hand, do not change with the bin size. For instance, if all events of z(t) have amplitude aj = 1, i.e. the amplitude distribution has a single peak at ξ = 1, the population consists of independent Poisson processes, irrespective of the bin size used for the analysis, or even the number of recorded neurons. In other words: patterns of complexity Z  2 can occur by chance, while events of amplitude ξ  2 do not occur by chance but imply correlations in the population. The order of the correlation, in turn, is determined by the amplitudes of the events in the corresponding CPP (see Appendix A for the proof):

Theorem 1 Letequation M31 be a compound Poisson process with amplitude distribution fA, and let equation M32 be the vector of counting variables obtained from the xi(t) with a bin width h. Then the components of X have correlations of order m if and only if fA assigns non-zero probabilities to amplitudes  m.

The above theorem confirms the intuitive conception that, within the framework of the CPP, correlations of a certain order m require injected coincidences into at least m processes (events of amplitudes ≥ m). Also, a population with only pairwise but no higher-order correlations has events of maximal amplitude m = 2. Importantly, correlations in the above theorem are defined strictly on the basis of the discretized counting variables Xi. As a consequence, they do not resolve (and do not depend on) the perfect temporal precision of the coincident events in the CPP: if the events of z(t) were copied into the individual processes with a temporal jitter that is small with respect to the bin size h, the correlations among the counting variables will hardly be affected (see Section 5.3.3 for a more detailed discussion). Note also that events of amplitude ξ induce not only correlations of order ξ, but also of orders < ξ.

Cumulant based inference of higher-order correlations (CuBIC)

This section describes our cumulant based inference procedure for higher-order correlations (CuBIC). The outcome of CuBIC is a lower bound equation M33 on the order of correlation in the spiking activity of large groups of simultaneously recorded neurons. As will be shown in Section 4, this lower bound can exceed the order of the cumulants that were estimated for the inference. Thereby, CuBIC can provide statistical evidence for large correlated groups without the discouraging requirements on sample size that direct tests for higher-order correlations have to meet. This is achieved by exploiting constraining relations among correlations of different orders. For illustration, consider as an example the extreme situation of 4 simultaneously recorded neurons, where all neuron pairs have a correlation coefficient of c = 1. As c = 1 implies identity for all pairs of spike trains, all four spike trains of this example must be identical (Fig. 2, left). In other words, a spike in one neuron implies joint spike events in all the other neurons. Thus, the data must have correlation of order 4.

Fig. 2
Strong pairwise correlations imply higher-order correlations. Raster plots show populations that have a single event amplitude that differs in the three panels (left: all events have amplitude 4, middle: amplitude 3, right: amplitude 2). Note that an ...

The key observation of this example is that the order of correlation we inferred (equation M34) exceeds the order m of the measured correlations (here m = 2). CuBIC formalizes this example in the framework of the CPP, and generalizes it in two aspects:

  1. Assume the correlation coefficients had been somewhat smaller than 1: do we still need correlation of order 4 to explain the measured pairwise correlation? Or would correlation of order 3 as in Fig. 2, middle, or even of order 2 as in Fig. 2, right, suffice to explain the measured correlations? In other words: What order of correlation is minimally required to explain the measured pairwise correlations?
  2. Can we formulate a similar reasoning for measured correlations of higher order? What order of correlation is minimally required to explain the measured correlations of third (fourth,...) order?

The next section details our approach to answer the first question. The results are then generalized to answer the second question (Section 3.2). Finally, we explain how CuBIC combines the results to infer a lower bound for the order of correlation in a given data set (Section 3.3).

Testing pairwise correlations (m = 2)

We approach the first aspect with a hierarchy of statistical hypothesis tests, labeled by the integer ξ. For fixed ξ, the null hypothesis equation M35 states that the pairwise correlations (thus the superscript “2”) in the population are compatible with the assumption that there is no correlation beyond order ξ; the alternative equation M36 states that correlation of order higher than ξ is necessary to explain the pairwise correlations. Rejection of e.g. equation M37 in favor of equation M38 thus implies the existence of correlation of at least order 5. The test statistics to decide between these alternatives is the second sample cumulant of the population spike count, and the compound Poisson process provides the framework to analytically derive the required confidence bounds.

The null hypothesis equation M39

The formulation of a null hypothesis requires a careful distinction between the random variable that describes the experiment, and the data sample that is tested against the hypothesis. For a fixed value of the test parameter ξ, this section formulates the null hypothesis equation M40, using the random variable Z′ that describes the population spike count of a neuronal population. Section 3.1.2 explains how a given data sample {Z1,...,ZL} is tested against equation M41.

The null hypotheses equation M42 is based on the CPP model whose population spike count has the maximal second cumulant under the constraints that (a) the first cumulant equals that of a given population spike count Z′, and (b) the maximal order of correlation in the model is ξ. Using the amplitude distribution fA to parametrize its correlation structure, this CPP model is the solution of the constrained maximization problem

equation M43

where ν and fA are the CPP parameters that determine the population spike count Z. We rewrite Eq. (4) by observing that the second constraint (equation M44) together with Eq. (3) implies equation M45. With νl = fA(lν (Section 2.3), we thus have

equation M46

where the dot denotes the standard scalar product. Using Eq. (5) and the vector notation equation M47 and equation M48, the maximization problem of Eq. (4) becomes

equation M49

In Eq. (6), both the function to maximize (equation M50) and the constraint (equation M51) depend linearly on the parameters equation M52. Problems of this type, so-called Linear Programming Problems, are uniquely solvable, e.g. using the Simplex Method (Press et al. 1992, Chapter 10.8). The solution yields the upper bound for the second cumulant equation M53 and the corresponding parameter vector equation M54. The carrier rate and amplitude distribution of the CPP that maximizes Eq. (6) are then given by equation M55 and equation M56.

Assume that the combination of firing rates and pairwise correlations in the population can be realized with correlations of order ≤ ξ. Then, the second cumulant of its population spike count Z′ must be smaller than the upper bound equation M57 computed in the previous section. We thus formulate the null hypothesis

equation M58

The alternative hypothesis

equation M59

states that, within the framework of the CPP model, the pairwise correlations in the population imply the presence of correlation of order > ξ.

Test statistics and their distributions

The derivations in the previous section require the cumulants κ1[Z′] and κ2[Z′] of the random variable Z′. Given only a data sample {Z1,...,ZL} of size L, we estimate these quantities by the standard sample mean and (unbiased) sample variance (Stuart and Ord 1987)

equation M60

The test of the data sample against equation M61 requires the distribution of the test statistics k2 under the null hypothesis. To derive this distribution, recall that equation M62 and equation M63, where the κi[Z′] are the unknown cumulants of the variable Z′ that underlies the sample (Stuart and Ord 1987). The trick of CuBIC is to assume in equation M64 that Z′ is the population spike count of the CPP that solves Eq. (6) after the unknown cumulant κ1[Z′] has been substituted with its estimate k1. Under this assumption, all cumulants of Z′ can be computed by inserting the model parameters ν* and equation M65 into Eq. (3), which yields

equation M66

Given sample sizes of L > 10.000 (roughly corresponding to 10–100 s of data, see Section 4.2 and Fig. 6(d)), the distribution of k2 under equation M67 is well approximated by a normal distribution. Taken together, under equation M68 the test statistics k2 is normal with mean equation M69 and variance given by Eq. (7), such that the p-value for the test against equation M70 is

equation M71
Fig. 6
Parameter dependence of CuBIC for m = 3. Shown are the percentiles ξ05 and ξ95 (lower and upper curves of same color, respectively; see Fig. 5(b)), as a function of (a) the population firing rate Λ, (b) ...

The rejection of the null hypothesis equation M72 for a given ξ implies that the pairwise correlations in the data are not compatible with the assumption that there is no correlation beyond order ξ. Thus, rejecting equation M73 implies that ξ + 1 is a lower bound for the order of correlation.

Testing higher-order correlations (m > 2)

We now generalize the tests with measured pairwise correlations (m = 2) to exploit also estimated correlations of higher order (m > 2). The difference to the tests with m = 2 is that the upper bound for the mth cumulant involves all cumulants up to order m  1. The generalization of the maximization problem (Eq. (4)) reads

equation M74

where, as before, the κi[Z′] are the cumulants of the random variable Z′. Using the ν-parametrization of Eq. (5), we can rewrite Eq. (8) as the Linear Programming Problem

equation M75

where Ξm  1 is a (ξ×m  1)-dimensional coefficient matrix with entries equation M76, and equation M77 is the column-vector of the first m  1 cumulants of Z′. Solving Eq. (9) yields both the maximal mth cumulant equation M78 and the rate vector equation M79, which again can be used to compute the corresponding model parameters ν* and equation M80.

In direct analogy to the case m = 2, the (m,ξ)-null hypothesis equation M81 states that correlation of orders 1,2,...,m are compatible with the assumption that there is no correlation beyond order ξ. This implies that κm[Z′] falls below the upper bound equation M82, and we define

equation M83

To test a sample {Z1,...,ZL} against equation M84, we estimate correlations of order m by the mth sample cumulant of the population spike count, the so-called the mth k-statistics km (Stuart and Ord 1987; the sample mean and unbiased sample variance defined above are the first two k-statistics). And again, we assume that Z′ is the population spike count of the CPP model that solves Eq. (4) after the unknown κi[Z′] have been replaced by the estimates ki (i = 1,...,m  1). Then, the test statistics km is normally distributed with mean equation M85 and variance equation M86, where expressions for equation M87 can be found in the literature (Stuart and Ord 1987). The p-value of equation M88 is thus given by

equation M89

Note that after substituting the true cumulant vector equation M90 by its estimate equation M91 in Eq. (9), the resulting constraint equation M92 consists of m  1 equations for the ξ positive parameters ν1,...,νξ. Unfortunately, there can be combinations of estimated cumulants for which these equations cannot be solved. In this case, Eq. (9) does not posses a solution, and the corresponding null hypothesis cannot be tested (see also Section 3.3).

Computing the lower bound

In the preceding section, CuBIC was presented as a collection of hypothesis tests equation M93, labeled by the indices m (the order of the estimated cumulant) and ξ (the maximal order of correlation in the null). We now combine these tests to infer a lower bound equation M94 for the order of correlation in a given data set. To do so, recall that rejecting equation M95 for given m and ξ implies that the combination of the first m cumulants requires correlation of order > ξ. As a consequence, every rejected hypothesis equation M96 implies that ξ + 1 is a lower bound for the highest order of correlation in the data. The question thus is: which of these bounds should we use?

Evidently, we aim to infer the highest order of correlation that is present in the data. Hence, the objective is to find the maximum of all the lower bounds that were obtained from the hierarchy of tests. We thus aim for the pair (m,ξ) with the highest value of ξ such that the corresponding null hypothesis equation M97 is rejected.

A conceptual algorithm for this search is presented in Fig. 3. It consists of two nested loops: for a fixed order of the estimated cumulant m, the inner loop searches the highest ξ for which equation M98 is rejected; the outer loop iterates over subsequent orders m. The free parameters of the algorithm are the test level α, and, to ensure termination of the loops, upper bounds for the cumulant (mmax) and the maximal order of correlation assumed in the null (ξmax; see Section 5 for their choice). After initializing the test variables (m = 2, ξ = 1, equation M99), we estimate the first m cumulants from the data by computing the corresponding k-statistics. Next, we check if Eq. (9) is solvable for the current values of m and ξ (yellow box). This check consists of two steps. The first step checks if the first m  1 estimated cumulants increase with order m, a requirement of the CPP model assumed to underlie the data (left rhomb in yellow box; compare last paragraph of Section 2.3). If this is not the case, then k1  k2  ...  km  1 is false for all m  m, which implies that Eq. (9) cannot be solved for m  m. If the solution, i.e. the maximal cumulant equation M100 and the model parameters ν* and equation M101, are not available, however, the corresponding p-values (Eq. (10)) cannot be computed. In this case, no hypothesis equation M102 with m  m can be tested, and the procedure is terminated. If the first m  1 estimated cumulants are in principle compatible with the CPP, we check in the second step if the current value of ξ can solve Eq. (9) (right rhomb in yellow box). If this is not the case, ξ is increased by 1 until either Eq. (9) can be solved, or ξ reaches its upper bound ξmax. In the former case, the data are tested against equation M103 (blue box), in the latter case the test is skipped. At rejection of equation M104 (rightward arrow at the box pm,ξ < α) we set equation M105 and, unless ξ reached the upper bound ξmax, repeat the inner loop with ξ set to ξ + 1. The loop terminates with equation M106, which is the largest lower bound for the current value of m. As we reject the “zeroth” hypothesis equation M107 by convention (every data set has correlations of order > 0), this holds even if the tests for the current value of m were skipped. The inner loop is repeated with new m set to m + 1, until m reaches the upper bound mmax. Finally, the lower bounds obtained for the different values of m are compared, and their maximum is returned as the absolute lower bound equation M108. A MatLab-implementation of the proposed algorithm is available upon request.

Fig. 3
Procedure to infer the highest lower bound equation M109 for the order of correlation in a given data set. The pair (m,ξ) with the highest ξ such that the corresponding null hypothesis equation M110 is rejected is found by two nested loops. For any given m, the ...

Simulation results


Before we present a systematic analysis of the proposed procedure in Section 4.2, we illustrate CuBIC’s sensitivity by a comparison of three sample data sets (Fig. 4). All three data sets are based on the same number of neurons (N = 100) with identical firing rates (λ = 10 Hz) and pairwise correlations (c = 0.01 in NC = 30 and c = 0 in the remaining 70 neurons). The difference between the data sets lies solely in their higher-order correlation structure, i.e. their maximal order of correlation (Fig. 4, top panels). All sets have amplitudes at ξ = 1 that generate independent “background spikes”. Correlations are induced by additional events of amplitudes ξsyn > 1 that differ across data sets: ξsyn = 2 in Set1 (left), ξsyn = 7 in Set2 (middle), and ξsyn = 15 in Set3 (right). Thus, Set1 has correlation of order 1 and 2, Set2 has correlations up to order 7, and Set3 has correlations up to order 15. The identical pairwise correlations in the three data sets are achieved by adjusting the rate equation M114 of the higher-order events (see Appendix B for the relationship between the model parameters ν,fA and the population statistics N,NC,λ and c). Test parameters were set to ξmax = 15 and mmax = 4.

The similarity of the raster displays and population spike counts (Fig. 4, 2nd and 3rd from top) of the three data sets may illustrate that such standard visualization methods cannot reveal the differences in the higher-order properties of the data sets. Larger probabilities for patterns of high complexity induced by the higher-order correlations present in Set2 and Set3 are visible in the complexity distributions (third row of panels) only when plotted on a logarithmic scale (green lines). The low rates of the higher-order events imposed by the small pairwise correlations (ν7 = 2.07 Hz in Set2, ν15 = 0.41 Hz in Set3) makes their detection in the three distributions (bars) very difficult on a linear scale.

The identical pairwise correlations in all three data sets are reflected in identical results for tests with m = 2. For all data sets, equation M117 is rejected (p < 0.05, indicated by red marks) only for ξ = 1, implying that the pairwise correlations are significant and events of amplitude 1 are not enough to explain the data. As equation M118 is retained for ξ > 1 in all data sets (p > 0.05), the lower bound for m = 2 is equation M119 (see Fig. 3), i.e. the tests with m = 2 do not imply correlation of order ≥ 2 in any of the data sets.

Test results with m = 3 differ strongly for the three data sets (Fig. 4, lower panels, right graphs). For Set1, only equation M120 is rejected and all equation M121 for ξ  2 are retained. Hence equation M122. For Set2, equation M123 is rejected for ξ = 1,...,6, yielding equation M124, while for Set3 the rejection of all null-hypotheses for ξ  12 yields equation M125. In all tests with m = 4 (data not shown), the smallest values for ξ that solved Eq. (9) yielded non-significant p values (for Set1 p4,2 = 0.5, for Set2 p4,9 = 0.8, for Set3 p4,16 = 0.67). Hence, no hypothesis with m = 4 was rejected except for equation M126 which is rejected by convention. Thus, equation M127 in all data sets. For the total lower bounds equation M128, we thus obtain equation M129 for Set1, equation M130 for Set2, and equation M131 for Set3.

The lower bound corresponds to the maximal order of correlation in Set1 (ξsyn = 2) and Set2 (ξsyn = 7). Only for Set3 the lower bound falls short of the true maximum by 2 (ξsyn = 15), but nevertheless indicates correlations of higher order than in Set2. Taken together, the three examples illustrate that CuBIC reliably detects differences in higher-order statistics, even if firing rates (here λ = 10 Hz) and very small pairwise correlations (here c = 0.01 in 30 out of 100 neurons) are identical.

Parameter dependence of test performance

Several parameters are likely to influence CuBIC’s performance. For instance, the number of recorded neurons, N, the collection of firing rates, λ1,...,λN, the duration of the experiment, T, and the bin size, h, influence the amount of available data and can thus be expected to effect test results. Also the correlation structure, i.e. the rate and the order of injected coincidences, parametrized by the N  1 probabilities of the amplitude distribution fA, is likely to affect the performance of the method. However, some of the parameters mentioned above are in fact redundant with respect to the test statistics used here, i.e. the k-statistics of the population spike count Z. For instance, Z depends only on the summed firing rate equation M132, but not on the precise combination of the rates of individual neurons. To analyze the parameter dependence of CuBIC (Section 4.2.3), we focus on parameters that are relevant in experimental data and avoid the above-mentioned redundancies. Furthermore, we restrict our analysis to tests with the third cumulant, thus study the inner loop in Fig. 3 for m = 3. To keep notation simple, we drop the subscript 3 and write equation M133 in the remainder of this section. The upper bound for ξ was ξmax = 30 in all simulations.


We investigate the performance of CuBIC with respect to five parameters: the duration of the experiment, T, the total firing rate of the population, equation M134, the bin width used to discretize the data, h, and two parameters, ξsyn and ρ, that characterize the correlation structure (amplitude distribution, fA) of the data. We here consider amplitude distributions that consist of two isolated peaks only: a “background-peak” at ξ = 1 and a second peak at ξ = ξsyn that determines the maximal order of correlation in the data (see for examples Fig. 4, top row). The relative height of the two peaks is parametrized by the moments ratio

equation M135

where equation M136 is the ith moment of A.

Substituting the moments μi[A] in Eq. (11) by the cumulants of Z via Eq. (3) shows that ρ equals the Fano Factor of the population spike count, i.e. equation M137 (see e.g., Kumar et al. 2008; Kriener et al. 2008). This interpretation is valid for arbitrary amplitude distributions, and provides a method to estimate ρ from data. For a more intuitive interpretation of ρ, consider a population of N neurons with identical firing rates, where a sub-population of NC neurons forms a homogeneously correlated subgroup, while the remaining N  NC neurons fire independently (compare rasters displays in Fig. 4). In this case ρ relates to the (average) pairwise count correlation coefficient c of the correlated subgroup

equation M138

For a derivation of this relation see Appendix B, Eq. (29). In case of a completely homogeneous population with N = NC, this simplifies to equation M139 (compare e.g., Kumar et al. 2008; Kriener et al. 2008). In either case, the correlation parameter ρ determines the effect of higher-order events on the pairwise correlations in the population.

Quantification of test performance

We asses CuBIC’s performance by Monte-Carlo techniques. That is, we compute the lower bound for the order of correlation equation M140 (Figs. 3 and and5(a))5(a)) for each of 1000 data sets that are simulated from the CPP model with identical parameter combinations. The resulting (estimated) distribution equation M141 then indicates the range of lower bounds for this particular parameter combination. We thus study the parameter dependence of CuBIC’s performance by the parameter dependence of the distribution of equation M142 (Fig. 5(b)). For the specific parameter combination of Fig 5(a), for instance, lower bounds fall almost exclusively between equation M143 and equation M144 (Fig. 5(b), bottom panel). Increasing the correlation parameter ρ gradually sharpens and shifts equation M145 to higher values of equation M146 (Fig. 5(b)). CuBIC performs optimally if the lower bound corresponds to the maximal order of correlation in all Monte-Carlo simulations, i.e. if equation M147 for all data sets. In this case, equation M148 is a delta-peak located at ξsyn (Fig. 5(b), top panel).

Fig. 5
Quantification of test performance. (a) p-values for tests against equation M149 and the resulting lower bound for the order of correlation equation M150 (here equation M151) for one data set (parameters: ξsyn = 15, Λ = 103 Hz, T = 100 s, ...

To further reduce complexity, we will discuss the parameter dependence of equation M156 by means of the two percentiles (black triangles in Fig. 5(b))

equation M157

CuBIC performs optimally if equation M158 is a delta peak as in Fig. 5(b), top panel, i.e. if ξ05 + 1 = ξsyn = ξ95. More generally, a value of ξ05 [dbl greater-than sign] 0 indicates that CuBIC reliably infers the existence of higher-order correlations (equation M159 in 95% of all cases), while a value of ξ95  ξsyn guarantees that the actual order is not overestimated (equation M160 in less than 5% of all cases, no false positives).

As mentioned before (Section 3.3), certain combinations of sample cumulants, e.g. k2 < k1, represent untestable data sets and yield lower bounds of equation M161. To avoid border effects for extremely small pairwise correlations (ρ < 1.05), we set equation M162 not only if k2 < k1 but also if k2 does not significantly exceed k1, i.e. if equation M163 is retained (marked as dotted lines in Fig 6(c)).


Figure 6 shows the dependence of CuBIC’s performance on the population firing rate Λ (Fig. 6(a)), the duration of the experiment T (Fig. 6(b)), the correlation strength (population Fano Factor) ρ (Fig. 6(c)) and the bin width h (Fig. 6(d)) for different values of ξsyn (crosses in Fig. 6(a–d) denote the default parameters that are fixed while only one of them is varied).

Changes in the population firing rate Λ, the duration of the experiment T, and the correlation strength ρ have qualitatively very similar effects on the performance (Fig. 6(a), (b), and (c), respectively). For all analyzed values of ξsyn (dark blue, green, and yellow lines) and small values of either of the parameters, i.e. Λ = 10 Hz (Fig. 6(a)), T = 1 s (Fig. 6(b)) or ρ = 1.0087 (Fig. 6(c)), the distributions of lower bounds equation M165 extend from ξ05 = 0 (lower lines) to ξ95 < ξsyn (upper lines below colored triangles). Hence, lower bounds might not exceed the trivial value of equation M166 for such small parameter values. Increasing either of the parameters gradually shifts equation M167, i.e. ξ05 and ξ95, to higher values of equation M168, thus indicating improved performance. The exact quality of the performance, captured by the width of equation M169, i.e. the distance between ξ05 and ξ95, and in particular its difference from ξsyn (colored triangles, compare also Fig. 5(b)), depends crucially on the order of correlation ξsyn present in the data.

For ξsyn = 7 (dark blue lines), CuBIC performs optimally (ξ05 + 1 = 7 = ξ95, compare Fig. 5(b)) for experiment durations of T > 100 s (Fig. 6(b)) and a wide range population Fano Factors (ρ  1.17, Fig. 6(c)). Such optimality is not guaranteed for the analyzed range of the population firing rate Λ, although population rates above Λ = 400 Hz yield a lower percentile of ξ05 = 5, indicating close-to-optimal performance (equation M170 in 95% of all cases, Fig. 6(a)). Smaller values of either of the parameters shift equation M171 to lower values, indicating impaired performance. Nevertheless, as the lower percentile ξ05 is considerably larger than zero even for population rates as low as Λ = 20 Hz, experiment durations of T = 6 s and population Fano Factors of ρ = 1.02, CuBIC reliably detects the presence of higher-order correlations for a wide range of parameter combinations.

In the other extreme, where the same pairwise correlations (identical value of ρ) are realized by correlations of order ξsyn = 30 (yellow lines), equation M172 falls short of its optimum (ξ05 + 1 = 30 = ξ95, see Fig. 5(b), top panel) even for high values of all investigated parameters, i.e. population rates of Λ = 104 Hz, experiment durations of T = 103 s and correlation strengths of ρ = 5.35 (Fig. 6(a), (b) and (c), respectively). Note, however, that also here ξ05 [dbl greater-than sign] 0 even for moderate values of Λ, T and ρ, implying that the presence of higher-order correlations is reliably detected, albeit not their actual order. For the default values Λ = 103 Hz, T = 100 s and ρ = 1.087 (yellow crosses in Fig. 6(a–d)), for instance, lower bounds for ξsyn = 30 fall between equation M173 and equation M174 in 90% of the analyzed data sets (ξ05 = 19, ξ95 = 24).

Test performance for ξsyn = 15 (green lines) is somewhere in between the often optimal results for ξsyn = 7 and the more indicative results for ξsyn = 30. In short, optimal performance is achieved only for high population Fano Factors (ρ  3.75, rightmost black triangle in Fig. 6(c)), while the presence of higher-order correlations is reliably detected (equation M175) for a wide range of parameters (Λ > 70 Hz, Fig. 6(a), T  80 s, Fig. 6(b), and ρ  1.02 Hz, Fig. 6(c)).

The qualitative similarity between Fig. 6(a) and (b) result from the identical effect of the varied parameters, Λ and T, on the expected spike count n = ΛT (Fig. 6(a), (b), top x axis). As n affects the quality of the estimators ki , both Λ and T influence test results in a similar manner. The similar ρ-dependence results from the fact that increasing ρ increases the probability for high-amplitude events fA(ξsyn), which simplifies their detection (Fig. 6(c)).

Changes in the bin width h (Fig. 6(d)) affect CuBIC’s performance for all values of ξsyn in a very similar way. Increasing h decreases the lower percentile from ξ05 [dbl greater-than sign] 0 at h = 0.1 ms to ξ05 = 0 for h = 100 ms, irrespective of ξsyn. The upper percentiles ξ95 (upper lines) can even show non-monotonic behavior: for ξsyn = 15 (green) and ξsyn = 30 (yellow), ξ95 increases for values between 1  h  60 ms and decreases for h > 60 ms. The median of equation M176, which lies roughly in the middle of ξ05 and ξ95, however, remains approximately constant for h  60 ms, at least for ξsyn = 15 and ξsyn = 30 (green and yellow, respectively). Therefore, the average lower bound remains approximately constant as long as h  10 (L  104, top axis), but only the variability of equation M177, i.e. the distance between ξ05 and ξ95, increases for this range of h.

We wish to stress that the larger variability of equation M178 for increased bin sizes is not so much due to the increased average spike count per bin μ = Λh, but rather the resulting decreased sample size equation M179 (Fig. 6(d), top x-axis). The expected spike count μ = Λh changes with the population rate Λ and the bin width h in the exact same manner. However, while increasing μ via Λ (Fig. 6(a)) from the default value Λ = 103 Hz (crosses) to Λ = 104 Hz does not alter test performance for ξsyn = 7 (dark blue lines), increasing μ via the bin width from h = 1 ms to h = 10 ms (Fig. 6(d)) drastically reduces the lower percentile ξ05 (lower dark blue line in Fig. 6(d)). The decrease in sample size increases the variance of our test statistic (compare Eq. (7)) and therefore decreases test power.

Summary of results

The simulation results of the previous section show that CuBIC reliably detects the presence of higher-order correlations for a wide range of parameter values. As CuBIC is not designed to directly estimate the order of correlation, the inferred lower bound might not always correspond to the maximal order of correlation present in a given data set. The most important parameter is the sample size L, which is required to be ≥ 104 to reliably detect higher-order correlations (Fig. 6(d)), corresponding to e.g. an experiment of T = 50 s duration analyzed with a bin size of h = 5 ms. Furthermore, total spike counts of n = 104 (e.g. N = 100 neurons firing at 1 Hz for T = 100 s; Fig. 6(a) and (b)) and a correlation parameter (population Fano Factor) in the order of ρ = 1.087 suffice for the test to reliably detect the presence of higher-order correlations (lower percentiles ξ05 > 0, green and dark blue lower lines Fig. 6(a–c)) if correlations in the data are of moderate order (ξsyn = 7 and ξsyn = 15, green and dark blue lines).


The identification of active cell-assemblies in simultaneously recorded spike trains has been the topic of extensive research (Gerstein et al. 1985, 1989; Sakurai 1998; Martignon et al. 2000; Grün et al. 2002a, b; Harris 2005; Schrader et al. 2008). Analysis tools that go beyond pairs of neurons, however, face severe limitations due to limited sample sizes of typical electrophysiological recordings (Martignon et al. 1995, 2000; Brown et al. 2004). The cumulant based inference of higher-order correlations (CuBIC) presented in this study avoids the need for extensive sample sizes. This is achieved by a combination of three ingredients: (1) pooling, i.e. using the population spike count and its cumulants to measure correlations in the population; (2) using the compound Poisson process (CPP) as a parametric model for correlated spiking populations; and (3) exploiting the constraining relations among correlations of lower and higher order in large neuronal pools.

Initially, we presented CuBIC as a hierarchy of hypothesis tests equation M180, indexed by the order of the estimated cumulants m and the maximal order of correlation ξ allowed in the null. The algorithm of Fig. 3 combines these tests to infer a lower bound equation M181 on the order of correlation in a given data set (see also Section 3.3). As free parameters, the algorithm requires upper bounds for m (mmax) and ξ (ξmax) to be chosen by the experimenter. To chose mmax, recall that CuBIC is extremely sensitive for higher-order correlations already with m = 3 (Fig. 6). For the parameters analyzed here, tests with m = 4 increased equation M182 only occasionally, as the smallest ξ that solved Eq. (9) typically yielded non-significant p-values (data not shown). Furthermore, the reliable estimation of high cumulants is known to require vast sample sizes. Therefore, although tests with m > 3 can in principle yield additional information, we do not expect that choosing the bound for m larger mmax = 4 improves test results in situations with sample sizes compatible with electrophysiological experiments. To chose the upper bound for ξ, note that ξmax is also an upper bound for equation M183. As the maximal order of correlation in the data has to be smaller than the number of recorded neurons N, we know that equation M184 must hold. This suggests to chose ξmax = N, as higher values of ξmax cannot increase equation M185. If the analysis was based on multi-unit activity and the number of recorded neurons is unknown, an a posteriori criterion on whether the chosen value for ξmax was high enough can be obtained by inspecting tests with ξ = ξmax. That is, if equation M186 was rejected for some m, hypothesis with ξ > ξmax might also have been rejected, potentially resulting in higher equation M187. Thus, if tests against equation M188 are retained for all m, the chosen ξmax was high enough.

The procedure was calibrated with rather restricted model classes, namely amplitude distributions with two isolated peaks. Thus, the question arises how differences in the correlation structure would influence CuBICs performance. In the numerical examples of this study, CuBIC used only the first three cumulants of the data, which depend on the first three moments of the amplitude distribution (compare Eq. (3)). As a consequence, if the data under consideration had a different amplitude distribution, however with similar first three moments, CuBIC will produce similar results. If, for instance, the additional isolated peak at ξsyn is exchanged by a somewhat distributed peak centered at ξsyn, the inferred lower bound will be similar.

Importantly, CuBIC is not designed to directly estimate the “true” maximal order of correlation in the data. Furthermore, the neuron IDs that realize the correlations are currently not resolved. It is important to keep in mind, however, that approaches that pose such specific question to the data face severe limitations with respect to the size of the analyzed populations (e.g., Shlens et al. 2006; see also Section 5.1). Thus, we regard the less specific lower bound equation M189 obtained by CuBIC as the price to pay when analyzing large populations. Renouncing the specific questions is precisely what allows CuBIC to detect higher-order correlations in data sets with previously unreached population size (N~100) and sample sizes compatible with physiological experiments (see Figs. 4 and and66).

Model dependence of higher-order correlations

The role of correlations for cortical information processing has been lively debated (Abeles 1991; Shadlen and Newsome 1994; Riehle et al. 1997; Singer 1999; Abbott and Dayan 1999; Kohn and Smith 2005; Schneidman et al. 2006; Shlens et al. 2006; Laurent 2006; Pillow et al. 2008). While in the meantime pairwise correlations are widely accepted as a relevant feature of parallel spike trains, the relevance of higher-order correlations is controversial. In fact, several recent studies raised doubts about the importance of higher-order interactions. In two studies that explore multi-neuron firing patterns in primate (Shlens et al. 2006) and salamander (Schneidman et al. 2006) retina, data were collected for fairly long times (~1 hour), and distributions of the binary patterns PN for sub-populations with N  10 were estimated. Based on this, a surrogate distribution was generated that reproduces the pairwise correlations in the population, but disregards all higher-order parameters: the “maximum entropy distribution” P2. A comparison of PN with P2 provides estimates on the amount of information contained in the higher-order parameters. The results of both studies indicate that only a negligible fraction of the total information is contained in the higher-order parameters (less than 2% in primate retina according to Shlens et al. 2006, and less than 10% in salamander retina according to Schneidman et al. 2006), and a recent study reports similar values for neo-cortical networks (Tang et al. 2008). In other words: almost the entire information available in the distribution of multi-neuron firing patterns can already be inferred from pairwise correlations. So why bother estimating higher-order parameters?

The consequences of the maximum entropy results for the relevance of CuBIC are difficult to judge. First of all, the maximum entropy approach does not appear to be very sensitive for higher-order correlations. As shown by Shlens et al. (2006), additional higher-order coincidences with a rate of ~0.3Hz among the N = 7 analyzed neurons would not have altered their results. This “coincidence rate”, however, lies in the parameter range of Set3 of Fig. 4 (here ν15 = 0.41Hz). In this range CuBIC detected correlations of order ≥ 13, illustrating the comparatively low sensitivity of the maximum entropy approach for existing higher-order correlations. Second, the combinatorial explosion of the parameter space limits maximum entropy approaches to only relatively small populations (equation M190), unless additional assumptions are imposed on the data (Montani et al. 2009; Shlens et al. 2009). Theoretical evidence suggests, however, that the relevance of higher-order correlations increases with the size of the neuronal population (Roudi et al. 2009). In fact, using N = 24 unsorted multi-unit spike trains of the rodent whisker system, Montani et al. (2009) report a significant increase in information if additional to the pairwise interactions also triplet interactions are taken into account. Taken together, we firmly believe that the analysis of large groups with methods that are sensitive to higher-order correlations is inevitable to uncover their role in cortical computation.

Our third remark concerns the conceptual differences between cumulant correlation as measured here, and the higher-order parameters used in the maximum entropy framework. The latter is based on the prominent binary exponential family (Martignon et al. 1995; Nakahara and Amari 2002; Martignon et al. 2000; Gütig et al. 2003; Del Prete et al. 2004; Shimazaki et al. 2009). In this setting, spike trains are represented as binary sequences, such that Xi(s) = 1 if the ith neuron has one or more spikes in the sth bin, and Xi(s) = 0 otherwise. Assuming stationarity and no memory, the multivariate distribution function of the binary pattern vector X admits the representation

equation M191

In this framework, higher-order correlations are reflected in nonzero “higher-order parameters” θ. For example, θijk quantifies third-order interactions of the triplet (Xi,Xj,Xk). Importantly, these higher-order parameters are neither equivalent, nor related in any straightforward way to the cumulant correlations employed by CuBIC (see Section 2.1.1 and Appendix A). A major difference, for instance, results from the equivalence of cumulant-based correlations with the concept of coincident events addressed by the CPP (Theorem 1). Given that these coincidences are perfectly synchronized, changes in the bin width h do not affect cumulant-based correlations. In contrast, as changes in the bin width alter the shape of pX, the interaction parameters θ inevitably change with h. Observe that the bin size dependence of θ does not result from the fact that larger bins measure correlations on a broader time scale as the time scale of correlation was assumed to be instantaneous, but is rather inherent to the exponential family. More generally, cumulant correlations determine to what extent expectation values of products factorize (e.g., Stratonovich 1967, and our Section 2.1.1), while the higher-order parameters in the exponential family measure to what extent the probability of certain binary patterns can be explained by the probabilities of its sub-patterns (e.g., Gütig et al. 2003). Taken together, the two frameworks are designed for different purposes, and can hence be expected to be sensitive to different features of the data. A better insight into their relationship presents an important path for future research (see also Staude et al. 2010, for a more detailed discussion).

Our main motivation to use cumulant correlations is that they enable us to analyze large populations of simultaneously recorded spike trains. As opposed to the θ parameters, cumulant correlations can be inferred from the population spike count Z,. Thereby, the multivariate problem to estimate correlations for all 2N  1 subgroups is transferred into a parsimoniously parametrized univariate problem with N parameters. The small sample size required by CuBIC and its ability to analyze populations of N > 100 neurons is a direct consequence of this property. Furthermore, the compound Poisson process presents a very intuitive interpretation of cumulant correlations (Theorem 1) that directly implements the concept of temporally coordinated spikes (see also Staude et al. 2008). And finally, the exponential family captures imprecise correlations only at the cost of ignoring spikes by clipping, as it requires the variables Xi to be binary (e.g., Del Prete et al. 2004), whereas cumulant correlations can be measured among arbitrary spike counts. We conclude that the adequate framework to analyze a given data set should be carefully chosen in view of the specific scientific question. For the detection of temporally coordinated spikes in large neuronal populations, we believe that cumulant-based correlations as exploited by CuBIC are well-suited.

Relation to other tools

While algorithms for the generation of correlated spike trains are becoming available (e.g., Kuhn et al. 2003; Niebur 2007; Staude et al. 2008; Brette 2009; Macke et al. 2008; Krumin and Shoham 2009), only few of these models can be exploited for data analysis. In an approach similar to ours, Ehm et al. (2007) derive estimators of the continuous time parameters (analogous to our carrier rate and amplitude distribution) from the Poissonized bin counts (here: “population spike count”) via Fourier inversion of the empirical characteristic function. Instead of inferring lower bounds as in CuBIC, “empirical de-Poissonization” presented by Ehm et al. (2007) thus directly estimates the orders of correlation in the data. However, the more explicit parameters appear to impose constraints on the data, such that particular care is required if the expected number of spikes per bin exceeds a critical value. A comparison of CuBIC and the empirical de-Poissonization to resolve their respective advantages and disadvantages is currently under way (Reimer et al. 2009).

Besides the higher-order approaches described above, the majority of available analysis tools for massively parallel spike trains either exploit pairwise statistics to infer inter-correlated groups by using some sort of clustering (Gerstein and Aertsen 1985; Gerstein et al. 1985; Berger et al. 2007; Grün et al. 2008b; Fujisawa et al. 2008), or test the occurrence of specific patterns against the assumption of complete (Grün et al. 2002a; Pipa et al. 2008) or partial independence (Berger et al. 2010). The three examples depicted in Fig. (4), however, illustrate the limitations of both approaches. First, the three data sets have identical second order statistics, making it impossible for pairwise approaches to differentiate between them. Second, the combinatorial explosion prevents the analysis of populations sizes of N~100 neurons with pattern-based approaches. And third, it is unclear whether either approach is sensitive enough to detect correlated groups when pairwise correlation are as low as discussed here.

A novel approach that is recently gaining attention starts from the analytically tractable generalized linear model (GLM), and estimates its parameters by maximum likelihood techniques (Paninski 2004; Truccolo et al. 2005; Okatan et al. 2005; Pillow et al. 2008). The power of this framework is that estimated parameters can be interpreted in terms of network properties like e.g. synaptic interactions between nerve cells, stimulus-response properties and the history dependence of individual spike trains. To understand how such anatomical and biophysical attributes of neuronal populations relate to higher-order correlations in the joint spiking of groups of neurons as detected by CuBIC is subject of current research.

Constraining relations among correlations of different order have also been studied by Amari et al. (2003), Johnson and Goodman (unpublished manuscript) and Benucci et al. (2004). The latter study defined a combinatorial algorithm to obtain upper bounds for certain observables. Closer to our approach is the study of Johnson and Goodman (unpublished manuscript), who derive explicit constraints among correlation coefficients of different orders within the CPP framework. However, both studies do not provide confidence intervals, hence are not applicable to data analysis in a straightforward way. Amari et al. (2003) investigated the relation between population variance and higher-order correlations in the information geometric framework (Amari and Nagaoka 2000) for the limit of large population sizes (N ∞). The result is similar to our findings in the sense that a high population variance implies higher-order correlations. However, results from large N provide only qualitative suggestions for data analysis.

Consequences of underlying model

An important ingredient of the present study is to use the CPP as a model for correlated spike trains. A consequence of this model-driven approach is that the reliability of obtained test results depends crucially on how well the CPP is suited to describe populations of spiking neurons in general, and assembly activity in particular.


We presented the CPP as a model of a stationary population, i.e. where individual spike trains have constant firing rates. However, cortical spike trains typically undergo both stimulus-driven and internally generated changes in firing rate. To analyze also such non-stationary data sets, the relatively small required sample size enables us to use a “sliding-window” approach (see e.g., Grün et al. 2002b), where data are cut into quasi-stationary windows of e.g. 100 ms and the statistics is therein obtained over trials. For the case of 100 trials analyzed with a bin width of h = 1 ms, this would yield L = 104 samples, which is perfectly in the required range for CuBIC to yield lower bounds close to the actual order of correlation (compare Fig. 6(d)). An alternative could be to directly include the non-stationarities in the model (Staude and Rotter 2009). Both approaches present important aspects of future research.

Single-process properties

Recent studies of cortical spike trains criticize the use of Poisson processes as models for single spike trains (Amarasingham et al. 2006; Nawrot et al. 2008). Without doubt, inclusion of process types deviating from Poisson in the generating model presents a major challenge for the presented procedure. However, criticism against the Poisson model for cortical spike trains is typically motivated by either non-exponential interval statistics (e.g., Maimon and Assad 2009), or non-Poissonian spike count distributions for large bin sizes (h > 50 ms) and high firing rates under ’optimal’ stimulation (λ > 50 Hz; e.g. Amarasingham et al. 2006). Both of these aspects are of little impact if the counting statistics is sampled with small (h  10 ms) bin sizes. In such a scenario, and in particular for low firing rates (e.g. Lennie 2003; Lee et al. 2006; Olshausen and Field 2006; Maldonado et al. 2008), neurons will have at most one spike per bin, leading to almost identical counting statistics for Poisson processes and other process types (e.g. gamma processes; van Vreeswijk 2006; Nawrot et al. 2008). Note also that the CPP does not directly model individual processes as Poisson processes, but only assumes the carrier processes to be Poisson. Importantly, the latter can be a good approximation even if individual processes deviate from strict Poisson statistics. For instance, if background spikes are Poisson and high-amplitude events are very regular but have a very low rate (as in Fig. 4), the carrier process is still well approximated by a Poisson process. More generally, the sum of increasingly sparse point processes converge to a Poisson process (Daley and Vere-Jones 2005, but see Lindner 2006; Câteau and Reyes 2006). Although these aspects need to be considered and will be subject of further research, we therefore expect CuBIC to yield reliable results even if single processes deviate from the Poisson assumption (Staude et al. 2007).

Precision of coincidences

Cortical correlations are reported on a variety of time scales, from precise correlations within a few few ms (Grün et al. 1999; Kohn and Smith 2005; Pazienti et al. 2008; Desbordes et al. 2008) to slow covariations of firing rates on a time scale of > 50ms (e.g., Kohn and Smith 2005; Smith and Kohn 2008). The perfect temporal precision of the joint spikes generated by the CPP may thus be questioned as a biologically plausible model of experimental spike trains. In the present study, the CPP was not intended to realistically model spike trains in continuous time, however. We presented the CPP strictly to interpret the correlated counting variables Xi, obtained with a given bin width h. Allowing a temporal jitter in the synchronous events, e.g. by jittering the coincident events during the inserting process (see e.g. Staude et al. 2008; Brette 2009), results in comparable counting variables Xi if this jitter is small compared to the applied bin size (Pazienti et al. 2008). Although its sensitivity for such imprecise coincidences remains to be assessed systematically, we expect CuBIC to tolerate imprecise coincidences with small but realistic jitter.

A somewhat different issue is the detection of correlations with systematic temporal offset, e.g. due to functional or “effective” connectivity (Aertsen et al. 1989; Fujisawa et al. 2008). As an extension of CuBIC to detect also non-zero-lag higher-order correlations, one may shift individual spike trains against each other and then analyze this shifted population (cmp. Grün et al. 1999). As tests for different shifts will generally not be independent, however, appropriate corrections that account for multiple testing have to be applied. Furthermore, the combinatorial explosion that emerges from shifting all spike train against each other with different delays requires smart algorithms to reduce the resulting computational cost.

Positivity of the correlations

A limitation of the CPP is that the insertion of synchronous spike events models spike coordination, not spike avoidance. As a consequence, correlations in the CPP can only be positive (Staude et al. 2008; Brette 2009). The CPP cannot model negatively correlated spike counts, that is, troughs in the cross-correlation function. Whether or not this is relevant for CuBICs applicability depends crucially on the properties of electrophysiological data. The rare reports of negative correlations, especially among cortical spike trains (e.g., Nevet et al. 2007), however, makes us confident that this is not an important limitation.


Based on Hebb’s original ideas (Hebb 1949), the coordinated spiking of large neuronal groups is the conceptual foundation of various “brain theories” (e.g., Eggermont 1990; Abeles 1991). However, limitations of available analysis tools has rendered the exploitation of the increasing number of simultaneously observable neurons extremely difficult. The CuBIC method presented in this study is a novel procedure to detect higher-order correlations in massively parallel spike trains, which overcomes these limitations. In particular, the sample size required by our method lies in a range that is compatible with electrophysiological experiments. CuBIC assumes data to be stationary and single processes the sum of which is appropriately modeled by the CPP, and its extension to account for dynamically changing, non-Poissonian data is the most important aspect of future work. Taken together, we regard CuBIC as a unique opportunity to analyze large populations of neurons for coordinated spike timing.


We would like to thank Abigail Morrison, Imke Reimer, and Stefano Cardanobile for helpful comments on earlier versions of this manuscript. Funded by NaFöG Berlin, the German Ministry for Education and Research (BMBF grants 01GQ01413 and 01GQ0420 to the BCCN Berlin and Freiburg, respectively), the IGPP Freiburg, the Stifterverband für die Deutsche Wissenschaft, the Initiative and Networking Fund of the Helmholtz Association within the Helmholtz Alliance on Systems Biology, and the German Research Foundation (DFG SFB 780).

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.


A Higher-order correlations and the CPP

A pair of random variables X1 and X2 is uncorrelated, if the expectation of their product factorizes

equation M192

and the covariance equation M193 measures the degree of correlation between them. As the covariance is bounded by the square root of the individual variances, one quantifies pairwise correlations by the dimensionless correlation coefficient

equation M194

Cumulant-correlations are a generalization of the (co-)variance to measure correlations among more than two r.v.s (e.g., Stratonovich 1967; Stuart and Ord 1987; Streitberg 1990; Papoulis and Pillai 2002; Gardiner 2003; Staude et al. 2010). The cumulants κm[Z] (m [set membership] N) of a r.v. Z are defined in terms of the coefficients of the power series expansion of the logarithm of its characteristic function

equation M195

where ΘZ(u) is called the “cumulant generating function” (c.g.f.) of Z. To relate the cumulants to the more familiar (raw) moments equation M196, recall that the latter are defined in terms of the coefficients of the power series expansion of the characteristic function

equation M197

Expressions relating cumulants to moments are obtained by writing the logarithm of the right hand side of Eq. (15) as a power series in u, and collecting the coefficients. For the first two cumulants, this yields

equation M198

Similar expressions for higher m are increasingly complicated, but algorithms for their computation are available (e.g. Stuart and Ord 1987; Di Nardo et al. 2008).

Given two independent r.v.s X1 and X2, we have

equation M199

Writing the c.g.f.s as power series and comparing coefficients yields the central property of cumulants: the cumulant of the sum of independent random variables is the sum of the individual cumulants

equation M200

If X1 and X2 are not independent, the cumulants of X1 + X2 contain also “mixed”, or “connected” cumulants. For m = 2, we obtain the well-known equation

equation M201
equation M202

Equation (16) implies that the second cumulant of X1 + X2 is additive-linear, if and only if equation M203, i.e. if X1 and X2 are uncorrelated. Higher-order cumulant correlations generalize the covariance in exactly this sense: they determine the “degree of linearity” of the higher cumulants (e.g., Stratonovich 1967; Gardiner 2003). Formally, higher-order cumulants arise as coefficients of the c.g.f. of vector-valued variables X = (X1,...,XN)

equation M204
equation M205

Higher-order cumulant correlations exist, if the corresponding mixed cumulant κJ[X] is non-zero. To be precise,

Definition 1 Let X be a random vector. Denote the number of non-zero elements of a multi-index equation M206 by nJ. Then the components of X are said to have correlations of order m if there exists J [set membership] NN with nJ = m such that the corresponding mixed cumulant is non-zero, i.e. κJ[X]  0.

For multi-indices J with only one non-zero element, e.g. J = (2,0,0...,0), the mixed cumulants κJ in Eq. (18) are simply the univariate cumulants of Eq. (14): equation M207. Furthermore, a generalization of equation M208 is that mixed cumulants with identical arguments are the univariate cumulants of the appropriate order

equation M209

Evidently, the c.g.f. of a sum of r.v.s equation M210 is the multivariate c.g.f. with all arguments set to u

equation M211

Inserting into Eq. (18) and comparing coefficients with Eq. (14) yields

equation M212

showing in particular that the mth cumulant of Z depends on correlations of maximal order m.

Proof of Theorem 1 To prove Theorem 1 of Section 2.4, we present a multivariate version of the compound Poisson process (CPP) defined in Section 2.2 (compare Ehm et al. 2007; Staude et al. 2010). Consider N neurons {xi(t)}i  N and let equation M213 be the set of all nonempty subsets of {1,...,N}. For each equation M214, we assign a stationary Poisson process yM(t) with rate νM that determines the points in time where the neurons {xi(t)}i [set membership] M have a joint spike event.Assuming all the yM(t)s to be independent, it is straightforward to verify that equation M215 is a CPP as defined in Section 2.2, where |M| denotes the number of elements of M. The processes yl(t) of Section 2.3 are given by yl(t) =  |M| = lyM(t), and we obtain the parameters of z(t) (carrier rate ν and amplitude distribution fA) by equation M216 and equation M217. Evidently, fA assigns non-zero probabilities to events of amplitude ≥ m if and only if at least one of the rates νM with |M|  m is positive. Using the above definition, Theorem 1 is thus equivalent to the equivalence of the following two statements:

  1. at least one mth order mixed cumulant is non-zero, i.e. there exists a multi-index J [set membership] NN with nJ  m such that κJ[X]  0
  2. at least one of the processes yM(t) that determines coincidences of at least m of the spike trains xi(t) has non-zero rate, i.e. there exists equation M218 with |M|  m and νM > 0

The proof is not difficult, yet notationally cumbersome for arbitrary N. The idea is to decompose Z into the counting variables YM of the yM′(t) and show that their mixed cumulants κJ[YM] vanish if J has more non-zero elements than the multi-index associated to M.With the variables YM, the argument of κJ[X], the pattern vector X = (X1,...,,XN), can be written as

equation M219

where equation M220 is the binary vector whose ith component is 1 if i [set membership] M and 0 otherwise. For given J [set membership] NN, we thus have

equation M221

where the second equality follows from the independence of the YMs. Clearly κJ[δM·YM] = 0 for all M with nJ > |M|, as in this case the order of the measured correlation nJ exceeds the number of non-zero components of the argument. In fact, κJ[δM·YM]  0 if and only if supp(J) [subset or is implied by] M, as these are the cases where we do not correlate with the variable 0 (the support supp(J) is the set of non-zero elements of J). In this case, κJ[δM·YM] is a cumulant of the Poisson variable YM which satisfies equation M222 (Eq. (20)). Hence,

equation M223

This implies

equation M224

Now we can prove the equivalence of statements (a) and (b).“(a)[implies](b)”: Let J +  [set membership] NN be such that equation M225 and equation M226. Then by Eq. (23), equation M227, and as all νM  0 there must be at least one equation M228 with supp(J) [subset or is implied by] M + and equation M229. As equation M230, this M + must have at least m elements, hence |M + |  m.“(b)[implies](a)”: Let equation M231 with |M + |  m and equation M232 and let J + be the associated multi-index. Then equation M233 and

equation M234

B Model parameters and population statistics

This appendix relates the parameters of the compound Poisson process, i.e. the carrier rate ν and the amplitude distribution fA, to parameters used in the context of physiological data, i.e. firing rates λ and pairwise count correlation coefficients c, for a population of N neurons.

Firing rates To compute the average rate λ, observe that equation M235. With Eq. (3), we thus have

equation M236

Pairwise Correlations To compute the count correlation coefficient c for a fixed bin width h (Eq. (13)), recall that by definition of the variance

equation M237

Denoting the population average variance and covariance by

equation M238

the population Fano-Factor equation M239 reads

equation M240

The Poisson property of the individual processes implies equation M241, which yields equation M242. Using the identity of the population Fano-Factor with the moments ratio equation M243, we finally obtain

equation M244

Note that by Eq. (25), the pairwise correlations, contained in r, do not depend on the details of the amplitude distribution fA, but are completely determined by the ratio of its first two moments ρ.If the population is homogeneous (all neurons have the same firing rate, all pairs are equally correlated, and so on), the average quantities v and r are the actual variance and covariance of the processes. In this case, the count correlation coefficient is given by equation M245, which implies

equation M246

Conversely, the moments ratio is given by

equation M247

Correlated subgroups Now consider the case where NC of the N neurons form a homogeneously correlated subgroup while the remaining NI = N  NC neurons are independent (see rasters in Fig. 4 for examples). In this case, we can decompose Z into the contribution of the correlated, ZC, and the uncorrelated, ZI, sub-populations: Z = ZC + ZI. Now equation M248 and equation M249. We can thus write

equation M250

Dividing by equation M251 yields

equation M252

For a homogeneous population as the correlated subgroup, the Fano Factor equation M253 is related to the pairwise correlation coefficient via Eq. (27), i.e. equation M254. Inserting into Eq. (28) and solving for c yields the correlation coefficient of the correlated subgroup

equation M255

Single interaction process The amplitude distributions used in Section 4 consist of two isolated peaks

equation M256

where δi,k is the Kronecker-Delta. The “single interaction process” (Kuhn et al. 2003) is the special case with ξin = N. Because fA(1) + fA(ξsyn) = 1, this model has three free parameters ν, fA(1) and ξsyn, which are related to N, λ and ρ (or c) by Eqs. (24) and (29). Straightforward computation yields that homogeneous populations with mean firing rate λ and correlation parameter ρ are achieved by setting

equation M257


equation M258

Given λ and ρ are constrained, the event amplitude ξsyn that realizes the correlation can be freely chosen as long as ν > 0 and η [set membership] [0,1]. This is guaranteed if

equation M259

Contributor Information

Benjamin Staude, Phone: +49-761-2039324, Fax: +49-761-2039559, ed.grubierf-inu.nccb@eduats.

Sonja Grün, pj.nekir.niarb@neurg.


  • Abbott LF, Dayan P. The effect of correlated variability on the accuracy of a population code. Neural Computation. 1999;11:91–101. doi: 10.1162/089976699300016827. [PubMed] [Cross Ref]
  • Abeles M. Role of cortical neuron: Integrator or coincidence detector? Israel Joural of Medical Sciences. 1982;18:83–92. [PubMed]
  • Abeles M. Corticonics: Neural circuits of the cerebral cortex. 1. Cambridge: Cambridge University Press; 1991.
  • Aertsen A, Gerstein G, Habib M, Palm G. Dynamics of neuronal firing correlation: Modulation of “effective connectivity” Journal of Neurophysiology. 1989;61(5):900–917. [PubMed]
  • Amarasingham A, Chen T-L, Harrison MT, Sheinberg DL. Spike count reliability and the poisson hypothesis. Journal of Neuroscience. 2006;26(4):801–809. doi: 10.1523/JNEUROSCI.2948-05.2006. [PubMed] [Cross Ref]
  • Amari S-i, Nagaoka H. Information geometry. New York: MS and Oxford University Press; 2000.
  • Amari S-i, Nakahara H, Wu S, Sakai Y. Synchronous firing and higher-order interactions in neuron pool. Neural Computation. 2003;15:127–142. doi: 10.1162/089976603321043720. [PubMed] [Cross Ref]
  • Benucci, A., Verschure, P. F. M. J., & König, P. (2004). High-order events in cortical networks: A lower bound. Physical Review E, 70(051909). [PubMed]
  • Berger, D., Borgelt, C., Louis, S., Morrison, A., & Grün, S. (2010). Efficient identification of assembly neurons within massively parallel spike trains. Computational Intelligence and Neuroscience, 21. doi:10.1155/2010/439648. [PMC free article] [PubMed]
  • Berger D, Warren D, Normann R, Arieli A, Grün S. Spatially organized spike correlation in cat visual cortex. Neurocomputing. 2007;70(10–12):2112–2116. doi: 10.1016/j.neucom.2006.10.141. [Cross Ref]
  • Bohte SM, Spekreijse H, Roelfsema PR. The effects of pair-wise and higher-order correlations on the firing rate of a postsynaptic neuron. Neural Computation. 2000;12:153–179. doi: 10.1162/089976600300015934. [PubMed] [Cross Ref]
  • Brette R. Generation of correlated spike trains. Neural Computation. 2009;21:188–215. doi: 10.1162/neco.2009.12-07-657. [PubMed] [Cross Ref]
  • Brown EN, Kaas RE, Mitra PP. Multiple neural spike train data analysis: State-of-the-art and future challenges. Nature Neuroscience. 2004;7(5):456–461. doi: 10.1038/nn1228. [PubMed] [Cross Ref]
  • Câteau H, Reyes A. Relation between single neuron and population spiking statistics and effects on network activity. Physical Review Letters. 2006;96:058101. doi: 10.1103/PhysRevLett.96.058101. [PubMed] [Cross Ref]
  • Daley DJ, Vere-Jones D. An introduction to the theory of point processes, Vol. 1: Elementary theory and methods. 2. New York: Springer; 2005.
  • Prete V, Martignon L, Villa AEP. Detection of syntonies between multiple spike trains using a coarse-grain binarization of spike count distributions. Network: Computation in Neural Systems. 2004;15:13–28. doi: 10.1088/0954-898X/15/1/002. [PubMed] [Cross Ref]
  • Desbordes G, Jin J, Weng CLNA, Stanley GB, Alonso J-M. Timing precision in population coding of natural scenes in the early visual system. PLoS Biology. 2008;6:e324. doi: 10.1371/journal.pbio.0060324. [PMC free article] [PubMed] [Cross Ref]
  • Nardo E, Guarino G, Senato D. A unifying framework for k-statistics, polykays and their multivariate generalizations. Bernoulli. 2008;14(2):440–468. doi: 10.3150/07-BEJ6163. [Cross Ref]
  • Diesmann M, Gewaltig M-O, Aertsen A. Stable propagation of synchronous spiking in cortical neural networks. Nature. 1999;402(6761):529–533. doi: 10.1038/990101. [PubMed] [Cross Ref]
  • Eggermont JJ. Studies of brain function (Vol. 16) Berlin: Springer; 1990. The correlative brain.
  • Ehm W, Staude B, Rotter S. Decomposition of neuronal assembly activity via empirical de-poissonization. Electronic Journal of Statistics. 2007;1:473–495. doi: 10.1214/07-EJS095. [Cross Ref]
  • Fujisawa S, Amarasingham A, Harrison M, Buzsaki G. Behavior-dependent short-term assembly dynamics in the medial prefrontal cortex. Nature Neuroscience. 2008;11:823–833. doi: 10.1038/nn.2134. [PMC free article] [PubMed] [Cross Ref]
  • Gardiner CW. Handbook of stochastic methods for physics, chemistry and the natural sciences. Springer Series in Synergetics. 3. New York: Springer; 2003.
  • Gerstein GL, Aertsen AMHJ. Representation of cooperative firing activity among simultaneously recorded neurons. Journal of Neurophysiology. 1985;54(6):1513–1528. [PubMed]
  • Gerstein GL, Bedenbaugh P, Aertsen A. Neuronal assemblies. IEEE Transactions on Biomedical Engineering. 1989;36:4–14. doi: 10.1109/10.16444. [PubMed] [Cross Ref]
  • Gerstein GL, Perkel DH, Dayhoff JE. Cooperative firing activity in simultaneously recorded populations of neurons: Detection and measurement. Journal of Neuroscience. 1985;5(4):881–889. [PubMed]
  • Grün, S., Abeles, M., & Diesmann, M. (2008a). Impact of higher-order correlations on coincidence distributions of massively parallel data. In Lecture notes in computer science. ‘Dynamic brain—from neural spikes to behaviors’ (Vol. 5286, pp. 96–114).
  • Grün S, Abeles M, Diesmann M. Lecture notes in computer science. ‘Dynamic brain—from neural spikes to behaviors’ (Vol. 5286) New York: Springer; 2008. Impact of higher-order correlations on coincidence distributions of massively parallel data.
  • Grün S, Diesmann M, Aertsen A. ‘Unitary events’ in multiple single-neuron spiking activity. Neural Computation. 2002;14(1):43–80. doi: 10.1162/089976602753284455. [PubMed] [Cross Ref]
  • Grün S, Diesmann M, Aertsen A. ‘Unitary Events’ in multiple single-neuron spiking activity. II. Non-stationary data. Neural Computation. 2002;14(1):81–119. doi: 10.1162/089976602753284464. [PubMed] [Cross Ref]
  • Grün S, Diesmann M, Grammont F, Riehle A, Aertsen A. Detecting unitary events without discretization of time. Journal of Neuroscience Methods. 1999;94(1):67–79. doi: 10.1016/S0165-0270(99)00126-0. [PubMed] [Cross Ref]
  • Gütig R, Aertsen A, Rotter S. Analysis of higher-order neuronal interactions based on conditional inference. Biological Cybernetics. 2003;88(5):352–359. doi: 10.1007/s00422-002-0388-0. [PubMed] [Cross Ref]
  • Harris K. Neural signatures of cell assembly organization. Nature Reviews. Neuroscience. 2005;5(6):339–407. [PubMed]
  • Hebb DO. The organization of behavior: A neuropsychological theory. New York: Wiley; 1949.
  • Holgate P. Estimation for the bivariate poisson distribution. Biometrika. 1964;51:241–245.
  • Johnson D, Goodman I. Inferring the capacity of the vector Poisson channel with a bernoulli model. Network: Computation in Neural Systems. 2008;19:13–33. doi: 10.1080/09548980701656798. [PubMed] [Cross Ref]
  • Kohn A, Smith MA. Stimulus dependence of neuronal correlations in primary visual cortex of the Macaque. Journal of Neuroscience. 2005;25(14):3661–3673. doi: 10.1523/JNEUROSCI.5106-04.2005. [PubMed] [Cross Ref]
  • König P, Engel AK, Roelfsema PR, Singer W. How precise is neuronal synchronization. Neural Computation. 1995;7:469–485. doi: 10.1162/neco.1995.7.3.469. [PubMed] [Cross Ref]
  • Kreiter AK, Singer W. Stimulus-dependent synchronization of neuronal responses in the visual cortex of awake macaque monkey. Journal of Neuroscience. 1996;16(7):2381–2396. [PubMed]
  • Kriener B, Tetzlaff T, Aertsen A, Diesmann M, Rotter S. Correlations and population dynamics in cortical networks. Neural Computation. 2008;20:2185–2226. doi: 10.1162/neco.2008.02-07-474. [PubMed] [Cross Ref]
  • Krumin M, Shoham S. Generation of spike trains with controlled auto-and cross-correlation functions. Neural Computation. 2009;21:1642–1664. doi: 10.1162/neco.2009.08-08-847. [PubMed] [Cross Ref]
  • Kuhn A, Aertsen A, Rotter S. Higher-order statistics of input ensembles and the response of simple model neurons. Neural Computation. 2003;1(15):67–101. doi: 10.1162/089976603321043702. [PubMed] [Cross Ref]
  • Kumar A, Rotter S, Aertsen A. Conditions for propagating synchronous spiking and asynchronous firing rates in a cortical network model. Journal of Neuroscience. 2008;28(20):5268–5280. doi: 10.1523/JNEUROSCI.2542-07.2008. [PubMed] [Cross Ref]
  • Laurent G. Shall we even understand the fly’s brain? In: Hemmen L, Sejnowski T, editors. 23 problems in systems neuroscience (Chap. 8) Oxford: Oxford University Press; 2006. p. 21.
  • Lee A, Manns I, Sakmann B, Brecht M. Whole-cell recordings in freely moving rats. Neuron. 2006;51:399–407. doi: 10.1016/j.neuron.2006.07.004. [PubMed] [Cross Ref]
  • Lennie P. The cost of cortical computation. Current Biology. 2003;13:493–497. doi: 10.1016/S0960-9822(03)00135-0. [PubMed] [Cross Ref]
  • Lindner B. Superposition of many independent spike trains is generally not a Poisson process. Physical Review E. 2006;73:022901. doi: 10.1103/PhysRevE.73.022901. [PubMed] [Cross Ref]
  • Louis, S., & Grün, S. (2009). Estimating the temporal precision and size of correlated groups of neurons from population activity. In BMC Neuroscience (Vol. 10, pp. P253).
  • Macke J, Berens P, Ecker AS, Tolias AS, Bethge M. Generating spike trains with specified correlation coefficients. Neural Computation. 2008;21(2):397–423. doi: 10.1162/neco.2008.02-08-713. [PubMed] [Cross Ref]
  • Maimon, G., & Assad, J. (2009). Beyond poisson: Increased spike-time regularity across primate parietal cortex. Neuron, 426–440. [PMC free article] [PubMed]
  • Maldonado P, Babul C, Singer W, Rodriguez E, Berger D, Grün S. Synchronization of neuronal responses in primary visual cortex of monkeys viewing natural images. Journal of Neurophysiology. 2008;100:1523–1532. doi: 10.1152/jn.00076.2008. [PubMed] [Cross Ref]
  • Martignon L, Deco G, Laskey K, Diamond M, Freiwald W, Vaadia E. Neural coding: Higher-order temporal patterns in the neurostatistics of cell assemblies. Neural Computation. 2000;12:2621–2653. doi: 10.1162/089976600300014872. [PubMed] [Cross Ref]
  • Martignon L, Hasseln H, Grün S, Aertsen A, Palm G. Detecting higher-order interactions among the spiking events in a group of neurons. Biological Cybernetics. 1995;73:69–81. doi: 10.1007/BF00199057. [PubMed] [Cross Ref]
  • Mattner L. What are cumulants? Documenta Mathematica. 1999;4:601–622.
  • Montani F, Ince RAA, Senatore R, Arabzadeh E, Diamond ME, Panzeri S. The impact of high-order interactions on the rate of synchronous discharge and information transmission in somatosensory cortex. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2009;367(1901):3297–3310. doi: 10.1098/rsta.2009.0082. [PubMed] [Cross Ref]
  • Nakahara H, Amari S. Information-geometric measure for neural spikes. Neural Computation. 2002;14:2269–2316. doi: 10.1162/08997660260293238. [PubMed] [Cross Ref]
  • Nawrot MP, Boucsein C, Rodriguez Molina V, Riehle A, Aertsen A, Rotter S. Measurement of variability dynamics in cortical spike trains. Journal of Neuroscience Methods. 2008;169:374–390. doi: 10.1016/j.jneumeth.2007.10.013. [PubMed] [Cross Ref]
  • Nevet A, Morris G, Saban G, Arkadir D, Bergman H. Lack of spike-count and spike-time correlations in the substantia nigra reticulata despite overlap of neural responses. Journal of Neurophysiology. 2007;98:2232–2243. doi: 10.1152/jn.00190.2007. [PubMed] [Cross Ref]
  • Niebur E. Generation of synthetic spike trains with defined pairwise correlations. NeuralComput. 2007;19:1720–1738. [PMC free article] [PubMed]
  • Okatan M, Wilson MA, Brown EN. Analyzing functional connectivity using a network likelihood model of ensemble neural spiking activity. Neural Computation. 2005;17(9):1927–1961. doi: 10.1162/0899766054322973. [PubMed] [Cross Ref]
  • Olshausen BA, Field DJ. What is the other 85 percent of V1 doing. In: Hemmen L, Sejnowski T, editors. 23 problems in systems neuroscience (Chap. 8) Oxford: Oxford University Press; 2006. pp. 182–211.
  • Palm, G. (1982). Neural assemblies. An alternative approach to artificial intelligence. Studies of Brain Function (Vol. 7). Berlin: Springer.
  • Paninski L. Maximum likelihood estimation of cascade point-process neural encoding models. Network: Computation in Neural Systems. 2004;15:243–262. doi: 10.1088/0954-898X/15/4/002. [PubMed] [Cross Ref]
  • Papoulis A, Pillai SU. Probability, random variables, and stochastic processes. 4. Boston: McGraw-Hill; 2002.
  • Pastalkova E, Itskov V, Amarasingham A, Buzsaki G. Internally generated cell assembly sequences in the rat hippocampus. Science. 2008;321:1322–1327. doi: 10.1126/science.1159775. [PMC free article] [PubMed] [Cross Ref]
  • Pazienti A, Maldonado PE, Diesmann M, Grün S. Effectiveness of systematic spike dithering depends on the precision of cortical synchronization. Brain Research. 2008;1225:39–46. doi: 10.1016/j.brainres.2008.04.073. [PubMed] [Cross Ref]
  • Pillow JW, Shlens J, Paninski L, Sher A, Litke AM, Chichilnisky EJ, et al. Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature. 2008;454:995–999. doi: 10.1038/nature07140. [PMC free article] [PubMed] [Cross Ref]
  • Pipa G, Wheeler D, Singer W, Nikolić D. Neuroxidence: A non-parametric test on excess or deficiency of joint-spike events. Journal of Computational Neuroscience. 2008;25:64–88. doi: 10.1007/s10827-007-0065-3. [PMC free article] [PubMed] [Cross Ref]
  • Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C. 2. Cambridge: Cambridge University Press; 1992.
  • Reimer, I. C. G., Staude, B., & Rotter, S. (2009). Detecting assembly activity in massively parallel spike trains. In H. Bähr, & I. Zerr (Eds.), Proceedings of the 8th meeting of the German neuroscience society/30th Göttingen neurobiology conference Neuroforum, Supplement (Vol. 1).
  • Riehle A, Grün S, Diesmann M, Aertsen A. Spike synchronization and rate modulation differentially involved in motor cortical function. Science. 1997;278(5345):1950–1953. doi: 10.1126/science.278.5345.1950. [PubMed] [Cross Ref]
  • Roudi Y, Nirenberg S, Latham PE. Pairwise maximum entropy models for studying large biological systems: When they can work and when they can’t. PLoS Computational Biology. 2009;5(5):e1000380+. doi: 10.1371/journal.pcbi.1000380. [PMC free article] [PubMed] [Cross Ref]
  • Sakurai Y. The search for cell assemblies in the working brain. Behavioural Brain Research. 1998;91:1–13. doi: 10.1016/S0166-4328(97)00106-X. [PubMed] [Cross Ref]
  • Sakurai Y, Takahashi S. Dynamic synchrony of firing in the monkey prefrontal cortex during working-memory tasks. Journal of Neuroscience. 2006;6(40):10141–10153. doi: 10.1523/JNEUROSCI.2423-06.2006. [PubMed] [Cross Ref]
  • Schneider G, Grün S. Analysis of higher-order correlations in multiple parallel processes. Neurocomputing. 2003;52–54:771–777. doi: 10.1016/S0925-2312(02)00772-5. [Cross Ref]
  • Schneidman E, Berry MJ, Segev R, Bialek W. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature. 2006;440:1007–1012. doi: 10.1038/nature04701. [PMC free article] [PubMed] [Cross Ref]
  • Schrader S, Grün S, Diesmann M, Gerstein G. Detecting synfire chain activity using massively parallel spike train recording. Journal of Neurophysiology. 2008;100:2165–2176. doi: 10.1152/jn.01245.2007. [PubMed] [Cross Ref]
  • Shadlen MN, Movshon AJ. Synchrony unbound: A critical evaluation of the temporal binding hypothesis. Neuron. 1999;24:67–77. doi: 10.1016/S0896-6273(00)80822-3. [PubMed] [Cross Ref]
  • Shadlen MN, Newsome WT. Noise, neural codes and cortical organization. Current Opinion in Neurobiology. 1994;4(4):569–579. doi: 10.1016/0959-4388(94)90059-0. [PubMed] [Cross Ref]
  • Shimazaki, H., Amari, S., Brown, E. N., & Grün, S. (2009). State-space analysis on time-varying correlations in parallel spike sequences. Proc. IEEE international conference on acoustics, speech, and signal processing (ICASSP), 3501–3504.
  • Shlens J, Field GD, Gauthier JL, Greschner M, Sher A, Litke AM, et al. The structure of large-scale synchronized firing in primate retina. Journal of Neuroscience. 2009;29(15):5022–5031. doi: 10.1523/JNEUROSCI.5187-08.2009. [PMC free article] [PubMed] [Cross Ref]
  • Shlens J, Field GD, Gauthier JL, Grivich MI, Petrusca D, Sher A, Litke AM, Chichilnisky E. The structure of multi-neuron firing patterns in primate retina. Journal of Neuroscience. 2006;26(32):8254–8266. doi: 10.1523/JNEUROSCI.1282-06.2006. [PubMed] [Cross Ref]
  • Singer W. Neuronal synchrony: A versatile code for the definition of relations? Neuron. 1999;24(1):49–65. doi: 10.1016/S0896-6273(00)80821-1. [PubMed] [Cross Ref]
  • Singer W, Engel AK, Kreiter AK, Munk MHJ, Neuenschwander S, Roelfsema PR. Neuronal assemblies: Necessity, signature and detectability. Trends in Cognitive Sciences. 1997;1(7):252–261. doi: 10.1016/S1364-6613(97)01079-6. [PubMed] [Cross Ref]
  • Singer W, Gray C. Visual feature integration and the temporal correlation hypothesis. Annual Review of Neuroscience. 1995;18:555–586. doi: 10.1146/ [PubMed] [Cross Ref]
  • Smith MA, Kohn A. Spatial and temporal scales of neuronal correlation in primary visual cortex. Journal of Neuroscience. 2008;28(48):12591–12603. doi: 10.1523/JNEUROSCI.2929-08.2008. [PMC free article] [PubMed] [Cross Ref]
  • Snyder DL, Miller MI. Random point processes in time and space. New York: Springer; 1991.
  • Staude B, Grün S, Rotter S. Higher order correlations. In: Grün S, Rotter S, editors. Analysis of parallel spike trains. New York: Springer; 2010. [PMC free article] [PubMed]
  • Staude, B., & Rotter, S. (2009). Higher-order correlations in non-stationary parallel spike trains: Statistical modeling and inference. In BMC Neuroscience (Vol. 10, pp. P108).
  • Staude, B., Rotter, S., & Grün, S. (2007). Detecting the existence of higher-order correlations in multiple single-unit spike trains. In Society for neuroscienceabstract viewer/itinerary planner (Vol. 103.9/AAA18). Washington, DC.
  • Staude B, Rotter S, Grün S. Can spike coordination be differentiated from rate covariation? Neural Computation. 2008;20:1973–1999. doi: 10.1162/neco.2008.06-07-550. [PubMed] [Cross Ref]
  • Stratonovich RL. Topics in the theory of random noise. New York: Gordon & Breach Science; 1967.
  • Streitberg B. Lancaster interactions revisited. The Annals of Statistics. 1990;18(4):1878–1885. doi: 10.1214/aos/1176347885. [Cross Ref]
  • Stuart A, Ord JK. Kendall’s advanced theory of statistics. 5. London: Griffin; 1987.
  • Tang A, Jackson D, Hobbs J, Chen W, Smith JL, Patel H, et al. A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. Journal of Neuroscience. 2008;28:505–518. doi: 10.1523/JNEUROSCI.3359-07.2008. [PubMed] [Cross Ref]
  • Truccolo W, Eden UT, Fellows MR, Donoghue JP, Brown EN. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. Journal of Neurophysiology. 2005;93:1074–1089. doi: 10.1152/jn.00697.2004. [PubMed] [Cross Ref]
  • Vaadia E, Aertsen A, Nelken I. ’dynamics of neuronal interactions’ cannot be explained by ’neuronal transients’ Proc Biol Sci. 1995;261(1362):407–410. doi: 10.1098/rspb.1995.0167. [PubMed] [Cross Ref]
  • Vreeswijk C. What is the neural code? In: Hemmen L, Sejnowski T, editors. 23 problems in systems neuroscience (Chap. 8) Oxford: Oxford University Press; 2006. pp. 143–159.
  • Wennekers T, Garagnani M, Pulvermüller F. Language models based on hebbian cell assemblies. Journal of Physiology (Paris) 2006;100(1–3):16–30. doi: 10.1016/j.jphysparis.2006.09.007. [PubMed] [Cross Ref]
  • Wennekers T, Sommer F, Aertsen A. Cell assemblies - editorial. Theory in Biosciences. 2003;122:1–4.

Articles from Springer Open Choice are provided here courtesy of Springer