## Abstract

Correlations among neurons are supposed to play an important role in computation and information coding in the nervous system. Empirically, functional interactions between neurons are most commonly assessed by cross-correlation functions. Recent studies have suggested that pairwise correlations may indeed be sufficient to capture most of the information present in neural interactions. Many applications of correlation functions, however, implicitly tend to assume that the underlying processes are stationary. This assumption will usually fail for real neurons recorded in vivo since their activity during behavioral tasks is heavily influenced by stimulus-, movement-, or cognition-related processes as well as by more general processes like slow oscillations or changes in state of alertness. To address the problem of nonstationarity, we introduce a method for assessing stationarity empirically and then “slicing” spike trains into stationary segments according to the statistical definition of weak-sense stationarity. We examine pairwise Pearson cross-correlations (PCCs) under both stationary and nonstationary conditions and identify another source of covariance that can be differentiated from the covariance of the spike times and emerges as a consequence of residual nonstationarities after the slicing process: the covariance of the firing rates defined on each segment. Based on this, a correction of the PCC is introduced that accounts for the effect of segmentation. We probe these methods both on simulated data sets and on in vivo recordings from the prefrontal cortex of behaving rats. Rather than for removing nonstationarities, the present method may also be used for detecting significant events in spike trains.

- weak-sense stationarity
- neural coding
- in vivo electrophysiology
- neural dynamics
- information
- statistics

mounting evidence suggests that an important part of the code that neural populations entertain in processing information about sensory inputs and behavioral outputs is reflected in the structure of neural correlations, i.e., the tendency of neurons to generate temporally coordinated spiking patterns. Recent studies (Cohen and Maunsell 2009; Pillow et al. 2008; Schneidman et al. 2006) found that in a variety of brain areas pairwise interactions not only dominate over higher-order interactions but may actually be sufficient to account for the joint distributional properties in the spiking activity of populations of neurons. For pairs of neurons, functional interactions are most commonly assessed by cross-correlation functions that chart the spiking probability (or frequency) in one neuron as a function of the time lag to the spiking times of the other neuron. In particular, precise phase-lag-zero synchronous firing has been hypothesized to be of special computational relevance, for instance, as a solution to the problem of binding multiple stimulus features into a single percept (Gray 1999; Singer and Gray 1995). Synchronous spiking has been proposed to enhance the strength and reliability of the response of postsynaptic neurons (Abeles et al. 1995; Diesmann et al. 1999; Herrmann et al. 1995) and to reflect functional subunits with strong interconnections in a network, so-called cell assemblies (Abeles 1991; Hebb 1949). More generally, correlated neural activity (at zero or nonzero phase lags) has been linked to diverse computational functions such as encoding and discrimination of stimulus features (Averbeck et al. 2006; Decharms and Merzenich 1996; Stopfer et al. 1997), orientation tuning and contrast sensitivity (Kohn and Smith 2005), recognition of visual objects (Hirabayashi and Miyashita 2005), encoding of places and spatial trajectories in hippocampus (Harris 2005), time perception (Hass and Herrmann 2012; Haß et al. 2008), motor behavior (Riehle et al. 1997; Vaadia et al. 1995), and attention (Cohen and Maunsell 2009; Steinmetz et al. 2000). On the other hand, correlations have also been considered to be a limiting factor in neural information processing as they introduce redundancy (Abbott and Dayan 1999; Averbeck et al. 2006; Zohary et al. 1994) and have been found by some studies to be attenuated by attention (Cohen and Maunsell 2009). In addition to their presumed role in neural coding, cross-correlations are often used to estimate synaptic connectivity among neurons (Barthó et al. 2004; Binder and Powers 2001; Constantinidis et al. 2001; Fujisawa et al. 2008; Sears and Stagg 1976; Snider et al. 1998; Türker and Powers 2001). Thus they convey important information about the structural and functional architecture of neural networks, although extracting this information may be a difficult problem in its own right (Alonso and Martinez 1998; Melssen and Epping 1987; Trong and Rieke 2008), and a number of recent computational studies have attempted to specify exactly how correlations depend on synaptic connectivity and other factors (Kriener et al. 2008; Litwin-Kumar et al. 2011; Ostojic et al. 2009; Renart et al. 2010; Tchumatchenko et al. 2011).

