Home | About | Journals | Submit | Contact Us | Français |

**|**Front Comput Neurosci**|**v.4; 2010**|**PMC2906200

Formats

Article sections

- Abstract
- Introduction
- The Stationary Case
- The Non-Stationary Case
- Case Studies
- Discussion
- Conflict of Interest Statement
- References

Authors

Related links

Front Comput Neurosci. 2010; 4: 16.

PMCID: PMC2906200

Edited by: Jakob H. Macke, Max Planck Institute for Biological Cybernetics, Germany

Reviewed by: Yasser Roudi, NORDITA, Sweden; Don H. Johnson, Rice University, USA; Jonathan D. Victor, Weill Cornell Medical College, USA

*Correspondence: Stefan Rotter, Bernstein Center Freiburg and Faculty of Biology, Albert-Ludwig University, Hansastrasse 9a, 79104 Freiburg, Germany. e-mail: ed.grubierf-inu.eigoloib@rettor.nafets

Received 2009 December 2; Accepted 2010 May 11.

Copyright © 2010 Staude, Grün and Rotter.

This is an open-access article subject to an exclusive license agreement between the authors and the Frontiers Research Foundation, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are credited.

This article has been cited by other articles in PMC.

The extent to which groups of neurons exhibit higher-order correlations in their spiking activity is a controversial issue in current brain research. A major difficulty is that currently available tools for the analysis of massively parallel spike trains (*N*>10) for higher-order correlations typically require vast sample sizes. While multiple single-cell recordings become increasingly available, experimental approaches to investigate the role of higher-order correlations suffer from the limitations of available analysis techniques. We have recently presented a novel method for cumulant-based inference of higher-order correlations (CuBIC) that detects correlations of higher order even from relatively short data stretches of length *T*=10–100s. CuBIC employs the compound Poisson process (CPP) as a statistical model for the population spike counts, and assumes spike trains to be stationary in the analyzed data stretch. In the present study, we describe a non-stationary version of the CPP by decoupling the correlation structure from the spiking intensity of the population. This allows us to adapt CuBIC to time-varying firing rates. Numerical simulations reveal that the adaptation corrects for false positive inference of correlations in data with pure rate co-variation, while allowing for temporal variations of the firing rates has a surprisingly small effect on CuBICs sensitivity for correlations.

It has long been suggested that fundamental insight into the nature of neuronal computation requires the understanding of the cooperative dynamics of populations of neurons (Hebb, 1949). A controversial issue in this debate is the role of correlations among nerve cells. On the one hand, an increasing body of both experimental (e.g., Gray and Singer, 1989; Vaadia et al., 1995; Riehle et al., 1997; Bair et al., 2001; Kohn and Smith, 2005; Shlens et al., 2006; Fujisawa et al., 2008; Pillow et al., 2008) and theoretical (Abeles, 1991; Diesmann et al., 1999; Kuhn et al., 2003) literature supports the concept of cooperative computation on various temporal and spatial scales. On the other hand, the mostly detrimental effect of correlations on rate-based information transmission and processing (Abbott and Dayan, 1999; Averbeck and Lee, 2006; Josić et al., 2009) has generated a strong opposition toward correlation-based concepts of cortical coding (Shadlen and Newsome, 1998; Averbeck et al., 2006; Schneidman et al., 2006; Ecker et al., 2010). Evidently, a thorough description of the correlation structure of neuronal populations is an indispensable prerequisite to resolve these opposing theoretical viewpoints (Brown et al., 2004).

Experimental reports on coordinated activity at the level of spike trains resort almost exclusively to correlations between pairs of nerve cells (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; Ecker et al., 2010). Such pairwise correlations cannot, as a matter of principle, resolve the cooperative activity of neuronal populations to the extent required for rigorous hypothesis testing (Gerstein et al., 1989; Martignon et al., 2000; Brown et al., 2004). In particular, whether or not coincident spikes of pairs of neurons participate in synchronized “cluster-events” cannot be decided on measurements of pairwise correlation alone; this can only be achieved by the systematic assessment of higher-order correlations, i.e., statistical couplings among triplets, quadruplets, and larger groups (Martignon et al., 1995; Staude et al., 2010). Importantly, the nonlinear dynamics of spike generation makes neurons extremely sensitive for synchrony in their input pools (Softky, 1995; König et al., 1996). Ignoring these higher-order correlations in the statistical description of spiking populations is therefore hardly advisable (Bohte et al., 2000; Kuhn et al., 2003).

Initially, the main obstacle for assessing the higher-order structure of neuronal populations were limitations in experimental methodology, as until recently state-of-the-art electrophysiological setups allowed to record only few neurons simultaneously. The advent of multi-electrode arrays and optical imaging techniques, however, now reveals fundamental shortcomings of available analysis tools (Brown et al., 2004). Mathematical frameworks to model and estimate higher-order correlations typically assign one “interaction parameter” for every subgroup of the population, leading to a 2* ^{N}*−1 dimensional model for a population comprising

We have recently presented a novel method for a cumulant-based inference for the presence of higher-order correlations (CuBIC) that avoids the need for extensive sample sizes (Staude et al., 2007, 2009). Instead of directly estimating correlation parameters from all subgroups, CuBIC aims only at population-average correlations, estimated via the cumulants of the pooled and discretely sampled spiking activity of all recorded neurons (population spike counts). The presence of higher-order correlations is then inferred from measured cumulants of low order by exploiting certain constraining relations among correlations of different orders in a statistical model of correlated spiking. CuBIC avoids the direct estimation of higher-order correlations, but decides whether or not lower order cumulants require the presence of higher-order correlations. Focusing on such less specific questions drastically reduces the requirements with respect to the sample size: when applied to artificial data, CuBIC reliably infers higher-order correlations from large (*N*100), even weakly correlated populations (pairwise correlation coefficient *c* ~ 0.01) that were generated with reasonable sample sizes (*T*<100s, Staude et al., 2009).

As a statistical model, CuBIC employs the compound Poisson process (CPP), where correlations are induced by the insertion of coincident events in continuous time, i.e., before binning is applied (Ehm et al., 2007; Johnson and Goodman, 2007; Brette, 2009; Staude et al., 2010). Interestingly, this model of correlation fits perfectly to measuring (higher-order) correlations via connected cumulants of the binned spike trains (Staude et al., 2010), a common framework for (higher-order) correlation measures. The simple relationship of the unknown model parameters, i.e., the orders of correlation present in the data, and the observable cumulants of the population spike count allows to devise null-hypotheses concerning the orders of correlation in the data (the details of the CPP and CuBIC are explained in Section “The Stationary Case”). Combining tests against different null-hypothesis yields a lower bound $\widehat{\xi}$ for the maximal order of correlation in the data.

A central assumption in the original presentation of CuBIC (Staude et al., 2009) was that the statistics of spiking in the population does not change over time (stationarity). As both experimental cues and/or internal processes often induce transients or fluctuations of firing rates, this central assumption is frequently violated in electrophysiological data.

In the present study, we describe a non-stationary version of the CPP by decoupling the correlation structure from the spike intensity of the population (see Section “The Non-stationary Case”). Using the “law of total cumulance” we are able to incorporate non-stationarities in firing rate into the computation of the cumulants of the population spike counts. These rate-adjusted cumulants are then used to adapt CuBIC to infer higher-order correlations also from non-stationary data. This adaptation requires a specification of the kind of non-stationarity in terms of a parametric family of distributions for the bin-wise mean firing rates (the “carrier distribution”). Allowing for uniform rate fluctuations, for instance, yields as a result that the data must have correlations of some minimal order $\widehat{\xi}$ even if firing rates fluctuated uniformly from bin to bin. In this sense, the choice of a family for the carrier distribution implies a demarcation line between “genuine” correlation and “artifacts” due to rate (co-)variation (Staude et al., 2008). Numerical simulations reveal that the adaptation corrects for false positive inference of correlations in data with pure rate co-variation, while allowing for potential variations in firing rates has a surprisingly small effect on CuBIC's sensitivity for correlations (see Case Studies). Furthermore, we find that a perfect match between the true carrier family and the family allowed in the adapted CuBIC does not seem to be fundamentally important to guarantee reliable test performance.

