Emerging evidence indicates that information processing, as well as learning and memory processes, in both the network and single-neuron levels are highly dependent on the correlation structure of multiple spike trains. Contemporary experimental as well as theoretical studies that involve quasi-realistic neuronal stimulation thus require a method for controlling spike train correlations. This letter introduces a general new strategy for generating multiple spike trains with exactly controlled mean firing rates and correlation structure (defined in terms of auto- and cross-correlation functions). Our approach nonlinearly transforms random gaussian-distributed processes with a predistorted correlation structure into nonnegative rate processes, which are then used to generate doubly stochastic Poisson point processes with the required correlation structure. We show how this approach can be used to generate stationary or nonstationary spike trains from small or large groups of neurons with diverse auto- and cross-correlation structures. We analyze and derive analytical formulas for the high-order correlation structure of generated spike trains and discuss the limitations of this approach.
The principal representation of neural information is in sequences of action potentials (APs). Within the mammalian central nervous system, each neuron receives thousands of such input sequences (spike trains) from presynaptic cells. A neuron's output is believed to be dependent on not only the mean firing rate of presynaptic cells but also the degree of spatial and temporal correlation between spike trains from different presynaptic neurons (for reviews, see, e.g., Abeles, 1991; Mel, 1994; Segev & London, 2000; Hausser & Mel, 2003; London & Hausser, 2005). Furthermore, synaptic plasticity processes are also known to be highly dependent on the correlation structure between pre- and postsynaptic neurons, a process termed spike-timing-dependent plasticity (STDP; Dayan & Abbott, 2001). The multicorrelation structure of spike trains is also believed to be crucially important to the representation and processing of network information (Abeles, 1991), where even weak pair-wise correlations can lead to highly organized network states (Schneidman, Berry, Segev, & Bialek, 2006).
Methodological advances in recent years led to an experimental ability to rapidly stimulate dendrites with a high temporal resolution and at a resolution level approaching a single synapse (Matsuzaki et al., 2001; Shoham, O'Connor, Sarkisov, & Wang, 2005; Araya, Eisenthal, & Yuste, 2006; Losonczy & Magee, 2006). Similar optical tools can also be used to individually stimulate groups of presynaptic neurons (Boucsein, Nawrot, Rotter, Aertsen, & Heck, 2005; Shoham, O'Connor et al., 2005; Farah, Reutsky, & Shoham, 2007; Nikolenko, Poskanzer, & Yuste, 2007). These tools enable the study of responses to controlled yet realistically complex input patterns. Of particular interest among the large space of possible stimulation pattern is the study of responses to varying degrees of correlations between spike trains. This, however, requires a standard method for generating spike trains with a specific degree of correlation between them in order to conduct systematic studies.
The central measures of correlations among and between random processes are the auto- and cross-correlation functions (estimated by the auto- and cross-correlograms). Spike train correlation functions describe not only the coincidence of spikes, but also probabilities for spikes to appear at different relative delays, whether in two different sequences (cross-correlation function) or in the same one (autocorrelation function). Currently, there is no method that enables spike train generation with these functions predefined. Recent work by Niebur (2007) proposes a partial solution that allows the generation of spike trains with a specific degree of coincident spikes (with predefined time domain binning) between different sequences. Here we present a new approach to spike train generation (stationary or nonstationary) in continuous time (no prebinning of the time domain) with predefined auto- and cross-correlation functions. Our approach is based on analytical solutions of the nonlinear distortion of correlation structure in a linear-nonlinear-Poisson (LNP) generative model with gaussian white noise inputs.
We consider the problem of generating spike trains as doubly stochastic Poisson point processes (Snyder & Miller, 1991; also known as Cox processes: Cox, 1955; Cox & Lewis, 1966) with stochastic intensity functions (instantaneous rates) having specific auto- and cross-correlation functions. A Poisson process is called doubly stochastic if its rate profile is itself a random process. Interestingly, point processes generated this way preserve the full correlation structure of the corresponding stochastic intensity processes (see Knox, 1970 and proof below). This important property reduces the problem of generating correlated point processes to the problem of generating correlated continuous stochastic nonnegative processes. We next reduce the complexity of the problem by generating these nonnegative processes through a nonlinear transformation applied to correlated gaussian processes, an approach reviewed in Johnson (1994). This transformation affects not only the marginal distributions of the processes but also their correlation functions, so we create correlated gaussian processes whose correlation structure is appropriately predistorted to yield the desired final correlations in the nonnegative processes. Table 1 and Figure 1 outline the full algorithm.
2.2. Correlation Distortions.
Our method is based on predistorting the correlation structure of random gaussian processes in an exact way such that a nonnegative transformation applied to them will yield the desired auto- and cross-correlation functions (Johnson, 1994). We derive below the distortions resulting from three basic nonnegative transformations: exponential, square, and absolute value. We begin by restricting the discussion to correlation coefficients between draws of random variables Xi and Λi. In section 2.4, we expand this treatment to auto- and cross-correlation functions between stochastic processes xi(t) and λi(t) by appropriately arranging and filtering uncorrelated processes.
2.2.1. Exponential Transformation.
2.2.2. Square Transformation.
2.2.3. Absolute Value Transformation.
In the appendix, we derive a numerical solution for this transformation, which is inherently computationally more expensive compared to the analytical solutions of the two other transformations above.
2.3. Doubly Stochastic Poisson Processes Possess the Same Correlations as Their Stochastic Intensity Functions.
The instantaneous cross-correlation function RPij(t, τ) between two doubly stochastic Poisson processes is equal to the cross-correlation function between their respective intensity functions Rij(t, τ). This property also holds for the autocorrelation function of a single point process, except for τ = 0, where it behaves as a delta function.
2.4. Generating the Correlated Processes.
In practice, is highly redundant as a consequence of the replication procedure applied in equations 2.19 and 2.20, and we select only N of its rows (the central row in each block) to perform the filtering. In the nonstationary case, the matrices Wi and rij do not have a Toeplitz structure, and the matrix is not redundant in this case.
We note that in addition to this FIR-based approach, an alterative method for constructing multivariate correlated gaussian processes using multivariate autoregressive models is also possible here.
Once the correlated gaussian processes have been generated, they are nonlinearly transformed into rate processes using one of the transformations discussed above. Generating the Poisson spike trains from these inhomogeneous rate processes λ(t) proceeds using the standard procedure described by Johnson (1996, eq. 7). In this procedure, a sequence of independent exponentially distributed random variables ζi ∼ exp(1) is generated, and then the event times tk are calculated by solving .
Using our approach, it is possible to generate spike trains with different correlation structures that are typical for real neural activity (see Figures 2C–2F for examples). Cross-correlation function between two signals changes while propagating through the steps of the algorithm and corresponds to its theoretical values at every step (see Figures 2A–2C).
The generated point processes have intrinsically much more variability than the underlying intensity processes. This can be seen reflected in the intrinsic variance of the correlation function estimates (see Figures 2C–2F).
In principle, the approach developed in this letter is general enough to also generate nonstationary point processes with a predefined correlation structure. We tested the algorithm's ability to generate spike trains with time-varying rates and instantaneous, time-varying correlation structures. Correlation structures in this general scenario are typically represented in neuroscience research using the joint peristimulus time histogram (JPSTH). Figure 3 presents a predefined and a data-estimated JPSTH for processes whose mean has a gaussian temporal profile and displays a complex biphasic time-varying interaction. As can be seen, the algorithm can jointly capture the full temporal structure of both time-varying rates and the interaction.
Next, we tested our algorithm's ability to simulate spike trains resulting from networks of pair-interacting neurons. Figure 4 shows typical activity patterns in a network of 100 units simulated by controlling pair-wise cross-correlation functions. Varying this correlation structure leads to a large diversity of network activity patterns that can be simulated.
4. High-Order Correlation Structure
Nonlinearly transforming gaussian processes, where the entire statistical structure is determined by the first- and second-order moments, introduces higher-order moments, which then carry over to the spike train statistics. These moments depend exclusively on the nonlinear function applied. In Figures 5A and 5B, we present the complexity distribution histogram and quantile-quantile plots for the distribution of synchronous events in spike trains with the same first- and second-order statistics generated by the three different nonlinear transformations. The complexity distribution is a useful way of characterizing the zero-lag high-order correlation structure that is defined in Grün, Abeles, and Diesmann (2008) as the probability p(ξ) to observe a particular number of synchronous spikes ξ (complexity). These results support the notion that the three nonlinear transformations lead to a different higher-order correlation structure. While the mean level of coincidences is the same (here E(ξ) = 5), the complexity distribution resulting from the exponential nonlinear transformation has significantly wider tails (corresponding to a higher proportion of highly synchronous events). The square and absolute value transformations appear to lead to a fairly similar high-order correlation structure.
The generation of continuous gaussian processes with a predefined correlation structure is a classical problem with well-known solutions and numerous applications in all fields of signal and image processing. However, while the empirical correlation structure of neural spike trains has been measured many thousands of times using well-known analysis tools introduced in the 1960s (Perkel, Gerstein, & Moore, 1967a, 1967b; Gerstein & Perkel, 1969), the development of a complementary simulation tool was left largely ignored for four decades and has only recently begun to attract attention. The recent interest in controlled-correlation spike trains is mostly motivated by the emergence of experimental techniques for rapid patterned stimulation of neurons with high temporal and spatial resolution (Shoham, O'Connor et al., 2005; Gasparini & Magee, 2006; Losonczy & Magee, 2006; Nikolenko, et al., 2007), as well as the recent suggestion in neural coding studies that second-order statistics dominate the observed neuroretinal spiking pattern and can be used to explain essentially the full structure (Schneidman et al., 2006; Shlens et al., 2006). In addition, these developments could potentially find applications in econometric time series, earthquake data, network traffic, reliability research, optical communication, and heartbeat variability.
It is interesting to compare our approach with other methods developed recently for controlling the correlation structure of spike trains. The method by Niebur (2007) probabilistically distributes spikes from a “mother” spike train (Kuhn, Aertsen, & Rotter, 2003) to tune synchrony in stationary homogeneous Poisson processes. For addressing the more general problem of temporal correlation functions among different spike trains, Brette (2009, method II) extended this approach by adding a controlled temporal jitter to the spike timing. However, this general strategy was solved only for stationary cases, introduces only positive correlations, and generally leads to strong coupling between the auto- and cross-correlation structure (Tetzlaff et al., 2008). It is also not well connected to central concepts in the theory of random point processes, such as stochastic intensity, and could potentially introduce spurious high-order interdependencies.
In contrast, our strategy was based on using correlated gaussian processes as an underlying structure for generating the spike trains. This general concept was also developed independently in two very recent studies (Brette, 2009; Macke, Berens, Ecker, Tolias, & Bethge, 2009), in both cases using a threshold nonlinearity. The study by Brette (2009, method I) generates rate processes by applying a rectifying nonlinearity to gaussian processes with a specific exponential autocorrelation function. The rectifying nonlinearity is conceptually close to the absolute-value transformation studied here, and formulas for the correlation distortion it creates can be derived using the procedures we introduce in the appendix (see also above). However, as noted by Brette (2009), it has the undesirable property that when the correlations are strong, the rates are zero most of the time, and it is therefore of limited use for simulating realistic spike trains. In addition, Brette limits his discussion to the case of low correlations, without the concept of correlation distortion, which is central to the generality of our approach. In an alternative “dichotomized gaussian” approach introduced by Macke et al. (2009), gaussian processes are thresholded to directly generate the binary spike trains in the time domain, without going through the intermediate step of a rate process (see also Bethge & Berens, 2008, and Johnson, 1994, section IIID, for related analysis). Because the second step is deterministic, the structure is single stochastic, in contrast to the doubly stochastic Poisson structure explored here. The dichotomized gaussian approach requires numerical solutions of nonlinear equations and was found to lead to spike trains with nearly maximal entropy.
Our methods can generate a flexible and very broad range of auto- and cross-correlation structures. For example, the method allows the generation of processes with the same autocorrelation and very different cross-correlation structures (see Figure 2), and either one can practically be very different from an exponential function (in contrast to Brette, 2009). The method is practically limited to correlation structures that lead to positive-definite matrices (see equation 2.20). More generally, the method is limited to correlation structures that are realizable by pure rate covariation, without any precise spike coordination. (Staude, Rotter, & Grün, 2008). The analysis by Staude et al. demonstrated that rate-covariation models have peak cross-correlation between spike trains that are weaker than the respective peak autocorrelations (the smaller of the two). Moreover, when two processes have slow dynamics (broad autocorrelations), they have a limited ability to sustain a narrow cross-correlation.
Appendix: Correlation Distortion for an Absolute Value Transformation
We consider the transformation Λi = |μi + σiXi|, where Xi ∼ N(0, 1) with and rii = 1 that results in folded normal distribution of Λi. The solution is divided into two steps. First, we derive a solution for μi and σi as a function of a desired Rii = E[Λ2i] and E[Λi] and then proceed to derive rij as a function of the calculated μi, σi and the desired Rij. Expressions for the first and second moments of truncated gaussian random variables (both univariate and bivariate) were derived earlier in Rosenbaum (1961). To derive moments for the absolute value transformation, we treat it as a superposition of two (univariate) or four (bivariate) rectified gaussian random variables.
A.1. Step 1.
Bounds on Rii can be derived from the above expressions: .
A.2. Step 2.
In summary, we calculate the transformation parameters (μ and σ) from the desired mean rate and the autocorrelation function at τ = 0 of the rate process, and then proceed to calculate the predistorted correlations rij for the gaussian processes from the desired correlations Rij of the nonnegative processes.
This work was supported by Israeli Science Foundation grant 1248/06 and European Research Council starting grant 211055. We thank the two anonymous reviewers for their comments and suggestions.