Most applications of cross-correlation functions and related measures implicitly assume that the observed system is stationary, that is, that the stochastic properties (statistical moments) of the neurons' spiking activities stay constant in time. Violations of this assumption can lead to serious artifacts in measured cross-correlations or their interpretation (Baker and Gerstein 2001; Ben-Shaul et al. 2001; Brody 1999a, 1999b; Grün 2009; Grün et al. 2002a, 2003). For instance, it has been shown that the strength of correlations between two neurons is related to their firing rates when nonlinearities are present in the input-output properties of neurons (de la Rocha et al. 2007). Nonstationarities are common in real neurons in which cofluctuations in the responses of a pair of neurons can arise over a range of timescales (Cohen and Kohn 2011). Variations that occur at short (<1 s) and medium (1–10^{2} s) timescales are usually related to stimulus conditions, motor events, or attentional or motivational modulations, while phenomena like synaptic plasticity and learning are typically associated with longer-term variations (>10^{2} s) (Aertsen et al. 1989). Such fluctuations modulate or interfere with the pairwise spike interactions among cells, making correlations difficult to interpret. Moreover, it is almost impossible to determine the contribution that different internally generated nonstationarities will have on correlations since these cannot be or can only partially be experimentally controlled or assessed. For these reasons, special care must be taken when analyzing or interpreting spike interaction data from behaving animals, in which the spiking activity of neurons often represents a wide variety of stimulus-, movement-, cognition-, motivation-, or emotion-related processes.

A number of methods have been proposed to correct for the effects of nonstationarities. Most of these methods focus on nonstationarities across trials, such as changes in excitability of the neurons (Ben-Shaul et al. 2001; Brody 1999b; Grün et al. 2002a, 2003; Ventura et al. 2005) or the latency in response to a stimulus (Baker and Gerstein 2001; Grün et al. 2003; Ventura et al. 2005). However, as outlined above, the type and origin of nonstationarities can be diverse and hard to determine. Furthermore, there are situations in which collapsing information across trials is not possible or sensible. In sensory domains, comparability among trials is usually assumed when presenting identical stimuli. However, this option is not available or is less applicable when investigating network activity without stimuli (Beggs and Plenz 2003; Ikegaya et al. 2004) or within “higher-order” cognitive areas such as the prefrontal cortex (Lapish et al. 2008) where the neural dynamics are largely internally driven and inputs can be less well directly experimentally controlled. Moreover, single-trial analysis of brain activity can convey information that would be lost by combining trials (Durstewitz et al. 2010). Therefore, it is desirable to have a method for filtering experimental spike trains, within or across trials, such that they become stationary, independent of the source and nature of the nonstationarities, and a way to correctly compute the cross-correlation from these trains.

In this report, we address this problem by first introducing a method for assessing nonstationarity of spike trains empirically by comparing running local interspike interval-based statistics to distributional assumptions. The method then “slices” spike trains into segments that collectively satisfy (approximately) the statistical definition of weak-sense stationarity. The spike train slicing process itself introduces another source of covariance that can be differentiated from the covariance of the spike times and is due to residual nonstationarities, namely, the remaining firing rate covariations across the sampled segments. Removing the corresponding terms from the Pearson cross-correlation (PCC) yields the stationarity-corrected PCC (scPCC). We probe these methods on simulated data sets that mimic spike trains from neurons perturbed by two different types of nonstationarity, namely, brief and randomly distributed comodulations of firing rates by external stimuli and slow comodulation of firing rates by a harmonic oscillation. We also probe our methods on multiple single-unit in vivo recordings from the prefrontal cortex of behaving rats and compare the results with the simulated cases.

## MATERIALS AND METHODS

#### Defining stationary conditions.

Given a spike train {*t*_{i}}, *i* = 1...*N* + 1, with interspike intervals ISI([*t*_{i+1} + *t*_{i}]/2) = ISI_{i} = *t*_{i+1} − *t*_{i}, *i* = 1...*N*, weak-sense stationarity requires that the expectancy *E*[ISI] and the variance of the ISI series do not change with time (i.e., are constant) and that the autocovariance depends only on the time lag τ but not on time itself, i.e., Cov[ISI(*t*),ISI(*t* + τ)] = Cov(τ) for all *t* (Gilgen 2006). For simplicity and practical reasons, we only impose the first two conditions in our method (autocorrelations in ISI series in our experience decay away very quickly anyway, after a few ISIs). Let
*k*, centered halfway between the last and the first spike of the sliding window. To assess whether the ISI series is, in the weak sense, stationary in time, we construct the variables
_{i} were all independently and normally distributed random variables, assuming stationarity and *N* >> *k* (>10-fold), *T*_{m},_{k} for each (fixed) *m* and *k* would be approximately normally distributed with mean 0 and variance 1, while *Q*_{m},_{k} is approximately χ^{2}-distributed with *k* degrees of freedom (df) (alternatively, *T*_{m,k} may be based on the difference between the moving window mean μ̂_{m,k} and the grand mean excluding the moving window, with pooled variance, giving rise to a *t*-distributed quantity with df = *N* − 2; similarly, *Q*_{m,k} may be replaced by an *F*-distributed quantity relating “in-window” and “out-of-window” variation). On the basis of these variables, we regard the collection of segments {*m*} of length *k* of the ISI series as stationary for which *T*_{m,k} falls into the [*T*_{min},*T*_{max}] confidence interval of the normal distribution and at the same time *Q*_{m,k} falls into the [*Q*_{min},*Q*_{max}] confidence interval of the χ^{2}-distribution. Note that according to the central limit theorem, for large enough *k* (given *N* >> *k*) the distributional assumptions about *T*_{m},_{k} and *Q*_{m},_{k} would approximately hold even if the ISI_{i} themselves were not normally (but still independently) distributed. Nevertheless, as ISIs are usually not normally distributed, we transform the ISI distributions before applying the above procedure (Box and Cox 1964),
*q* ∈ [0,1] that is determined through a maximum-likelihood criterion (Box and Cox 1964). Finally, we standardize the ξ distribution (such that

#### Pearson cross-correlation.

The Pearson cross-correlation PCC_{xy}(Δ*i*) of two binary random variables *x* and *y* separated by a time lag of Δ*i* is given by
*p _{xy}* =

*pr*{

*x*(

*i*) = 1 ∧

*y*(

*i*+ Δ

*i*) = 1} is the joint probability and

*p*

_{x}=

*pr*{

*x*= 1} and

*p*

_{y}=

*pr*{

*y*= 1} are the marginal probabilities associated with each of the random variables. In the following we assume that the two spike trains

*x*and

*y*have already been shifted by lag Δ

*i*with respect to each other (such that Δ

*i*can be dropped for the derivations below). After shifting, the simultaneously stationary time intervals at that lag are determined as illustrated in Fig. 1. In this example we have two jointly stationary segments from each spike train with spiking probabilities

*p*

_{x1}and

*p*

_{x2}for neuron

*x*and

*p*

_{y1}and

*p*

_{y2}for neuron

*y*, respectively, where we estimate these probabilities as the number (#) of spikes divided by the number of bins (total of

*n*

_{1}+

*n*

_{2}=

*N*bins). The total probability

*p*

_{x}(and likewise

*p*

_{y}and

*p*

_{xy}) can now be expressed as

*x*

_{l}refers to random variable

*x*within segment

*l*. More generally, for

*L*segments of length

*n*

_{l}one could write

*f*(

*x*

_{l}) and

*f*(

*y*

_{l}) the firing rates of processes

*x*and

*y*in segment

*l*and

*b*the bin width (assuming a maximum of 1 spike can occur per bin). Hence, the total covariance splits up into a weighted sum of the covariances across the

*L*segments and a weighted (by relative segment length) form of the covariance between segment-specific firing rates. The variances of the two spiking processes can be split up in a similar fashion, leading to the following expression for the PCC:

*i*) indicates that process

*y*had been shifted by Δ

*i*bins with respect to process

*x*. In the limit, if the segmented processes were truly stationary, Var(

*p*

_{xl}) → 0, Var(

*p*

_{yl}) → 0, and Cov(

*p*

_{xl},

*p*

_{yl}) → 0. Hence, these terms reflecting “cross talk” among segments appear as a consequence of nonstationarities left after the segmentation process, and we drop them from the final definition to have a stationarity-corrected Pearson cross-correlation (scPCC) that captures solely precise spike time coincidences:

*p*

_{xl},

*p*

_{yl}) or Var(

*p*

_{xl}), Var(

*p*

_{yl}) were significant in size, this would imply that there were considerable nonstationarities left, and either the confidence limits will have to be tightened or

*k*has to be increased to diminish variation among segments. Specifically, we suggest

*p*

_{xl}) and Var(

*p*

_{yl}) have to be weighted by relative segment length

*w*

_{l}]. Also note that if the processes were truly stationary, Var(

*x*

_{l}) and Var(

*y*

_{l}

^{(Δi)}) would be constants and

*Eq. 9*could be simplified further.

Finally, in order to have a robust estimate of the cross-correlation we apply a smoothing spline (Garcia 2010), based on cosine functions and penalized least-squares estimation, to the scPCC:
**S** is the spline operator (which combines values across lags for the smooth function estimate).

## RESULTS

We proposed a method for assessing and enforcing weak stationarity empirically by restricting the spike train data to those segments for which the local estimate of the mean transformed ISI is confined within the [*T*_{min},*T*_{max}] confidence interval of the normal distribution and the local second moments within the [*Q*_{min},*Q*_{max}] confidence interval of the χ^{2}-distribution (see materials and methods). For these stationarity-segmented spike train data we then derived an expression for the PCC that is largely free from nonstationarity artifacts, within the confidence bounds defined, as demonstrated in the following. Before investigating its application to different simulated and physiological data sets, the next section first demonstrates the basic operation of the stationarity-enforcing segmentation technique.

#### Segmenting spike trains into stationary periods.

Spike trains were generated by a stationary Poisson process at a rate of 5 spikes/s, yielding an exponential ISI distribution (Fig. 2*A*). As shown in Fig. 2*B*, an exponent *q* ∼ 0.2 yields the highest likelihood for the Box-Cox transform, converting the exponential ISI distribution and thus quantity *T*_{m,k} (*Eq. 3*) approximately into a Gaussian (Fig. 2*C*), while *Q*_{m},_{k} (*Eq. 3*) will be approximately χ^{2}-distributed (Fig. 2*D*). A sliding window of 10 ISIs was chosen to obtain local estimates of first and second moments of the standardized transformed ISI series (see materials and methods). Collectively, stationary segments are now defined as those periods of the transformed ISI series that have both *T*_{m},_{k} (Fig. 3*A*) and *Q*_{m},_{k} (Fig. 3*B*) statistics ranging within the [2–98%] confidence limits of the normal and χ^{2}-distributions, respectively. These segments will be cut out from the original ISI series for further processing, as illustrated in Fig. 3*C*. Of course, the final stationary process is not continuous in time anymore since we have eliminated the nonstationary intervals from it, leaving a set of disconnected stationary segments.

#### Pulse distortion.

In this section we probe our method for determining cross-correlations among pairs of neurons from stationary segments (stationarity-corrected Pearson cross-correlation, scPCC) on simulated pairs of independent Poisson spiking processes at a rate of 5 Hz, jointly perturbed by randomly distributed (at a rate of 0.1 Hz) “stimulus” pulses. These pulses consist of joint increases in the spiking rate to 50 Hz for periods of 200 ms in both “neurons” (Poisson processes) simultaneously (Fig. 4*A*). We then apply our stationarity criterion, which results in the removal of periods from the spike trains for which one of the statistics *T*_{m},_{k} or *Q*_{m},_{k} falls outside of a 96% confidence interval ([0.02, 0.98]), thereby eliminating exactly those periods where the pulse perturbations occurred, as shown in Fig. 4*B* (for these confidence settings and *k* = 10, one obtains η ≈ 3.5 × 10^{−3} << 1, indicating that the remaining nonstationarity compared with the variation in the spiking processes is negligible).

Figure 4*C* gives one example of a PCC function for a pair of spike trains, both for the original spike trains and with stationarity enforced. As described in materials and methods, a spline smoothing was applied to the scPCC, in order to obtain statistically more reliable estimates. Figure 4*D* gives the average PCC functions (across *n* = 100 pairs) for these simulated spike train data both before and after enforcing stationarity. While the cross-correlation function for the raw, nonstationary spike trains peaks at ∼0.033 (well within the range of empirically observed significant cross-correlations; see, e.g., Cohen and Kohn 2011) and exhibits a triangular shape that falls off to both sides within the 200-ms window caused by the pulse distortion, the stationarity-corrected cross-correlation function is completely flat, as it should be for two unrelated Poisson processes.

Note that our stationarity method has two free parameters, the width of the confidence interval for statistics *T*_{m},_{k} and *Q*_{m},_{k}, which controls the “size” of the nonstationarities tolerated (e.g., the pulse amplitude in the present example), and the number *k* of ISIs used to determine *T*_{m,k} and *Q*_{m,k}. Parameter *k* controls a kind of “bias-variance trade-off” (see, e.g., Hastie et al. 2009): The larger *k*, the more robust (less variable) the estimates of *T*_{m,k} and *Q*_{m,k} will become, but the lower the temporal precision will be (e.g., with regard to pulse duration in the present example; in this sense introducing more of a “stationarity bias”). That is, with larger *k* briefer nonstationarities may go unnoticed (Fig. 4, *E* and *F*). This ambiguity in the size of the window has also been brought up by Grün et al. (2002a). Ideally, the size of the window (and in our case also the confidence limits) has to be adjusted to the timescale of the nonstationarities, which in principle is unknown. Therefore, we recommend analyzing the data with a few different combinations of settings for *k* and the width of the confidence band, such that there is convergence in the scPCC estimate while η is sufficiently small (e.g., η < α, assuming a [α/2, 1 − α/2] confidence band).

#### Slow oscillations with synchronized spike times.

In this second example for the application of the scPCC method, pairs of initially independent Poisson spiking processes are generated, with a common slow harmonic oscillatory modulation of firing rates, *v*(*t*) = *B* + *C*cos(2π*ft* + φ) (Staude et al. 2010b), with *B* the mean rate, *C* the amplitude, and *f* the frequency of the modulation (Fig. 5*A*; here we chose *B* = 27.5, *C* = 22.5, *f* = 0.5 Hz). In addition, 500 randomly selected spike times (representing <5% of all spike times) were synchronized among pairs with a phase lag of 25 ms and temporal jitter of σ = 2 ms. Spike trains were dissected with confidence bounds set to [0.02, 0.98] and *k* = 10 (yielding η ≈ 2.5 × 10^{−2}) as illustrated in Fig. 5*B*. (Note that if phase φ is treated as a random variable and observation time is long enough, a harmonic periodic process is stationary with *E*[ν] = *B* in the case above, independent of *t*. Practically, as noted above, the relevant timescale for treating a process as nonstationary is defined here by the window size *k*.)

Figure 5*C* gives one example of a PCC function for a pair of spike trains, both for the original spike trains and with stationarity enforced. As described in materials and methods, spline smoothing was applied to the scPCC. Figure 5*D* presents results averaged across a sample of 100 independent cross-correlation pairs. As evident from Fig. 5, *C* and *D*, enforcing stationarity removes the slow oscillatory component from the cross-correlation functions while leaving the precise (within 2 ms) spike coincidences at lag ± 25 ms intact. Thus the method successfully separates the temporally precise spiking patterns (which will always be preserved regardless of the parameter settings of the method) from the slow comodulation of firing rates. To further confirm the validity of the method, bootstrap data sets were constructed by permuting blocks of 10 consecutive ISIs (retaining autocorrelations up to that length, as induced in this example by the oscillation). As evident in Fig. 5*D*, the precise spike patterns are gone within these bootstrap data.

#### Application to experimental data.

In this final section the stationarity detection and scPCC methods are illustrated on an experimental data set consisting of in vivo multiple single-unit recordings from the prefrontal cortex of behaving rats (Balaguer-Ballester et al. 2011; Lapish et al. 2008). These data contain a number of stimulus-, memory-, and behavior-related nonstationarities. Figure 6 gives an example for the translation of an experimental ISI series and distribution into corresponding test statistics and their empirical distributions (analog to the theoretical example in Fig. 2). Although transformed distributions will generally be approximately normal (as in the example illustrated), this will not always be the case. Larger sliding window sizes *k* may have to be chosen for ISI distributions that violate normality too much (but see discussion), as Gaussianity is ensured at least in the limit of large *k* (and *N*) because of the central limit theorem (assuming approximate independence of ISIs). Figure 7 illustrates the slicing of an experimentally obtained spike train according to the stationarity criteria.

The whole experimental data set analyzed here consists of recordings with 6 to 39 simultaneously recorded neurons from 12 different animals with 1 to 6 trials from each. Figure 8*A* gives some example spike trains. We fixed a 96% confidence interval ([0.02, 0.98]) for these analyses and computed the PCC on the original (not stationarity corrected) data sets as well as scPCC for sliding window sizes of *k* = 10 (Fig. 8*B*; η ≈ 8.5 × 10^{−3}) and *k* = 5 (Fig. 8*D*; η ≈ 1.05 × 10^{−2}). For control, block permutation bootstraps with a block length of 10 ISIs, as introduced in the previous section, were processed in the very same way as the original data. As Fig. 8, *B* and *D*, illustrate, the stationarity correction removes a kind of slowly decaying background present in the original data, potentially reflecting a very slow oscillation. Side peaks at ±25 ms are visible in the grand average across all pairs for *k* = 5. However, if cross-correlation pairs are split into those with a positive value at zero lag and those with a negative value, it becomes apparent that these two humps or very flat flanks in the total average were caused by the overlay of wider positive and narrower negative peak correlations (Fig. 8*C*).

Figure 8, *B–D*, illustrate that our method of stationarity correction may actually bring out true phase-coupling or temporally precise spiking patterns more clearly, once a “background of noise” due to nonstationarity has been removed. They also show that the cross-correlation peaks hardly change in size when window size *k* is reduced from 10 to 5, further confirming that they reflect true precise phase relationships in the spike data, rather than being an artifact of nonstationarity.

## DISCUSSION

Nonstationarities in neuronal data can severely complicate or compromise the interpretation of statistical measures such as cross-correlations. We therefore devised a method for removing nonstationarities by “slicing” the data into intervals that together approximately satisfy the statistical definition of weak-sense stationarity. These data segments can then be used to compute statistical measures that are sensitive to violations of stationarity or to extract a form of “unperturbed” baseline activity from in vivo data.

Here we focused on the PCC between spike trains. We showed that for the stationarity-segmented data the correlation splits up into a part that reflects precise spike time correlations within the segments and a part that reflects covariation of firing rates across segments. Thus if the interest is in precise spike time coincidences, the part due to firing rate (co)variations should be neglected if it is sufficiently small. If it is not, this would indicate that stationarity conditions were not enforced rigorously enough and tighter confidence limits or larger moving window sizes *k* should be used. A similar partitioning of the cross-correlation has been introduced by Staude et al. (2008), although with a slightly different focus.

In two simulations, we mimicked typical nonstationarities that can be found in real experiments, namely, distortions by brief pulses as they may arise from noncontrolled sensory inputs or internal events and slow oscillations as dominant especially during awake resting activity in vivo (Steriade et al. 1993). Both types of nonstationarities produced a pronounced background of cross-correlation that could easily obscure peaks due to precise spike time coordination, and our method reliably removed these effects. In the case of slow oscillations, these were completely filtered after stationary constraints were enforced, while precise spike time correlations among the neurons remained. Similar results were obtained with experimental recordings, revealing pronounced peaks in the cross-correlation.

#### Removing nonstationarities from cross-correlations and other applications.

Several studies have shown that the cross-correlation is particularly sensitive to violations of stationarity (Baker and Gerstein 2001; Ben-Shaul et al. 2001; Brody 1999a, 1999b; Grün 2009; Grün et al. 2002a, 2002b, 2003). Computing cross-correlations often requires averaging over a number of identically prepared trials. Consequently, the detection and removal of nonstationarities across trials has gained most attention so far, particularly those that are due to covariations of the firing rate of both neurons across trials (Ben-Shaul et al. 2001; Brody 1999b; Grün et al. 2002a, 2003; Ventura et al. 2005) and to the latency in response to a stimulus (Baker and Gerstein 2001; Grün et al. 2003; Ventura et al. 2005). All of these methods focus on a specific source of nonstationarity and often require considerable effort in terms of calibration and computational demand. In contrast, the present method is based on the definition of weak-sense stationarity itself and thus removes all kinds of nonstationarities (but see below) by means of a relatively simple procedure in which only two parameters, the length of the sliding window and the width of the confidence interval, need to be specified. We explicitly tested the procedure on two types of nonstationarities that are likely to occur in real data, short pulse distortions and slow oscillations. Both of these may occur within a single trial and in that case would not be addressed by the methods mentioned above (see Grün et al. 2003). Thus our method addresses the issue of nonstationarity more generally and would remove them either within or across trials.

Using simulated data, we have shown that our method completely filters out slow oscillations in the cross-correlogram, while precise spike time correlations among the neurons remained. Likewise in the analyzed experimental recordings, an elevated correlation baseline was filtered out by the stationarity correction method, leaving clear narrow peaks centered around zero when separating pairs with positive versus negative contributions. Our method suggests that these peaks reflect real precise spike time synchrony, rather than being produced by nonstationarity-related artifacts. Subtracting off pairs with negative from those with positive correlations tended to leave two side peaks at approximately ±25 ms. Whether these reflect true side peaks obscured by the center peak, potentially indicating interactions within the gamma frequency range (30–60 Hz), or are just an artifact of averaging remains an interesting question to explore. It is important to point out that small and hard to detect pairwise correlations can still have large effects at the population level (Zohary et al. 1994).

Note that our method could also be used to detect interesting nonstationarities, rather than just filtering them out. In certain behavioral or physiological contexts exactly these nonstationary events may be most interesting. The methods introduced here could reveal such events in an unsupervised fashion (i.e., independent from knowledge about stimulus conditions for instance) and provide a measure of their significance through the definition of the confidence bands (and possible correction for multiple testing).

#### Relation to other methods of spike train segmentation.

The approach of segmenting a time series of spike times and computing correlations for the individual segments has also been used recently by Nikolić et al. (2012), mainly to remove correlations induced by slow oscillations. To attenuate slower components, the authors fixed the length of the segments according to the longest timescale to be considered while retaining correlations induced by fast components, but the issue of stationarity itself was not addressed. Thus intervals with strong firing rate covariations that do not follow a (regular) oscillatory cycle would not be detected by that method, unlike the method presented here, which will remove all sources of nonstationarity (but see below). Furthermore, with the present method the length of the segments is not preset but generically defined through the stationarity conditions, i.e., the statistical features of the time series itself define the appropriate timescale for spike time correlations.

#### Limitations and extensions of the present method.

Our stationarity-segmentation method works from the assumption that there is one underlying process with constant mean and variance, perturbed by comparably brief nonstationary excursions. A different situation arises when there are systematic changes on a much slower timescale, e.g., a constant drift in the mean rate, or switches between different mean rates that are each maintained for considerable and comparable periods of time (seconds). Such scenarios may be indicated by significant deviations of the Box-Cox-transformed ISI distributions from Gaussians, e.g., by multimodality in the distribution if window size *k* is chosen large enough (large *k* in this case would average out local fluctuations, acting as a low-pass filter, and thus bring out features due to long-term nonstationarity more clearly). Under these conditions, the present method, without further modifications, may not yield optimal or useful results. “Classical” remedies like detrending the time series, available also for point processes (see, e.g., Pruscha 1994), or adjusting for changing mean ISIs, imply modifying the spiking times in a way that may corrupt precise spike time coordination between spike trains as well, and may thus not be advisable. Rather, spike trains may have to be first dissected by eye at obvious change points, in the case of different stable mean rates or variances, or, in the case of a systematic trend, may have to be cut into sufficiently short segments of roughly constant mean rate (and variance). Explicit stationarity-segmentation may then be run on each of these pieces separately, and either all segments may be recombined in the calculation of the scPCC or separate scPCCs may be computed if it is suspected that the changes in the underlying processes also affect their cross-correlation structure. Alternatively, although not explicitly designed for this purpose, the present segmentation method may also be run iteratively on the time series using large *k* and strict confidence bounds, i.e., by reapplying it to the set of cut-out segments, until a satisfactory dissection of the whole process has been obtained. In the case of multimodality, another possibility is to extract the different means first from the *T*_{m,k} distribution and perform different segmentation runs around each of these different means.

#### Application to other statistical measures.

In the present work, we focused on pairwise cross-correlations between spike trains, but of course the present method of enforcing stationarity can be used for computing any other statistical measure that may be confounded by nonstationarities. In particular, higher-order correlations between three or more spike trains has recently gained a lot of interest (Balaguer-Ballester et al. 2011; Ohiorhenuan et al. 2010; Onken et al. 2012; Pipa and Munk 2011; Staude et al. 2010a), although there are also studies that suggest that pairwise correlations may not only dominate over higher-order interactions but may actually be sufficient to account for the joint distributional properties in the spiking activity of populations of neurons (Cohen and Maunsell 2009; Pillow et al. 2008; Schneidman et al. 2006). Generally, our method may help to extract many other spike train statistics free from nonstationarity confounds, like the coefficient of variation (*C*_{V}) or other measures of spike train irregularity (e.g., Shinomoto et al. 2003), measures of intrinsic bursting propensity of neurons (e.g., Gourévitch and Eggermont 2007), or simply estimates of the spontaneous firing rate or ISI distribution not confounded by stimulus conditions or other strongly time-varying components. Also, many time series models and frequency domain methods would usually assume stationarity across the studied periods (e.g., Fan and Yao 2003).

Another method that is frequently utilized to analyze the interdependence of neurons is unitary event analysis (Grün et al. 2002a, 2002b). Similar to the cross-correlogram, it counts the co-occurrences of spikes in several neurons, but—at least practically—is limited to a time lag of zero. This count is then compared to the number of coincidences that would be expected from independent Poisson processes with the rate of the actual neurons, and an event is labeled as unitary if the actual number of coincidences significantly exceeds the predicted number. In its original form, unitary event analysis strongly depends on stationarity of the neurons' spike trains. It has been extended to nonstationary spike trains (Grün et al. 2002a) by considering short sliding windows within which stationarity can be assumed to hold at least approximately and including firing rate nonstationarities in the Poisson predictor. The present method may contribute to unitary event analysis by enforcing stationarity as defined here within all recorded spike trains and then applying unitary event analysis on the segmented trains using a single stationary firing rate for each neuron.

## GRANTS

This work was funded by grants from the German ministry of education and research through the Bernstein Center for Computational Neuroscience (BMBF, 01GQ1003B) and the Deutsche Forschungsgemeinschaft to D. Durstewitz (Du354/6-1 and 7-2).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: C.Q.L. and D.D. conception and design of research; C.Q.L. analyzed data; C.Q.L., J.H., and D.D. interpreted results of experiments; C.Q.L. prepared figures; C.Q.L., J.H., and D.D. drafted manuscript; C.Q.L., J.H., and D.D. edited and revised manuscript; C.Q.L., J.H., and D.D. approved final version of manuscript.

## ACKNOWLEDGMENTS

We are grateful to Drs. Chris Lapish and Jeremy Seamans for providing the in vivo data sets for probing our methods.

- Copyright © 2013 the American Physiological Society