## Abstract

Experimental advances allowing for the simultaneous recording of activity at multiple sites have significantly increased our understanding of the spatiotemporal patterns in neural activity. The impact of such patterns on neural coding is a fundamental question in neuroscience. The simulation of spike trains with predetermined activity patterns is therefore an important ingredient in the study of potential neural codes. Such artificially generated spike trains could also be used to manipulate cortical neurons in vitro and in vivo. Here, we propose a method to generate spike trains with given mean firing rates and cross-correlations. To capture this statistical structure we generate a point process by thresholding a stochastic process that is continuous in space and discrete in time. This stochastic process is obtained by filtering Gaussian noise through a multivariate autoregressive (AR) model. The parameters of the AR model are obtained by a nonlinear transformation of the point-process correlations to the continuous-process correlations. The proposed method is very efficient and allows for the simulation of large neural populations. It can be optimized to the structure of spatiotemporal correlations and generalized to nonstationary processes and spatiotemporal patterns of local field potentials and spike trains.

## INTRODUCTION

The generation of artificial spike trains with predetermined spatiotemporal characteristics is of considerable interest in both theoretical and experimental neuroscience. In recent years there have been significant advances in simultaneous recording of neural activity at multiple sites. The resulting data provide strong evidence of spatial and temporal correlations in neural tissue (Bair et al. 2001; Gray et al. 1989; Heck et al. 2002; Kohn and Smith 2005; Zohary et al. 1994). Complex spatiotemporal patterns, such as synchronous and phase-locked oscillatory activity, have been observed experimentally and shown to play an important role in information processing (Engel et al. 2001; Fries et al. 2001; Gray et al. 1989; Kohn and Smith 2005; MacLeod and Laurent 1996; Riehle et al. 1997; Steinmetz et al. 2000; Vaadia et al. 1995; Wehr and Laurent 1996; Womelsdorf et al. 2006). However, the functional role of correlated activity is not fully understood. Theoretical tools can guide our efforts to understand the impact of spatiotemporal activity structure on neural coding (Averbeck and Lee 2006; de la Rocha et al. 2007; Gutnisky and Dragoi 2008; Markowitz et al. 2008; Moreno et al. 2002; Poort and Roelfsema 2009; Romo et al. 2003; Salinas and Sejnowski 2000, 2001; Shea-Brown et al. 2008).

Artificially generated spike trains with predetermined statistical properties could be used to study and manipulate cortical neurons in silico, in vitro, and in vivo *(*Netoff et al. 2005). For instance, optical recording and controlled stimulation of single neurons at many different locations along their dendritic tree is now a possibility (Callaway and Katz 1993; Duemani Reddy et al. 2008; Gasparini and Magee 2006; Huber et al. 2008). These techniques are dramatically increasing our understanding of the computational capabilities of single neurons. Advances in two-photon stimulation (Nikolenko et al. 2007, 2008) allow neurotransmitter uncaging at different synaptic locations in a precisely determined spatiotemporal pattern. Optical stimulation of neuroengineered microbial opsins can be used to activate and inactivate specific neural populations (Boyden et al. 2005; Han and Boyden 2007; Zhang et al. 2006, 2007a,b). These approaches will offer unprecedented insights into the input–output properties of single cells and neuronal networks. However, because the number of potential activation patterns is practically infinite, it will be possible to characterize them only statistically. The generation of neural activity patterns with a given statistical structure is therefore essential in the design of such experiments.

We propose a method to generate spike trains with predetermined mean firing rates and cross-correlations. These measures are typically used to characterize data from multiple electrode recording experiments. The cross-correlogram is widely used to describe the degree of interdependence between the activity of two neurons (Perkel et al. 1967) and study aspects of cortical processing. For instance, the presence and strength of monosynaptic connections can be usually inferred from the cross-correlogram (Alonso and Martinez 1998; Alonso et al. 2001; Fetz and Gustafsson 1983; Usrey et al. 1999; Veredas et al. 2005), whereas periodically recurring peaks in the cross-correlogram indicate the presence of oscillatory neural activity (Muresan et al. 2008). Furthermore, cross-correlograms have been used to study how functional connectivity among neurons depends on stimulus selectivity (Bair et al. 2001; Gray et al. 1989; Gutnisky and Dragoi 2008; Kohn and Smith 2005; Smith and Kohn 2008), task demands (Cohen and Newsome 2008; Poort and Roelfsema 2009; Riehle et al. 1997; Steinmetz et al. 2000; Vaadia et al. 1995), previous stimulation (Gutnisky and Dragoi 2008), and cognitive processing (Engel et al. 2001; Fries et al. 2001; Womelsdorf et al. 2006). Despite some efforts (Bohte et al. 2000; Kuhn et al. 2003; Niebur 2007) there are no general methods to efficiently generate a large number of spike trains with predetermined firing rates and arbitrary spatiotemporal correlations [but see the discussion section and similar methods proposed independently in Brette (2009), Krumin and Shoham (2009), and Macke et al. (2009)].

The proposed algorithm generalizes the following simple idea—a Poisson spike train can be approximated by discretizing time into bins of length Δ*t* and assuming a spike in each bin occurs with probability *r*Δ*t* independently of all other events. More generally, for a collection of spike trains with complex spatial and temporal structure, this probability becomes dependent on past and current spikes in all the spike trains. To capture this structure we generate a continuous stochastic process^{1} by filtering Gaussian noise through a multivariate autoregressive model (AR). We threshold this process to generate the desired collection of spike trains. The parameters of the AR are obtained by transforming the spike train correlations to continuous-process correlations. This approach allows us to optimize our algorithm to efficiently simulate spatiotemporal activity patterns in large groups of neurons. The ideas are easily extended to nonstationary statistics and simulations of neural populations with predetermined spike-triggered averages (STAs) of the local field potential (Fries et al. 2001). The algorithms described are implemented in Matlab (freely available at: http://www.math.uh.edu/∼josic/code).

We first illustrate these ideas in the case of a single spike train. We next extend the approach to generate correlated spike train pairs. Several examples are provided that illustrate the method: single spike trains with almost absolute refractory period, populations with oscillatory cross-correlations, and experimentally constrained spike trains. We demonstrate how the algorithm can be extended to generate spike trains with time-varying statistics and spike-field coherence. Finally, we show how the level of synchrony in a presynaptic neural population can affect the firing rate of a conductance-based neural model.

## METHODS

### Setup

The dependencies of a collection of spike trains are frequently characterized by covariance density functions (cross-correlation function). If *N _{i}*(

*a*,

*b*) denotes the number of spikes fired by neuron

*i*in the time interval (

*a*,

*b*), then the

*covariance density*is defined by (Cox and Isham 1980)

*o*(Δ

*t*

^{2}) represents terms that decay faster than Δ

*t*

^{2}.

In trying to statistically describe a collection of spike trains, one is faced with the problem of estimating *c*_{i,j}(*t*) (Bair et al. 2001; Brody 1999a,b; Perkel et al. 1967). Our goal is the opposite—we prescribe *c*_{i,j}(*t*) and aim to simulate a collection of spike trains with these statistics.

We start by discussing the second-order characterization of a collection of spike trains. Next we review the theory of autoregressive (AR) processes and discuss how to transform an AR process to a spike train using thresholding. Finally, we describe how to find an AR process that produces a collection of spike trains with the desired statistics after thresholding.

### Covariance matrix of univariate and multivariate point process

A spike train *x _{i}*(

*t*) is a sequence of 1s and 0s and the occurrence of a spike between the times

*t*

_{0}and

*t*

_{0}+ Δ

*t*is denoted by

*x*(

_{i}*t*

_{0}) = 1. The magnitude of Δ

*t*determines the resolution at which the simulation is performed.

We assume that each spike train *x _{i}*(

*t*) is stationary, with fixed rate

*r*(generalizations are discussed in a following subsection). The probability of firing during a time interval Δ

_{i}*t*is therefore Pr [

*x*(

_{i}*t*) = 1] =

*r*Δ

_{i}*t*and 〈

*x*(

_{i}*t*)〉 =

*r*Δ

_{i}*t*, with the expectation 〈·〉, taken over all time bins. The covariance between the number of spikes in time bins that are separated by τ units of time is given by

*c*

_{i,i}(0) = 〈

*x*(

_{i}*t*)

^{2}〉 − 〈

*x*(

_{i}*t*)〉

^{2}=

*r*Δ

_{i}*t*+

*O*(Δ

*t*

^{2}) is the variance of

*x*(

_{i}*t*).

To simplify notation, in the following we measure time in discrete increments of Δ*t*, so that *t* = *n* refers to the time bin of width Δ*t*, starting at time *n*Δ*t*. Since *c*_{i,j}(τ) is the covariance of *x _{i}*(

*t*) and

*x*(

_{i}*t*+ τ), we see that the Toeplitz matrix

*x*(

_{i}*t*),

*x*(

_{i}*t*+ 1), … ,

*x*(

_{i}*t*+

*L*− 1). Since

*C*must be positive definite, the covariances

_{i}*c*

_{i,i}cannot be arbitrary.

We denote the binary vector of events at time *t* by *x→*(*t*) = [*x*_{1}(*t*), *x*_{2}(*t*), … , *x _{N}*(

*t*)], where

*N*is the number of spike trains. The matrices

*C*(τ) are the covariance matrix of the random variables

*x→*(

*t*) and

*x→*(

*t*+ τ),

The covariance matrices obtained from the collection of random variables *x→*(*t*), *x→*(*t* + 1), … , *x→*(*t* + *L* − 1) are *C*(0), *C*(1), … , *C*(*L* − 1) (as defined by *Eq. 4*). We can define *C* as the (*NL*) × (*NL*) block–Toeplitz matrix
*N* spike trains are described by *C* and *r→* = (*r*_{1}, *r*_{2}, … , *r _{N}*). Such matrices completely describe the second-order statistics of spike counts over windows of arbitrary size (Cox and Isham 1980).

If the cross-covariances are equal for every neuron pair, that is *c _{cross}*(τ) =

*c*

_{i,j}(τ) for all (

*i*,

*j*), where

*i*≠

*j*, and autocovariances are equal for every neuron,

*c*(τ) =

_{auto}*c*

_{i,i}(τ) for all

*i*, we can write

*C*(τ) as the Kronecker tensor product (⊗)

*I*is the identity matrix of size

_{N}*N*, 1

_{N}is a matrix of 1s, and

*C*and

_{off}*C*are

_{diag}*L*×

*L*Toeplitz matrices, defined as

We next review the theory of autoregressive processes.

### Univariate autoregressive processes

A univariate autoregressive (AR) process of order *L* is defined as
*y*(*t*) is a collection of continuous random variables, each with 0 mean and variance *R*(0), *n*(*t*) are independent univariate Gaussian random variables with 0 mean and variance σ^{2}, and the constants *a _{k}* characterize the AR process. Using arbitrary initial values

*y*(0) to

*y*(

*L*), and iterating

*Eq. 6*starting at time

*t*=

*L*, we see that for subsequent times each

*y*(

*t*) is a sum of Gaussian random variables and thus is itself Gaussian. Multiplying by

*y*(

*t*− τ) and taking the expected value of both sides in

*Eq. 6*, we obtain

Defining *R*(τ) = 〈*y*(*t*)·*y*(*t* − τ)〉 we obtain

*Equations 7* and *8* are known as the Yule–Walker equations (Shanmugan and Breipohl 1988) and can be expressed in a matrix form as

The parameters needed to obtain an AR process with a given autocovariance *R*(τ) can therefore be obtained by

To define a Gaussian process in discrete time with prescribed autocovariance we simply follow the steps described earlier in reverse. Given an autocovariance function *R*(τ) defined for τ ∈ {0, 1, … , *L*}, we find coefficients *a _{k}* from

*Eq. 10*and the variance of the Gaussian noise

*n*(

*t*) from

*Eq. 8*. The AR process obtained by iterating

*Eq. 6*has the desired statistics.

### Thresholding the AR process to obtain a univariate spike train

The simplest way of mapping continuous to binary variables is by assigning the value “1” whenever a threshold is exceeded and “0” otherwise. This operation can be used to convert the AR process described in the previous section to a binary sequence, i.e., a spike train. In particular, if *y*(*t*) is defined using an AR process (*Eq. 6*), we define the spike train *x*(*t*) as
*x*(*t*) = *H*_{θ}[*y*(*t*)], where *H*_{θ} is the Heaviside function. For each *t*, the random variable *y*(*t*) is Gaussian with zero mean and variance *R*(0). Therefore the probability that *x*(*t*) = 1—i.e., the firing rate of the spike train obtained by thresholding—is

The threshold θ divides the *y*(*t*)–*y*(*t* + τ) plane into four regions (Fig. 1) corresponding to a spike at time *t* and time *t* + τ (red), no spikes (black), or a spike either at *t* or *t* + τ (blue). The probability that both *y*(*t*) and *y*(*t* + τ) jointly exceed the threshold is given by the integral of their joint density over the upper right quadrant (red region) in Fig. 1. The dependencies between *y*(*t*) and *y*(*t* + τ) [i.e., *R*(τ)] imply dependencies between *x*(*t*) and *x*(*t* + τ). We define *p*(τ) as the mass of the bivariate Gaussian distribution above the threshold θ. Defining^{2} *a* = *y*(*t*) and *b* = *y*(*t* + τ)

For the spike train *x*(*t*), the covariance is determined by the probability of having two spikes τ time units apart, since *c*(τ) = 〈*x*(*t*)*x*(*t* + τ)〉 − *r*^{2}. Thresholding therefore results in a spike train with autocovariance function related to *p*(τ) by

### Generating a univariate spike train

We now return to the problem of generating a univariate spike train with a prescribed covariance function. This can be achieved by inverting the steps in the previous section—we start with the firing rate *r* of a neuron. Since we assumed that *R*(0)—the variance of *y*(*t*)—is unity, we can invert *Eq. 12* to obtain θ in terms of *r*. Similarly, the autocovariance function *c*(τ) is related to *p*(τ) via *Eq. 14*. We can therefore use *Eq. 13* to determine *R*(τ) from the threshold θ and *c*(τ).

In particular, fixing θ and solving for *R*(τ) in terms of *p*(τ) using numerical integration gives a monotonic function that can be used to convert the autocovariance function of the spike train into the autocovariance of a continuous stochastic process (Fig. 2).

The resulting autocovariance *R*(τ) determines the coefficients *a _{k}* of the AR process (

*Eqs. 9*and

*10*). The continuous variable

*y*(

*t*) obtained by filtering the Gaussian noise

*n*(

*t*) using the filter defined by the coefficients

*a*has autocorrelation

_{k}*R*(τ), which equals the autocovariance when the variance of

*y*(

*t*) is unity. By construction, applying the Heaviside function with threshold θ to

*y*(

*t*) generates a spike train

*x*(

*t*) with mean firing rate

*r*and autocovariance

*c*(τ).

### Extension to a multivariate process

In this section we give a brief overview of the extension to populations of spike trains with firing rates *r→* = (*r*_{1}, *r*_{2}, … , *r _{N}*) and cross-correlation functions

*c*

_{i,j}(τ). The multivariate point process is thus defined through the (

*N*×

*L*) × (

*N*×

*L*) block–Toeplitz matrix

*C*(

*Eq. 5*). A vector of thresholds θ→ = {θ

_{1}, θ

_{2}, … , θ

_{N}} is calculated using

*Eq. 12*based on the vector of firing rates

*r→*. The covariance matrix of the continuous stochastic process is obtained by applying

*Eq. 13*to each pair (

*i*,

*j*) for all time lags τ→ = (τ

_{1}, τ

_{2}, … , τ

_{L−1}) to get the covariances

*R*

_{i,j}(τ). We obtain the following (

*N*×

*L*) × (

*N*×

*L*) block–Toeplitz covariance matrix

*R*(τ) is an

*N*×

*N*matrix [note that

*R*

_{i,j}(τ) =

*R*

_{j,i}(−τ)]

A multivariate autoregressive process of order *L* is defined as
*A _{k}* is a symmetric

*N*×

*N*matrix. We again have to solve

*A*and

*R*given by

The covariance matrix of the multivariate Gaussian random variable *n→*(*t*) is
*x→*(*t*) the vector variable *y→*(*t*) is thresholded using the vector of thresholds θ→.

### Summary of the algorithm to generate a multivariate process

We summarize the steps used to generate spike trains with given spatiotemporal structure. The inputs to the algorithm are the desired cross-correlations *c*_{i,j}(τ) for every pair of neurons (*i*, *j*) for a finite number of time lags τ ∈ {0, 1, … , *L*} and the firing rates of each neuron *r _{i}*. An outline schematic of the algorithm is shown in Fig. 3.

In the simulations, the spike-train covariances were normalized by the geometric mean of the firing rates
*R*(0) = 1 (i.e., the continuous variables will be Gaussian with unit variance). Thus the covariances equal the correlations and are measured in units of coincidences/spike. Using this normalization, the probability of having coincident spikes in two bins τ time units apart for two neurons *i, j* is *t* = 1 in all simulations). We chose this normalization because it is one of the most commonly used in experimental work (e.g., Bair et al. 2001; Kohn and Smith 2005). Using *Eqs. 12* and *13* we calculate the thresholds θ_{i}. The continuous-process correlations *R*_{i,j}(τ) are obtained by finding bivariate Gaussian distributions with mass above θ equal to *p*(τ) (*Eq. 13*).

To obtain the parameters of the multivariate autoregressive process, *A*_{i,j}(*k*) and Σ, we solve the linear system in *Eqs. 16* and *18*. To obtain the continuous variables *y _{i}*(

*t*), we filter the time-independent Gaussian random noise

*n*(

_{i}*t*) with covariance matrix Σ, through the AR with parameters

*A*

_{i,j}(

*k*) using

*Eq. 15*(Fig. 4). Finally, the population spike trains

*x*(

_{i}*t*) are obtained by thresholding

*y*(

_{i}*t*) with θ

_{i}for each neuron

*i*.

### Notes on the implementation: complexity and memory requirements

We next discuss the computational costs and optimization strategies available in certain special cases. After converting the block–Toeplitz covariance matrix *C* into the continuous-process covariance matrix *T*, we have to solve a linear system *A* = *T*^{−1}·*R*. The matrix *T* is *NL* × *NL*, whereas *A* and *R* are *NL* × *L* (*N* is the number of neurons and *L* is the order of the AR). The Cholesky factorization of *T* is *O*(*N*^{3}*L*^{3}), but the subsequent operations needed to solve the linear system are only *O*(*N*^{2}*L*^{2}). More efficient implementations are available by taking advantage of the block–Toeplitz structure of *T* (it has only *N*^{2}*L* different entries) and many algorithms are applicable to positive-definite block–Toeplitz matrices (Akaike 1973; Ammar 1996; Ammar and Gragg 1988; Kamm and Nagy 2000) and the special case of Toeplitz–block–Toeplitz linear systems (Wax and Kailath 1983).

After obtaining the parameter matrix *A* we have to generate Gaussian noise with covariance Σ and apply the autoregressive equations (*Eq. 15*). Defining *M* as the number of time steps that will be computed, the number of computations needed to obtain *y→*(*t*) is *O*(*N*^{2}*LM*). Therefore the computational time is linear in the spike-train length but quadratic in the number of neurons. Alternatively, the AR process can be written as an infinite-impulse response (IIR) filter (Shanmugan and Breipohl 1988). In this case, the computation of *y→*(*t*) is *O*(*N*^{2}*M* log *M*) for *L* ≪ *M*. The last equation indicates that the computational time is independent of the maximum correlation length *L* when using the frequency-domain approach. With programs such as Matlab the use of loops may be less efficient than implementing the filter approach, especially when *L* ≫ log (*M*). It is also possible to obtain *y→*(*t*) using a moving-average process (FIR filter) instead of an AR process (Krumin and Shoham 2009).

### Optimizations

In many practical cases the spatiotemporal covariances have additional structure. An important case is when the covariance is the product of a function of the neuron number and a function of the time lag

Note that *R* will be symmetric since *R _{ij}*(τ) =

*R*(−τ) implies

_{ji}*g*(τ) =

*g*(−τ). The last equation allows us to have a different time dependence for the cross-correlations and the autocorrelations. Using the Kronecker tensor product we can write

*T*compactly

*G*and

*D*are

*L*×

*L*Toeplitz matrices with row entries

*g*(τ) and

*d*(τ), respectively (τ = 0, 1, … ,

*L*− 1), and

*F*is an

*N*×

*N*symmetric matrix with values [

*F*]

_{ij}=

*f*(

*i*,

*j*). We also define the matrices

*R*and

_{off}*R*as the

_{diag}*L*×

*L*Toeplitz matrices, with row entries

*g*(τ) and

*d*(τ), respectively, with τ = 1, 2, … ,

*L*.

We can factorize *F* = *P*Λ*P*^{t} with the columns of *P* having the eigenvectors of *F* and Λ the eigenvalues (λ_{i}) on its diagonal. We define the *N* × *N* matrices Λ_{i} to have a single nonzero entry [Λ_{i}]_{ii} = λ_{i}. The linear system *TA* = *R* takes the form
*A _{off}*(

*i*) and

*A*(

_{diag}*i*) are the parameter matrices to be determined. Using the fact that

*PP*

^{t}=

*P*

^{t}

*P*=

*I*and Λ

_{N}_{i}Λ

_{j}= 0 for

*i*≠

*j*

Instead of having to solve one *NL* × *NL* linear system we have to solve *N* + 1 Toeplitz-structured systems of size *L* × *L*. For such systems the calculations are linear in the number of neurons and the memory requirement for the AR parameters drop to (*N* + 1)*L*.

We also consider the special case of a population consisting of *N _{sub}* subpopulations, where subpopulation

*i*contains

*n*neurons. The cross-correlations for neurons within a subpopulation are identical, whereas neurons in different subpopulations are uncorrelated. The matrix

_{i}*F*can then be written as

*F*are 0 except for the

*N*eigenvalues (λ

_{sub}_{1}, λ

_{2}, … , λ

*) = (*

_{Nsub}*n*

_{1},

*n*

_{2}, … ,

*n*). In this way,

_{Nsub}*Eq. 20*reduces to

*N*+ 1 equations to solve and thus the number of linear systems depends on the number of subpopulations but does not depend on the number of cells within each subpopulation.

_{sub}### Generalization to nonstationary spike trains

In the previous sections we assumed that the firing rate and the correlations do not change in time. A spike train with time-dependent firing rate *r*(*t*) can be obtained by using a time-dependent threshold and solving *Eq. 12* at each time step. The continuous-process autocovariance *R*(*t*_{1}, *t*_{2}) = cov [*y*(*t*_{1}), *y*(*t*_{2})] = 〈*y*(*t*_{1})·*y*(*t*_{2})〉 is obtained by solving a time-dependent version of *Eq. 13*
*R*(*t*_{1}, *t*_{2}) = *R*(*t*_{2} − *t*_{1}) = *R*(τ). Note that *R*(*t*_{1}, *t*_{2}) may be time dependent, even if the point-process autocovariance *c*_{i,i}(τ) is not because *R*(*t*_{1}, *t*_{2}) depends on the threshold θ(*t*).

The next step is to generate a continuous process with zero mean and covariance *R*(*t*_{1}, *t*_{2}). We generate the nonstationary continuous process by means of a time-varying autoregressive process
*a _{k}*(

*t*) are time dependent. For a given, nonstationary cross-covariance

*R*(

*t*

_{1},

*t*

_{2}), these parameters can be obtained by solving the nonstationary Yule–Walker equations (Hallin and Ingenbleek 1983)

*Eq. 23*to obtain

*A*(

*t*) and then calculate

*a*

_{0}(

*t*) in

*Eq. 24*to get the parameters

*a*

_{1}(

*t*),

*a*

_{2}(

*t*), … ,

*a*(

_{L}*t*). The matrices

*T*(

*t*) and

*V*(

*t*) are defined as

Note that in the stationary case we used *a*_{0} = 1 and found the SD using *Eq. 8*. For consistency with Hallin and Ingenbleek (1983), here we let *a*_{0}(*t*) be a normalizing factor and let the SD of *n*(*t*) equal 1.

### Generalization to STAs of the LFP

The generation of continuous random variables *y→*(*t*) with given spatiotemporal correlations is an intermediate step in our algorithm to create spike trains. However, the vector of continuous variables can also be used to directly simulate a collection of local field potentials (LFPs) with specific spectral components. In this section we show how we can extend our algorithm to simulate given spike-triggered averages (STAs) of the LFP.

Let *y _{LFP}*(

*t*) be the LFP signal and

*x*(

*t*) be the spike train associated with the continuous variable

*y*(

_{spike}*t*). As before,

*x*(

*t*) is obtained by thresholding

*y*(

_{spike}*t*) at θ (

*Eq. 11*). The STA is defined as

*t*is the time of each spike in

_{k}*x*(

*t*). We can find the correlation

*R*(τ) between the bivariate Gaussian distribution [

*y*(

_{LFP}*t*),

*y*(

_{spike}*t*)] to obtain a desired value of

*STA*(τ). We have to find the distribution of

*y*(

_{LFP}*t*+ τ) given that a spike occurred at time

*t*, i.e.,

*y*(

_{spike}*t*) ≥ θ (Cox and Wermuth 1999). The corresponding density is given by

*f*,

_{STA}*f*, and

_{yspike}*f*are Gaussian probability density functions (pdfs) of the corresponding random variables. The conditional pdf

_{yLFP|yspike≥θ}*f*has mean μ

_{yLFP|yspike≥θ}*and variance σ*

_{yLFP|yspike≥θ}

_{yLFP|yspike≥θ}^{2}We use σ* _{yspike}* = 1. Calculating the integral in

*Eq. 27*we obtain

In this way, we can generate the STA of the LFP (*Eq. 30*) by choosing an appropriate σ* _{yLFP}* and finding the correlation function

*R*(τ).

### Conductance-based model

We implemented a standard conductance-based integrate-and-fire model to study the impact of correlated input on the firing properties of downstream neurons. The membrane voltage *V*(*t*) is described as follows
*C* is the membrane capacitance; *G _{l}* is the leakage conductance;

*V*is the resting potential; and

_{r}*V*and

_{e}*V*are the excitatory and inhibitory reversal potentials, respectively. The excitatory and inhibitory conductances

_{i}*G*(

_{e}*t*) and

*G*(

_{i}*t*) are sums of unitary synaptic events,

*g*(

_{e}*t*) and

*g*(

_{i}*t*)

*t*and

_{j}*t*.

_{k}The unitary synaptic events were modeled by alpha-functions with peak conductance *B _{e}* and

*B*and time constants τ

_{i}_{e}and τ

_{i}

*H*(

*t*) is the Heaviside step function. The parameters and the model are the same as those in Kuhn et al. (2004):

*C*= 250 pF;

*G*= 1/60 ms;

_{l}*V*= −70 mV;

_{r}*V*= 0 mV;

_{e}*V*= −75 mV;

_{i}*B*= 7.1 ns;

_{e}*B*= 3.7 ns; τ

_{i}_{e}= 0.2 ms; τ

_{i}= 2 ms. An action potential is elicited whenever the membrane voltage

*V*(

*t*) crosses the threshold

*V*= −50 mV. An absolute refractory period is obtained by holding

_{th}*V*(

*t*) at qj

*V*= −60 mV for τ

_{reset}_{refr}= 2 ms after a spike. We used a time step of 0.1 ms for the numerical simulations.

## RESULTS

To test our algorithm we generated a variety of spike trains with different spatiotemporal structures. We first present single spike trains with prescribed autocorrelations. These are followed by examples of multivariate spike trains that illustrate the impact of pairwise correlation on the behavior of the neural population.

### Single-cell examples

We briefly summarize the steps used to generate a single spike train with a given cross-correlation (for details see methods). In the simplest case, a Poisson spike train with mean firing rate *r* can be obtained by discretizing time into bins of length Δ*t*. A spike is generated with probability *r*Δ*t* in each bin independently of all other events. We can implement a Poisson spike train by thresholding a continuous-in-space, discrete-in-time Gaussian stochastic process *y*(*t*). A spike is elicited whenever the Gaussian random variable is higher than a threshold θ (*Eq. 11*). Using *Eq. 12* we can obtain the appropriate θ for a spike train with mean firing rate *r*.

More generally, for spike trains with complex temporal structure, the probability of having a spike in a given time bin *r*Δ*t* becomes dependent on past spikes. To capture this structure we have to generate a stochastic process *y*(*t*) with specific autocovariance function. After finding the threshold (θ) we have to map the desired spike-train autocovariance, *c*(τ), into a continuous process covariance function *R*(τ). This transformation is computed for every time lag in such a way that the application of the threshold θ to a pair of continuous variables [*y*(*t*), *y*(*t* + τ)] with a bivariate Gaussian distribution gives the desired point-process autocovariance *c*(τ) (Figs. 1 and 2; *Eq. 13*).

The second step consists of generating a continuous random variable *y*(*t*) with prescribed correlation *R*(τ). This random variable follows a Gaussian distribution with mean 0 and is thus fully described by *R*(τ). We have to design a filter to transform independent random Gaussian noise *n*(*t*) to *y*(*t*). This is accomplished by generating an AR process (*Eq. 6*). The parameters of the AR process are calculated using *Eqs. 8* and *10*. Finally, *y*(*t*) is obtained by applying the AR process equations (*Eq. 6*).

The last step involves the only nonlinear operation in the algorithm: we assume a spike occurs at time *t*, so that *x*(*t*) = 1, whenever *y*(*t*) > θ. In this way, we obtain a spike train *x*(*t*) with the specified firing rate *r* and correlation *c*(τ). The explicit representation using a continuous function *y*(*t*) and a threshold θ is similar to the ubiquitous integrate-and-fire model: *y*(*t*) plays the role of the scaled membrane potential, which elicits a spike when it crosses the threshold θ.

To demonstrate the algorithm we provide several illustrative examples. First, we generated a spike train with vanishing autocorrelation, except at time *T* = 10 ms (for all simulations Δ*t* = 1 ms; see Fig. 5*A*). In Fig. 5*B* we present a collection of simulated spike trains aligned to the spike occurring at the origin. Trials were sorted by the time of the second spike (red dots) and an increase in spiking probability at *T* = 10 ms is clearly visible. This is precisely what we expect, in that the autocorrelation is directly related to the conditional probability of a spike at time *T* = 10 ms, given a spike at *time 0*.

In the second example we examined how the cross-correlation estimated from the simulated spike trains converges to the predetermined cross-correlation with an increase in the length of the simulation and the amount of data generated. We simulated a spike train with exponentially decaying autocorrelation and set the probability of coincident spikes at *t* ± 1 ms to 10% of chance (see following text for an example of absolute refractory period). In Fig. 5*C*, the red and black curves correspond to the autocorrelation function obtained from simulations (6.4 × 10^{6} bins) and theory, respectively. To assess the precision of the algorithm we calculated the mean squared error (RMS) for spike trains of different lengths (Fig. 5*D*). The error decayed as 1/

In the third example, we approximated a spike train with absolute refractory period of 1 ms. Our algorithm does not include any mechanism to remember when a spike was elicited to prevent the firing of another action potential during the absolute refractory period. However, under certain conditions, such as short refractory periods or low firing rates, we can obtain a good approximation to absolute refractory periods. To ensure the absence of pairs of spikes in adjacent time bins, it is necessary that *R*(±1) = −1 for the corresponding AR process. However, unless *R*(±1) > −0.5 the covariance matrix of *y*(*t*), *y*(*t* + 1), … , *y*(*t* + *L* − 1) may not be positive definite and thus the AR process cannot be generated. However, it is possible to generate processes with large relative refractory periods. The degree of refractoriness that can be achieved depends on the firing rate of the neuron and the discretization Δ*t*. For a Poisson spike train with a firing rate of 25 Hz the probability of having two consecutive spikes is (0.025)^{2} = 6.25 × 10^{−4} (Fig. 6*A*). We generated spike trains where this probability is reduced to about 4 × 10^{−6}, a decrease of nearly two orders of magnitude compared with chance. We show that the distribution of the ISI of an inhomogeneous Poisson spike train with absolute refractory period of 1 ms fits the empirical histogram well (Fig. 6*B*; the *inset* is a zoom of the histogram in the main figure).

We explored how well the algorithm could generate longer refractory periods. We derived the minimum correlation coefficient of the continuous process that could be achieved as a function of the length of the refractory period (Fig. 6*C*; see the appendix). The minimum correlation coefficient *R*_{min} is inversely related to the length of the refractory period *m*, as *R*_{min} = −0.5/*m*. The probability of having a spike during the refractory period depends not only on *R*_{min}, but also on the firing rate of the neuron. We defined the refractoriness index as the ratio of the probability of having a spike during the refractory period to the probability of the same event for a homogeneous Poisson process. In Fig. 6*D* we plot the refractoriness index for different firing rates and refractory period lengths. The algorithm is better at generating refractory periods for low firing rates and short refractory periods.

Cortical responses are typically nonstationary. In sensory cortices, neurons respond to external stimulation with a fast increase in the firing rate followed by a slower adaptation response (e.g., Muller et al. 1999). To demonstrate that our algorithm can generate such responses we simulated a neuron with a nonstationary firing rate (Fig. 7*A*), but stationary cross-correlation function. The simulation was started at a baseline firing rate *r _{baseline}* = 5 Hz. After 50 ms, the firing rate was first set to increase exponentially with rate τ

_{rise}= 10 ms and then set to exponentially decay back to baseline with rate τ

_{decay}= 30 ms. In between, the peak firing rate was set at

*r*= 30 Hz. The firing rate was therefore described as

_{peak}*z*is a normalization constant. This behavior emulates a typical neural response evoked by a stimulus at 50 ms: a fast phasic increase in firing rate followed by a slower adapted response.

The autocorrelation function was chosen to be time independent, with a peak of 0.02 coincidence/spike and an exponential decay rate of 30 ms (Fig. 7*B*). To generate nonstationary spike trains we used a time-dependent threshold θ(*t*) so that the underlying continuous Gaussian process, *R*(*t*_{1}, *t*_{2}), was also time dependent. The main difference with the stationary process algorithm is the determination of time-dependent parameters of the autoregressive process. In Fig. 7, *C* and *D* we show that we can generate a nonstationary continuous Gaussian process using *Eqs. 22–26* (see methods). As in the case of the stationary process, we obtained spike trains by thresholding the continuous process at θ(*t*). Because we chose the autocorrelation to be stationary, after thresholding we obtained the prescribed spike-train autocorrelations *C*(*t*_{1}, *t*_{2}) (defined as *C*(*t*_{1}, *t*_{2}) = cov [*x*(*t*_{1}), *x*(*t*_{2})]; Fig. 7, *E* and *F*).

### Multiple-cell examples

We next provide examples of both homogeneous and nonhomogeneous population activity. In the first two examples we focus on the temporal aspects of the correlation structure. In the last example, we illustrate the impact of spatial correlations in the activity of a presynaptic population on the firing of an afferent neuron. In all simulations we used an optimized version of the algorithm applicable when the spatial and temporal correlations are “decoupled” (see methods).

We first illustrate how correlations between cells can affect population activity in ways that are difficult to detect at the level of small groups of neurons, but may be significantly amplified in larger populations. To this effect, we simulated the activity of homogeneous neural populations of 200 neurons with exponentially decaying autocorrelations (Fig. 8*A*). In contrast, we chose the cross-correlations to be oscillatory, with frequency varying between 20 and 150 Hz (Fig. 8, *B–D*). As a consequence, no periodic patterns were observable in the activity of individual cells, whereas the pooled activity of many neurons clearly exhibited rhythmic activity at the predetermined frequency. As expected, the activity of individual cells displayed approximately constant power at all frequencies. However, the population activity (i.e., the summed response of all neurons) displayed a distinguishable peak at the chosen frequencies (Fig. 9*A*). We computed STAs of the pooled activity to determine whether single neurons were synchronized to the population activity. We illustrate for an example of an oscillatory frequency of 50 Hz that the STA exhibits oscillations at the frequency of the cross-correlation functions. To further characterize this behavior we computed the power spectrum of the STA to find peaks at the specified frequencies (Fig. 9*C*).

We note that LFPs are believed to represent the pooled activity of large populations. LFPs show strong oscillations in different frequency bands depending on brain area, behavioral state, and stimulus (e.g., Engel et al. 2001; Fries et al. 2007; Schroeder and Lakatos 2009). On the other hand, there are many electrophysiological reports (e.g., Bair et al. 2001; Kohn and Smith 2005) of single units that failed to find oscillatory cross-correlations. As this example illustrates, it is possible that oscillations in spike-train cross-correlograms are simply difficult to detect with limited data. However, even weak, oscillatory features that are swamped by noise and undetectable at the single-cell or small population level could be amplified at the level of larger populations and are thus detectable in the LFP.

We also studied the impact of temporal correlation in the activity of a presynaptic population on the firing rate of a postsynaptic neuron. In this simulation the autocorrelations and cross-correlations of spike trains in the simulated presynaptic population exhibited exponential decay in time (auto- and cross-correlations had the same parameters). We kept the area of the auto- and cross-correlogram fixed at 0.2 while varying the exponential decay constant between 1 and 100 ms (we used Δ*t* = 0.1 ms). The area under the cross-correlogram is proportional to the Pearson spike-count correlation between two neurons (Bair et al. 2001; Cox and Isham 1980). Therefore the Pearson correlation coefficient was held constant, whereas the precision of synchrony between the cells was varied. The postsynaptic neuron also received inhibitory Poisson input at 1,000 Hz. Figure 10*A* shows some of the cross-correlation profiles of single cells in the presynaptic population. The pooled population activity (*n* = 100) displayed autocorrelation similar to that of individual cells (Fig. 10*B*), but increased magnitude. This increase was not due to an increased signal strength (i.e., increase in firing rates), but an amplification caused by effective averaging over correlated variables. The amplification of the autocorrelation function is explained by the fact that a random variable that is the sum of small correlated variables exhibits an autocorrelation approximately *Nc*, with *N* the number of variables and *c* the correlation magnitude.

Although the mean population firing rates were equal for all decay constants, the instantaneous population firing rate variability was greatly increased when the decay rate was smaller (i.e., the population was more synchronized; Fig. 10*C*). As expected, this caused an increase in the firing rate of the postsynaptic neuron (Fig. 10*D*; Salinas and Sejnowski 2000). Moreover, changes in the decay rate of input correlations influenced the variability of the output of the afferent neuron. As shown in Fig. 10*E* slower decay rates led to outputs that are more bursty. Indeed, an increased correlation timescale in the input led to broader synchronous events in the input and an increased probability that a threshold crossing will be closely followed by others in the afferent cell. We quantified the variability in the ISIs by computing their coefficient of variation (CV). We found that the output CV was close to 1 for tightly synchronous inputs. Longer decay rates in input correlations led to increased output burstiness and thus a higher CV. Although the impact of changes in decay rate on input CV (i.e., single-cell input) was negligible, the effect on output CV was considerable (Fig. 10*F*).

In the next example we generated different examples of heterogeneous populations to demonstrate the impact of variations in spatial correlation profile on the output of a postsynaptic neuron. We simulated several populations of 128 neurons, comprised of a different number of homogeneous subpopulations. The neurons within a subpopulation shared the same auto- and cross-correlations but were uncorrelated with cells in other subpopulations. The number of subpopulations ranged between 1 (a single subpopulation with statistically identical neurons and neuron pairs) and 128 (all neurons uncorrelated). The mean firing rate of each neuron was 20 Hz, whereas the auto- and cross-correlations were exponentially decaying functions with τ_{decay} = 10 ms, with the area of the cross-correlogram equal to 0.4. The structure of these populations allowed for the use of very efficient versions of the algorithm (see methods).

As shown previously (Fig. 10), neurons are sensitive to highly synchronized inputs (Kuhn et al. 2003; Salinas and Sejnowski 2000). However, if the input contains a high number of synchronous spikes, adding further synchronous spikes does not cause an increase in the firing rate of the output neuron (Kuhn et al. 2003). The distance of the resting membrane voltage to the threshold determines the optimal synchronization level at which the number of spikes in a synchronous volley is exactly sufficient to cross threshold and no spikes are “wasted.” Another factor that affects the firing rate is the frequency of highly synchronous events. In homogeneous populations, the level of synchronization depends on the level of correlation between neurons. To show the interplay between synchronized neural populations and the frequency of synchronized events we generated eight different populations by modifying the number of subpopulations (from 1 to 128 in steps determined by powers of 2) and used their output to drive a conductance-based neural model. Increasing the number of subpopulations caused a higher rate of synchronized events because the synchronous events in each subpopulation were independent. On the other hand, there were also fewer neurons in each subpopulation causing a reduction in the impact of the individual synchronous events.

In this way, we kept the same mean population firing rate, but altered the probability and the number of neurons involved in synchronous events. In addition, we manipulated the background excitation of the output neuron. We used positive rates to indicate excitatory inputs and negative ones to indicate inhibitory inputs. We found that the optimal number of subpopulations depended on the level of inhibition to the cell (Fig. 11). For a highly inhibited cell, a critical factor to make it fire was to have highly synchronized events. This is because to drive the neuron to fire, a large number of cells have to be active at the same time to cross the threshold. For less-inhibited cells, less synchronization was needed (i.e., it is more efficient to have more synchronous events of smaller amplitudes). Populations with fewer subpopulations were less efficient in driving the neuron because the synchronous events involved a larger number of neurons than necessary. Thus the optimal number of subpopulations depended on the background inhibitory input to the output neuron.

Finally, we tested whether our algorithm could generate spatiotemporally correlated spike trains that follow experimentally observed statistics. We used spike trains from a publicly available database (http://crcns.org/data-sets/pvc/pvc-3/). The data were recorded from cat area 17 (primary visual cortex) using multielectrode silicon probes. Anesthetized cats were stimulated with drifting gratings for 722.8 s. We computed the mean firing rate and covariance functions of four selected neurons that had a relatively large number of spikes. In this example we did not normalize the covariance function to display the raw number of spike coincidences. Next, we generated spike trains with the measured first- and second-order statistics obtained from experimental data. Figure 12 shows an excellent agreement between our simulations (red curves) and experimental data (black curves). Our algorithm was able to generate spike trains with irregular peaks in the auto- and cross-covariance functions. Interestingly, the algorithm was able to capture the strongly asymmetric covariance functions of cell 1 and cell 2.

### LFPs and STA of the LFP

In Fig. 13 we show how the present algorithm can be used to simulate LFPs and spike trains synchronized to the LFP. We simulated a particular case where the spike train had a constant firing rate with no oscillatory autocorrelation, although it was synchronized to one of the frequency components of the LFP. The spike train had a constant firing rate of 20 Hz and an exponentially decaying autocorrelation (Fig. 13*A*; τ_{decay} = 30 ms). The LFP had a longer time constant (τ_{decay} = 70 ms) and spectral peaks at 30 and 80 Hz (Fig. 13*B*). The nonnormalized spike-field coherence (Fig. 13*C*, *inset*) shows that the spike train was synchronized only to the 80-Hz component of the LFP.

Here we showed an example of a single neuron spike train synchronized to an LFP. In the general case of multiple neurons the LFPs are not required to be associated one to one to each of the spike trains. For instance, we can define fewer LFPs than neurons.

## DISCUSSION

We have presented a conceptually simple and computationally efficient method that uses thresholding to translate a multivariate AR process to a multivariate point process with predetermined spatiotemporal statistics. The parameters of the AR process are determined by the desired statistical properties of the spike trains. The algorithm is particularly efficient when simulating homogeneous populations or when the spatial and temporal dependencies of the correlations can be decoupled. Moreover, the Yule–Walker equations, used to determine the parameters of the AR, can be naturally extended to nonstationary first- and second-order statistics and to spike-field coherence simulations.

We next compare our algorithm to previously proposed methods and discuss the limitations, statistical properties of the generated spike trains, and possible applications of the algorithm.

### Relation to alternative methods

In most previous efforts to devise algorithms to generate spike trains with specified statistics the proposed methods did not allow a full control of mean firing rates and spatiotemporal correlations (Bohte et al. 2000; Kuhn et al. 2003; Mazurek and Shadlen 2002; Niebur 2007; Oram et al. 1999; Song et al. 2000). However, approaches similar to the present method were very recently developed in parallel (Brette 2009; Krumin and Shoham 2009; Macke et al. 2009). We first describe some generalities of the previously developed algorithms and then relate these new algorithms to our own approach.

Of particular interest is the “thinning and shifting” process described in Bauerle and Grubel (2005). In a strong mathematical sense, this algorithm can be used to generate any stationary multivariate point process in which the constituent point processes are marginally Poisson. Since this algorithm contains a number of others as special examples (method II of Brette 2009; Kuhn et al. 2003; Niebur 2007), we provide a brief description and discussion.

To generate *N* spike trains start with defining a probability distribution on the 2^{N} − 1 nonempty subsets of the set of numbers from 1 through *N*. Generate a “mother” Poisson spike train and assign to each spike one of the 2^{N} − 1 subsets with the given probability. Now generate *N* daughter processes by creating train *i* from spikes that have been marked by subsets containing the number *i*. For instance, when *N* = 3, the set {1, 3} could have probability *p*_{1,3} = 0.1, so that spikes in the mother spike train are marked with 1,3 with probability 0.1. A spike at time *t* marked by {1, 3} is copied to spike trains 1 and 3.

The final step of the algorithm is to choose a vector *v→ ^{i}* from a distribution

*Q*on ℜ

^{N}for each spike

*t*in the mother spike train. Now shift (or jitter) the spike at time

_{i}*t*in daughter train

_{i}*j*by the amount

*v*. Note that the spike

_{j}^{i}*t*appears only in some daughter trains and only those will be affected.

_{i}Although this procedure results in Poisson spike trains with correlations of all orders easily determined from *Q*, it also has several shortcomings. First, only Poisson spike trains preserve their characteristics under thinning and jittering. Thus the procedure is difficult to generalize to cases when the individual spike trains are not marginally Poisson. Second, the cross-correlation functions are given in terms of double integrals of the density of *Q*. Thus to obtain the density *Q* that produces a desired cross-correlation one would have to solve a set of nontrivial integral equations.

The most recent methods described in Brette (2009), Krumin and Shoham (2009), and Macke et al. (2009) are more closely related to our approach. The method proposed by Macke et al. (2009) also uses the thresholding of Gaussian processes to generate spatial correlations (see also Qaqish 2005). However, in contrast to the present algorithm, this method samples a multivariate Gaussian distribution conditioned on the last *K* samples to generate temporal correlations. We instead used a multivariate AR that we generalized to the nonstationary case and provided extensive simulations of large neural populations exploiting the structure of the covariance matrix.

Krumin and Shoham (2009) proposed generating a doubly stochastic Poisson process (i.e., a Cox process) with the prescribed mean firing rates and cross-correlations by nonlinearly transforming a Gaussian processes into a nonnegative rate process. They start by choosing an appropriate nonlinear transformation (exponential, square, or absolute value) that is applied to a correlated Gaussian process to yield a continuous nonnegative process that represents the instantaneous firing rate of each neuron. The final step is to use the Gaussian process to determine the intensity of a doubly stochastic Poisson process. The correlated Gaussian process is generated by applying a finite-impulse response (FIR) filter to an uncorrelated Gaussian process (in contrast to the use of AR in our case). The parameters of the FIR are derived from the nonlinear transformation of the desired firing rates and the desired correlation structure. One advantage of the algorithm is that the transformation of the desired firing rates and correlations can be directly calculated, whereas in our case this requires numerical evaluation. Their implementation is slower because, first, a nonnegative process needs to be generated by filtering Gaussian noise (similarly to our case), but in addition the resulting rates need to be used to generate a point process.

Brette (2009) proposed two algorithms. One was an extension of the method of “thinning and shifting” point processes discussed earlier, whereas another involved doubly stochastic processes, with the analysis restricted to exponential correlations. The stochastic rate is obtained by generating an Ornstein–Ulhenbeck process. The rate process has to be nonnegative, although the Ornstein–Ulhenbeck process is not. To generate a nonnegative rate the realization of the Ornstein–Ulhenbeck process is rectified with a Heaviside step function. Brette (2009) notes that this transformation implies that spike trains with strong correlations will have vanishing rates most of the time. Thus this method is not suitable for the generation of strongly correlated spike trains on fine timescales.

Krumin and Shoham (2009) explained how to generate nonstationary firing rates and cross-correlations explicitly, whereas Brette (2009) described time-varying firing rates. In addition, Brette's method had stronger restrictions on the possible cross-correlation functions. As far as we know, the present approach is the first to exploit the covariance matrix structure for efficient simulations of large neural populations with a natural extension to generate LFPs and STA of the LFP.

### Statistical properties of the generated spike trains

The statistical structure of a collection of spike trains is not fully characterized by the first two moments. The structure of higher-order spatiotemporal correlations may have a significant impact on the neural code (Kuhn et al. 2003). It is therefore desirable that the generated spike trains have the least possible amount of structure beyond the first two moments that are specified. In other words, the generated spike trains should have the maximal entropy among all the spike trains with the same first- and second-order statistics (Schneidman et al. 2006). As shown in Macke et al. (2009), point processes generated by thresholding multivariate Gaussian distributions appear to have close to maximal entropy in some parameter regimes (Bethge and Berens 2008).

### Limitations

One requirement of our method is that the continuous-process covariance matrix has to be positive definite. In general, a positive-definite point-process covariance matrix *C* does not imply a correspondent positive-definite continuous-process covariance matrix *R*. This implies that point processes with covariance matrices with strong negative correlations will not be able to be generated with our method. This could potentially be circumvented by using an AR process with other nonlinearities or alternative spike-generation mechanisms (Krumin and Shoham 2009). In the appendix, we derive the conditions under which we can obtain a spike train with a prescribed refractory period.

### Applications

Multielectrode recording technology, calcium imaging (Greenberg et al. 2008), and voltage-sensitive dyes (Chen et al. 2006) offer unprecedented insights into the structure of correlated activity in neuronal populations. Because the role of correlations in the neural code is difficult to intuit, theoretical tools that aid in the exploration of their impact will be of significant importance. The techniques we have presented will aid this study in at least two ways: *1*) the algorithm offers a way to generate a hypothesis about the role of spatiotemporal structure in spiking activity; and, similarly, *2*) it can be used to test the validity of reduced models that have been proposed to explain this role.

The algorithm also has direct experimental applications. Since there are thousands of synapses per cell, the space of input patterns can be best described using statistical measures. Spatial and temporal correlations between synaptic inputs on the dendritic tree are examples of such measures. New methods developed for random access microscopy (Duemani Reddy et al. 2008) can be used to experimentally simulate almost arbitrary spatiotemporal patterns of synaptic input to a single neuron. The present algorithm provides a natural way of generating such patterns from their statistical description. With the possibility of generating controlled and complex spatiotemporal input patterns we could in principle study more in depth what makes a more realistic, multicompartment neuron model to fire and to assess whether the same conclusions would hold for simpler neuron models.

Simulated spike trains could also be useful in conjunction with more traditional electrophysiological techniques. For instance, the spatiotemporal structure of the inputs plays a significant role in the response of ubiquitous cortical microcircuits (Cruikshank et al. 2001; Pouille and Scanziani 2001). Synthetic spike trains could be used to probe the impact of spatiotemporal structure in such microcircuits and to experimentally study the propagation of coherent behavior through layered neural networks (Reyes 2003).

Local field potential oscillations in particular frequency bands are implicated in a variety of sensory, motor, and cognitive processes (Engel and Singer 2001; Engel et al. 2001; Fries 2005; Fries et al. 1997, 2001; Gray et al. 1989; Rickert et al. 2005; Roelfsema et al. 1997; Schoffelen et al. 2005; Siegel et al. 2008; Womelsdorf and Fries 2007; Womelsdorf et al. 2006, 2007). However, the sources of the LFP activity are usually unknown. In this way, it is usually difficult to understand the role of spike–LFP synchronization in information processing. We believe that progress can be made by simulating the experimentally observed statistical relationship between spikes and LFP. For instance, we could test specific conditions for efficient information transmission by parametrically manipulating the synchronization of spikes and LFPs and studying the impact in the firing patterns of neural networks. Specific hypotheses about the role of spike–LFP synchronization in different frequency bands could be obtained by comparing the simulations results with experimental data.

## GRANTS

This work was supported by a Texas Advanced Research Program/Advanced Technology Program grant to K. Josić and National Science Foundation Grants DMS-0604429 and DMS-0817649.

## ACKNOWLEDGMENTS

We thank R. Rosenbaum, H. Rey, and P. Berens for comments on the manuscript.

## Footnotes

↵1 All the processes that we use in simulations are discrete in time, although some are continuous and some discrete in space. Rather than using the phrases “continuous in space,” “discrete in time,” or “discrete in space and time,” we will refer to them as continuous and discrete processes, respectively.

↵2 For simplicity we fix the variance of

*y*to be 1 because it is just a scaling factor.

- Copyright © 2010 the American Physiological Society

## APPENDIX

Our method of generating spike trains requires that the continuous process covariance matrix be positive definite. Spike trains with absolute refractory periods cannot be generated with our algorithm because the associated continuous-process covariance matrix is nonpositive definite. However, the algorithm can generate relative refractory periods, close to absolute ones, for short refractory periods (see Fig. 6). Here, we derive the maximum possible refractoriness (i.e., how much the probability of firing a spike is reduced after a spike was elicited in the past) as a function of the length of the refractory period. The maximum possible refractoriness for a fixed refractory period length is obtained when the continuous-process covariance matrix is positive definite but it is as close as possible to be singular. This means that we have to minimize the correlation of the continuous variable *y*(*t*) and *y*(*t* + τ), with τ shorter than or equal to the refractory period length.

Let *m* be the length of the refractory period and *a* the minimum correlation value. The Toeplitz continuous-process covariance matrix *R* of size *N* × *N* is a banded matrix, where

The definition of a positive definite matrix requires that *x→Rx→ ^{t}* > 0 for every

*x*. If we take

_{i}*x*=

_{i}*x*for ∀

*i*

*N*→ ∞, the condition for the matrix to be positive definite is

*a*> −1/(2

*m*). This indicates that for longer refractory periods, the restriction of possible covariance matrices becomes tighter. Next, we demonstrate that

*a*> −1/(2

*m*) is also a sufficient condition when

*N*→ ∞. We can rewrite

*x→Rx→*> 0 as

^{t}The factor *k* incorporates the terms corresponding to the boundaries of the Toeplitz matrix, becoming negligible when *N* → ∞. If every term in the sum is >0, then *x→Rx→ ^{t}* > 0

*a*> −1/(2

*m*) ensures that this inequality holds for every

*x*,

_{i}*x*.

_{j}