A generative spike train model with time-structured higher order correlations

Emerging technologies are revealing the spiking activity in ever larger neural ensembles. Frequently, this spiking is far from independent, with correlations in the spike times of different cells. Understanding how such correlations impact the dynamics and function of neural ensembles remains an important open problem. Here we describe a new, generative model for correlated spike trains that can exhibit many of the features observed in data. Extending prior work in mathematical finance, this generalized thinning and shift (GTaS) model creates marginally Poisson spike trains with diverse temporal correlation structures. We give several examples which highlight the model's flexibility and utility. For instance, we use it to examine how a neural network responds to highly structured patterns of inputs. We then show that the GTaS model is analytically tractable, and derive cumulant densities of all orders in terms of model parameters. The GTaS framework can therefore be an important tool in the experimental and theoretical exploration of neural dynamics.


INTRODUCTION
Recordings across the brain suggest that neural populations spike collectively-the statistics of their activity as a group are distinct from that expected in assembling the spikes from one cell at a time (Bair et al., 2001;Salinas and Sejnowski, 2001;Harris, 2005;Averbeck et al., 2006;Schneidman et al., 2006;Shlens et al., 2006;Pillow et al., 2008;Ganmor et al., 2011;Bathellier et al., 2012;Hansen et al., 2012;Luczak et al., 2013). Advances in electrode and imaging technology allow us to explore the dynamics of neural populations by simultaneously recording the activity of hundreds of cells. This is revealing patterns of collective spiking that extend across multiple cells. The underlying structure is intriguing: For example, higher-order interactions among cell groups have been observed widely (Amari et al., 2003;Schneidman et al., 2006;Shlens et al., 2006Shlens et al., , 2009Ohiorhenuan et al., 2010;Ganmor et al., 2011;Vasquez et al., 2012;Luczak et al., 2013). A number of recent studies point to mechanisms that generate such higher-order correlations from common input processes, including unobserved neurons. This suggests that, in a given recording or given set of neurons projecting downstream, higher-order correlations may be quite ubiquitous (Barreiro et al., 2010;Macke et al., 2011;Yu et al., 2011;Köster et al., 2013). Moreover, these higher-order correlations may impact the firing statistics of downstream neurons (Kuhn et al., 2003), the information capacity of their output (Ganmor et al., 2011;Cain and Shea-Brown, 2013;Montani et al., 2013), and could be essential in learning through spike-time dependent synaptic plasticity (Pfister and Gerstner, 2006;Gjorgjieva et al., 2011).
What exactly is the impact of such collective spiking on the encoding and transmission of information in the brain? This question has been studied extensively, but much remains unknown. Results to date show that the answers will be varied and rich. Patterned spiking can impact responses at the level of single cells (Salinas and Sejnowski, 2001;Kuhn et al., 2003;Xu et al., 2012) and neural populations (Amjad et al., 1997;Tetzlaff et al., 2003;Rosenbaum et al., 2010Rosenbaum et al., , 2011. Neurons with even the simplest of non-linearities can be highly sensitive to correlations in their inputs. Moreover, such non-linearities are sufficient to accurately decode signals from the input to correlated neural populations (Shamir and Sompolinsky, 2004).
An essential tool in understanding the impact of collective spiking is the ability to generate artificial spike trains with a predetermined structure across cells and across time (Brette, 2009;Gutnisky and Josić, 2009;Krumin and Shoham, 2009;Macke et al., 2009). Such synthetic spike trains are the grist for testing hypotheses about spatiotemporal patterns in coding and dynamics. In experimental studies, such spike trains can be used to provide structured stimulation of single cells across their dendritic trees via glutamate uncaging (Gasparini and Magee, 2006;Reddy et al., 2008;Branco et al., 2010;Branco and Häusser, 2011). In addition, entire populations of neurons can be activated via optical stimulation of microbial opsins (Han and Boyden, 2007;Chow et al., 2010). Computationally, they are used to examine the response of non-linear models of downstream cells (Carr et al., 1998;Salinas and Sejnowski, 2001;Kuhn et al., 2003). Therefore, much effort has been devoted to developing statistical models of population activity. A number of flexible, yet tractable probabilistic models of joint neuronal activity have been proposed. Pairwise correlations are the most common type of interactions obtained from multi-unit recordings. Therefore, many earlier models were designed to generate samples of neural activity patterns with predetermined first and second order statistics (Brette, 2009;Gutnisky and Josić, 2009;Krumin and Shoham, 2009;Macke et al., 2009). In these models, higher-order correlations are not explicitly and separately controlled.
A number of different models have been used to analyze higher-order interactions. However, most of these models assume that interactions between different cells are instantaneous (or near-instantaneous) (Kuhn et al., 2003;Johnson and Goodman, 2009;Staude et al., 2010;Shimazaki et al., 2012). A notable exception is the work of Bäuerle and Grübel (2005), which developed such methods for use in financial applications. In these previous efforts, correlations at all orders were characterized by the increase, or decrease, in the probability that groups of cells spike together at the same time, or have a common temporal correlation structure regardless of the group.
The aim of the present work is to provide a statistical method for generating spike trains with more general correlation structures across cells and time. Specifically, we allow distinct temporal structure for correlations at pairwise, triplet, and all higher orders, and do so separately for different groups of cells in the neural population. Our aim to describe a model that can be applied in neuroscience, and can potentially be fit to emerging datasets.
A sample realization of a multivariate generalized thinning and shift (GTaS) process is shown in Figure 1. The multivariate spike train consists of six marginally Poisson processes. Each event was either uncorrelated with all other events across the population, or correlated in time with an event in all other spike trains. This model was configured to exhibit activity that cascades through a sequence of neurons. Specifically, neurons with larger index tend to fire later in a population wide event (this is similar to a synfire chain (Abeles, 1991), but with variable timing of spikes within the cascade). In Figure 1B, we plot the "population cross-cumulant density" for three chosen neurons-the summed activity of the population triggered by a spike in a chosen cell. The center of mass of this function measures the average latency by which spikes of the neuron in question precede those of the rest of the population (Luczak et al., 2013). Finally, Figure 1C shows the third-order cross-cumulant density for the three neurons. The triangular support of this function is a reflection of a synfirelike cascade structure of the spiking shown in the raster plot of panel (A): when firing events are correlated between trains, they tend to proceed in order of increasing index. We demonstrate the impact of such structured activity on a downstream network in section 2.2.3.

RESULTS
Our aim is to describe a flexible multivariate point process capable of generating a range of high order correlation structures. To do so, we extend the TaS (thinning and shift) model of temporally-and spatially-correlated, marginally Poisson counting processes (Bäuerle and Grübel, 2005). The TaS model itself generalizes the SIP and MIP models (Kuhn et al., 2003) which have been used in theoretical neuroscience (Tetzlaff et al., 2008;Rosenbaum et al., 2010;Cain and Shea-Brown, 2013). However, the TaS model has not been used as widely. The original TaS model is too rigid to generate a number of interesting activity patterns observed in multi-unit recordings (Ikegaya et al., 2004;Luczak et al., 2007Luczak et al., , 2013. We therefore developed the GTaS which allows for a more diverse temporal correlation structure. We begin by describing the algorithm for sampling from the GTaS model. This constructive approach provides an intuitive understanding of the model's properties. We then present a pair of examples, the first of which highlights the utility of the GTaS framework. The second example demonstrates how sample point processes from the GTaS model can be used to study population dynamics. Next, we present the analysis which yields the explicit forms for the cross-cumulant densities derived in the context of the examples. We do so by first establishing a useful distributional representation for the GTaS process, paralleling Bäuerle and Grübel (2005). Using this representation, we derive cross-cumulants of a GTaS counting process, as well as explicit expressions for the cross-cumulant densities. After explaining the derivation at lower orders, we present a theorem which describes cross-cumulant densities at all orders.

GTaS MODEL SIMULATION
The GTaS model is parameterized first by a rate λ which determines the intensity of a "mother process"-a Poisson process on R. The events of the mother process are marked, and the markings determine how each event is distributed among a collection of N daughter processes. The daughter processes are indexed by the set D = {1, . . . , N}, and the set of possible markings is the power set 2 D , the set of all subsets of D. We define a probability distribution p = (p D ) D ⊂ D , assigning a probability to each possible marking, D. As we will see, p D determines the probability of a joint event in all daughter processes with indices in the set D. Finally, to each marking, D, we assign a probability distribution Q D , giving a family of shift (jitter) The rate λ, the distribution p over the markings, and the family of jitter distributions (Q D ) D ⊂ D , define a vector X = (X 1 , . . . , X N ) of dependent daughter Poisson processes described by the following algorithm, which yields a single realization (see Figure 2): 1. Simulate the mother Poisson process of rate λ on R, generating a sequence of event times {t j }. (Figure 2A).

FIGURE 2 | An illustration of a GTaS simulation. (A)
Step 1: Simulate the mother process-a time-homogeneous Poisson process with event times Step 2: For each t j in step 1, select a set D j ⊂ D according to the distribution p D , and project the event at time t j to the subsets with indices in D j . The legend indicates the colors assigned to three possible markings in this example. (C) Step 3: For each pair (t j , D j ) generated in the previous two steps, draw Y j from Q D j , and shift the event times in the daughter processes by the corresponding values Y j i .
2. With probability p D j assign the subset D j ⊂ D to the event of the mother process at time t j . This event will be assigned only to processes with indices in D j . (Figure 2B).

Sample a vector
For each i ∈ D, the time t j + Y j i is set as an event time for the marginal counting process X i . (Figure 2C).
Hence copies of each point of the mother process are placed into daughter processes after a shift in time. A primary difference between the GTaS model and the TaS model presented in Bäuerle and Grübel (2005) is the dependence of the shift distributions Q D on the chosen marking. This allows for greater flexibility in setting the temporal cumulant structure.

Relation to SIP/MIP processes
Two simple models of correlated, jointly Poisson processes were defined in Kuhn et al. (2003). The resulting spike trains exhibit spatial correlations, but only instantaneous temporal dependencies. Each model was constructed by starting with independent Poisson processes, and applying one of two elementary point process operations: superposition and thinning (Cox and Isham, 1980). We show that both models are special cases of the GTaS model.
In the single interaction process (SIP), each marginal process X i is obtained by merging an independent Poisson process with a common, global Poisson process. That is, where Z c and each Z i are independent Poisson counting processes on R with rates λ c , λ i , respectively. An SIP model is equivalent to a GTaS model with mother process rate λ = λ c + N i = 1 λ i , and marking probabilities Note that if λ c = 0, each spike will be assigned to a different process X i , resulting in N independent Poisson processes. Lastly, each shift distribution is equal to a delta distribution at zero in every coordinate (i.e., q D (y 1 , . . . , y N ) ≡ N i = 1 δ(y i ) for every D ⊂ D). Thus, all joint cumulants (among distinct marginal processes) of orders 2 through N are delta functions of equal magnitude, λp D .
The multiple interaction process (MIP) consists of N Poisson processes obtained from a common mother process with rate λ m by thinning (Cox and Isham, 1980). The ith daughter process is formed by independent (across coordinates and events) deletion of events from the mother process with probability p = (1 − ). Hence, an event is common to k daughter processes with probability k . Therefore, if we take the perspective of retaining, rather than deleting events, the MIP model is equivalent to a GTaS process with λ = λ m , and p D = |D| (1 − ) d−|D| . As in the SIP case, the shift distributions are singular in every coordinate. Below, we present a general result (Theorem 1.1) which immediately yields as a corollary that the MIP model has cross-cumulant functions which are δ functions in all dimensions, scaled by k , where k is the order of the cross-cumulant.

Generation of synfire-like cascade activity
The GTaS framework provides a simple, tractable way of generating cascading activity where cells fire in a preferred order across the population-as in a synfire chain, but (in general) with variable timing of spikes (Abeles, 1991;Abeles and Prut, 1996;Aertsen et al., 1996;Aviel et al., 2002;Ikegaya et al., 2004). More generally, it can be used to simulate the activity of cell assemblies (Hebb, 1949;Harris, 2005;Buzsáki, 2010;Bathellier et al., 2012), in which the firing of groups of neurons is likely to follow a particular order.
In the Introduction, we briefly presented one example in which the GTaS framework was used to generate synfire-like cascade activity (see Figure 1), and we present another in Figure 3. In what follows, we will present the explicit definition of this second model, and then derive explicit expressions for its cumulant structure. Our aim is to illustrate the diverse range of possible correlation structures that can be generated using the GTaS model.
Consider an N-dimensional counting process X = (X 1 , . . . , X N ) of GTaS type, where N ≥ 4. We restrict the marking distribution so that p D ≡ 0 unless |D| ≤ 2 or D = D. That is, events are either assigned to a single, a pair, or all daughter processes. For sets D with |D| = 2, we set Q D ∼ N (0, )-a Gaussian distributions of zero mean and some specified covariance. The choice of the precise pairwise shift distributions is not important. Shifts of events attributed to a single process have no effect on the statistics of the multivariate process-this will become clear in section 2.3, when we exhibit that a GTaS process is equivalent in distribution to a sum of independent Poisson processes. In that context, the shifts of marginal events are applied to the event times of only one of these Poisson processes, which does not impact its statistics.
It remains to define the jitter distribution for events common to the entire population of daughter processes, i.e., events marked by D. We will show that we can generate cascading activity, and analytically describe the resulting correlation structure. We will say that a random variable T follows the exponential distribution Exp(α) if it has probability density where (t) is the Heaviside step function. We generate random vectors Y ∼ Q D according to the following rule, for each i = 1, . . . , N: In particular, note that these shift times satisfy Y N ≥ . . . ≥ Y 2 ≥ Y 1 ≥ 0, indicating the chain-like structure of these joint events.
From the definition of the model and our general result (Theorem 1.1) below, we immediately have that κ X ij (τ), the second Red spikes indicate population-wide events which have shift-times given by cumulative sums of independent exponentials, as described in the text. Arrows indicate the location of the first spike in the cascade. (B) A second-order cross-cumulant κ X 13 (black line) of this model is composed of contributions from two sources: correlations due to second-order markings, which have Gaussian shifts (c 2 13 -dashed red line), and correlations due to the occurrence of population wide events (c N 13 -dashed blue line). (C) Density plots of the third-order cross-cumulant density for triplets (i) (1, 2, 3) and (ii) (1, 2, 4)-the latter is given explicitly in Equation (6). System parameters are given in the Appendix.

Frontiers in Computational Neuroscience
www.frontiersin.org July 2013 | Volume 7 | Article 84 | 4 order cross-cumulant density for the process (i, j), is given by where define the contributions to the second order cross-cumulant density by the second-order, Gaussian-jittered events and the population-level events, respectively. Therefore, correlations between spike trains in this case reflect distinct contributions from second order and higher order events. The functions q D D indicate the densities associated with the distribution Q D , projected to the dimensions of D . All statistical quantities are precisely defined in the methods. By exploiting the hierarchical construction of the shift times, we can find an expression for the joint density q D , necessary to explicitly evaluate Equation (1). For a general N-dimensional distribution, Substituting this in to the identity Equation (3), we have Using Theorem 1.1 (Equation A8) we obtain the Nth order crosscumulant density (see the Methods), where, for notational convenience, we define τ 0 = 0. A raster plot of a realization of this model is shown in Figure 3A. We note that the cross-cumulant densities of arbitrary subcollections of the counting processes X can be obtained by finding the appropriate marginalization of q D via integration of Equation (4). In the case that common distributions are used to define the shifts, symbolic calculation environments (i.e., Mathematica) can quickly yield explicit formulas for cross-cumulant densities. Mathematica notebooks for Figure 1 available upon request. As a particular example, we consider the cross-cumulant density of the marginal processes X 1 , X 3 . Using Equations (2, 4), we find An expression for c 2 13 (τ) may be obtained similarly using Equation (2) and recalling that Q {i, j} ≡ N (0, ) for all i, j. In Figure 3B, we plot these contributions, as well as the full covariance density.
Similar calculations at third order yield, as an example, where the cross-cumulant density κ X 124 (τ 1 , τ 2 ) is supported only on τ 2 ≥ τ 1 ≥ 0. Plots of the third-order cross-cumulants for triplets (1, 2, 3) and (1, 2, 4) in this model are shown in Figure 3C. Note that, for the specified parameters, the conditional distribution of Y 4 -the shift applied to the events of X 4 in a joint population event-given Y 2 follows a gamma distribution, whereas Y 3 |Y 2 follows an exponential distribution, explaining the differences in the shapes of these two cross-cumulant densities.
General cross-cumulant densities of at least third order for the cascading model will have a form similar to that given in Equation (6), and will contain no signature of the correlation of strictly second order events. This highlights a key benefit of cumulants as a measure of dependence: although they agree with central moments up to third order, we know from Equation (23) below [or Equation (22) in the general case] that central moments necessarily exhibit a dependence on lower order statistics. On the other hand, cumulants are "pure" and quantify only dependencies at the given order which cannot be inferred from lower order statistics .
One useful statistic for analyzing population activity through correlations is the population cumulant density (Luczak et al., 2013). The second order population cumulant density for cell i is defined by (see the Methods) This function is linearly related to the spike-triggered average of the population activity conditioned on that of cell i. In Figure 4 Frontiers in Computational Neuroscience www.frontiersin.org July 2013 | Volume 7 | Article 84 | 5 we show three different second-order population-cumulant functions for the cascading GTaS model of Figure 3A. When the second order population cumulant for a neuron is skewed to the right of τ = 0 (as is κ X 1, pop -blue line), a neuron tends to precede other cells in the population in pairwise spiking events. Similarly, skewness to the left of τ = 0 (κ X 6, pop -orange line) indicates a neuron which tends to trail other cells in the population in such events. A symmetric population cumulant density indicates a neuron is a follower and a leader. Taken together, these three second order population cumulants hint at the chain structure of the process.
Greater understanding of the joint temporal statistics in a multivariate counting process can be obtained by considering higherorder population cumulant densities. We define the third-order population cumulant density for the pair (i, j) to be The third-order population cumulant density is linearly related to the spike-triggered population activity, conditioned on spikes in cells i and j separated by a delay τ 1 . In Figures 4B-D, we present three distinct third-order population cumulant densities.
Examining κ X 12, pop (τ 1 , τ 2 ) (panel B), we see only contributions in the region τ 2 > τ 1 > 0, indicating that the pairwise event 1 → 2 often precedes a third spike elsewhere in the population (i.e., with a probability above chance). The population cumulant κ X 34, pop (τ 1 , τ 2 ) has contributions in two sections of the plane (panel C). Contributions in the region τ 2 > τ 1 > 0 can be understood following the preceding example, while contributions in the region τ 2 < 0 < τ 1 imply that the firing of other neurons tends to precede the joint firing event 3 → 4. Lastly, contributions to κ X 16, pop (τ 1 , τ 2 ) (panel D) are limited to 0 < τ 2 < τ 1 , indicating an above chance probability of joint firing events of the form 1 → i → 6, where i indicates a distinct neuron within the population.
A distinct advantage of the study of population cumulant densities as opposed to individual cross-cumulant functions in Second order population cumulant densities for processes 1, 3, and 6. Greater mass to the right (resp. left) of τ = 0 indicates that a cell tends to lead (resp. follow) in pairwise-correlated events. (B) Third order population cumulant for processes X 1 , X 2 in the cascading GTaS process. Concentration of the mass in different regions of the plane indicates temporal structure of events correlated between X 1 , X 2 relative to the remainder of the population (see the text). (C) Same as (B), but for processes X 3 , X 4 . (D) Same as (B), but for processes X 1 , X 6 . System parameters are given in the Appendix.

Frontiers in Computational Neuroscience
www.frontiersin.org July 2013 | Volume 7 | Article 84 | 6 practical applications is related to data (i.e., sample size) limitations. In many practical applications, where the temporal structure of a collection of observed point processes is of interest, we often deal with a small, noisy samples. It may therefore be difficult to estimate third-or higher-order cumulants. Population cumulants partially circumvent this issue by pooling (Tetzlaff et al., 2003;Rosenbaum et al., 2010Rosenbaum et al., , 2011) (or summing) responses, to amplify existing correlations and average out the noise in measurements. We conclude this section by noting that even cascading GTaS examples can be much more general. For instance, we can include more complex shift patterns, overlapping subassemblies within the population, different temporal processions of the cascade, and more.

Timing-selective network
The responses of single neurons and neuronal networks in experimental (Meister and Berry II, 1999;Singer, 1999;Bathellier et al., 2012) and theoretical studies (Jeffress, 1948;Hopfield, 1995;Joris et al., 1998;Thorpe et al., 2001;Gütig and Sompolinsky, 2006) can reflect the temporal structure of their inputs. Here, we present a simple example that shows how a network can be selective to fine temporal features of its input, and how the GTaS model can be used to explore such examples.
As a general network model, we consider N leaky integrateand-fire (LIF) neurons with membrane potentials V i obeying When the membrane potential of cell i reaches a threshold Vth, an output spike is recorded and the membrane potential is reset to zero, after which evolution of V i resumes the dynamics in Equation (7). Here w ij is the synaptic weight of the connection from cell j to i, w in is the input weight, and we assume time to be measured in units of membrane time constants. The function is a delayed, unit-area exponential synaptic kernel with time-constant τ syn and delay τ d . The output of the ith neuron is where the event times {s j i } correspond to those of a GTaS counting process X. Thus, each input spike results in a jump in the membrane potential of the corresponding LIF neuron of amplitude w in . The particular network we consider will have a ring topology (nearest neighbor-only connectivity)-specifically, for i, j = 1, . . . , N, we let We further assume that all neurons are excitatory, so that w syn > 0. A network of LIF neurons with synaptic delay is a minimal model which can exhibit fine-scale discrimination of temporal patterns of inputs without precise tuning (Izhikevich, 2006) (that is, without being carefully designed to do so, with great sensitivity to modification of network parameters). To exhibit this dependence we generate inputs from two GTaS processes. The first (the cascading model) was described in the preceding example. To independently control the mean and variance of relative shifts we replace the sum of exponential shifts with sums of gamma variates. We also consider a model featuring population-level events without shifts (the synchronous model), where the distribution Q D is a δ distribution at zero in all coordinates.
The only difference between the two input models is in the temporal structure of joint events. In particular, the rates, and all long timescale spike count cross-cumulants (equivalent to the total "area" under the cross-cumulant density, see the Methods) of order two and higher are identical for the two processes. We focus on the sensitivity of the network to the temporal cumulant structure of its inputs.
In Figures 5A,B, we present two example rasters of the nearestneighbor LIF network receiving synchronous (left) and cascading (right) input. In the second case, there is an obvious pattern in the outputs, but the firing rate is also increased. This is quantified in Figure 5C, where we compare the number of output spikes fired by a network receiving synchronous input (horizontal axis) with the same for a network receiving cascading input (vertical axis), over a large number of trials. On average, the cascading input increases the output rate by a factor of 1.5 over the synchronous inputs-we refer to this quantity as the cascade amplification factor (CAF).
Finally, in Figure 5D, we illustrate how the cascade amplification factor depends on the parameters that define the timing of spikes for the cascading inputs. First, we study the dependence on the standard deviation σ shift of the gamma variates determining the shift distribution. We note that amplification factors above 1.5 hold robustly (i.e., for a range of shift σ shift values). The amplification factors decrease with shift variance. In the inset to panel (D), we show how the gain depends on the mean of the shift distribution μ shift . On an individual trial, the response intensity will depend strongly on the total number of input spikes. Thus, in order to enforce a fair comparison, the mother process and markings used were identical in each trial of every panel of Figure 5. We note that network properties, such as the membrane properties of individual cells or synaptic timescales, may have an equally large impact on the cascade amplification factor-indeed, as we explain below, the observed behavior of the CAF is a result of synergy between the timescales of input and interactions within the network. These observations have simple explanations in terms of the network dynamics and input statistics. Neglecting, for a moment, population-level events, the network is configured so that correlations in activity decrease with topographic distance. Accordingly, the probability of finding neurons that are simultaneously close to threshold also decreases with distance. Under the synchronous input model, a population-level event results in a simultaneous increase of the membrane potentials of all neurons by an amount w in , but unless the input is very strong (in which case every, or almost every, neuron will fire regardless of finescale input structure), the set of neurons sufficiently close to threshold to "capitalize" on the input and fire will typically be restricted to a topographically adjacent subset. Neurons which do not fire almost immediately will soon have forgotten about this population-level input. As a result, the output does not significantly reflect the chain-like structure of the inputs (Figure 5A, right).
On the other hand, in the case of the cascading input, the temporal structure of the input and the timescale of synapses can operate synergistically. Consider a pair of adjacent neurons in the ring network, called cells 1 and 2, arranged so that cell 2 is downstream from cell 1 in the direction of the population-level chain events. When cell 1 spikes, it is likely that cell 2 will also have an elevated membrane potential. The potential is further elevated by the delayed synaptic input from cell 1. If cell 1 spikes in response to a population-level chain event, then cell 2 imminently receives an input spike as well. If the synaptic filter and time-shift of the input spikes to each cell align, then the firing probability of cell 2 will be large relative to chance. This reasoning can be carried on across the network. Hence synergy between the temporal structure of inputs and network architecture allows the network to selectively respond to the temporal structure of the inputs (Figure 5B, right).
In Kuhn et al. (2003), the effect of higher order correlations on the firing rate gain of an integrate-and-fire neuron was studied by driving single cells using sums of SIP or MIP processes with equivalent firing rates (first order cumulants) and pairwise correlations (second order cumulants). In contrast, in the preceding example, the two inputs have equal long time spike count cumulants, and differ only in temporal correlation structure. An increase in firing rate was due to network interactions, and is therefore a population level effect. We return to this comparison in the Discussion.
These examples demonstrate how the GTaS model can be used to explore the impact of spatio-temporal structure in population activity on network dynamics. We next proceed with a formal derivation of the cumulant structure for a general GTaS process.

CUMULANT STRUCTURE OF A GTaS PROCESS
The GTaS model defines an N-dimensional counting process. Following the standard description for a counting process, X = (X 1 , . . . , X N ) on R N , given a collection of Borel subsets A i ∈ B(R), i = 1, . . . , N, then X(A 1 × · · · × A N ) = (X 1 (A 1 ), . . . , X N (A N )) ∈ N N is a random vector where the value of each coordinate i indicates the (random) number of points which fall inside the set A i . Note that the GTaS model defines processes that are marginally Poisson. All GTaS model parameters and related quantities are defined in Table 1.

Frontiers in Computational Neuroscience
www.frontiersin.org July 2013 | Volume 7 | Article 84 | 8 Population cumulant density defined in Equation (25) For each D ⊂ D = {1, . . . , N}, define the tail probabilityp D bȳ Since p D is the probability that exactly the processes in D are marked,p D is the probability that all processes in D, as well as possibly other processes, are marked. An event from the mother process is assigned to daughter process X i with probabilityp {i} . As noted above, an event attributed to process i following a marking D i will be marginally shifted by a random amount determined by the distribution Q {i} D which represents the projection of Q D onto dimension i. Thus, the events in the marginal process X i are shifted in an independent and identically distributed (IID) manner according to the mixture distribution Q i given by Note that IID shifting of the event times of a Poisson process generates another Poisson process of identical rate. Thus, the process X i is marginally Poisson with rate λp {i} (Ross, 1995).
In deriving the statistics of the GTaS counting process X, it will be useful to express the distribution of X as Here, each ξ(D; A 1 , . . . , A N ) is an independent Poisson process, and the notation = distr indicates that the two random vectors are equal in distribution. This process counts the number of points which are marked by a set D ⊃ D, but (after shifting) only the points with indices i ∈ D lie in the corresponding set A i . Precise definitions of the processes ξ and a proof of Equation (9) may be found in the Appendix. We emphasize that the Poisson processes ξ(D) do not directly count points marked for the set D, but instead points which are marked for a set containing D that, after shifting, only have their D-components lying in the "relevant" sets A i . Suppose we are interested in calculating dependencies among a subset of daughter processes, {X i j } i j ∈D for some setD ⊂ D, consisting of |D| = k distinct members of the collection of counting processes X. Then the following alternative representation will be useful: ⎛ ⎜ ⎝ We illustrate this decomposition in the cases k = 2, 3 in Figure 6. The sums in Equation (10) run over all sets D ⊂ D containing the indicated indices i j and contained withinD. The processes ζ D are comprised of a sum of all of the processes ξ(D ) (defined below Equation 9) such that D contains all of the indices D, but no other indices which are part of the subsetD under consideration. These sums are non-overlapping, implying that the ζ D are also independent and Poisson.
The following examples elucidate the meaning and significance of Equation (10). We emphasize that the GTaS process is a completely characterized, joint Poisson process, and we use Equation (10) to calculate cumulants of a GTaS process. In principle, any other statistics can be obtained similarly.

Second order cumulants (covariance)
We first generalize a well-known result about the dependence structure of temporally jittered pairs of Poisson processes, X 1 , X 2 . Assume that events from a mother process with rate λ, are assigned to two daughter processes with probability p. Each event time is subsequently shifted independently according to a univariate distribution f . The cross-cumulant density (or crosscovariance function; see the Methods for cumulant definitions) then has the form (Brette, 2009) We generalize this result within the GTaS framework. At second order, Equation (10) has a particularly nice form. Following Bäuerle and Grübel (2005) we write for i = j (see Figure 6A) The process ζ {i, j} sums all ξ(D ) for which {1, 2} ⊂ D , while the process ζ {i} sums all ξ(D ) such that i ∈ D , j / ∈ D , and ζ {j} is defined likewise.
Using the representation in Equation (11), we can derive the second order cumulant (covariance) structure of a GTaS process.
The third equality follows from the construction of the processes ζ D : if D = D , then the processes ζ D , ζ D are independent. The final equality follows from the observation that every cumulant of a Poisson random variable equals its mean.
The covariance may be further expressed in terms of model parameters (see Theorem 1.1 for a generalization of this result to arbitrary cumulant orders): In other words, the covariance of the counting processes is given by the weighted sum of the probabilities that the (i, j) marginal of the shift distributions yield values in the appropriate sets. The weights are the intensities of each corresponding component processes ξ(D) which contribute events to both of the processes i and j.
In the case that Q D ≡ Q, Equation (12) reduces to the solution given in Bäuerle and Grübel (2005). Using the tail probabilities defined in Equation (8), if Q D ≡ Q for all D, the integral in Equation (12) no longer depends on the subset D , and the equation may be written as Using Equation (12), we may also compute the second crosscumulant density (also called the covariance density) of the processes. From the definition of the cross-cumulant density [Equation (24) in the Methods], this is given by Before continuing, we note that given a random vector Y = (Y 1 , . . . , Y N ) ∼ Q, where Q has density q(y 1 , . . . , y N ), the vector Z = (Y 2 − Y 1 , . . . , Y N − Y 1 ) has density q Z given by Assuming that the distributions Q D have densities q D , and denoting by q According to Equation (14), the integrals present in Equation (15) are simply the densities of the variables Thus κ X ij (τ), which captures the additional probability for events in the marginal processes X i and X j separated by τ units of time beyond what can be predicted from lower order statistics is given by a weighted sum (in this case, the lower order statistics are marginal intensities-see the discussion around Equation (24) of the Methods). The weights are the "marking rates" λp D for markings contributing events to both component processes, while the summands are the probabilities that the corresponding shift distributions yield a pair of shifts in the proper arrangementspecifically, the shift applied to the event as attributed to X i precedes that applied to the event mapped to X j by τ units of time. This interpretation of the cross-cumulant density is quite natural, and will carry over to higher order cross-cumulants of a GTaS process. However, as we show next, this extension is not trivial at higher cumulant orders.

Third order cumulants
To determine the higher order cumulants for a GTaS process, one can again use the representation given in Equation (10). The distribution of a subset of three processes may be expressed in the form (see Figure 6B) where, for simplicity, we suppressed the arguments of the different ζ D on the right hand side. Again, the processes in the representation are independent and Poisson distributed. The variable ζ {i, j, k} is the sum of all random variables ξ(D) (see Equation 9) with D ⊃ {i, j, k}, while the variable ζ {i, j} is now the sum of all ξ(D) with D ⊃ {i, j}, but k / ∈ D. The rest of the variables are defined likewise. Using properties (C1) and (C2) of cumulants given in the Methods, and assuming that i, j, k are distinct indices, we have The second equality follows from the fact that all cumulants of a Poisson distributed random variable equal its mean. Similar to Equation (12), we may write The third cross-cumulant density is then given similarly to the second order function by Here, we have again assumed the existence of densities q D , and denoted by q {i, j, k} D the joint marginal density of the variables Y i , Y j , Y k under q D . The integrals appearing in the expression for the third order cross-cumulant density are the probability densities of the vectors

General cumulants
Finally, consider a general subset of k distinct members of the vector counting process X as in Equation (10). The following theorem provides expressions for the cross-cumulants of the counting processes, as well as the cross-cumulant densities, in terms of model parameters in this general case. The proof of Theorem 1.1 is given in the Appendix.
where YD represents the projection of the random vector Y onto the dimensions indicated by the members of the setD. Furthermore, assuming that the shift distributions possess densities (q D ) D ⊂ D , the cross-cumulant density is given by where qD D indicates the kth order joint marginal density of q D in the dimensions ofD.
An immediate corollary of Theorem 1.1 is a simple expression for the infinite-time-window cumulants, obtained by integrating the cumulant density across all time lags τ i . From Equation (A8), we have This shows that the infinite time window cumulants for a GTaS process are non-increasing with respect to the ordering of sets, i.e., We conclude this section with a short technical remark: Until this point, we have considered only the cumulant structure of sets of unique processes. However occasionally, one may wish to calculate a cumulant for a set of processes including repeats. Take, for example, a cumulant κ(X 1 (A 1 ), X 1 (A 2 ), X 3 (A 3 )). Owing to the marginally Poisson nature of the GTaS process, we would have (referring to the Methods for cumulant definitions) For a general counting process X, it may be shown that κ X 113 (τ 1 , τ 2 ) = δ(τ 1 )κ X 13 (τ 2 ) + "non-singular contributions". (21) In addition, the second order auto-cumulant density may be written (Cox and Isham, 1980) κ X ii (τ) = r i δ(τ) + "non-singular contributions", where r i is the stationary rate. The singular contribution shown in Equation (21) at third order is in analogy to the delta contribution proportional to the firing rate which appears in the second-order auto-cumulant density. For a GTaS process, the non-singular contributions in Equation (21) are identically zero, following directly from Equation (20). Expressions similar to Equations (20, 21) hold for general cases.

DISCUSSION
We have introduced a general method of generating spike trains with flexible spatiotemporal structure. The GTaS model is completely analytically tractable: all statistics of interest can be obtained directly from the distributions used to define it. It is based on an intuitive method of selecting and shifting point processes from a "mother" train. Moreover, the GTaS model can be used to easily generate partially synchronous states, cluster firing, cascading chains, and other spatiotemporal patterns of neural activity. Processes generated by the GTaS model are naturally described by cumulant densities of pairwise and higher orders. This raises the question of whether such statistics are readily computable from data, so that realistic classes of GTaS models can be defined in the first place. One approach is to fit mechanistic models to data, and to use the higher order structure that is generated by the underlying mechanisms (Yu et al., 2011). A synergistic blend of other methods with the GTaS framework may also be fruitful-for example, the CuBIC framework of Staude et al. (2010) could be used to determine relevant marking orders, and the parametrically-described GTaS process could then be fit to allow generation of surrogate data after selection of appropriate classes of shift distributions. When it is necessary to infer higher order structure in the face of data limitations, population cumulants are an option to increase statistical power (albeit at the cost of spatial resolution; see Figure 4).
While the GTaS model has flexible higher order structure, it is always marginally Poisson. While throughout the cortex spiking is significantly irregular (Holt et al., 1996;Shadlen and Newsome, 1998), the level of variability differs across cells, with Fano factors ranging from below 0.5 to above 1.5-in comparison with the Poisson value of 1 (Churchland et al., 2010). Changes in variability may reflect cortical states and computation (Litwin-Kumar and Doiron, 2012;White et al., 2012). A model that would allow flexible marginal variability would therefore be very useful. Unfortunately, the tractability of the GTaS model is closely related to the fact that the marginal processes are Poisson. Therefore, an immediate generalization does not seem possible.
A number of other models have been used to describe population activity. Maximum entropy (ME) approaches also result in models with varied spatial activity; these are defined based on moments or other averaged features of multivariate spiking activity (Schneidman et al., 2006;Roudi et al., 2009). Such models are often used to fit purely spatial patterns of activity, though (Tang et al, 2008;Marre et al., 2009) have extended the techniques to treat temporal correlations as well. Generalized linear models (GLMs) have been used successfully to describe spatiotemporal patterns at second (Pillow et al., 2008), and third order (Ohiorhenuan et al., 2010). In comparison to the present GTaS method, both GLMs and ME models are more flexible. They feature well-defined approaches for fitting to data, including likelihood-based methods with well-behaved convexity properties. What the GTaS method contributes is an explicit way to generate population activity with explicitly specified high order spatio-temporal structure. Moreover, the lower order cumulant structure of a GTaS process can be modified independently of the higher order structure, though the reverse is not true.
There are a number of possible implications of such spatiotemporal structure for communication within neural networks. In section 2.2.3, we showed that these temporal correlations can play a role similar to that of spatial correlations established in Kuhn et al. (2003) for determining network input-output transfer. Our model allowed us to examine that impact of such temporal correlations on the network-level gain of a downstream population (cascade amplification factor). Even in a very simple network it was clear that the strength of the response is determined jointly by the temporal structure of the input to the network, and the connectivity within the network. Kuhn et al.

Frontiers in Computational Neuroscience
www.frontiersin.org July 2013 | Volume 7 | Article 84 | 12 examined the effect of higher order structure on the firing rate gain of an integrate-and-fire neuron by driving it with a mixture of SIP or MIP processes (Kuhn et al., 2003). However, in these studies, only the spatial structure of higher order activity was varied. The GTaS model allows us to concurrently change the temporal structure of correlations. In addition, the precise control of the cumulants allows us to derive models which are equivalent up to a certain cross-cumulant order, when the configuration of marking probabilities and shift distributions allow it (as for the SIP and MIP processes of Kuhn et al. (2003), which are equivalent at second order). Such patterns of activity may be useful when experimentally probing dendritic information processing (Gasparini and Magee, 2006), synaptic plasticity (Pfister and Gerstner, 2006;Gjorgjieva et al., 2011), or investigating the response of neuronal networks to complex patterns of input (Kahn et al., 2013). Spatiotemporal patterns may also be generated by cell assemblies (Bathellier et al., 2012). The firing in such assemblies can be spatially structured, and this structure may not be reflected in the activity of participating cells. Assemblies can exhibit persistent patterns of firing, sometimes with millisecond precision (Harris et al., 2002). The GTaS framework is well suited to describe exactly such activity patterns. The examples we presented can be easily extended to generate more complex patterns of activity with overlapping cell assemblies, different cells leading the activity, and other variations.
Understanding impact of spatiotemporal patterns on neural computations remains an open and exciting problem. Progress will require coordination of computational, theoretical, and experimental work-the latter taking advantage of novel stimulation techniques. We hope that the GTaS model, as a practical and flexible method for generating high-dimensional, correlated spike trains, will play a significant role along the way.

CUMULANTS AS A MEASURE OF DEPENDENCE
We first define cross-cumulants (also called joint cumulants) (Stratonovich and Silverman, 1967;Kendall et al., 1969;Gardiner, 2009) and review some important properties of these quantities. Define the cumulant generating function g of a random vector X = (X 1 , . . . , X N ) by The r-cross-cumulant of the vector X is given by where r = (r 1 , . . . , r N ) is a N-vector of positive integers, and |r| = N i = 1 r i . We will generally deal with cumulants where all variables are considered at first order, without excluding the possibility that some variables are duplicated. In this case, we define the cross-cumulant κ(X), of the variables in the random vector X = (X 1 , . . . , X N ) as κ(X) := κ 1 (X) = ∂ N ∂t 1 · · · ∂t N g(t 1 , . . . , t N ) (1, . . . , 1).
This relationship may be expressed in combinatorial form: where π runs through all partitions of D = {1, . . . , N}, and B runs over all blocks in a partition π. More generally, the r-crosscumulant may be expressed in terms of moments by expanding the cumulant generating function as a Taylor series, noting that similarly expanding the moment generating function M(t) = e g (t) , and matching the polynomial coefficients. Note that the nth cumulant κ n of a random variable X may be expressed as a joint cumulant via κ n (X) = κ(X, . . . , X) n copies of X .
We will utilize the following two principal properties of cumulants (Brillinger, 1965;Stratonovich and Silverman, 1967;Mendel, 1991;Staude et al., 2010): (C1) Multilinearity -for any random variables X, Y, {Z i } N i = 2 , we have κ(aX + bY, Z 2 , . . . , Z N ) = aκ(X, Z 2 , . . . , Z N ) This holds regardless of dependencies amongst the random variables. (C2) If any subset of the random variables in the cumulant argument is independent from the remaining, the crosscumulant is zero-i.e., if {X 1 , . . . , X N 1 } and {Y 1 , . . . , Y N 2 } are sets of random variables such that each X i is independent from each Y j , then To exhibit another key property of cumulants, consider a 4vector X = (X 1 , X 2 , X 3 , X 4 ) with non-zero fourth cumulant and a random variable Z independent of each X i . Define Y = (X 1 + Z, X 2 + Z, X 3 + Z, X 4 ). Using properties (C1), (C2) above, it follows that On the other hand, it is also true that that is, adding the variable Z to only a subset of the variables in X results in changes to cumulants involving only that subset, but not to the joint cumulant of the entire vector. In this sense, an rth order cross-cumulant of a collection of random variables captures exclusively dependencies amongst the collection which cannot be described by cumulants of lower order. In the example above, only the joint statistical properties of a subset of X were changed. As a result, the total cumulant κ(X) remained fixed. From Equation (22), it is apparent that κ(X i ) = E[X i ], and κ(X i , X j ) = cov X i , X j . In addition, the third cumulant, like the second, is equal to the corresponding central moment: As cumulants and central moments agree up to third order, central moments up to third order inherit the properties discussed above at these orders. On the other hand, the fourth cumulant is not equal to the fourth central moment. Rather: Higher cumulants have similar (but more complicated) expansions in terms of central moments. Accordingly, central moments of fourth and higher order do not inherit properties (C1), (C2).

TEMPORAL STATISTICS OF POINT PROCESSES
In the Results, we present an extension of previous work (Bäuerle and Grübel, 2005) in which we construct and analyze multivariate counting processes X = (X 1 , . . . , X N ) where each X i is marginally Poisson.
Formally, a counting process X is an integer-valued random measure on B(R N ). Evaluated on subset A 1 × · · · × A N of B(R N ), the random vector (X 1 (A 1 ), . . . , X N (A N )) counts events in d distinct categories whose times of occurrence fall in to the sets A i . A good general reference on the properties of counting processes (marginally Poisson and otherwise) is Daley and Vere-Jones (2002).
The assumption of Poisson marginals implies that for a set A i ∈ B(R), the random variable X i (A i ) follows a Poisson distribution with mean λ i (A i ), where is the Lebesgue measure on R, and λ i is the (constant) rate for the ith process. The processes under consideration will further satisfy a joint stationarity condition, namely that the distribution of the vector (X 1 (A 1 + t), . . . , X N (A N + t)) does not depend on t, where A i + t denotes the translated set {a + t : a ∈ A i }.
We now consider some common measures of temporal dependence for jointly stationary vector counting processes. We will refer to the quantity X i [0, T] as the spike count of process i over [0, T]. The quantity γ X i 1 ···i k (T) (which we will refer to as a spike count cumulant) is given by measures kth order correlations amongst spike counts for the listed processes which occur over windows of length T. At second order, γ X ij (T) measures the covariance of the spike counts of processes i, j over a common window of length T. The infinite window spike count cumulant quantifies dependencies in the spike counts of point processes over arbitrarily long windows, and is given by A related measure is the kth order cross-cumulant density The cross-cumulant density should be interpreted as a measure of the likelihood-above what may be expected from knowledge of the lower order cumulant structure-of seeing events in processes i 2 , . . . , i k at times τ 1 + t, . . . , τ k − 1 + t, conditioned on event in process i 1 at time t. The infinite window spike count cumulant is equal to the total integral under the cross-cumulant density, As an example, we again consider the familiar second-order cross-cumulant density κ X ij (τ)-often referred to as the crosscovariance density or cross-correlation function. Defining the conditional intensity h ij (τ) of process j, conditioned on process i to be that is, the intensity of j conditioned on an event in process i which occurred τ units of time in the past, then it is not difficult to show that That is, the second order cross-cumulant density supplies the probability of chance of observing an event attributed to process i, followed by one attributed to process j, τ units of time later, above what would be expected from knowledge of first order statistics (given by the product of the marginal intensities, λ i λ j ). More generally, at higher orders, the cross-cumulant density should be interpreted as a measure of the likelihood (above what may be expected from knowledge of the lower order correlation structure) of seeing events attribute to processes i 2 , . . . , i k at times τ 1 + t, . . . , τ k − 1 + t, conditioned on an event in process i 1 at time t.
Another statistic useful in the study of a correlated vector counting process X is the population cumulant density. At second-order, the population cumulant density for X i takes the form (Luczak et al., 2013) More generally, the kth order population cumulant density corresponding to the processes X i 1 , . . . , X i k − 1 is given by

APPENDIX PROOF OF THE DISTRIBUTIONAL REPRESENTATION OF THE GTaS MODEL IN EQUATION (9)
The construction of the GTaS model allows us to provide a useful distributional representation of the process. We describe this representation in a theorem that generalizes Theorem 1 in Bäuerle and Grübel (2005). This theorem also immediately implies that the GTaS process is marginally Poisson. Some definitions are required: first, for subsets A 1 , . . . , In addition, setting 1 = (1, . . . , 1) to be the N-dimensional vector with all components equal to unity, and if Q D is a measure on R N , then we define the measure ν(Q D ) by The measure ν(Q D ) can be interpreted as giving the expected Lebesgue measure of the subset L of R for which uniform shifts by the elements of L translate a random vector Y ∼ Q D in to A. Heuristically, one may imagine sliding the vector Y over the whole real line, and counting the number of times every coordinate ends up in the "right" set-the projection of A onto that dimension. In equation form, this means where the subscript Y indicates that we take the average over the distribution of Y ∼ Q D . A short proof of this representation is presented below. We now present the theorem, with a proof indicating adjustments necessary to that of Bäuerle and Grübel (2005).
Theorem 0 Let X be an N-dimensional counting process of GTaS type with base rate λ, thinning mechanism p = (p D ) D ⊂ D , and family of shift distributions (Q D ) D ⊂ D . Then, for any Borel subsets A 1 , . . . , A N of the real line, we have the following distributional representation: where the random variables ξ(D; A 1 , . . . , A N ), ∅ = D ⊂ D, are independent and Poisson distributed with Proof. For each marking D ⊂ D, define X D to be an independent TaS (Bäuerle and Grübel, 2005) counting process with mother process rate λp D , shift distribution Q D , and markings (p D D ) D ⊂ D where p D D = 1 if D = D and is zero otherwise (i.e., the only possible marking for X D is D ). We first claim that To see this, note that spikes in the mother process of the GTaS process of X marked for a set D occur at a rate λp D , which is the rate of the process X D . In addition, these event times are then shifted by Q D , exactly as they are for X D . Thus, the distribution of event times (and hence the counting process distributions) are equivalent.
Let A 1 , . . . , A N be any Borel subsets of the real line. Applying Theorem 1 of Bäuerle and Grübel (2005) to each X D gives the following distributional representation: ⎛ where the random variables ξ D (D; , A 1 , . . . , A N ) are taken to be identically zero unless D ⊂ D . In the latter case, they are independent and Poisson distributed with The second equality above follows from the fact that p D D = 1 if D = D and is zero otherwise.
Next, define As the sum of independent Poisson variables is again Poisson with rate equal to the sum of the rates, we have that ξ (D; A 1   Finally, combining Equations (A4, A5), we may write which, along with Equation (A6), establishes the theorem.
A short note: The variable ξ(D; A 1 , . . . , A N ) counts the number of points which are marked by a set D ⊃ D, but after shifting, only the points attributed to the processes with indices i ∈ D remain in the corresponding subsets A i . Thus, to determine the number of points attributed to the ith process which lie in A i (X i (A i )), one simply sums the variables ξ for all D containing i, as in Equation (A3). Thus, the intensity of ξ(D; A 1 , . . . , A N ), is simply the expected number of such points. Keeping in mind these natural interpretations of terms, Theorem 1 is easier to digest, and the result is not surprising.

PROOF OF EQUATION (27)
In Equation (A2), we gave a more intuitive representation of the measure ν(Q D ) than the one first defined in Bäuerle and Grübel (2005), which we prove here. Suppose that Q is a measure on B(R d ), and A ∈ B(R d ). Then we have where YD represents the projection of the random vector Y onto the dimensions indicated by the members of the setD. Furthermore, assuming that the shift distributions possess densities (q D ) D ⊂ 2 D , the cross-cumulant density is given by where qD D indicates the kth order joint marginal density of q D in the dimensions ofD.
Proof. First, as noted in the text, we may rewrite the distributional representation of Theorem 0 (Equation A3) as . . . , A N ) . . .
Repeating the description from the main text, the processes ζ D are comprised of a sum of all of the processes ξ(D ) (defined above, in Theorem 0) such that D contains all of the indices D, but no other indices which are part of the subsetD under consideration. These sums are non-overlapping, implying that the ζ D are also independent and Poisson. Using the representation of Equation (A9), we first find that where we suppressed the dependence of the variables ζ D on the subsets A i . The first equality in the previous equation is simply the representation defined in Equation (A10), and the second is from the multilinear property of cumulants (property (C1) in the Methods). Note that the sums are over the sets D 1 , . . . , D k satisfying the given conditions. Recall that, by construction, the Poisson processes ζ D (see Equation A10) are independent for distinct marking sets. Accordingly, the cumulant κ[ζ D 1 , . . . , ζ D k ] is zero unless D 1 = . . . = D k , by property (C2) of cumulants-that is, κ[ζ D 1 (A 1 , . . . , A N  This completes the proof of Equation (A7), and (A8) follows from the definition of the cross-cumulant density in Equation (24) of the Methods. Figure 1. For Figure 1, the GTaS process of size N = 6 consisted of only first order and population-level events which were assigned marking probabilities . The rate of the mother process was λ = 0.5 kHz, and the shift times for population level events were generated as in section 2.2.2 with

Parameters for figures in the text
where the Gamma distribution has density Figures 3, 4. For Figures 3, 4, the GTaS process of size N = 6 consisted of first and second order as well as population-level events. These events had marking probabilities . The rate of the mother process was λ = 0.5 kHz, and the shift times for population level events were generated as in section 2.2.2 with T i ∼ Exp(0.5), i = 1, . . . , 6.
The shift times of the second order events were drawn from an independent Gaussian distribution with each coordinate having standard deviation 5 ms. The rate of the mother process was λ = 1.5 kHz, and the shift times for population level events were generated as in section 2.2.2 with Frontiers in Computational Neuroscience www.frontiersin.org July 2013 | Volume 7 | Article 84 | 20 T i ∼ (k, θ), i = 1, . . . , 6.
The shift parameters k, θ (representing shape and scale) were determined by the given shift mean μ shift and standard deviation σ shift as μ shift = kθ, σ shift = kθ 2 .
The shift times of the second order events were drawn from an independent Gaussian distribution with each coordinate having standard deviation 0.3 ms.