The basic observable of this study is the pattern vector **X**(*s*):=(*X*_{1}(*s*),…,*X _{N}*(

$$Z(s)={\displaystyle \sum _{i=1}^{N}{X}_{i}(s).}$$

In the case where the *X _{i}* are binary (“1” for one or more spikes in the bin, “0” for no spike),

We here assume that *Z*(*s*) and *Z*(*s*+*k*) are independent for *k*≠0 (zero memory). Furthermore, let us for now assume that the distribution of *Z*(*s*) does not depend on the time bin *s* (stationarity). This critical assumption will be relaxed in Section “The Non-stationary Case”.

In the present framework, correlations among the variables *X _{i}* are measured by mixed or “connected” cumulants. Like the more familiar (raw) moments E[

Multivariate, or “connected”, cumulants arise when the variable under consideration is a sum of correlated variables. For *m*=2 and $Z={\Sigma}_{i=1}^{N}{X}_{i}$, for instance, we have the well-known formula

$$\begin{array}{c}{\text{\kappa}}_{\text{2}}\left[Z\right]=\text{Var}\left[Z\right]=\text{Var}\left[{\displaystyle \sum _{i=1}^{N}{X}_{i}}\right]={\displaystyle \sum _{i=1}^{N}\text{Var}\left[{X}_{i}\right]}+{\displaystyle \sum _{i\ne j}\text{Cov}\left[{X}_{i},{X}_{j}\right]}\\ =:{\displaystyle \sum _{i=1}^{N}{\kappa}_{2}\left[{X}_{i}\right]+}{\displaystyle \sum _{i\ne j}{\kappa}_{1,1}\left[{X}_{i},{X}_{j}\right]}\end{array}$$

(1)

Hence, the second-order cumulant correlations Cov[*X _{i}*,

**Definition 1** *Let* **X**=(*X*_{1},…,*X _{N}*)

The following generalization of Eq. 1 is a straightforward consequence of the construction of connected cumulants (Staude et al., 2010).

**Theorem 1** *The mth cumulant* κ* _{m}*[

By the above theorem, κ* _{m}*[

Before a discussion of our model, we wish to stress that the cumulant correlations presented here do not comply with the interaction parameters of the more familiar log-linear model. In particular, data sets can have higher-order log-linear interactions without having higher-order cumulant correlations and vice versa (see Staude et al., 2010 for concrete examples, and e.g., Darroch and Speed, 1983; Streitberg, 1990; Staude et al., 2009 for more general discussions).

As opposed to the *discretized*, 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 $z(t)={\Sigma}_{i=1}^{N}{x}_{i}(t)$, where ${x}_{i}(t)={\Sigma}_{j}\delta (t-{t}_{j}^{i})$ denotes the *i*th unbinned, continuous-time spike train with spike event times ${t}_{j}^{i}\left(i=1,\dots ,N,j\in \mathbb{N}\right)$. The model we propose for *z*(*t*) is that of a CPP

$$z\left(t\right)={\displaystyle \sum _{j}\text{\delta}\left(t-{t}_{j}\right){a}_{j}},$$

(2)

where the event times *t _{j}* constitute a Poisson process, and the marks

With the above model, the generation of a population of spike trains proceeds in two steps. First, realize a Poissonian carrier processes *m*(*t*)= ∑* _{j}*δ(

**Theorem 2** *Let* $z(t)={\Sigma}_{i=1}^{N}{x}_{i}(t)$ *be a CPP with amplitude distribution f _{A} and carrier rate*ν,

*The components of***X***have correlations of order m**(in the sense of Definition 1)**if and only if f*_{A}*assigns non-zero probabilities to amplitudes*≥*m.**With*${\mu}_{m}:=\text{E}[{A}^{m}]={\Sigma}_{k=1}^{N}{k}^{m}{f}_{A}(\xi )$,*the cumulants of Z**are given as*

$${\kappa}_{m}[Z]={\mu}_{m}\nu h.$$

(3)

Note that correlations in the above theorem are measured strictly on the basis of the discretized counting variables *X _{i}*. As a consequence, they do not resolve (and do not depend on) the perfect temporal precision of the coincident events in the CPP. That is, if the events of

Now let ξ≤ *N* be the maximal order of correlation in the model, i.e., *f _{A}*(

$${\kappa}_{m}[Z]={\overrightarrow{\xi}}_{m}\cdot {\overrightarrow{\nu}}_{\xi}h,$$

(4)

where ${\overrightarrow{\xi}}_{m}=({1}^{m},\dots ,{\xi}^{m})$, ${\overrightarrow{\nu}}_{\xi}:=({\nu}_{1},\dots ,{\nu}_{\xi})$, and ${\overrightarrow{\xi}}_{m}\cdot {\overrightarrow{\nu}}_{\xi}:={\Sigma}_{i=1}^{\xi}{i}^{m}{\nu}_{i}$ is the vector dot product.

This section summarizes the stationary version of CuBIC to the extent that is needed to understand its adaptation to non-stationary populations (see Staude et al., 2009 for details). In brief, CuBIC quantifies the following thought experiment. Consider the situation of four 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 must in fact be identical. In the framework of the CPP, this translates to the existence of events of amplitude *a _{j}*=4, and hence correlation of order 4 (Theorem 1). This illustrates that it is possible, in principle, to infer the existence of fourth order correlations from estimated pairwise correlations. In Staude et al. (2009), this inference was generalized by a hierarchy of statistical hypothesis tests ${H}_{0}^{m,\xi}$, labeled by the order

Assume one is given the first three cumulants of a population spike count variable *Z*′. Then, for a fixed value of the test parameter ξ, consider the following constrained maximization problem

$${\kappa}_{3,\xi}^{\ast}:={\text{max}}_{\nu ,{f}_{A}}\left\{{\kappa}_{3}[Z]\right\}$$

(5)

subject to κ_{2}[*Z*′]= κ_{2}[*Z*]

κ_{1}[*Z*′]= κ_{1}[*Z*]

*f _{A}*(

where *Z* is the population spike count of a model with parameters ν and *f _{A}*. The model that solves Eq. 5 has the maximal third cumulant, i.e., triplet correlations, among all models that do not have any correlations beyond order ξ, and have the same population-averaged first- and second-order properties, i.e., firing rates and pairwise correlations, as the given spike count variable

To solve Eq. 5, we use Eq. 4 and obtain the equivalent problem

$${\kappa}_{3,\xi}^{*}=\underset{{\overrightarrow{\nu}}_{\xi}}{\text{max}}\left\{{\overrightarrow{\xi}}_{3}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h\right\}$$

(6)

$$\begin{array}{c}\text{subject}\hspace{0.17em}\text{to}\hspace{0.17em}{\kappa}_{2}\left[Z\prime \right]={\overrightarrow{\xi}}_{2}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h\\ {\kappa}_{1}\left[Z\prime \right]={\overrightarrow{\xi}}_{1}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h.\end{array}$$

In Eq. 6, the objective function and the constraints depend linearly on the model parameters ${\overrightarrow{\nu}}_{\xi}$. Problems of this type, so-called Linear Programming Problems, are uniquely solvable, e.g., by the Simplex Method or any of its variants (Press et al., 1992, Chapter 10.8). The solution yields an upper bound for the third cumulant ${\kappa}_{3,\xi}^{\ast}$ and the corresponding parameter vector ${\overrightarrow{\nu}}_{\xi}^{\ast}$. As it turns out, the only non-zero components of the solver ${\overrightarrow{\nu}}_{\xi}^{\ast}=({\nu}_{1}^{\ast},\dots ,{\nu}_{\xi}^{\ast})$ are ${\nu}_{1}^{\ast}$ and ${\nu}_{\xi}^{\ast}$, i.e., ${\nu}_{k}^{\ast}=0$ for all *k*{1,ξ} (see The Solution of the Maximization in Appendix). The carrier rate and amplitude distribution of the CPP that maximizes Eq. 6 are then given by ${\nu}^{\ast}={\nu}_{1}^{\ast}+{\nu}_{\xi}^{\ast}$ and ${f}_{{A}^{*}}(l)={\nu}_{l}^{\ast}/{\nu}^{\ast}$.

With the solution of Eq. 5, the null hypothesis ${H}_{0}^{3,\xi}$ is

$${H}_{0}^{3,\xi}:{\kappa}_{3}[Z\prime ]\le {\kappa}_{3,\xi}^{\ast}$$

To test a sample $\{Z\prime \}=\{{Z}_{1}^{\prime},\dots ,{Z}_{L}^{\prime}\}$ against ${H}_{0}^{3,\xi}$, we estimate its cumulants by the so-called *k*-statistics *k _{m}* (Stuart and Ord, 1987; the well-known sample mean and unbiased sample variance are the first two

$$\text{Var}[{k}_{3}]=\frac{{\kappa}_{6}[Z\prime ]}{L}+9\frac{{\kappa}_{2}[Z\prime ]{\kappa}_{4}[Z\prime ]}{L-1}+6\frac{{\kappa}_{2}{[Z\prime ]}^{3}}{(L-1)(L-2)},$$

(7)

where κ* _{i}*[

$${p}_{3,\text{\xi}}={\displaystyle {\int}_{{k}_{m}}^{\infty}\frac{1}{\sqrt{2\text{\pi Var}\left[{k}_{3}\right]}}}\mathrm{exp}\left(-\frac{{\left(t-{\text{\kappa}}_{\text{3,\xi}}^{\text{*}}\right)}^{2}}{2\text{Var}\left[{k}_{3}\right]}\right)dt.$$

(8)

As mentioned above, the rejection of a specific hypothesis ${H}_{0}^{3,\xi}$ implies that the data have correlations beyond order ξ or, in other words, that ξ +1 is a lower bound for the order of correlation. The final result of CuBIC is the maximum of these lower bounds

$$\widehat{\xi}:=\text{max}\{\xi |{p}_{3,\xi}<\alpha \}+1,$$

(9)

where α is a predefined test level (see Staude et al., 2009 for details).

In the previous section, the CPP was presented as a model for populations with constant firing rates. However, (co-)variations in firing rate are a common feature of neuronal populations. To incorporate potential non-stationarities into the CPP, recall its central ingredients: the intensity of the population is described by the carrier rate ν, the population-averaged correlation structure is determined by the amplitude distribution *f _{A}*, and the precise composition of spikes and correlations within the population is determined by the assignment distribution. Given this parametrization, non-stationarities can in principle be included into the CPP by allowing a time-dependent carrier rate, amplitude distribution, assignment distribution, or combinations thereof.

Rather than presenting a general model for non-stationary populations, the focus of this study is to adapt the analysis technique CuBIC for potential variations in the firing rates. CuBIC, however, aims only at the population-average correlation structure *f _{A}*, and inference is based on the population spike count

To relate the model parameters to the cumulants of *Z* for time varying ν(*t*), observe that the value of *Z*(*s*) in the window [*sh*,(*s*+1)*h*) does not depend on the precise time course of the carrier rate ν(*t*) in this window, but only on the integral ${R}_{s}:=1/h{\int}_{sh}^{(s+1)h}\nu (t)dt.$ Substituting the carrier rate ν(*t*) with a piecewise constant function whose value in the interval [*sh*,(*s*+1)*h*) is *R _{s}* thus results in an identical population spike count

The above setting characterizes the population spike count *Z* as a parameter-dependent random variable, where the outcome of *Z* in the *s*th bin depends not only on the outcome of the CPP realization, but also on the (random) value of the rate variable *R _{s}*. For such “doubly stochastic” variables, the raw moments are given as

$${\text{\mu}}_{m}\left[Z\right]:=\text{E}\left[{Z}^{m}\right]={\text{E}}^{R}\left[\text{E}\left[{Z}^{m}|R\right]\right]\left(m\in \mathbb{N}\right),$$

where the inner expectation is the expectation value of *Z ^{m}* for a given value of the rate

$${\varphi}_{Z}(s):=E\left[{e}^{isZ}\right]$$

(10)

$$={\displaystyle \sum _{m}{i}^{m}\frac{{s}^{m}}{m!}{\text{\mu}}_{m}\left[Z\right]},$$

(11)

such that the moments can be obtained from ϕ* _{Z}* via

$${\mu}_{m}[Z]={\frac{1}{{i}^{m}}\frac{{\partial}^{m}{\varphi}_{Z}(s)}{\partial {s}^{m}}|}_{s=0}.$$

Analogously, cumulants are the coefficients of the logarithm of the characteristic function

$${\psi}_{Z}(s)=\text{log}E[{e}^{isZ}]$$

(12)

$$={\displaystyle \sum _{m}{i}^{m}\frac{{s}^{m}}{m!}{\text{\kappa}}_{m}\left[Z\right]},$$

(13)

such that

$${\kappa}_{m}[Z]=\frac{1}{{i}^{m}}{\frac{{\partial}^{m}{\psi}_{Z}(s)}{\partial {s}^{m}}|}_{s=0}$$

(14)

A straightforward strategy to relate cumulants to moments is to insert Eq. 11 into the right hand side of Eq. 12, writing the logarithm as a power series, and collecting coefficients for identical powers of *s*. This procedure illustrates in particular that the *m*th cumulant is a function of the first *m* raw moments only, and, reversely, that the *m*th moment can be expressed as a function of the first *m* cumulants. We denote the maps relating cumulants to moments and moments to cumulants by *F _{m}* and

$${\mu}_{m}[Z]={F}_{m}\left({\kappa}_{1}[Z],\dots ,{\kappa}_{m}[Z]\right)$$

(15)

$${\kappa}_{m}[Z]={G}_{m}\left({\mu}_{1}[Z],\dots ,{\mu}_{m}[Z]\right).$$

(16)

With these maps at hand, the *m*th cumulant of the non- stationary CPP can be computed by the following procedure.

- For
*i*=1,…,*m*express the conditional moment μ_{i}[*Z*|*R*] in terms of the first*i*cumulants, i.e., write μ_{i}[*Z*|*R*] =*F*_{i}(κ_{1}[*Z*|*R*],…,κ[_{i}*Z*|*R*]). - Apply Eq. 3 to the individual cumulants, such that μ
[_{i}*Z*|*R*] =*F*_{i}(μ_{1}*Rh*,…,μ), where, as before, μ_{i}Rh=μ_{i}[_{i}*A*] are the moments of the amplitude distribution*f*._{A} - Compute the
*m*th cumulant of*Z*by applying*G*_{m}

$${\kappa}_{m}[Z]={G}_{m}\left({F}_{1}({\mu}_{1}Rh),\dots ,{F}_{m}({\mu}_{1}Rh,\dots ,{\mu}_{m}Rh)\right).$$

(17)

The results are summarized as the “law of total cumulance”. The first three orders read

$${\kappa}_{1}[Z]={\mu}_{1}{\kappa}_{1}[R]h$$

(18)

$${\kappa}_{2}[Z]={\mu}_{2}{\kappa}_{1}[R]h+{\mu}_{1}^{2}{\kappa}_{2}[R]{h}^{2}$$

(19)

$${\kappa}_{3}[Z]={\mu}_{3}{\kappa}_{1}[R]h+{\mu}_{1}^{3}{\kappa}_{3}[R]{h}^{3}+3{\mu}_{1}{\mu}_{2}{\kappa}_{2}[R]{h}^{2}.$$

(20)

A similar parameter transformation that leads from Eq. 3 to Eq. 4 simplifies Eqs. 18–20. In slight abuse of notation, we use the same symbol for the *expected* compound rates as for the constant compound rates of Eq. (4), i.e., write ν* _{k}*: =

$${\kappa}_{1}[Z]={\overrightarrow{\xi}}_{1}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h$$

(21)

$${\kappa}_{2}[Z]={\overrightarrow{\xi}}_{2}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h+{\left({\overrightarrow{\xi}}_{1}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}\right)}^{2}{h}^{2}{\beta}_{2}$$

(22)

$${\kappa}_{3}[Z]={\overrightarrow{\xi}}_{3}\cdot {\overrightarrow{\nu}}_{\xi}h+{\left({\overrightarrow{\xi}}_{1}\cdot {\overrightarrow{\nu}}_{\xi}\right)}^{3}{h}^{3}{\beta}_{3}+3\left({\overrightarrow{\xi}}_{1}\cdot {\overrightarrow{\nu}}_{\xi}\right)\left({\overrightarrow{\xi}}_{2}\cdot {\overrightarrow{\nu}}_{\xi}\right){h}^{2}{\beta}_{2}$$

(23)

To adapt CuBIC to non-stationary populations, we need to formulate the general maximization problem (Eq. 5) for the case of time- dependent carrier rates. Using Eqs. 21–23 instead of Eq. 4, we obtain

$${\text{\kappa}}_{3,\text{\xi}}^{*,ns}=\underset{{\overrightarrow{\nu}}_{\xi},{\text{\beta}}_{2},{\text{\beta}}_{3}}{\mathrm{max}}\left\{{\overrightarrow{\xi}}_{\text{3}}\cdot {\overrightarrow{\nu}}_{\text{\xi}}h+{\left({\overrightarrow{\xi}}_{\text{3}}\cdot {\overrightarrow{\nu}}_{\text{\xi}}\right)}^{3}{h}^{3}{\text{\beta}}_{3}+3\left({\overrightarrow{\xi}}_{\text{1}}\cdot {\overrightarrow{\nu}}_{\text{\xi}}\right)\left({\overrightarrow{\xi}}_{\text{12}}\cdot {\overrightarrow{\nu}}_{\text{\xi}}\right){h}^{2}{\text{\beta}}_{\text{2}}\right\}$$

(24)

$$\begin{array}{c}\text{subject}\hspace{0.17em}\text{to}\hspace{0.17em}{\kappa}_{2}[Z\prime ]={\overrightarrow{\xi}}_{2}\cdot {\overrightarrow{\nu}}_{\xi}h+{\left({\overrightarrow{\xi}}_{1}\cdot {\overrightarrow{\nu}}_{\xi}\right)}^{2}{h}^{2}{\beta}_{2}\\ {\kappa}_{1}[Z\prime ]={\overrightarrow{\xi}}_{1}\cdot {\overrightarrow{\nu}}_{\xi}h,\end{array}$$

with ${\overrightarrow{\nu}}_{\xi}\in {[0,\infty )}^{\xi}$, β_{2}[0,∞) and β_{3}(−∞,∞). After some algebra and substitution of the unknown cumulants κ* _{i}*[Z′] by their estimates

$${\kappa}_{3,\xi}^{\ast ,ns}=\underset{{\overrightarrow{\nu}}_{\xi},{\beta}_{2},{\beta}_{3}}{\text{max}}\left\{{\overrightarrow{\xi}}_{3}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h+{k}_{1}^{3}{\beta}_{3}-3{k}_{1}^{3}{\beta}_{2}^{2}+3{k}_{1}{k}_{2}{\beta}_{2}\right\}$$

(25)

$$\begin{array}{c}\text{subject}\hspace{0.17em}\text{to}\hspace{0.17em}{k}_{2}={\overrightarrow{\xi}}_{2}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h+{k}_{1}^{2}{\beta}_{2}\\ {k}_{1}={\overrightarrow{\xi}}_{1}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h.\end{array}$$

As opposed to the Linear Programming Problem of the stationary case (Eq. 6), the constraints in Eq. 25 do not apply to all free variables: the third standardized cumulant of the rate variable β_{3} can, in principle, be arbitrarily large. As a consequence, also the objective function in Eq. 25 is unbounded. We therefore have to impose additional constraints on the carrier distribution *f _{R}* in order to ensure convergence of the maximization. The approach taken here is to prescribe a two-dimensional parametric family for the distribution of the rate variable, such that its third (standardized) cumulant can be expressed in terms of its first two cumulants. The choice for the family of carrier distributions determines the form of the objective function. Let us present two specific cases (see Section “Discussion” for more details on the role of this choice).

If the carrier distribution is symmetric about its mean, like e.g., the uniform distribution, we can exploit the fact that the third cumulant κ_{3}[*R*] of symmetric distributions vanishes^{1}. As a consequence, also β_{3} =κ_{3}[*R*]/κ_{1}[*R*]^{3} =0. The objective function of the maximization problem (Eq. 25) for symmetric carrier distributions thus becomes

$$F\left({\nu}_{1},\dots ,{\nu}_{\xi},{\beta}_{2}\right)={\overrightarrow{\xi}}_{3}\cdot {\overrightarrow{\nu}}_{\xi}h+3{\beta}_{2}{k}_{1}{k}_{2}-3{k}_{1}^{3}{\beta}_{2}^{2}.$$

(26)

This objective function is linear in the (expected) compound rates ν* _{k}* (

If the carrier variable *R* follows a Gamma-distribution, the third cumulant can be expressed in terms of the first two as ${\kappa}_{3}[R]=2{\kappa}_{2}{[R]}^{2}\text{\hspace{0.05em}}/{\kappa}_{1}[R]$. For the normalized third cumulant β_{3} we thus have ${\beta}_{3}=2{\kappa}_{2}{[R]}^{2}/{\kappa}_{1}{[R]}^{4}=2{\beta}_{2}^{2}$. The objective function then reads

$$F\left({\nu}_{1},\dots ,{\nu}_{\xi},{\beta}_{2}\right)={\overrightarrow{\xi}}_{3}\cdot {\overrightarrow{\nu}}_{\xi}h+3{\beta}_{2}{k}_{1}{k}_{2}-{k}_{1}^{3}{\beta}_{2}^{2}.$$

(27)

As for symmetric carrier distributions, the objective function in Eq. 25 is linear in the expected component rates ν* _{k}* (

We wish to stress once again that the choice made here concerns only the family of the carrier distribution, not its particular shape. That is, if we chose e.g., a uniform distribution, we only determine that κ_{3}[*R*] =0 but do not fix the support of the distribution. Finally, we note that the only non-zero components of the solution ${\overrightarrow{\nu}}_{\xi}^{\ast}=\left({\nu}_{1}^{\ast},\dots ,{\nu}_{\xi}^{\ast}\right)$ to the maximization problem are ${\nu}_{1}^{\ast}$ and ${\nu}_{\xi}^{\ast}$, i.e., ${\nu}_{k}^{\ast}=0$ for all *k* {1,ξ}, as in the stationary case (see The Solution of the Maximization in Appendix).

As a final ingredient, CuBIC requires the variance of the test statistic *k*_{3} in order to be able to compute the *p*-values. Assuming that the solution of Eq. 25 is available, we compute the second, fourth and sixth cumulant of the solution by inserting the solving ${\overrightarrow{\nu}}_{\xi}^{\ast}$ and ${\beta}_{2}^{\ast}$ into the algorithm leading to Eq. 17. Plugging these cumulants into Eq. 7 yields Var[*k*_{3}], the variance of the test statistic *k*_{3} under the non- stationary null- hypothesis. Insertion into Eq. 8 yields the corresponding *p*-value.

To illustrate the application of the adapted CuBIC, we consider here two typical experimental situations. In the first situation, data are recorded in early sensory systems, where characteristic firing rate profiles closely follow the stimulus. In such a scenario, information about the rate distribution can be obtained from the type of the stimulus. In the second situation, there is no direct relation between the experimental paradigm and non-stationarities in firing rates. In this case, ad hoc assumptions of the family of carrier rate are the only option. We illustrate both scenarios and the steps required when applying the non-stationary version of CuBIC by analyzing simulated spike trains.

Our first example mimics a recording in visual cortex with a oriented sinusoidal grating as stimulus. We model the population response in such an experiment by a CPP model with sinusoidal carrier rate ν(*t*), i.e.,

$$v(t)=B+C\mathrm{cos}(2\text{\pi}\omega t-d),$$

(28)

where *B* is the offset, *C*≤ *B* is the amplitude of the modulation, *ω* is the temporal frequency of the stimulus and *d* is the phase, i.e., the sum of the stimulus phase and the delay it takes for the stimulus-driven activity to propagate to the recording electrodes (Figure (Figure3,3, top panels of left and right columns). The first data set of this example models the case of pure rate non-stationarity (Figure (Figure3,3, left column). The amplitude distribution is a delta-peak at ξ=1, such that all events of the marked process *z*(*t*) have amplitude *a _{j}*=1. As a consequence,

To analyze the data, we chose a bin size of *h*=5ms, computed the population spike count *Z*(*s*), and estimated its distribution *f _{Z}*(

Application of the adapted CuBIC requires the specification of a parametric family of carrier distributions *f _{R}* (“carrier family”). We put ourselves in the position of an experimenter, to whom the carrier rate is unknown and cannot be estimated directly; the only observable quantity is the population activity, which is a combination of the carrier rate and potential events of high amplitude. However, as the stimulus was a cosine, we assume also the carrier rate to be a cosine as in Eq. 28, only that the parameters

In the first data set, the stationary CuBIC assigned the rate-generated correlation among the counting variables *X _{i}* to events of amplitude ≥2. Allowing for potential (co-)variations of firing rates in the adapted CuBIC corrected for this faulty inference of correlations. Mathematically, allowing for rate (co-)variations allows non-zero values of the parameter β

Compared to the size of the population (*N*=50), the rate of the high-amplitude events (ν_{7}= κ_{1}[*I*]·*f _{A}*(7)≈ 6Hz) is relatively small. As a consequence, these are hardly visible in the raster displays and population spike counts (Figure (Figure3,3, second and third panels from top in the middle column). In the distributions of the population spike count, they lead to a slight increase for the frequency of patterns of size

In this data set, the stationary CuBIC rejects all hypotheses up to a value of ξ =6. Hence, CuBIC performs optimally in this data set, as the resulting lower bound $\widehat{\xi}=7$ corresponds to the maximal order of correlation ξ_{syn} =7 in this data set.

To our big surprise and satisfaction, *p*-values for the adapted version of CuBIC were very similar to those of the stationary CuBIC (Figure (Figure3,3, bottom panel in middle column, red and green bars respectively). In particular, the adapted CuBIC rejected the same hypotheses, and, as a consequence, also yielded the same lower bound $\widehat{\xi}=7.$ Contrary to our assumption, the proposed adaptation thus did not compromise CuBIC's sensitivity in this stationary data set.

Finally, we generated a third data set that combined the properties of the first two examples (Figure (Figure3,3, right column). The amplitude distribution was the same as in the example of pure correlation, i.e., with an additional entry at ξ_{syn} =7, while the carrier rate was the same cosine as in the data with pure non-stationarity (Figure (Figure3,3, top panels of right and left columns, respectively).

For this data set, we expected the stationary CuBIC to overestimate the order of correlation, i.e., to yield $\widehat{\xi}>{\xi}_{\text{syn}}=7$, as the considerable rate co-variation produces correlation among the *X _{i}* in addition to the events of high amplitude (Staude et al., 2008). To the contrary, however,

To investigate the robustness of the proposed adaptation with respect to faulty choices for the carrier family, we also analyzed the data under the assumption that the carrier rate switches between a state of low rate, ν_{min}, and of high rate, ν_{max}. This can be described by a bimodal carrier distribution

*f _{R}*(ν;ν

where δ(ν) is the Dirac-delta function and η[0,1] parametrizes the relative proportion of bins where ν(*t*) is in the low (high) rate state. Inspecting the time course of *Z*(*s*), we chose the mass of the two peaks identical (η =0.5), which leaves a two-parameter family for *f _{R}* (its cumulants are derived in Section “Cumulants of Carrier Distributions” in Appendix). For the data shown in Figure Figure3,3, allowing for such drastic rate dynamics hardly changed

In the examples of the previous section, we inferred the type of non-stationarity from the experimental paradigm, i.e., from the properties of the stimulus. In many experimental situations, such a priori assumptions on the rate variable cannot be justified, and changes in firing rates can have diverse statistical properties. General excitability of the tissue can change firing rates either in the form of slow drifts or abrupt transitions between so called up- and down-states. Furthermore, both local computations and additional cortical or sub-cortical inputs may change firing rates of the recorded population.

A further feature of the previous examples was that the carrier rate ν(*t*) changed rather slowly. As mentioned above, however, the temporal dynamics of the carrier rate does not influence the statistics required for CuBIC, i.e., the population spike count *Z*, as long as the distribution of the rate values per bin, *f _{R}*, is not altered. The carrier rate ν(

$$\nu (t)={\displaystyle \sum _{s=1}^{L}{R}_{s}{\Pi}_{s}(t),}$$

where *L*=*T*/*h* is the sample size (number of bins) Π* _{s}*(

The gamma-distributed rate variable generates strongly fluctuating rate profiles, with peak values of the carrier rate above 3000Hz. This leads to strong fluctuations in the population spike count even in the case of pure rate non-stationarity (Figure (Figure4,4, third panel from top in left column), such that its distribution has a fairly elongated tail (Figure (Figure4,4, fourth panel from top in left column).

Despite quantitative differences, stationary CuBIC performs qualitatively similarly for gamma- and cosine-like carrier rate. For pure rate non-stationarity, it wrongly interprets the rate-induced correlations among the *X _{i}* as events of high amplitudes. The null-hypotheses are rejected up to ξ=3, yielding a lower bound of $\widehat{\xi}=4$ (Figure (Figure4,4, green arrowheads in the bottom panel of the left column). Similarly, adding gamma-like non-stationarities to a data set with correlation decreases the inferred lower bound, here from $\widehat{\xi}=7$ in the stationary data (Figure (Figure4,4, green arrowheads in the bottom panel of the middle column) $\widehat{\xi}=6$ to in the non- stationary data (Figure (Figure4,4, green arrowheads in the bottom panel of the right column).

As opposed to Section “Stimulus-driven Non-stationarity” where properties of the carrier rate could be inferred from the stimulus, we here cannot make qualitative guesses about the type of non-stationarity. The fact that firing rates fluctuate on a bin-to-bin basis makes it very difficult to infer the type of non-stationarity from the raster plots, the population spike counts *Z* or its distribution *f _{Z}*. As a consequence, we can only make ad hoc assumptions on

For the data with pure rate non-stationarity (Figure (Figure4,4, left column), allowing a uniform rate variable rejects hypotheses up to ξ=3 and thus yields a lower bound of $\widehat{\xi}=4$. Allowing *f _{R}* to be gamma-distributed, on the other hand, produces

For data with correlation, the choice of the carrier distribution has only little influence on the resulting *p*-values (Figure (Figure4,4, bottom panels in middle and left column). For pure correlation, the lower bounds are identical for all three rate models ($\widehat{\xi}=7={\xi}_{\text{syn}}$). In the combined data set (Figure (Figure4,4, right column), the additional non-stationarity reduces the lower bound as compared to the stationary data with correlation to $\widehat{\xi}=6$, irrespective of the non-stationarity allowed in the null-hypotheses.

Besides the stationarity, CuBIC's second central assumption is that spike trains of single neurons can be described as Poisson processes, i.e., have exponential interspike interval (ISI) distributions. While tails of ISI distributions can often be relatively well described by exponentials, the high probability for short intervals of the Poisson process contradicts the absolute refractory period of a few milliseconds found in most nerve cells. We investigated the extent to which refractoriness influences test results of CuBIC by re-analyzing the data of the previous sections after short ISIs were removed. Specifically, for each data set we assigned the events of the carrier process randomly to the *N*=50 spike trains, removed spikes of all spike trains that occurred with a temporal difference of τ≤ 2 ms, and constructed the population spike counts of these thinned spike trains. The analysis of the refractory data showed that introducing refractoriness has a very limited effect on test results (outlined bars and arrowheads in lower panels of Figures Figures33 and and4).4). If *p*-values changed, they increased, which makes CuBIC more conservative but does not generate false positives. From all simulation, the increased *p*-value reduced the lower bound $\widehat{\xi}$ only in a single instance (cosine rate modulation with correlation, analyzed with stationary CuBIC at ξ=4; green arrowhead in bottom right panel of Figure Figure3);3); in all remaining cases the refractory period changed *p*-values where they were above the significance level of α=0.05. We thus conclude that CuBIC is reasonably robust if spike deviate from Poisson processes in terms of short refractory periods (here 2ms).

The analysis of massively parallel spike trains for higher-order correlations is a prerequisite for understanding the mechanisms of cooperative neuronal computation in the brain. However, analysis techniques to estimate higher-order correlations typically require enormous sample sizes, rendering the analysis of large (*N*>10) populations for higher-order effects extremely difficult. In Staude et al. (2009), we have presented a novel method (CuBIC) that avoids the need for extensive sample sizes. Numerical simulations illustrated that CuBIC reliably infers correlations of very high order (>10) from recordings of large (*N* ~ 100), even weakly correlated neuronal populations (average pairwise correlation coefficient *c*<0.01) with sample sizes corresponding to 10–100s recording time.

As a statistical model, CuBIC assumes the CPP, where correlation among the spike trains is modeled by the insertion of additional coincident events (Ehm et al., 2007; Johnson and Goodman, 2007; Brette, 2009; Staude et al., 2010). The linear relation between the model parameters, i.e., the orders of correlation present in the data, and the cumulants κ* _{m}*[

Before going into detail, we need to make a general remark. CuBIC is a parametric procedure, and assumes that the data, i.e., the population spike count, can be described sufficiently well by a discretized, potentially non-stationary CPP. If this model does not fit, results of CuBIC may be unreliable. The extent to which results are reliable then depends on the mismatch between the CPP and the data. In practice, where single spike trains deviate from Poisson processes, for example due to refractory properties, this mismatch may be evaluated by analyzing surrogate data (e.g., Grün, 2009). If, for instance, CuBIC returns $\widehat{\xi}=10$ in a data set of non-Poissonian spike trains, and the analysis of surrogate data with the same single process properties but without correlation yields $\widehat{\xi}=1$, we can conclude that the value of $\widehat{\xi}=10$ is really due to existing correlations. This kind of analysis, however, has to be performed specifically for a given data set. As a first step, we have here investigated the effect of a spike train's most obvious deviation from the Poisson process: absolute refractory periods (here 2ms, Figures Figures33 and and4).4). Its relatively small effect on *p*-values makes us confident that CuBIC is a promising analysis method even if single spike trains deviate from the Poisson assumption. The detailed analysis of CuBIC's robustness with respect to these violation is a central aspect of our current research (e.g., Staude et al., 2007; Reimer et al., 2009). We wish to stress that in the case of small bin sizes and hard refractory periods, assuming single processes to have Poisson statistics is essentially identical to the popular assumption of independence of subsequent bins in Bernoulli models (as in e.g., Martignon et al., 1995; Shlens et al., 2006). A more detailed discussion of this issue, together with an analysis of the parameter dependence of CuBIC can be found in Staude et al. (2009).

A central assumption in the original presentation of CuBIC (Staude et al., 2009) is that the statistics of the population does not change over time (stationarity). In the present study, we have adapted CuBIC to be able to analyze also non-stationary data. The basis of this adaptation is a non-stationary version of the CPP, where the intensity of the population, parametrized by the ν, is decoupled from the correlation structure, parametrized by the amplitude distribution *f _{A}*. Describing the population spike count as a doubly stochastic CPP, potential non-stationarities in the data can by integrated into the quantification of the null-hypotheses of CuBIC. We wish to stress once again, however, that non- stationarities in single neurons do not necessarily imply time varying carrier rates (see, e.g., Figure Figure2A),2A), such that not every non-stationary data set requires the application of the adapted CuBIC.

In this study, we presented the adaptation only for the third cumulant, i.e., *m*=3. Although the derivation of the mathematical equations is straightforward also for higher values of *m*, the resulting expressions become increasingly complicated. This may result in particular in strong non-linearities in the maximization problem, such that additional care is necessary to ensure convergence of numerical solvers at the global maximum. Furthermore, the estimation of cumulants of order >3 becomes unreliable and their estimators have large variance, which may compromise test power. As CuBIC proved to be very sensitive even for *m*=3 (Staude et al., 2009), we currently have not developed the theory for higher *m*. Nevertheless, this might be a fruitful direction of future research.

The main difference between the stationary CuBIC and its non-stationary adaptation lies in the maximization of the third cumulant of *Z*. Here, the adapted version requires that the third standardized cumulant β_{3} =κ_{3}[*R*]/κ_{1}[*R*]^{3} of the binned carrier rate *R* does not assume arbitrary large values. In this study, this is achieved by prescribing a two-parameter family for the carrier distribution *f _{R}* (the “carrier family”). The remainder of this section is therefore primarily concerned with elucidating the role of this particular choice.

Classically, correlations are measured on a small time scale, and subsequently corrected for effects from the firing rates. The estimation of the firing rate, in turn, proceeds along one of the following two lines. Either, rate variations are identified with artifacts that are locked to the stimulus or some other external cue. Then, firing rates are estimated by averaging over repeated presentation of the stimulus, or trials (e.g., Perkel et al., 1967a,b; Aertsen et al., 1989). Alternatively, changes in firing rates are assumed to fluctuate only on slow time scales; they are then estimated by averaging over time. Although combinations, refinements and optimizations of both strategies have been developed (e.g., Grün et al., 2002b; Ventura et al., 2005; Shimazaki et al., 2009), we wish to stress that all such corrections make strong a priori assumptions on the differences between “genuine” correlations and “artifacts” induced by non-stationary firing rates (see also Staude et al., 2008).

The correction of CuBIC, i.e. the choice of a carrier family, follows a fundamentally different strategy. Rather than assuming that firing rate profiles are identical over trials (first strategy) or that rate-fluctuations are band-limited (second strategy), the carrier family limits the extent to which correlations among the counting variables *X _{i}* are assigned to rate-variations. As discussed above, this choice ensures boundedness of the third standardized cumulant β

In Section “Case Studies”, we illustrated how CuBIC operates in two alternative scenarios. In Section “Stimulus-driven Non-stationarity”, we assumed that properties of the firing rates can be inferred from the stimulus. Here, the stimulus was a cosine, but the reasoning there can be generalized to a broad class of stimuli. The effect of oriented bars, for instance, might be modeled by a bimodal carrier distribution as in Section “Different Carrier Family”,

$${f}_{R}\left(\nu ;{\nu}_{\text{min}},{\nu}_{\text{max}},\eta \right)=(1-\eta )\delta \left(\nu -{\nu}_{\text{min}}\right)+\eta \delta \left(\nu -{\nu}_{\text{max}}\right),$$

(29)

where the mixing parameter η determines the relative duration of the respective stimulus phases (light/dark). The carrier family of Eq. 29 also describes experiments with a well-defined stimulus onset, as e.g., odor presentation or whisker stimulation. In such a scenario, η has to be chosen as the relative duration of the stimulus-off/stimulus-on epoch. Furthermore, properties of the carrier family in “free viewing” conditions might be estimated from the statistical properties of the visual scene.

In the second scenario (see Internally Generated Non-stationarities), the experimental paradigm did not provide information about the carrier family. Here, we argued that ad hoc assumptions are the only option. However, if the activity of individual neurons is available, their statistics can provide additional information that can be exploited. If, for instance, individual spike trains do not show evidence of time varying firing rates, or provide upper bounds for the cumulants of their firing rates, this information may be extrapolated to the level of the carrier rate. Together with general moment inequalities (e.g., Kumar, 2002), such information may help to dispose of the explicit choice of the carrier family by providing an explicit upper bound for β_{3}. Although the absence of a precise parametric model for the carrier distribution impedes the faithful computation of Var[*k*_{3}] under ${H}_{0}^{3,\xi}$, the close similarity of the error-bars in Figure Figure5A5A for the different methods makes us confident that Var[*k*_{3}] can be reasonably approximated by any carrier distribution as long as the upper bound ${\kappa}_{3,\xi}^{\ast ,ns}$ is faithful (see also Section ‘What is the “true” Carrier Family?’ for conceptual issues concerning the choice for *f _{R}*).

The analysis of artificial data (Figures (Figures33 and and4)4) can be summarized by four major points:

- For data with pure firing rate non-stationarity but no higher-order correlations, the stationary CuBIC misinterprets the common rate variations as events of amplitude ≥2, i.e., as correlations. The order of the inferred correlation depends on the kind of non-stationarity in the data ($\widehat{\xi}=2$ for cosine carrier rate, $\widehat{\xi}=4$ for gamma-distributed rate variable).
- Allowing for potential non-stationarities in the null-hypotheses can reduce this lower bound. Using the correct family for the carrier distribution, i.e., a cosine-like model for the data in Figure Figure33 and a gamma-distribution for the data in Figure Figure4,4, corrects entirely for the false positive inference of the stationary CuBIC, such that the lower bound becomes $\widehat{\xi}=1$. Choosing the wrong family, however, may not account properly for rate variations, especially if the family assigns smaller values of β
_{3}(e.g., allowing a uniform carrier distribution if the data had a gamma-distributed carrier rate as in Figure Figure44). - For stationary data with correlation, allowing for non- stationarities in the null-hypotheses has no effect on the inferred lower bound. This result holds irrespective of the type of non-stationarity allowed in the null-hypotheses.
- Non-stationarities in data with correlation reduce the inferred lower bound as compared to data with the same correlation structure but constant firing rates. The degree of reduction depends mostly on the kind of non-stationarity in the data. The family allowed for the carrier distribution did not affect the lower bound.

It is well-known that co-variations in firing rates induce correlations between spike counts. Thus, it comes as no surprise that analyzing non-stationary data with the stationary CuBIC generates false positives (point 1). In the adaptation, parts of the correlation can be assigned to (co)-variations in firing rates. Allowing for non-stationarities therefore corrects for this faulty inference (point 2).

To understand why the adaptation did not generally reduce the lower bound (point 3), we investigated the interplay between the constraints and the objective function in the maximization of the third cumulant in further detail. Due to the increased number of degrees of freedom, our intuition was that allowing for non- stationarities additional to the correlations of order ξ would simply increase the maximal value of the third cumulant ${\kappa}_{3,\xi}^{\ast}$ as compared to the stationary version. This, however, seems not to be the case. In the data sets with pure correlation (Figures (Figures33 and and4,4, middle panels), for example, the values of ${\kappa}_{3,\xi}^{\ast}$ and ${\kappa}_{3,\xi}^{\ast ,ns}$ are identical for ξ ≥4 (Figure (Figure5B5B shows results for stationary, uniform and gamma-distributed rate variables, results of combined data sets are very similar, data not shown). The reason is that for a given value of the first estimated cumulant *k*_{1}, the constraint ${k}_{2}={\overrightarrow{\xi}}_{2}\cdot {\overrightarrow{\nu}}_{\xi}h+{k}_{1}^{2}{\beta}_{2}$ penalizes stronger non-stationarity, i.e., larger values of the standardized variance β_{2}, with a reduction of the (expected) component rates ν_{1},…,ν_{ξ}, and vice versa. In the objective function (Eq. 25)

$$F\left({\nu}_{1},\dots ,{\nu}_{\xi},{\beta}_{2}\right)={\overrightarrow{\xi}}_{3}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}h+{k}_{1}^{3}{\beta}_{3}-3{k}_{1}^{3}{\beta}_{2}^{2}+3{k}_{1}{k}_{2}{\beta}_{2},$$

(30)

the component rates enter via ${\overrightarrow{\xi}}_{3}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}$, while the β_{2}-dependence is a parabola with negative curvature. As

$${\overrightarrow{\xi}}_{3}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi}={\displaystyle \sum _{k=1}^{\xi}{k}^{\text{3}}{\nu}_{k}}\gg {\displaystyle \sum _{k=1}^{\xi}{k}^{\text{2}}{\nu}_{k}}={\overrightarrow{\xi}}_{\text{2}}\text{\hspace{0.05em}}\cdot {\overrightarrow{\nu}}_{\xi},$$

(31)

especially for large values of ξ, the objective function profits more from high component rates than the constraint penalizes these. As a consequence, the maximization favors high component rates over strong rate fluctuations, especially for large ξ. The results of the maximization procedure for the data with pure correlation supports this interpretation, as the standardized variance of the model that maximizes κ_{3}[*Z*], ${\beta}_{2}^{\ast}$, decreases with ξ (Figure (Figure5A).5A). For ξ ≥4, we have ${\beta}_{2}^{\ast}=0$, and, as a consequence, the maximizing model of the adapted maximization problem is the same as that of the stationary maximization. Consequently, also the solutions ${\kappa}_{3,\xi}^{\ast}$ and the *p*-values are identical. Contrary to our initial intuition, the maximization thus does not generally favor strong rate variations. Evidently, however, the extent to which the inclusion of non- stationarities in the null-hypothesis alter test results depends crucially on the parameters of the data.

The results of the case studies summarized above suggest that including potential non-stationarities in the null-hypothesis is always a safe bet: it corrects for false positive inference if correlation originates from rate effects, but does not alter *p*-values if the stationary CuBIC did not overestimate the order of correlation. To sketch the parameter range where including non-stationarities reduces *p*-values only if necessary, recall that we have identified the reason for the unchanged maximal third cumulant in the interplay between the constraint and the objective function. Considering the general objective function (Eq. 25) however, we find that the influence of the non-stationarity (via β_{2} and β_{3}) depends on the value of the first sample cumulant *k*_{1}. For *k*_{1} <1, non-stationarities (positive β_{2} and β_{3}) hardly influence the objective function, hence the maximization can be assumed to favor high component rates over non-stationarities, yielding identical test results for the stationary and the adapted CuBIC. For *k*_{1}1 a small increase in β_{3} has a strong effect on the objective function, which may in turn favor strong non-stationarities, thereby producing different test results for the stationary and the adapted method. Thus, a crucial parameter for the performance of the adaptation is the first sample cumulant *k*_{1}. Now *k*_{1} is the estimator of ${\kappa}_{1}[Z]={\kappa}_{1}\left[{\Sigma}_{i=1}^{N}{X}_{i}\right]=h{\Sigma}_{i=1}^{N}{\lambda}_{i},$ where λ* _{i}* is the (average) firing rate of the

In Section “Correcting for Rate Effects” we have presented a few guidelines for the choice of the family of carrier rate distributions. There are cases, however, where this choice is not easily justifiable by resorting to observable quantities. To discuss the status of this problem in more depth, we wish to stress that measures of correlation and their various corrections are purely statistical in nature. To be of scientific value, statistical results have to be put in context and must be interpreted, e.g., in terms of biophysical mechanisms. The “firing rate” of a neuron, for example, is not a biophysical entity as such. The term arises only if one describes the variable behavior of a neuron using point processes. As a consequence, whether or not the choice of a chosen carrier family *f _{R}* is “valid” depends entirely on the intended biological interpretation. If, for example, the choice of

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

We are grateful to Stuart Baker for a fruitful discussion, Leona Schild for help in the initial phase of this project, and Imke Reimer for comments on an earlier version of the manuscript. Supported by the German Federal Ministry of Education and Research (BMBF grants 01GQ0420 and 01GQ01413 to the BCCN Freiburg and Berlin), the Helmholz Alliance on Systems Biology (Germany), and SFB 780.

###### Symbol

###### Meaning

- x
_{i}(*t*) *i*th spike train in continuous time*X*_{i}- Counting variable of
*i*th spike train *z*(*t*)- Carrier process, i.e., summed activity in continuous time
*Z*- population spike count
- κ
_{i}[*Z*] *i*th cumulant of*Z**k*_{i}*i*th sample cumulant of*Z*, i.e.,*k*-statistic*h*- Bin width
- λ
- Average firing rate of individual spike trains
*f*_{A}- Amplitude distribution, i.e., population-average correlation structure
- μ
_{i} *i*th raw moment of amplitude distribution- ν
- Carrier rate
- ν
_{k} - (Expected) compound rate of all events with amplitude
*k* - ${\overrightarrow{\nu}}_{\xi}$
- Vector of (expected) compound rates ν
_{1},…,ν_{ξ} - ${\overrightarrow{\xi}}_{i}$
- Vector of
*i*th powers of 1,…,ξ - ${H}_{0}^{3,\xi}$
- Null-hypothesis stating that the data has correlation of maximal order ξ
- ${\kappa}_{3,\xi}^{\ast}$
- Maximal value of κ
_{3}[*Z*] given correlations of orders ≤ξ *R*_{s}- Carrier variable, mean value of the carrier rate in the sth bin
*f*_{R}- Carrier distribution: distribution of the {
*R*}_{s}_{s} - β
_{k}*k* - th standardized cumulant of the carrier distribution
- ${\kappa}_{3,\xi}^{\ast ,ns}$
- Maximal value of κ
_{3}[*Z*] given correlations of orders ≤ξ and non-stationarity

This Appendix shows that the solution $({\nu}_{1}^{\ast},\dots ,{\nu}_{\xi}^{\ast})$ of the stationary (Eq. 6) and the non-stationary maximization problem (Eq. 25) fulfills ${\nu}_{k}^{\ast}=0$ for *k*=2,…,ξ −1.

The non-stationary problem (Eq. 25) has the objective function

$$F={\displaystyle \sum _{k=1}^{\xi}{k}^{3}{\nu}_{k}+{k}_{1}^{3}{\beta}_{3}-3{k}_{1}^{3}{\beta}_{2}^{2}+3{k}_{1}{k}_{2}{\beta}_{2}}$$

(32)

with constraints

$${k}_{1}={\displaystyle \sum _{k=1}^{\xi}k{\nu}_{k}}$$

(33)

$${k}_{2}={\displaystyle \sum _{k=1}^{\xi}{k}^{2}{\nu}_{k}+{k}_{1}^{2}{\beta}_{2}}$$

(34)

Simple computations starting from Eqs. 33 and 34 yield

$$\begin{array}{l}{\nu}_{1}=\frac{1}{\xi -1}\left(\xi {k}_{1}-{k}_{2}-{k}_{1}^{2}{\beta}_{2}-{\displaystyle \sum _{k=2}^{\xi -1}\left(\xi k-{k}^{2}\right){\nu}_{k}}\right)\\ {\nu}_{\xi}=\frac{1}{\xi (\xi -1)}\left({k}_{2}-{k}_{1}+{k}_{1}^{2}{\beta}_{2}-{\displaystyle \sum _{k=2}^{\xi -1}\left({k}^{2}-k\right){\nu}_{k}}\right),\end{array}$$

which, after insertion into Eq. 32 yield

$$F=\frac{1}{\xi -1}{\displaystyle \sum _{k=2}^{\xi -1}\left({k}^{3}+{k}^{2}\left(1-{\xi}^{2}\right)+k\left({\xi}^{2}-\xi \right)\right){\nu}_{k}+H,}$$

(35)

with

$$H={k}_{1}\left({k}_{1}{\beta}_{2}(\xi +1)-\xi \right)+{k}_{2}(\xi +1)+{k}_{1}^{3}{\beta}_{3}-3{k}_{1}^{3}{\beta}_{2}^{2}+3{k}_{1}{k}_{2}{\beta}_{2}.$$

Now for *k*=2,…,ξ −1 we have

$$\frac{\partial F}{\partial {\nu}_{k}}=\frac{k}{\xi -1}\left({k}^{2}+k\left(1-{\xi}^{2}\right)+{\xi}^{2}-\xi \right)$$

(36)

$$k\le \xi -1\hspace{0.17em}{k}^{2}+{\xi}^{2}-{\xi}^{2}k+k-\xi $$

(37)

$$={k}^{2}+{\xi}^{2}-{\xi}^{2}k+k-\xi $$

(38)

$$k<\xi \hspace{0.17em}{\xi}^{2}+{\xi}^{2}-{\xi}^{2}k+\xi -\xi $$

(39)

$$=2{\xi}^{2}-k{\xi}^{2}$$

(40)

$$2\le k\hspace{0.17em}0$$

(41)

The gradient of the objective function therefore points to negative values of ν* _{k}* for

For the computation of Var[*k*_{3}] (Eq. 7) from the solutions of the maximization procedure, i.e., the parameters ${\overrightarrow{\nu}}_{\xi}^{\ast}$ and ${\beta}_{2}^{\ast}$, we require explicit expressions for the cumulants of the carrier distribution *f _{R}* up to order

For ν(*t*) =*B*+*C*cos(2π*ω*t +*d*), we derive the distribution of the rate values ${R}_{s}=1/h{\int}_{sh}^{(s+1)h}\nu (t)dt$ under the assumption that subsequent values of *R _{s}* sample the cosine faithfully. In this case, we can express

$$\begin{array}{c}{f}_{R}\left(r\right)={f}_{T}\left({g}^{-1}\left(r\right)\right)\left|\left(\frac{d{g}^{-1}\left(r\right)}{dr}\right)\right|\\ =\frac{1}{C\sqrt{1-{\left(\frac{R-B}{C}\right)}^{2}}}\end{array}$$

The first six moments can be computed by solving the integrals

$$\text{E}\left[{R}^{m}\right]={\displaystyle \underset{B-C}{\overset{B+C}{\int}}{r}^{m}}\frac{1}{C\sqrt{1-{\left(\frac{r-B}{C}\right)}^{2}}}dr,$$

which yields

$$\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}E\left[R\right]=B\\ E\left[{R}^{2}\right]={B}^{2}+\frac{{C}^{2}}{2}\\ E\left[{R}^{3}\right]={B}^{3}+\frac{3B{C}^{2}}{2}\\ E\left[{R}^{4}\right]={B}^{4}+3{B}^{2}{C}^{2}+\frac{3{C}^{4}}{8}\\ E\left[{R}^{5}\right]={B}^{5}+5{B}^{3}{C}^{2}+\frac{15B{C}^{4}}{8}\\ E\left[{R}^{6}\right]={B}^{6}+\frac{15{B}^{4}{C}^{2}}{2}+\frac{45{B}^{2}{C}^{4}}{8}+\frac{5{C}^{6}}{16}\end{array}$$

Let *f _{R}*(ν;ν

$$E\left[{R}^{m}\right]=(1-\eta ){\nu}_{\text{min}}^{m}+\eta {\nu}_{\text{max}}^{m}.$$

If *R* is uniformly distributed between *a* and *b*, the raw moments of *R* are given as

$$E\left[{R}^{m}\right]=\frac{{b}^{m+1}-{a}^{m+1}}{(m+1)(b-a)}.$$

The moments of the gamma distribution $f(R;k,\theta )={R}^{k-1}\frac{{e}^{-R/\theta}}{{\theta}^{k}\Gamma (k)}$ with parameters *k* and θ are given as

$$E\left[{R}^{m}\right]=\frac{{\theta}^{m}\Gamma (k+m)}{\Gamma (k)},$$

where Γ denotes the gamma-function.

^{1}Generally, κ_{3}[*R*] is a measure for the skewness of *fR*, such that κ_{3}[*R*]<0 for left- skewed, κ_{3}[*R*]=0 for symmetric, and κ_{3}[*R*]>0 for right-skewed distributions.

- Abbott L. F., Dayan P. (1999). The effect of correlated variability on the accuracy of a population code. Neural Comput. 11, 91–10110.1162/089976699300016827 [PubMed] [Cross Ref]
- Abeles M. (1991). Corticonics: Neural Circuits of the Cerebral Cortex, 1st Edn.Cambridge: Cambridge University Press
- Aertsen A., Gerstein G., Habib M., Palm G. (1989). Dynamics of neuronal firing correlation: modulation of “effective connectivity”. J. Neurophysiol. 61, 900–917 [PubMed]
- Averbeck B. B., Latham P. E., Pouget A. (2006). Neural correlations, population coding and computation. Nat. Rev. Neurosci. 7, 358–36610.1038/nrn1888 [PubMed] [Cross Ref]
- Averbeck B. B., Lee D. (2006). Effects of noise correlations on information encoding and decoding. J. Neurophysiol. 95, 3633–364410.1152/jn.00919.2005 [PubMed] [Cross Ref]
- Bair W., Zohary E., Newsome W. (2001). Correlated firing in macaque visual area MT: time scales and relationship to behavior. J. Neurosci. 21, 1676–1697 [PubMed]
- Bohte S. M., Spekreijse H., Roelfsema P. R. (2000). The effects of pair-wise and higher-order correlations on the firing rate of a postsynaptic neuron. Neural Comput. 12, 153–17910.1162/089976600300015934 [PubMed] [Cross Ref]
- Brette R. (2009). Generation of correlated spike trains. Neural Comput. 21, 188–21510.1162/neco.2009.12-07-657 [PubMed] [Cross Ref]
- Brown E. N., Kaas R. E., Mitra P. P. (2004). Multiple neural spike train data analysis: state-of-the-art and future challenges. Nat. Neurosci. 7, 456–46110.1038/nn1228 [PubMed] [Cross Ref]
- Darroch J., Speed T. (1983). Additive and multiplicative models and interactions. Ann. Stat. 11, 724–73810.1214/aos/1176346240 [Cross Ref]
- Del Prete V., Martignon L., Villa A. E. P. (2004). Detection of syntonies between multiple spike trains using a coarse-grain binarization of spike count distributions. Network 15, 13–2810.1088/0954-898X/15/1/002 [PubMed] [Cross Ref]
- Diesmann M., Gewaltig M.-O., Aertsen A. (1999). Stable propagation of synchronous spiking in cortical neural networks. Nature 402, 529–53310.1038/990101 [PubMed] [Cross Ref]
- Di Nardo E., Guarino G., Senato D. (2008). A unifying framework for k-statistics, polykays and their multivariate generalizations. Bernoulli 14, 440–46810.3150/07-BEJ6163 [Cross Ref]
- Ecker A. S., Berens P., Keliris G. A., Bethge M., Logothetis N. K., Tolias A. S. (2010). Decorrelated neuronal firing in cortical microcircuits. Science 327, 584–58710.1126/science.1179867 [PubMed] [Cross Ref]
- Eggermont J. J. (1990). The Correlative Brain, Vol. 16,
*Studies of Brain Function*Berlin: Springer-Verlag - Ehm W., Staude B., Rotter S. (2007). Decomposition of neuronal assembly activity via empirical de- poissonization. Electron. J. Stat. 1, 473–49510.1214/07-EJS095 [Cross Ref]
- Fujisawa S., Amarasingham A., Harrison M., Buzsaki G. (2008). Behavior-dependent short-term assembly dynamics in the medial prefrontal cortex. Nat. Neurosci. 11, 823–83310.1038/nn.2134 [PMC free article] [PubMed] [Cross Ref]
- Gardiner C. W. (2003). Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 3 Edn.
*Springer Series in Synergetics*, Vol. 13 Berlin: Springer - Gerstein G. L., Bedenbaugh P., Aertsen A. (1989). Neuronal assemblies. IEEE Trans. Biomed. Eng. 36, 4–1410.1109/10.16444 [PubMed] [Cross Ref]
- Gray C. M., Singer W. (1989). Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex. Proc. Natl. Acad. Sci. U.S.A. 86, 1698–170210.1073/pnas.86.5.1698 [PubMed] [Cross Ref]
- Grün S. (2009). Data-driven significance estimation of precise spike correlation. J. Neurophysiol. 101, 1126–1140 (invited review).10.1152/jn.00093.2008 [PubMed] [Cross Ref]
- Grün S., Diesmann M., Aertsen A. (2002a). ‘Unitary Events’ in multiple single-neuron spiking activity. I. Detection and significance. Neural Comput. 14, 43–8010.1162/089976602753284455 [PubMed] [Cross Ref]
- Grün S., Diesmann M., Aertsen A. (2002b). ‘Unitary Events’ in multiple single-neuron spiking activity. II. Non-Stationary data. Neural Comput. 14, 81–11910.1162/089976602753284464 [PubMed] [Cross Ref]
- Hebb D. O. (1949). The organization of behavior: a neuropsychological theory. New York: John Wiley & Sons
- Johnson D., Goodman I. (2007). Jointly Poisson processes. arXiv:0911.2524.
- Josić K., Shea-Brown E., Doiron B., de la Rocha J. (2009). Stimulus-dependent correlations and population codes. Neural Comput. 21, 2774–280410.1162/neco.2009.10-08-879 [PubMed] [Cross Ref]
- Kohn A., Smith M. A. (2005). Stimulus dependence of neuronal correlations in primary visual cortex of the Macaque. J. Neurosci. 25, 3661–367310.1523/JNEUROSCI.5106-04.2005 [PubMed] [Cross Ref]
- König P., Engel A. K., Singer W. (1996). Integrator or coincidence detector? The role of the cortical neuron revisited. Trends Neurosci. 19, 130–13710.1016/S0166-2236(96)80019-1 [PubMed] [Cross Ref]
- Kreiter A. K., Singer W. (1996). Stimulus-dependent synchronization of neuronal responses in the visual cortex of awake macaque monkey. J. Neurosci. 16, 2381–2396 [PubMed]
- Kuhn A., Aertsen A., Rotter S. (2003). Higher-order statistics of input ensembles and the response of simple model neurons. Neural Comput. 1, 67–10110.1162/089976603321043702 [PubMed] [Cross Ref]
- Kumar P. (2002). Moments inequalities of a random variable defined over a finite interval. J. Ineq. Pure Appl. Math. 3, 41
- Martignon L., Deco G., Laskey K., Diamond M., Freiwald W., Vaadia E. (2000). Neural coding: higher-order temporal patterns in the neurostatistics of cell assemblies. Neural Comput. 12, 2621–265310.1162/089976600300014872 [PubMed] [Cross Ref]
- Martignon L., von Hasseln H., Grün S., Aertsen A., Palm G. (1995). Detecting higher-order interactions among the spiking events in a group of neurons. Biol. Cybern. 73, 69–8110.1007/BF00199057 [PubMed] [Cross Ref]
- Montani F., Ince R. A. A., Senatore R., Arabzadeh E., Diamond M. E., Panzeri S. (2009). The impact of high-order interactions on the rate of synchronous discharge and information transmission in somatosensory cortex. Philos. Trans. R. Soc. Lond. A 367, 3297–331010.1098/rsta.2009.0082 [PubMed] [Cross Ref]
- Nakahara H., Amari S. (2002). Information-geometric measure for neural spikes. Neural Comput. 14, 2269–231610.1162/08997660260293238 [PubMed] [Cross Ref]
- Perkel D. H., Gerstein G. L., Moore G. P. (1967a). Neuronal spike trains and stochastic point processes. I. The single spike train. Biophys. J. 7, 391–41810.1016/S0006-3495(67)86596-2 [PubMed] [Cross Ref]
- Perkel D. H., Gerstein G. L., Moore G. P. (1967b). Neuronal spike trains and stochastic point processes. II. Simultaneous spike trains. Biophys. J. 7, 419–44010.1016/S0006-3495(67)86597-4 [PubMed] [Cross Ref]
- Pillow J. W., Shlens J., Paninski L., Sher A., Litke A. M., Chichilnisky E. J., Simoncelli E. P. (2008). Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature 454, 995–99910.1038/nature07140 [PMC free article] [PubMed] [Cross Ref]
- Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. (1992). Numerical Recipes in C, 2 Edn.New York, NY: Cambridge University Press
- Reimer I. C. G., Staude B., Rotter S. (2009). “Detecting assembly activity in massively parallel spike trains,” in Proceedings of the 8th Meeting of the German Neuroscience Society/30th Göttingen Neurobiology Conference, Vol. 1, Neuroforum, Supplement, eds Bähr H., Zerr I., editors. Heidelberg: Spektrum Akademischer Verlag
- Riehle A., Grün S., Diesmann M., Aertsen A. (1997). Spike synchronization and rate modulation differentially involved in motor cortical function. Science 278, 1950–195310.1126/science.278.5345.1950 [PubMed] [Cross Ref]
- Roudi Y., Nirenberg S., Latham P. E. (2009). Pairwise maximum entropy models for studying large biological systems: when they can work and when they can't. PLoS Comput. Biol. 5, e1000380.10.1371/journal.pcbi.1000380 [PMC free article] [PubMed] [Cross Ref]
- Sakurai Y., Takahashi S. (2006). Dynamic synchrony of firing in the monkey prefrontal cortex during working-memory tasks. J. Neurosci. 6, 10141–1015310.1523/JNEUROSCI.2423-06.2006 [PubMed] [Cross Ref]
- Schneidman E., Berry M. J., Segev R., Bialek W. (2006). Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440, 1007–101210.1038/nature04701 [PMC free article] [PubMed] [Cross Ref]
- Shadlen M. N., Newsome W. T. (1998). The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J. Neurosci. 18, 3870–3896 [PubMed]
- Shimazaki H., Amari S., Brown E. N., Grün S. (2009). “State-space analysis on time-varying correlations in parallel spike sequences,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP) (Washington, DC: IEEE Computer Society), 3501–3504
- Shlens J., Field G. D., Gauthier J. L., Grivich M. I., Petrusca D., Sher A., Litke A. M., Chichilnisky E. J. (2006). The structure of multi- neuron firing patterns in primate retina. J. Neurosci. 26, 8254–826610.1523/JNEUROSCI.1282-06.2006 [PubMed] [Cross Ref]
- Snyder D. L., Miller M. I. (1991). Random Point Processes in Time and Space. New York: Springer
- Softky W. R. (1995). Simple codes versus efficient codes. Curr. Opin. Neurobiol. 5, 239–247 (commentary).10.1016/0959-4388(95)80032-8 [PubMed] [Cross Ref]
- Staude B., Grün S., Rotter S. (2010). “Higher order correlations and cumulants,” in Analysis of Parallel Spike Trains, eds Grün S., Rotter S., editors. Springer Series in Computational Neuroscience. Berlin: Springer-Verlag
- Staude B., Rotter S., Grün S. (2007). “Detecting the existence of higher-order correlations in multiple single-unit spike trains,” in Abstract Viewer/Itinerary Planner, Vol. 103.9/AAA18 (Washington, DC: Society for Neuroscience; ).
- Staude B., Rotter S., Grün S. (2008). Can spike coordination be differentiated from rate covariation? Neural Comput. 20, 1973–199910.1162/neco.2008.06-07-550 [PubMed] [Cross Ref]
- Staude B., Rotter S., Grün S. (2009). CuBIC: cumulant based inference of higher-order correlations. J. Comput. Neurosci. 10.1007/s10827-009-0195-x [PMC free article] [PubMed] [Cross Ref]
- Stratonovich R. L. (1967). Topics in the Theory of Random Noise. New York, Gordon & Breach Science
- Streitberg B. (1990). Lancaster interactions revisited. Ann. Stat. 18, 1878–188510.1214/aos/1176347885 [Cross Ref]
- Stuart A., Ord J. K. (1987). Kendall's Advanced Theory of Statistics, 5 Edn.London: Griffin and Co
- Tsodyks M., Kenet T., Grinvald A., Arieli A. (1999). Linking spontaneous activity of single cortical neurons and the underlying functional architecture. Science 286, 1943–194610.1126/science.286.5446.1943 [PubMed] [Cross Ref]
- Vaadia E., Aertsen A., Nelken I. (1995). ‘Dynamics of neuronal interactions’ cannot be explained by ’neuronal transients’. Proc. Biol. Sci. 261, 407–41010.1098/rspb.1995.0167 [PubMed] [Cross Ref]
- Vaadia E., Haalman I., Abeles M., Bergman H., Prut Y., Slovin H., Aertsen A. (1995). Dynamics of neuronal interactions in monkey cortex in relation to behavioural events. Nature 373, 515–51810.1038/373515a0 [PubMed] [Cross Ref]
- Ventura V., Cai C., Kass R. E. (2005). Trial-to-trial variability and its effect on time-varying dependency between two neurons. J. Neurophysiol. 94, 2928–293910.1152/jn.00644.2004 [PubMed] [Cross Ref]

Articles from Frontiers in Computational Neuroscience are provided here courtesy of **Frontiers Media SA**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |