where:

$$G(\mathbf{x}, \mathbf{y}) = \frac{1}{\sqrt{1 - c^2} \pi} e^{-\frac{1}{1 - c^2} (\mathbf{x}^2 + \mathbf{y}^2 - 2c\mathbf{x}\mathbf{y})} \tag{41}$$

and *R*(φ) remains to be determined. Note that *G*(*x*, *y*) is just the standard stationary solution to the multivariate OU equation. At this juncture, we remark that our main goal is to find *R*(φ), which is the marginal density of the phase differences between the two oscillators.

#### *Solving Equation 38*

Both Equations 38 and 39 have the form *L*0ρ = *b*(*x*, *y*, θ). *L*0 operates on the space of functions defined on *R*2 × *S*1 that are twice continuously differentiable in *x*, *y* and continuously differentiable in θ. In this space, *L*0 has a one-dimensional nullspace spanned by *G*(*x*, *y*) (constant in θ) and so *L*0 is not invertible. However, *L*0ρ(*x*, *y*, θ) = *b*(*x*, *y*, θ) does have a solution provided that *b*(*x*, *y*, θ) is orthogonal to the null space of *L*∗ 0, which is the adjoint linear operator of *L*0. Since *L*0 is a standard probability operator, its adjoint is always 1 (i.e., the function that is 1 everywhere).

Since:

$$b\_1(\mathbf{x}, \mathbf{y}, \boldsymbol{\theta}) = \frac{\mathbf{x} \mathbf{G}(\mathbf{x}, \mathbf{y})}{2\pi} [\Delta\_1'(\boldsymbol{\theta}) \mathbf{R}(\boldsymbol{\phi}) - \Delta\_1(\boldsymbol{\theta}) \mathbf{R}'(\boldsymbol{\phi})] + \boldsymbol{\mathcal{y}} \mathbf{G}(\mathbf{x}, \mathbf{y}) [\Delta\_2'(\boldsymbol{\theta} + \boldsymbol{\phi}) \mathbf{R}(\boldsymbol{\phi}) + \Delta\_2(\boldsymbol{\theta} + \boldsymbol{\phi}) \mathbf{R}'(\boldsymbol{\phi})] \tag{42}$$

we see that - *b*1(*x*, *y*, θ)*dxdyd*θ = 0. Thus, *L*0ρ1 = *b*1 has a solution. Since:

$$L\_0[\mathbf{x}G(\mathbf{x},\mathbf{y})] = -\mathbf{x}G(\mathbf{x},\mathbf{y})/\mathbf{r} \tag{43}$$

$$L\_0[\chi G(\mathbf{x}, \mathbf{y})] = -\chi G(\mathbf{x}, \mathbf{y})/\pi \tag{44}$$

we look for a solution of the form:

$$\rho\_1(\mathbf{x}, \mathbf{y}, \theta, \phi) = \frac{\mathbf{w}\_1(\theta, \phi)\mathbf{x}\mathbf{G}(\mathbf{x}, \mathbf{y}) + \mathbf{w}\_2(\theta, \phi)\mathbf{y}\mathbf{G}(\mathbf{x}, \mathbf{y})}{2\pi} \tag{45}$$

Inserting this into Equation 38, we find that *wj*(θ, φ) must satisfy:

$$\frac{\partial}{\partial\theta}\mathbb{w}\_1(\theta,\phi) + \frac{\mathbb{w}\_1(\theta,\phi)}{\mathfrak{r}} = -\Delta\_1'(\theta)R(\phi) + \Delta\_1(\theta)R'(\phi) \tag{46}$$

$$\frac{\partial}{\partial \theta} \mathcal{w}\_2(\theta, \phi) + \frac{\mathcal{w}\_2(\theta, \phi)}{\mathfrak{r}} = -\Delta\_2'(\theta + \phi)R(\phi) - \Delta\_2(\theta + \phi)R'(\phi) \tag{47}$$

*wj* must be periodic functions of θ; we defer their exact solution to later, but note that there is always a unique periodic solution to each of these equations.

#### *Solving Equation 39*

We now have:

$$b\_2(\mathbf{x}, \mathbf{y}, \boldsymbol{\theta}) = \frac{\partial}{\partial \boldsymbol{\theta}} [\Delta\_1(\boldsymbol{\theta}) \mathbf{x} \boldsymbol{\rho}] + \frac{\partial}{\partial \boldsymbol{\phi}} [(\Delta\_2(\boldsymbol{\theta} + \boldsymbol{\phi}) \mathbf{y} - \Delta\_1(\boldsymbol{\theta}) \mathbf{x}) \boldsymbol{\rho}] + a \frac{\partial}{\partial \boldsymbol{\phi}} \boldsymbol{\rho} \mathbf{o} \tag{48}$$

In order to solve Equation 39, for ρ2(*x*, *y*, θ, φ), we must have:

$$\begin{split} 0 &= \iint \int\_{-\infty}^{\infty} \int\_{0}^{2\pi} b\_{2}(\mathbf{x}, \mathbf{y}, \boldsymbol{\Phi}) d\mathbf{x} d\boldsymbol{\Phi} d\boldsymbol{\Phi} \\ &= 0 + \frac{\partial}{\partial \boldsymbol{\Phi}} \left\{ \iint\_{-\infty}^{\infty} \int\_{0}^{2\pi} \left\{ \frac{\Delta\_{1}(\boldsymbol{\Phi} + \boldsymbol{\Phi})}{2\pi} [w\_{2}(\boldsymbol{\Phi}, \boldsymbol{\Phi})\boldsymbol{\chi}^{2} + w\_{1}(\boldsymbol{\Phi}, \boldsymbol{\Phi})\boldsymbol{\chi}\boldsymbol{\chi}] G(\mathbf{x}, \boldsymbol{\eta}) \right. \right. \\ & \left. - \frac{\Delta\_{1}(\boldsymbol{\Phi})}{2\pi} [w\_{1}(\boldsymbol{\Phi}, \boldsymbol{\Phi})\boldsymbol{\chi}^{2} + w\_{2}(\boldsymbol{\Phi}, \boldsymbol{\Phi})\boldsymbol{\chi}\boldsymbol{\chi}] G(\mathbf{x}, \boldsymbol{\eta}) + \frac{a}{2\pi} R(\boldsymbol{\Phi}) G(\mathbf{x}, \boldsymbol{\eta}) \right\} d\mathbf{x} d\boldsymbol{y} d\boldsymbol{\theta} \\ &= \frac{1}{4\pi} \frac{\partial}{\partial \boldsymbol{\Phi}} \left\{ \int\_{0}^{2\pi} [\Delta\_{2}(\boldsymbol{\Phi} + \boldsymbol{\Phi})] \nu\_{2}(\boldsymbol{\Phi}, \boldsymbol{\Phi}) + c \cdot w\_{1}(\boldsymbol{\Phi}, \boldsymbol{\Phi})] - \Delta\_{1}(\boldsymbol{\Phi}) [w\_{1}(\boldsymbol{\Phi}, \boldsymbol{\Phi}) + c \cdot w\_{2}(\boldsymbol{\Phi}, \boldsymbol{\Phi})]] d\boldsymbol{\Phi} + 4\pi a R(\boldsymbol{\Phi}) \right\} \\ &= \frac{1}{4\pi} \frac{\partial}{\partial \boldsymbol{\Phi}} [f(\boldsymbol{\Phi}) + 4\pi a R(\boldsymbol{\Phi})] \end{split} \tag{49}$$

where:

$$f(\boldsymbol{\phi}) = \int\_0^{2\pi} [\Delta\_2(\boldsymbol{\theta} + \boldsymbol{\phi})\nu\_2(\boldsymbol{\theta}, \boldsymbol{\phi}) - \Delta\_1(\boldsymbol{\theta})\nu\_1(\boldsymbol{\theta}, \boldsymbol{\phi})]d\boldsymbol{\theta} \tag{50}$$

$$\nu\_1(\boldsymbol{\theta}, \boldsymbol{\phi}) = \nu\_1(\boldsymbol{\theta}, \boldsymbol{\phi}) + c\nu\_2(\boldsymbol{\theta}, \boldsymbol{\phi})$$

$$\nu\_2(\boldsymbol{\theta}, \boldsymbol{\phi}) = \nu\_2(\boldsymbol{\theta}, \boldsymbol{\phi}) + c\nu\_1(\boldsymbol{\theta}, \boldsymbol{\phi})$$

Given Equations 46 and 47, we see that *v*1(θ) and *v*2(θ) satisfy:

$$\begin{aligned} \nu\_1'(\theta) + \frac{\nu\_1(\theta)}{\pi} &= -[\Delta\_1'(\theta) + c \cdot \Delta\_2'(\theta + \phi)]R(\phi) + [\Delta\_1(\theta) - c \cdot \Delta\_2(\theta + \phi)]R'(\phi) \\ &:= q\_1(\theta) \\ \nu\_2'(\theta) + \frac{\nu\_2(\theta)}{\pi} &= -[c \cdot \Delta\_1'(\theta) + \Delta\_2'(\theta + \phi)]R(\phi) + [c \cdot \Delta\_1(\theta) - \Delta\_2(\theta + \phi)]R'(\phi) \\ &:= q\_2(\theta) \end{aligned} \tag{51}$$

For Equations 51 and 52, we can write down the solution of *v*1(θ) and *v*2(θ) in terms of *q*1(θ) and *q*2(θ) (see Appendix):

$$\nu\_n(\theta) = \int\_0^\infty e^{-\frac{t}{\epsilon}} q\_n(\theta - s) ds, \quad n = 1, 2 \tag{53}$$

**Frontiers in Computational Neuroscience www.frontiersin.org** August 2013 | Volume 7 | Article 113 |

We substitute *vn*(φ) into *f*(φ),

$$f(\boldsymbol{\theta}) = \int\_0^{2\pi} [\Delta\_2(\boldsymbol{\theta} + \boldsymbol{\phi})\nu\_2(\boldsymbol{\theta}) - \Delta\_1(\boldsymbol{\theta})\nu\_1(\boldsymbol{\theta})]d\boldsymbol{\theta}$$

$$= \int\_0^{\infty} e^{-\frac{t}{\tau}} ds \int\_0^{2\pi} [\Delta\_2(\boldsymbol{\theta} + \boldsymbol{\phi})q\_2(\boldsymbol{\theta} - \boldsymbol{s}) - \Delta\_1(\boldsymbol{\theta})q\_1(\boldsymbol{\theta} - \boldsymbol{s})]d\boldsymbol{\theta}$$

$$= \int\_0^{\infty} e^{-\frac{t}{\tau}} ds \int\_0^{2\pi} h(\boldsymbol{\theta}, \boldsymbol{\phi}, \boldsymbol{s})d\boldsymbol{\theta} \tag{54}$$

where:

$$h(\theta, \phi, s) = \Delta\_2(\theta + \phi)q\_2(\theta - s) - \Delta\_1(\theta)q\_1(\theta - s)$$

$$= \Delta\_1(\theta - s)[c \cdot \Delta\_2(\theta + \phi) - \Delta\_1(\theta)]R'(\phi)$$

$$+ \Delta\_2(\theta + \phi - s)[-\Delta\_2(\theta + \phi) + c \cdot \Delta\_1(\theta)]R'(\phi)$$

$$+ \Delta\_1'(\theta - s)[-c \cdot \Delta\_2(\theta + \phi) + \Delta\_1(\theta)]R(\phi)$$

$$+ \Delta\_2'(\theta + \phi - s)[-\Delta\_2(\theta + \phi) + c \cdot \Delta\_1(\theta)]R(\phi) \tag{55}$$

Define:

$$g\_{mn}(\phi) = \int\_0^\infty h\_{mn}(s+\phi)e^{-\frac{s}{\epsilon}}ds\tag{56}$$

$$h\_{mn}(s) = \int\_0^{2\pi} \Delta\_m(\theta)\Delta\_n(\theta+s)d\theta$$

Since -1(θ) and -2(θ) are periodic functions,

$$\begin{aligned} \int\_{0}^{2\pi} h(\theta,\phi,s)d\theta &= \int\_{0}^{2\pi} \{\Delta\_{1}(\Theta-s)[[c\Delta\_{2}(\Theta+\phi)-\Delta\_{1}(\Theta)]R'(\phi)+[c\Delta\_{2}'(\Theta+\phi)-\Delta\_{1}'(\Theta)]R(\phi)\} \\ &+ \Delta\_{2}(\Theta+\phi-s)[[c\Delta\_{1}(\Theta)-\Delta\_{2}(\Theta+\phi)]R'(\phi)-[c\Delta\_{1}'(\Theta)-\Delta\_{2}'(\Theta+\phi)]R(\phi)]\}d\theta \\ &= [c[h\_{12}(s+\phi)+h\_{21}(s-\phi)]-[h\_{11}(s+h\_{22}(s))]R'(\phi) \\ &+ \left\{c\frac{d}{d\phi}[h\_{12}(s+\phi)+h\_{21}(s-\phi)]-\frac{d}{d\phi}[h\_{11}(s+\phi)-h\_{22}(s+\phi)]\Big|\_{\phi=0}\right\}R(\phi) \\ f(\phi) &= \int\_{0}^{\infty} e^{-\frac{\phi}{\tau}}d\phi \int\_{0}^{2\pi} h(\theta,\phi,s)d\theta \\ &= [c[g\_{12}(\theta)+g\_{21}(-\phi)]-[g\_{11}(0)+g\_{22}(0)]R'(\phi) \\ &+ \left\{c\cdot\frac{d}{d\phi}[g\_{12}(\phi)+g\_{21}(-\phi)]-[g\_{11}'(0)-g\_{22}'(0)]\right\}R(\phi) \\ &= \frac{d}{d\phi}\left[l(\cdot,g(\theta)-C\_{1})R(\phi)\right]-C\_{2}R(\phi) \end{aligned} \tag{58}$$

where:

$$\mathbf{g}(\phi) = \mathbf{g}\_{12}(\phi) + \mathbf{g}\_{21}(-\phi) \tag{59}$$

$$C\_1 = \lg\_{11}(0) + \lg\_{22}(0)\tag{60}$$

$$C\_2 = g\_{11}'(0) - g\_{22}'(0) \tag{61}$$

Combined with Equations 49–58, we have a boundary value problem (BVP):

$$\frac{d}{d\phi}\left[\left(\mathbf{c}\cdot\mathbf{g}(\phi)-\mathbf{C}\_{1}\right)R(\phi)\right]+\left(4\pi\omega-\mathbf{C}\_{2}\right)R(\phi)=K\tag{62}$$

$$R(\phi)=R(\phi+2\pi)$$

$$g(\phi)=g(\phi+2\pi)$$

$$\int\_{-\pi}^{\pi}R(\phi)d\phi=1$$

$$K=2\omega-\frac{\mathbf{C}\_{2}}{2\pi}$$

To solve this BVP, we need to compute *g*(φ) for a given PRC. We use two forms of the PRC in this paper:

$$
\Delta(\theta) = \sin(a) - \sin(\theta + a) + b\sin(2\theta) \tag{63}
$$

and

$$
\Delta(\theta) = A[\sin(B) - \sin(\theta + B)]e^{C(\theta - 2\pi)} \quad (\theta \in (0, 2\pi), \Delta(\theta) = \Delta(\theta + 2\pi))\tag{64}
$$

The required integrals can be computed for both PRC forms. More generally, all PRCs can be written in Fourier form and, again, the integrals are readily computed to obtain *g*(φ) (see below).

#### *Small correlation approximation for R(*φ*)*

We use a BVP solver to get the numerical solution for *R*(φ), but we would like to better understand the form of *R*(φ) at low values of correlation, *c*, so we expand *R* as a series in *c*. As *K* is dependent on *c*, we must also expand *K* in *c*. Finally, we need to keep the normalization condition for *R*(φ), hence:

$$\begin{aligned} R(\phi) &= R\_0(\phi) + cR\_1(\phi) + \dots & K = K\_0 + cK\_1 + \dots \\ R\_0(\phi) &= R(\phi)|\_{\varepsilon=0}, \quad \int\_{-\pi}^{\pi} R\_0(\phi) d\phi = 1 \\ \int\_{-\pi}^{\pi} R\_n(\phi) &= 0, n \ge 1 \end{aligned} \tag{65}$$

We substitute these expressions into the BVP, Equation 62 and find after equating powers of *c*:

$$-C\_1 R\_0'(\phi) + (4\pi\omega - C\_2)R\_0(\phi) = K\_0 \tag{66}$$

$$-C\_1R\_1'(\phi) + (4\pi a - C\_2)R\_1(\phi) + [\![g(\phi)R\_0(\phi)]\!]' = K\_1 \tag{67}$$

We can integrate both left and right side over [−π,π], then use the assumptions and periodicity requirements above to get:

$$K\_0 = \frac{4\pi\alpha - C\_2}{2\pi}, \quad K\_1 = 0 \tag{68}$$

Rewriting these equations,

$$R\_0'(\phi) + DR\_0(\phi) = Q\_0 \tag{69}$$

$$D = \frac{C\_2 - 4\pi\omega}{C\_1}, \quad Q\_0 = -\frac{K\_0}{C\_1}$$

$$R\_1'(\phi) + DR\_1(\phi) = Q\_1(\phi) \tag{70}$$

$$Q\_1(\phi) = \frac{[g(\phi)R\_0(\phi)]'}{C\_1}$$

we see:

$$R\_0(\phi) = \frac{1}{2\pi}$$

$$[R\_1(\phi)e^{D\phi}]' = Q\_1(\phi)e^{D\phi}$$

We can use numerical methods to get the solution of *R*1(φ) given *R*0(φ) and for some choices of PRCs we can get exact expressions. For example, for the two term double sinusoidal form PRC (Equation 63) we get:

$$R\_1(\phi) = \frac{1}{C\_1} \left\{ \frac{\pi}{\pi^2 + 1} \frac{\cos(\phi + a\_2 - a\_1) - D \sin(\phi + a\_2 - a\_1)}{1 + D^2} + \frac{2b\_1 b\_2 \pi}{4\pi^2 + 1} \frac{2 \cos(2\phi) - D \sin(2\phi)}{4 + D^2} \right\} \tag{71}$$

#### **DIFFERENT PRCS**

*Double sines form PRC*

For the PRC:

$$\Delta\_m(\theta) = \sin(a\_m) - \sin(\theta + a\_m) + b\_m \sin(2\theta)$$

We have:

$$\begin{split} h\_{\text{nn}}(s) &= \int\_0^{2\pi} \Delta\_{\text{m}}(\theta) \Delta\_{\text{n}}(\theta + s) d\theta \\ &= 2\pi \sin(a\_{\text{m}}) \sin(a\_{\text{n}}) + \pi \cos(s - a\_{\text{m}} + a\_{\text{n}}) + b\_{\text{m}} b\_{\text{n}} \pi \cos(2s) \\ \text{where} \quad \epsilon\_{\text{n}}(\theta) &= \int\_0^{\infty} h\_{\text{nn}}(s \pm \lambda) e^{-\frac{s}{\lambda}} ds \end{split} \tag{72}$$

$$\begin{aligned} g\_{mn}(\phi) &= \int\_0^\pi h\_{mn}(s+\phi)e^{-\frac{\pi}{4}}ds \\ &= 2\pi\tau\sin(a\_m)\sin(a\_n) + \pi\tau \frac{\cos(\phi+a\_n-a\_m) - \tau\sin(\phi+a\_n-a\_m)}{\pi^2+1} \end{aligned}$$

$$\cos(2\phi) - 2\pi\sin(2\phi)$$

$$a\_n + b\_m b\_n \pi \tau \frac{\cos(2\phi) - 2\tau \sin(2\phi)}{4\pi^2 + 1} \tag{73}$$

$$g\_{mn}'(\phi) = -\pi \text{tr} \frac{\pi \cos(\phi + a\_n - a\_m) + \sin(\phi + a\_n - a\_m)}{\pi^2 + 1} - 2b\_m b\_n \pi \text{tr} \frac{2\pi \cos(2\phi) + \sin(2\phi)}{4\pi^2 + 1} \tag{74}$$

$$\begin{split} g(\phi) &= g\_{12}(\phi) + g\_{21}(-\phi) \\ &= 4\pi\tau \sin(a\_1)\sin(a\_2) + 2\pi\tau \frac{\cos(\phi + a\_2 - a\_1)}{\pi^2 + 1} + 2b\_1 b\_2 \pi \tau \frac{\cos(2\phi)}{4\pi^2 + 1} \end{split} \tag{75}$$

$$\begin{split} g'(\phi) &= g'\_{12}(\phi) - g'\_{21}(-\phi) \\ &= -2\pi \mathfrak{r} \frac{\sin(\phi + a\_2 - a\_1)}{\mathfrak{r}^2 + 1} - 4b\_1 b\_2 \,\pi \mathfrak{r} \frac{\sin(2\phi)}{4\tau^2 + 1} \end{split} \tag{76}$$

$$\text{tr}\,C\_{l} = \mathcal{g}\_{11}(0) + \mathcal{g}\_{22}(0) = 2\pi\tau(\sin^{2}(a\_{1}) + \sin^{2}(a\_{2})) + \frac{2\pi\tau}{\tau^{2} + 1} + \frac{(b\_{1}^{2} + b\_{2}^{2})\pi\tau}{4\tau^{2} + 1} \tag{77}$$

$$C\_2 = g\_{11}'(0) - g\_{22}'(0) = \frac{4\pi\mathfrak{r}^2(b\_2^2 - b\_1^2)}{4\mathfrak{r}^2 + 1} \tag{78}$$

#### *Exponential-sine form PRC*

For empirical PRCs:

We have:

$$h\_{mt}(\mathbf{s}) = B\_1 \cdot e^{-\epsilon\_{ns}s} \langle hc\_1, h\nu\_1(\mathbf{s}) \rangle + B\_0 \cdot e^{\epsilon\_{qs}} \langle hc\_0, h\nu\_0(\mathbf{s}) \rangle \tag{79}$$

$$g\_{mn}(\phi) = C\_0 \cdot e^{\frac{\phi}{\varepsilon}} - e^{-c\_0 \phi} \langle \g c\_1, g\nu\_1(\phi) \rangle - e^{\varepsilon\_0 \phi} \langle \g c\_0, g\nu\_0(\phi) \rangle \tag{80}$$

Where *s* ∈ [0, 2π) and φ ∈ [0, 2π), both have the period of 2π. Note also that **(,)** means the inner product of two vectors. The above quantities are defined as:

*B*1 = *am* · *an* · *e* −2π*cn* · [1 − *e* −2π*cm* ] *B*0 = *am* · *an* · [1 − *e* −2π*cn* ] *D*1 = *cm* + 1 τ *D*0 = *cn* − 1 τ *k*0 = sin(*bm*) · sin(*bn*) *cm* + *cn* − sin(*bm*) 1 + (*cm* + *cn*)2 [(*cm* + *cn*) · sin(*bn*) − cos(*bn*)] *k*1 = 1 2(*cm* + *cn*) , *k*2 = 1 4 + (*cm* + *cn*)2 , *k*3 = − *cm* + *cn* 2[4 + (*cm* + *cn*)2] *k*4 = (*cm* + *cn*) · sin(*bn*) 1 + (*cm* + *cn*)2 , *k*5 = sin(*bn*) 1 + (*cm* + *cn*)2 *j*0 = sin(*bm*) · sin(*bn*) *cm* + *cn* − sin(*bn*) 1 + (*cm* + *cn*)2 [(*cm* + *cn*) · sin(*bm*) − cos(*bm*)] *j*1 = 1 2(*cm* + *cn*) , *j*2 = − 1 4 + (*cm* + *cn*)2 , *j*3 = − *cm* + *cn* 2[4 + (*cm* + *cn*)2] *j*4 = −(*cm* + *cn*) · sin(*bm*) 1 + (*cm* + *cn*)2 , *j*5 = sin(*bm*) 1 + (*cm* + *cn*)2 *hc*1 = [*k*0, *k*1, *k*2, *k*3, *k*4, *k*5] *hc*0 = [*j*0, *j*1, *j*2,*j*3,*j*4, *j*5] *hv*1 = [1, cos(*s* + *bn* − *bm*),sin(*s* − *bm* − *bn*), cos(*s* − *bm* − *bn*),sin(*s* − *bm*), cos(*s* − *bm*)] *hv*0 = [1, cos(*s* + *bn* − *bm*),sin(*s* + *bm* + *bn*), cos(*s* + *bm* + *bn*),sin(*s* + *bn*), cos(*s* + *bn*)] *gc*1 = *B*1 1 + *D*2 1 [− 1 + *D*2 1 *D*1 *k*0, −*D*1*k*1, *k*1, *k*3 − *D*1*k*2, −*k*2 − *D*1*k*3, *k*5 − *D*1*k*4, −*k*4 − *D*1*k*5] *gc*0 = *B*0 1 + *D*2 0 [ 1 + *D*2 0 *D*0 *j*0, *D*0*j*1,*j*1,*j*3 + *D*0*j*2, −*j*2 + *D*0*j*3,*j*5 + *D*0*j*4, −*j*4 + *D*0*j*5] *gv*1(φ) = [1, cos(φ + *bn* − *bm*),sin(φ + *bn* − *bm*),sin(φ − *bm* − *bn*), cos(φ − *bm* − *bn*),sin(φ − *bm*), cos(*s* − *bm*)] *gv*0(φ) = [1, cos(φ + *bn* − *bm*),sin(φ + *bn* − *bm*),sin(φ + *bm* + *bn*), cos(φ + *bm* + *bn*),sin(φ + *bn*), cos(*s* + *bn*)] *C*0 = *e*−2πτ 1 − *e* − 2π τ [(*e* −2π*cm* − 1) · **(***gc*1, *gv*1(0)**)** + (*e* 2π*cn* − 1) · **(***gc*0, *gv*0(0)**)**]

*Fourier form PRC*

For the Fourier form of the PRC:

$$\Delta\_m(\theta) = \sum\_{k=-\infty}^{\infty} a\_{m,k} e^{i\theta \cdot \theta}$$

$$\begin{aligned} h\_{mn}(s) &= \int\_0^{2\pi} \Delta\_m(\theta) \Delta\_n(\theta) ds \\ &= \sum\_{k\_1, k\_2} a\_{m, k\_1} a\_{n, k\_2} \int\_0^{2\pi} e^{ik\_1 \theta} e^{ik\_2(\theta + s)} d\theta \end{aligned}$$

$$=2\pi \sum\_{k=-\infty}^{\infty} a\_{m,k} a\_{n,-k} e^{-iks} \tag{81}$$

$$g\_{mn}(\phi) = \int\_0^{\infty} h\_{mn}(s+\phi) e^{-\frac{t}{\tau}} ds$$

$$= 2\pi \sum\_{k=-\infty}^{\infty} a\_{m,k} a\_{n,-k} \int\_0^{\infty} e^{-ik(s+\phi)} e^{-\frac{t}{\tau}} ds$$

$$= 2\pi \sum\_{k=-\infty}^{\infty} \frac{a\_{m,k} a\_{n,-k}}{k^2 + \frac{1}{\tau^2}} (\frac{1}{\tau} - ik) e^{-ik\phi} \tag{82}$$

## Direct connections assist neurons to detect correlation in small amplitude noises

## *E. Bolhasani , Y. Azizi and A. Valizadeh\**

*Department of Physics, Institute for Advanced Studies in Basic Sciences, Zanjan, Iran*

#### *Edited by:*

*Ruben Moreno-Bote, Foundation Sant Joan de Deu, Spain*

#### *Reviewed by:*

*Germán Mato, Centro Atomico Bariloche, Argentina Zachary P. Kilpatrick, University of Houston, USA*

#### *\*Correspondence:*

*A. Valizadeh, Department of Physics, Institute for Advanced Studies in Basic Sciences, Gava zang, PO Box 45195-1159, Zanjan, Iran e-mail<EMAIL_ADDRESS>
## **1. INTRODUCTION**

The recent advent of novel recording techniques has made it easier to simultaneously record from a large number of neurons and has provided new possibilities to relate population activity to coding and information processing in the brain (Greenberg et al., 2008; Cohen and Kohn, 2011). Many researchers suggest that studying the correlated activity of neurons in a population is essential for understanding how information is coded in the brain (Zohary et al., 1994; Abbott and Dayan, 1999; Nirenberg and Latham, 2003; Averbeck et al., 2006; Biederlack et al., 2006; Schneidman et al., 2006; Pillow et al., 2008).  In addition to the functional effects of such correlations between populations of neurons on neural coding, understanding how different parameters such as biological, network or stimulus parameters tune them is eventually being revealed (Shadlen and Newsome, 1998; Binder and Powers, 2001; Moreno et al., 2002; Moreno-Bote and Parga, 2006; Tchumatchenko et al., 2010b; Rosenbaum and Josic, 2011b ´ ). Correlation between neuronal activities is measured frequently by pairwise correlation coefficients and spike count correlations, and the ability of a neuronal system to transfer correlation can be quantified by the correlation transfer function (CTF), which determines the relation between the output correlation of a system under stimulus and a specific input correlation (Doiron et al., 2006; Shea-Brown et al., 2008; Rosenbaum and Josic, 2011b ´ ).

A periodic common input on two (or more) uncoupled oscillators can cause coherent behavior when both oscillators lock to the external force (Pikovsky et al., 2003). A very common example is the control of circadian rhythms of humans/animals

We address a question on the effect of common stochastic inputs on the correlation of the spike trains of two neurons when they are coupled through direct connections. We show that the change in the correlation of small amplitude stochastic inputs can be better detected when the neurons are connected by direct excitatory couplings. Depending on whether intrinsic firing rate of the neurons is identical or slightly different, symmetric or asymmetric connections can increase the sensitivity of the system to the input correlation by changing the mean slope of the correlation transfer function over a given range of input correlation. In either case, there is also an optimum value for synaptic strength which maximizes the sensitivity of the system to the changes in input correlation.

**Keywords: correlation, correlation transfer, coupling, inhomogeneity, synchrony**

by the light-dark stimulation (Roberts, 2005). In case of noisy inputs the counterpart of the phenomena appears as stochastic synchronization (SS) which is a general topic that addresses the phenomenon of irregular phase locking between two noisy non-linear oscillators (Neiman et al., 1999). In nervous systems, cross-correlations can arise either from the presence of direct synaptic connections (Csicsvari et al., 1998; Barthó et al., 2004) or from shared inputs from the surrounding network or sensory layers (Binder and Powers, 2001; Türker and Powers, 2001, 2004). Effect of direct synaptic connections and common inputs have been widely studied, but these two sources of correlation can be present concurrently in many physical and biological systems and their interplay can result in quite interesting phenomena. Couplings can regulate the activity of noisy oscillators and less variability in neuronal dynamics emerges through synchronization in networks of coupled noisy oscillators (Ly and Ermentrout, 2010; Tabareau et al., 2010; Zilli and Hasselmo, 2010). Studies on the correlation of spike trains have reported increase and decrease of correlation due to the presence of excitatory and inhibitory synapses, respectively (Rosenbaum and Josic, 2011a; Ly et al., ´ 2012). When delay in communication and type of excitabilty of neurons are taken into account, the generality of these results can be debated since both excitatory and inhibitory synapses can be sources of synchrony and may increase correlation in different parameter ranges (Vreeswijk et al., 1994; Wang et al., 2012; Sadeghi and Valizadeh, 2013). Regarding the type of excitability and categorizing couplings as synchronizing and desynchronizing, it has been shown that shared inputs and direct couplings can show cooperative or disruptive effects on the correlation of noisy coupled oscillators (Ly and Ermentrout, 2009).

Possible differences between intrinsic parameters of neurons causes the message from the environment to the system to be decoded differently by the system components. Another aim of the current study is to investigate how the correlation is transferred by two neurons when the neurons are not identical. In such a heterogeneous system, the temporal symmetry of spike correlation is lost (Tchumatchenko et al., 2010b). We will show that with small amplitude stochastic inputs, even a slight inhomogeneity in the intrinsic parameters can lead to a large reduction of the pairwise correlation coefficient in the case of uncoupled neurons. As expected, the results depend on the time bins over which the correlation is calculated: spike count correlations over long time bins are less affected by the heterogeneity but synchrony—alignment of the action potential in small time bins—is tightly dependent on the homogeneity of the system.

We have shown that correlated inputs and direct connections can either show cooperative or disruptive effects in different ranges of parameters. For uncoupled neurons, correlation susceptibility increases by increasing the amplitude of noise for mildly correlated inputs (De La Rocha et al., 2007; Shea-Brown et al., 2008; Tchumatchenko et al., 2010b). We show that when direct connections are present between non-identical neurons, the mean susceptibility is not a monotonic function of the amplitude of the correlated noisy input anymore. Reminiscence of stochastic resonance phenomena, an intermediate noise amplitude in this case, leads to larger a sensitivity of the system to the changes in input correlation. We have also shown that with monosynaptic connections between two neurons, presence of inhomogeneity in the intrinsic firing rate of the neurons can enhance correlation of spike trains while for symmetric couplings, maximum correlation is seen for homogeneous system. Changing mismatch and synaptic strengths between two neurons, it is possible to change the functional form of the correlation transfer function to optimize the mean correlation susceptibility which is an indicator of the sensitivity of the system to the change of input correlation in different ranges. In this way, as the most important result of current study, we will show that with direct couplings it is possible to detect correlation in small amplitude noises by increasing the sensitivity of the system to the change of correlation in small amplitude noisy inputs.

## **2. MATERIALS AND METHODS**

The system under investigation consists of two coupled leaky integrate and fire (LIF) neurons (Knight, 1972), subjected to correlated stochastic inputs (see **Figure 1**). Subthreshold dynamics of the LIF neuron obeys the following first order equation:

$$
\pi\_m \frac{d\nu\_i}{dt} = V\_{\text{rest}} - \nu\_i + I\_i + I\_{ij}, \tag{1}
$$

in which *vi* is a voltage-like variable for each neuron labeled by *i* = 1, 2 with τ*m* = 20 ms and *V*rest = −70 mV. A severe nonlinearity is imposed on the model by considering a threshold value *vth* = −54 mV. Whenever this value is reached, the neuron *spikes* and the voltage resets to *v*reset = −60 mV. [Parameters taken from Troyer and Miller (1997)]. The spikes of the neurons are recorded as *xi*(*t*) = - *m* δ(*t* − *t m i* ) where *t m i* is the time of *m* th spike of the neuron *i*, and δ(*x*) is the Dirac delta function.

Each model neuron receives a synaptic current through the direct connection from the other neuron *Iij*, and an external current *Ii* representing the sensory input or the effect of the surrounding networks. In the model equations, external current to

the neuron *i* comprises a constant (dc) and a stochastic component with amplitude σ. The stochastic inputs are sum of a common component ξ*c*(*t*) and an individual component ξ*i*(*t*):

$$I\_i(t) = (1 \pm \delta)I + \sigma \left[ \sqrt{1 - c} \xi\_i(t) + \sqrt{c} \xi\_c(t) \right],\tag{2}$$

where ξ*c*(*t*) and ξ*i*(*t*) are mutually independent Gaussian stochastic processes with zero mean and unit variance ξ*i*(*t*)ξ*j*(*t* ) = δ*ij*δ(*t* − *t* ). The parameter *c* ∈ [0, 1] determines correlation of external currents which will be referred to as the input correlation. With the minimal model we used, inhomogeneity in the intrinsic activity rates is imposed by different constant currents which are chosen as *I*1 = (1 + δ)*I* and *I*2 = (1 − δ)*I*, where δ is referred to as the parameter of inhomogeneity. With non-zero δ the neurons 1 and 2 will be the high frequency (fast) and low frequency (slow) neurons, respectively. The currents are chosen suprathreshold (>14 mV) such that the neurons fire periodically at vanishing noise. Note that in this mean driven regime presence of small amplitude noise results in small jitters in firing times and a narrow distribution of interspike intervals.

Neurons are pulse coupled. The neuron *i*receives a pulse by the strength *ij* every time the neuron *j* fires, so the synaptic current in Equation 1 can be written as *Iij* = *ijxj*(*t*) where the synaptic strength *ij* can be positive (excitatory) or negative (inhibitory). For convenience, we call the connections 21 and 12, the forward and backward connections, respectively. Although the external and synaptic inputs appear as currents, they are actually measured in units of the membrane potential (mV) since a factor of the membrane resistance has been absorbed into their definition.

Co-fluctuations in the activity of neurons are measured over a range of timescales (for a review see Cohen and Kohn, 2011). Spike count correlation is usually measured over the time scales from tens of milliseconds to seconds, while synchrony, that is almost precise alignment of the spikes, is measured over the time scale of the typical width of an action potential. It has been shown that spike count correlation over the small bins, bins of the order of one millisecond, can be largely determined by zero-lag conditional firing rate which quantifies exact synchrony (Tchumatchenko et al., 2010a). In this study we focus on synchrony, by describing spike counts and correlation coefficients in discrete bins of duration *T* = 0.5 ms. Correlation coefficient of spike counts *ni*(*t*) = *t* +*T t xi*(*s*)*ds*, is defined as the zero lag cross-correlation between *n*1 and *n*2:

$$\rho\_T = \frac{\langle n\_1(t)n\_2(t)\rangle - \langle n\_1(t)\rangle \langle n\_2(t)\rangle}{\sqrt{\langle n\_1(t)^2\rangle - \langle n\_1(t)\rangle^2}\sqrt{\langle n\_2(t)^2\rangle - \langle n\_2(t)\rangle^2}}. \tag{3}$$

Dependence of the output correlation to the input correlation shows how correlation is transferred along neuronal layers in the nervous system (Rosenbaum and Josic, 2011a ´ ). With varying input correlation while other parameters are fixed, we compute ρ*T*(*c*), correlation of spike trains as a function of input correlation. To study sensitivity of correlation of output spike trains to the change of input correlation, we use *mean correlation susceptibility* (MCS), the mean slope of ρ*T*(*c*) in a given range of *c* ∈ [*c*1,*c*2]:

$$S\_T(c\_1, c\_2) = \frac{\Delta \rho\_T}{\Delta c}.\tag{4}$$

which shows ratio of the change of correlation of spike trains ρ*T* = ρ*T*(*c*2) − ρ*T*(*c*1) to the change of input correlation *c* = *c*2 − *c*1. For two identical neurons with no direct connection, this value is equal to one when it is evaluated over the full range of input correlation [0, 1].

## **3. RESULTS**

We first present the results for two uncoupled neurons. In **Figure 2A** we have shown the cross-correlation coefficient as a function of the mismatch between intrinsic firing rates of neurons for low noise amplitude and different values of the input correlation. When there is no direct connection between the neurons, highly correlated inputs lead to a large output correlation in case of identical neurons. Even a small mismatch decreases the output correlation considerably if the noise is small amplitude. In this case, even common noises lead to a relatively low output correlation in the presence of a slight inhomogeneity (e.g., δ = 0.01 in **Figure 2A**). For larger noise amplitudes, the output correlation is less sensitive to inhomogeneity (**Figure 2B**). The system is also less sensitive to inhomogeneity when the inputs are weakly correlated where both homogeneous and inhomogeneous systems show a small output correlation. In **Figures 2C,D** we have shown the correlation transfer function. It can be seen that while the slope of the correlation transfer function decreases with mismatch for all the values of input correlation, this dependence is only noticeable when inputs are highly (completely) correlated. Increasing the noise amplitude (while decreasing the constant input to avoid a change in the mean firing rate as explained below) makes the output correlation less sensitive to inhomogeneity, yet the maximum sensitivity to mismatch is observed for highly correlated inputs (**Figure 2D**).

To show how sensitive are the correlation of spike trains to the input correlation, in **Figure 2E** we have plotted MCS (mean slope

**FIGURE 2 | Correlation of spike trains for two uncoupled neurons. (A)** Correlation coefficient is plotted against inhomogeneity, the mismatch between input current of neurons, for different values of input correlation and low noise amplitude σ = 1 mV. In **(B)** the same results are shown for larger value of noise amplitude σ = 5 mV with the same mean firing rate as **(A)** (see materials and methods). **(C,D)** Correlation transfer function, which shows the dependence of correlation of spike trains to the input

correlation, is plotted for different values of inhomogeneity for the same noise amplitudes as **(A)** and **(B)**. **(E)** Mean correlation susceptibility (MCS) is plotted for homogeneous and slightly inhomogeneous systems, as a function of noise amplitude, which shows the mean sensitivity of the output correlation to the change of input correlation over the range [0, 0.5]. In **(F)** the geometric mean of the firing rate of the two neurons is shown when σ is varied.

of ρ*T*(*c*) as described in materials and methods) as a function of the amplitude of the stochastic input for two uncoupled neurons over the range *c* ∈ [0 − 0.5] for homogeneous (δ = 0) and slightly inhomogeneous (δ = 0.02) systems. The system shows low sensitivity to the change in input correlation for small amplitude noises and the sensitivity smoothly increases with noise amplitude. Also, the presence of inhomogeneity has negligible effect on the mean correlation suceptibility: as noted above, for uncoupled neurons effect of inhomogeneity is only significant when inputs are highly correlated and while MCS is calculated over a range of weakly correlated inputs, it is almost insensitive to small inhomogeneity. While increasing the amplitude of the fluctuations, we have decreased mean value of the input currents to keep the mean firing rate almost constant (∼64 Hz) as is shown in **Figure 2F**. In such a way the results observed in **Figure 2E** can not be attributed to the increase in firing rate which is known to increase the spike train correlation (De La Rocha et al., 2007; Shea-Brown et al., 2008). These results show that the correlation in small amplitude noises can not be suitably detected by a system of uncoupled neurons, whether the neurons have equal firing rates or their firing rates are different. To investigate the effect of direct couplings we have first considered a two neurons motif with just one unidirectional excitatory synapse. In many cases this configuration is favored when the synapses change through spike timing-dependent plasticity (Song et al., 2000). We considered an excitatory forward coupling from the high frequency neuron (as the presynaptic) to low frequency neuron (as the postsynaptic). In the absence of noise, any finite value of the forward coupling strength can lead to a zone of 1:1 synchrony, in which the dissimilar neurons fire in a causal master-slave fashion (Takahashi et al., 2009; Bayati and Valizadeh, 2012). In such causal limit the postsynaptic neuron fires immediately after receiving presynaptic stimulation (Woodman and Canavier, 2011; Wang et al., 2012). In our model delays in communication have been ignored, so in the causal 1:1 synchrony zones the postsynaptic neuron fires just one simulation time step after the firing of presynaptic neuron. Since the time bin on which the correlation is calculated contains several time steps (see materials and methods), such a causal master-slave firing leads to ρ = 1 (gray curves in **Figure 3**).

Stochastic inputs have non-trivial effects on the correlation of the spike trains of these two neurons. The output correlation is not a monotonically decreasing function of mismatch anymore, and in the presence of noise a small mismatch can increase the output correlation (**Figure 3A**). With zero mismatch, in the presence of one excitatory connection from neuron 1 to neuron 2 and in the absence of noise, the only stable state is the phase locked state in which neuron 2 fires one time step after neruon 1 (Bayati and Valizadeh, 2012). In the presence of noise this state looses stability as follows: because of the initial phase difference between the two neurons after master-slave firing (even though the phase difference is very small, just one time step), they respond slightly differently even to common noises. The different responses of the two neurons lead to a cumulative phase difference and if this phase difference results in the firing of neuron 2 before neuron 1

are shown when the neurons are bidirectionally coupled by symmetric connections. In **(C)** and **(D)** the results are presented for larger noise amplitude σ = 5 mV. Noise amplitude in **(A)** and **(B)** is σ = 1 mV. The gray curves correspond to autonomous case when no stochastic input is present.

reaches threshold, the excitatory pulse from neuron 1 would be desynchronizing and makes the next firing of the two neurons further apart. The probability of the advancement of the phase of neuron 2 decreases in the presence of inhomogeneity (with *I*1 > *I*2), and with larger inhomogeneity it is less likely that the firing of neuron 2 (low frequency neuron) exceeds the firing of neuron 1 (high frequency neuron). When neuron 2 fires, before neuron 1 has reached the threshold, the excitatory pulse to the low frequency neuron will be synchronizing and if the voltage of neuron 2 is in the range [*vth* − -21, *vth*] at the time of the firing of neuron 1, the neurons maintain causal master-slave firing. Further increasing the inhomogeneity lowers the probability of the voltage of the low frequency neuron reaching the range [*vth* − -21, *vth*] at the time of the firing of the high frequency neuron, which results in the reduction of the spike trains correlation. A similar argument can explain the other notable rise and fall of the correlation which is seen in 1:2 locking zone of the noiseless system.

With symmetric bidirectional couplings, maximum correlation is obtained when the neurons are of the same firing rate (**Figure 3B**). When the neurons have equal firing rates (with *I*1 = *I*2) and in the absence of noise, each of the neurons can play the role of the master in a causal master-slave firing: in this case the connection from the master is synchronizing and the other connection has a desynchronizing effect (Bayati and Valizadeh, 2012). In the presence of small amplitude noise, the system can maintain causal locking by interchanging the role of two connections as synchronizing and desynchronizing. Suppose the firing of neuron 1 (master) is followed by the firing of neuron 2 (slave). Firing of neuron 2 exerts an excitatory pulse on neuron 1 but the phase advance of neuron 1 is relatively small because of the weak response of the LIF neuron at the beginning of its cycle (Mirollo and Strogatz, 1990). So it is probable that neuron 2 fires before neuron 1 reaches the threshold, then the excitatory pulse to neuron 1 would be synchronizing and neuron 1 fires immediately at the time it receives the pulse if its voltage is within the range [*vth* − -12, *vth*](note that the argument holds also in the presence of an absolute refractory period where the desynchonizing pulse from the slave neuron is ineffective). In the presence of inhomogeneity, it is the high frequency neuron that more probably plays the role of the master in a locked causal firing in the absence of the noise. In this case, in the presence of noise, inhomogeneity increases the probability that the voltage of low frequency neuron takes a value outside the range [*vth* − -21, *vth*] at the time of the firing of the high frequency neuron, which reduces the correlation of spike trains as can be seen in **Figure 3B** for small values of inhomogeneity. For larger values of inhomogeneity, a bump can be seen again which belongs to the other main locking zone of the system in the absence of noise.

Intuitively, the relative amplitudes of noise and recurrent stimulations determine the behavior of the system and the most notable results can be expected when these two sources are of the same order, i.e., when neither the external noises nor recurrent stimulations are dominant. The results of **Figures 3A,B** are produced in this regime. For larger values of the noise amplitude, qualitative behavior of the system becomes more similar to the uncoupled system as shown in **Figures 3C,D**. For all partially correlated inputs, correlation of the spike trains is independent of the inhomogeneity and no signature of the locking zones is observed in the presence of large amplitude noises. It is only for common noise (γ = 1) that the effect of the unidirectional direct connection can be seen in the presence of strong noise in the region of the main locking zone.

In **Figure 4** we have plotted correlation of spike trains as a function of input correlation to inspect the effect of changing the correlation of the stochastic inputs on the correlation of the spike trains for a fixed value of the synaptic strength. When the noise amplitude is not large, depending on the mismatch, different dependencies of the output correlation to the input correlation can be observed (**Figures 4A,B**. Notably with changing mismatch it is possible to generate, for example, a system with higher sensitivity to the input correlation in different ranges of input correlation, or a negative slope ρ*T*(*c*). Comparing with the results of **Figure 3** it can be deduced that high sensitivities on the input correlation is seen on the main locking zone (where the neurons are causally locked in 1:1 zone in the absence of the noise), and a negative slope is seen between two main locking zones. Again, as can be seen in **Figures 4C,D**, strong noises wash the signature of the direct couplings, and ρ*T*(*c*) for large amplitude noises is qualitatively similar to the uncoupled neurons.

Impact of direct connections on the detection of the input correlation of low amplitude noisy inputs is more apparent in a plot of MCS. In **Figure 5A** we have plotted *ST*(0, 0.5) as a function of noise amplitude for several values of synaptic strength, for unidirectionally coupled neurons and in the presence of a small mismatch in the intrinsic firing rates. As shown in **Figure 5A**, a forward monosynaptic connection (from high frequency to low frequency neuron) can considerably change the performance of the heterogeneous system in detecting variable input correlation.

In an intermediate synaptic strength (-21 = 1) MCS shows a faster growth and a higher maximum in relatively small amplitude noise. Further increasing of the synaptic strength or the noise amplitude reduces the performance of the system in the detection of the input correlation. With very large noise amplitudes, the

effect of the direct connections is washed out and all the curves,

rate of the neurons and ρ*T* (*c*) are shown for the corresponding curves in **(A)**, respectively. Shadings in **(C)** are guide to eye for a comparison

including that of the uncoupled neurons, merge together and the MCS smoothly increases with noise amplitude. Overall increase of the correlation of the spike trains is an intuitive expectation when direct excitatory couplings are present in the systems (although this can be dependent on the type of excitability of the neurons). But how can direct connections increase the sensitivity to the changes in input correlation? In **Figure 5B** we have shown the geometric mean of the firing rate of the two neurons √ν1ν2 for the curves plotted in **Figure 5A**. Note that ν2 may be different from the intrinsic firing rate of neuron 2 because of the presence of an excitatory afferent synapse. The results show that the increase in the mean correlation susceptibility cannot be attributed to the increase of the mean firing rate of neurons, since then, larger coupling constants would lead to more sensitivity as they increase the mean firing rate of the

system. A simple explanation can be found in **Figure 5C**: the degree of amplification of the output correlation depends on the input correlation. A suitable choice of the synaptic strength would result in more amplification for higher input correlations and would increase the slope of ρ*T*(*c*). Increasing the synaptic strength further, decreases the sensitivity due to the saturation of the correlation of the spike trains for the upper bound of the input correlation. In calculating MCS we have considered the range [0, 0.5] for the input correlation. Reducing the upper bound of this range increases the synaptic strength which saturates the correlation of the spike trains, so the synaptic strength which gives the maximum sensitivity increases with decreasing the range over which the mean sensitivity is calculated.

coupling constant is increased. Vertical dotted lines are plotted to show

where the mean sensitivity is maximized.

The *best* synaptic strength, which maximizes sensitivity, depends also on the mismatch between the intrinsic firing rate of the neurons as can be implicitly deduced from the results shown in **Figures 3A,B**. In **Figure 5D** we have shown MCS as a function of the strength of the forward unidirectional coupling for three values of mismatch. Optimum value of synaptic strength is larger when the intrinsic firing rate of the neurons are more different. Plots of the spike train correlation ρ for upper limiting value of the input correlation *c* = 0.5 again shows that the maximum mean sensitivity in this range is obtained when the spike train correlation is not saturated for the upper bound of the range of *c* (**Figure 5E**).

All the results presented in this study have focused on the degree of zero-lag synchrony which is measured by the zero-lag cross correlation of the binned spike trains with small bin size. In the presence of an inhomogeneity and with asymmetric direct connections, it is possible that the maximum correlation of the spike trains appears in non-zero lag. In **Figure 6** we have shown the cross-correlation coefficient of spike trains as a function of the time lag for three values of noise strength and two values of the input correlation (*c* = 0 and *c* = 0.5). It can be seen that the maximum cross correlation for all the values appears in zero time lag (more precisely at a time lag equal to one simulation time step). Presence of other maximums is an indicator of almost periodic firing of the neurons which arises from the suprathreshold mean and the small amplitude stochastic fluctuations of the input current. Results in **Figure 6** are presented for one forward unidirectional coupling and sample values of inhomogeneity and synaptic strength. The results for other parameters are similar while the system is in the main locking zone in the absence of noise. This result shows a drawback of the simplified models we have used: LIF neurons with pulsatile instantaneous couplings can be synchronized with zero phase lag even in the presence of frequency mismatch, which is revealed as a maximum in correlation at zero lag (one simulation time step) when a small amplitude noise is added. Both mismatch and delay (synaptic and axonal) can be source of phase lag, when the neurons are modeled by limit cycle oscillators and more realistic models are used for synaptic currents. Our results are still valid when such phase lags are small, of the order of the time bins in the calculation of the correlation.

Above results were obtained for bidirectional symmetric couplings or for one unidirectional coupling. To find the *best* configuration through which direct couplings can improve the performance of the system in the detection of a variable input correlation, we have tested mutual couplings with different ratios of forward -21 and backward -12 connections. While the synaptic cost (sum of two synaptic strengths) is kept constant, different configurations can be designed by changing the ratio of the coupling constants *r* = -21/-12 (**Figures 7A,B**). In the absence of mismatch, the best configuration is that which preserves symmetry, i.e., the best performance results with equal forward and backward couplings. On the other hand, in the presence of mismatch, an asymmetric arrangement of couplings in which the forward coupling (from the high frequency neuron) is larger, improves the performance of the system. Interestingly, asymmetric excitatory couplings in favor of backward coupling (from the low frequency neuron), significantly decreases the sensitivity of the system since it plays the role of desynchronizing coupling as discussed above.

## **4. DISCUSSION**

Both direct connections and common inputs can be sources of the correlated activity of neurons in the nervous system. Effect of direct connections is widely studied as a general problem in dynamical systems and in particalur in nervous systems (Kuramoto, 1991; Strogatz and Mirollo, 1991; Abbott and van Vreeswijk, 1993). Stochastic inputs are usually a source of temporal disorder but spatial order can be induced in a neuronal pool when the neurons share stochastic inputs from common sources (Binder and Powers, 2001; Türker and Powers, 2001, 2004). Because of the possible cooperative/competitive effects of common inputs and direct connections, interesting results can be

correlation *c* = 0 and *c* = 0.5. The results are shown for three different noise amplitudes (shown above each panel) and two different values of

the high-frequency to the low-frequency neuron. The difference between two curves at zero lag gives MCS which has been shown in the plots by arrows.

**(A,B)** MCS is plotted against sum of synaptic weights for **(A)** homogeneous system δ = 0, and **(B)** inhomogenous system δ = 0.1. Different curves are plotted for different ratios *r* of forward and backward couplings indicated in the legends. When the neurons have equal intrinsic firing rates, symmetric configuration *r* = 1 shows the best performance with a suitable choice of synaptic strengths. For inhomogeneous case when the imbalance of couplings is in favor of forward coupling (from the high frequency to low frequency neuron) the sensitivity is considerably improved. When the backward coupling is larger *r* < 1, the system performance is quite poor. As is shown in axes labels, MCS is calculated over the range [0, 0.5] of input correlation.

expected when they are concurrently present in a system (Ostojic et al., 2009; Ly and Ermentrout, 2010; Tabareau et al., 2010; Zilli and Hasselmo, 2010; Rosenbaum and Josic, 2011a; Ly et al., ´ 2012). In this study we have numerically inspected the effect of correlated stochastic inputs on the correlation of spike trains of two coupled LIF neurons. We have mainly focused on the correlation of spike trains when correlated small amplitude noises were imposed on a system of two coupled neurons, and the neurons were regularly and synchronously firing in the absence of noise. We have shown that such a system shows high sensitivity to the changes of input correlation, and therefore can be a suitable detector of the correlation in small amplitude noises. To study the system in a more general framework, we have considered neurons with different intrinsic firing rates. We have assumed neurons have equal membrane time constants, and inhomogeneity is imposed on the system by feeding the neurons with unequal suprathreshold constant currents. The inhomogeneity, determined by the difference in the mean input currents, along with synaptic strengths are the key-parameters that specify the response of the system to stochastic inputs.

While for uncoupled neurons the output correlation is a monotonically decreasing function of inhomogeneity, for coupled neurons with low noise amplitudes, spike trains correlation can be increased by increasing inhomogeneity in some ranges. This result holds for sufficiently small noise amplitudes and the system inherits this property from *n*:*m* locking zones for the autonomous system when there is no stochastic input present. This introduces inhomogeneity as an important parameter with non-trivial impact on the correlation of spike trains in coupled systems.

Another feature of the system is that the two sources of correlation, correlated inputs and direct excitatory connections, do not necessarily cooperate in the formation of correlated spike trains. For uncoupled neurons output correlation is a monotonically increasing function of input correlation and for weakly correlated inputs, the slope decreases with lowering noise amplitude (De La Rocha et al., 2007; Shea-Brown et al., 2008) and with increasing mismatch. With different choices of the synaptic strengths and the inhomogeneity, it is possible to change functional form of correlation transfer (dependence of output correlation to the input correlation) and design a system with different sensitivity to the input correlation. In particular, it is possible to design a system with negative mean slope of correlation transfer, showing a case with destructive effect of common noises on the correlation of spike trains, or a system with maximum sensitivity to the changes in input correlation in a given range by maximizing the slope of correlation transfer. The latter proposes that direct connections can increase the sensitivity of the system to the correlation of the neuron's stochastic inputs, especially when the noises are small amplitude. We have further shown that for a homogeneous system (where the neurons have equal intrinsic firing rates), the best configuration of the couplings which maximizes the mean sensitivity of the system in a given range, is a symmetric configuration with equal coupling constants. On the other hand, in the presence of inhomogeneity, an asymmetric configuration in which the synaptic constant from the high frequency neuron to the low frequency neuron is larger, improves the sensitivity. In either case, there is an optimum value of the synaptic constant which maximizes the sensitivity.

Competitive learning through conventional spike timingdependent plasticity (STDP) in feed-forward networks leads to the potentiation of the synapses which convey correlated data and depression of those with uncorrelated activity (Babadi and Abbott, 2010). How does STDP change the lateral connections transverse to the path of data flow? It has been shown that in the recurrent networks, asymmetric connections arise through STDP and in the presence of inhomogeneity, such an asymmetric change is in favor of the connection from the high frequency to the low frequency neuron (Takahashi et al., 2009; Bayati and Valizadeh, 2012). Our results show that asymmetric connections can enhance the performance of inhomogeneous systems in the detection of input correlation, and interestingly such an optimum configuration of connections emerges through STDP (with asymmetric profile) in inhomogeneous neuronal pools (Bayati and Valizadeh, 2012).

Type of neuronal excitability can also affect the correlation transfer in neuronal pools (Galán et al., 2008; Abouzeid and Ermentrout, 2009; Barreiro et al., 2010). Phase resetting curve characterizes how small perturbations influence the oscillator's subsequent timing or phase. It has been recently shown that uncoupled type-II neurons with both negative and positive regions in their PRC transfer correlations more faithfully when the correlation is calculated over short time bins (Abouzeid and Ermentrout, 2011). Since the phase of a LIF neuron always advances in response to the external pulses, the results for LIF neurons are likely to apply for type-I neurons.

Correlation of spike trains over such small time bins that we have used *T* = 0.5 ms, is a measure of (almost) precise alignment of the action potentials. Similar results were obtained

## **REFERENCES**

and extracellular features. *J. Neurophysiol.* 92, 600–608. doi: 10.1152/jn.01170.2003

when we repeated the experiments with *T* = 1 ms but we expect qualitatively different results when the correlation of the spike counts is measured over the time scales comparable, or larger than the mean inter-spike interval. Less sensitivity to the inhomogeneity is expected when the correlation is evaluated over large time bins, but the effect of direct couplings warrants further studies to find out if correlation in small amplitude stochastic inputs can be revealed in co-variation of spike trains of coupled neurons over large time scales.

## **ACKNOWLEDGMENTS**

Authors gratefully acknowledge Bahman Farnoudi for proof reading of the manuscript, and the reviewers Germán Mato and Zachary P. Kilpatrick, for careful reading of the manuscript and giving valuable comments.

with the spatial correlation of the stochastic input stimulus. *Phys. Rev. Lett.* 93:48101. doi: 10.1103/PhysRevLett.93.048101

J. J. (1999). Synchronization of noisy systems by stochastic signals. *Phys. Rev. E* 60:284. doi: 10.1103/PhysRevE.60.284

population. *Nature* 440, 1007–1012. doi: 10.1038/nature04701

protects from noise. *PLoS Comput. Biol.* 6:e1000637. doi: 10.1371/journal.pcbi.1000637

*PLoS Compu. Biol.* 8:e1002306. doi: 10.1371/journal.pcbi.1002306

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

*Received: 29 April 2013; accepted: 24 July 2013; published online: 14 August 2013. Citation: Bolhasani E, Azizi Y and Valizadeh A (2013) Direct connections assist neurons to detect correlation in small amplitude noises. Front. Comput. Neurosci. 7:108. doi: 10.3389/fncom. 2013.00108*

*Copyright © 2013 Bolhasani, Azizi and Valizadeh. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

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

#### *James Trousdale1 \*, Yu Hu2, Eric Shea-Brown2,3 and Krešimir Josic´ 1,4*

*1 Department of Mathematics, University of Houston, Houston, TX, USA*

*2 Department of Applied Mathematics, University of Washington, Seattle, WA, USA*

*3 Program in Neurobiology and Behavior, University of Washington, Seattle, WA, USA*

*4 Department of Biology and Biochemistry, University of Houston, Houston, TX, USA*

#### *Edited by:*

*Robert Rosenbaum, University of Pittsburgh, USA*

#### *Reviewed by:*

*John A. Hertz, Niels Bohr Institute, Denmark Sonja Grün, Research Center Juelich, Germany*

#### *\*Correspondence:*

*James Trousdale, Department of Mathematics, University of Houston, 641 PGH Building, Houston, TX 77204-3008, USA e-mail<EMAIL_ADDRESS>
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.

**Keywords: correlations, spiking neurons, neuronal networks, cumulant, neuronal modeling, neuronal network model, point processes**

## **1. 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., 2006, 2009; Ohiorhenuan 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., 2010, 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 Josic, 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 Josic, 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.

## **2. 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., 2007, 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

**FIGURE 1 | (A)** Raster plot of event times for an example multivariate Poisson process **X** = (*X*1,..., *X*6) generated using the methods presented below. This model exhibits independent marginal events (blue) and population-level, chain-like events (red). **(B)** Some second order population cumulant densities (i.e., second order correlation between individual unit activities and population activity) for this model (Luczak et al., 2013). Greater mass to the right (resp. left) of τ = 0 indicates that the cell tends to lead (resp. follow) in pairwise-correlated events. **(C)** Third-order cross-cumulant density for processes *X*1, *X*2, *X*3. The quantity κ**X** 123(τ1, τ2) yields the probability of observing spikes in cells 2 and 3 at an offset τ1, τ2 from a spike in cell 1, respectively, in excess of what would be predicted from the first and second order cumulant structure. All quantities are precisely defined in the Methods. Note: system parameters necessary to reproduce results are given in the Appendix for all figures.

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.

## **2.1. 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 2D, the set of all subsets of D. We define a probability distribution *p* = (*pD*)*D* ⊂ D, assigning a probability to each possible marking, *D*. As we will see, *pD* 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 *QD*, giving a family of shift (jitter) distributions (*QD*)*D* ⊂ D. Each (*QD*) is a distribution over R*N*.

The rate λ, the distribution *p* over the markings, and the family of jitter distributions (*QD*)*D* ⊂ D, define a vector **X** = (*X*1,..., *XN*) 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**).

distribution *pD*, and project the event at time *tj* to the subsets with indices in *Dj* . The legend indicates the colors assigned to three possible markings in this example. **(C)** Step 3: For each pair (*tj* , *Dj* ) generated in the previous two steps, draw **Y***j* from *QDj* , and shift the event times in the daughter processes by the corresponding values *Y j i* .

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 *QD* on the chosen marking. This allows for greater flexibility in setting the temporal cumulant structure.

## **2.2. EXAMPLES**

## *2.2.1. 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 *Xi* is obtained by merging an independent Poisson process with a common, global Poisson process. That is,

$$X\_i(\cdot) = Z\_i(\cdot) + Z\_c(\cdot), \quad i = 1, \ldots, N,$$

where *Zc* and each *Zi* 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

$$p\_D = \begin{cases} \frac{\lambda\_i}{\lambda} & D = \{i\} \\ \frac{\lambda\_\zeta}{\lambda} & D = \mathbb{D} \\ 0 & \text{otherwise} \end{cases}$$

Note that if λ*c* = 0, each spike will be assigned to a different process *Xi*, resulting in *N* independent Poisson processes. Lastly, each shift distribution is equal to a delta distribution at zero in every coordinate (i.e., *qD*(*y*1,..., *yN*) ≡ *N i* = 1 δ(*yi*) 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 *i*th 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 *pD* = |*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.

#### *2.2.2. 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,..., *XN*) of GTaS type, where *N* ≥ 4. We restrict the marking distribution so that *pD* ≡ 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 *QD* ∼ *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

$$f(t|\alpha) = \alpha e^{-\alpha t} \Theta(t),$$

where (*t*) is the Heaviside step function. We generate random vectors **Y** ∼ *Q*D according to the following rule, for each *i* = 1,..., *N*:

1. Generate independent random variables *Ti* ∼ Exp(α*i*) where α*i* > 0. 2. Set *Yi* = *i j* = 1 *Tj*.

$$\text{In particular, note that these shift times satisfy } Y\_N \ge \dots \ge Y\_2 \ge 1$$

*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

**FIGURE 3 | An example of a six dimensional GTaS model exhibiting synfire-like cascading firing patterns. (A)** A raster-plot of spiking activity over a 100 ms window. Blue spikes indicate either marginal or pairwise events (i.e., events corresponding to markings for sets *D* ⊂ D with |*D*| ≤ 2). 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 (*cN* 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.

order cross-cumulant density for the process (*i*, *j*), is given by

$$\kappa\_{ij}^{\mathbf{X}}(\mathfrak{r}) = c\_{ij}^2(\mathfrak{r}) + c\_{ij}^N(\mathfrak{r}),\tag{1}$$

where

$$c\_{ij}^2(\mathbf{r}) = \lambda p\_{\{i,j\}} \int q\_{\{i,j\}}^{\{i,j\}}(t, t + \mathbf{r}) dt,$$

$$c\_{ij}^N(\mathbf{r}) = \lambda p\_{\mathbb{D}} \int q\_{\mathbb{D}}^{\{i,j\}}(t, t + \mathbf{r}) dt \tag{2}$$

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 *qD D* indicate the densities associated with the distribution *QD*, 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,

$$f(\boldsymbol{\jmath}\_1, \ldots, \boldsymbol{\jmath}\_N) = f(\boldsymbol{\jmath}\_N | \boldsymbol{\jmath}\_1, \ldots, \boldsymbol{\jmath}\_{N-1}) f(\boldsymbol{\jmath}\_{N-1} | \boldsymbol{\jmath}\_1, \ldots, \boldsymbol{\jmath}\_{N-2}) \cdots$$

$$\cdot f(\boldsymbol{\jmath}\_2 | \boldsymbol{\jmath}\_1) f(\boldsymbol{\jmath}\_1) . \tag{3}$$

Since *Y*1 ∼ Exp(α1), we have *f*(*y*1) = exp −α1*y*1 (*y*1), where (*y*) is the Heaviside step function. Further, as (*Yi* − *Yi* − 1)|(*Y*1,..., *Yi* − 1) ∼ Exp(α*i*) for *i* ≥ 2, the conditional densities of the *yi*'s take the form

$$f(\boldsymbol{y}\_{i}|\boldsymbol{y}\_{1},\ldots,\boldsymbol{y}\_{i-1}) = f(\boldsymbol{y}\_{i}|\boldsymbol{y}\_{i-1}) = \alpha\_{i} \exp\left[-\alpha\_{i}(\boldsymbol{y}\_{i} - \boldsymbol{y}\_{i-1})\right],$$

$$\cdot \oplus (\boldsymbol{y}\_{i} - \boldsymbol{y}\_{i-1}), \quad i \ge 2.$$

Substituting this in to the identity Equation (3), we have

$$q\_{\mathbb{D}}(\boldsymbol{\chi}\_{1},\ldots,\boldsymbol{\chi}\_{N}) = \begin{cases} \boldsymbol{\alpha}\_{1}\exp\left[-\boldsymbol{\alpha}\_{1}\boldsymbol{\chi}\_{1}\right] \prod\_{i=2}^{N} \boldsymbol{\alpha}\_{i} & \boldsymbol{\chi}\_{N} \ge \ldots \ge \\\ \cdot \ \exp\left[-\boldsymbol{\alpha}\_{i}(\boldsymbol{\chi}\_{i} - \boldsymbol{\chi}\_{i-1})\right] & \boldsymbol{\chi}\_{2} \ge \boldsymbol{\chi}\_{1} \ge \boldsymbol{0} \\\ 0 & \text{otherwise} \end{cases} (4)$$

Using Theorem 1.1 (Equation A8) we obtain the *N*th order crosscumulant density (see the Methods),

$$\begin{aligned} \mathsf{x}\_{1}^{\mathsf{X}} & \quad \mathsf{x}\_{1} \mathsf{x}\_{1} \left( \mathsf{r}\_{1}, \dots, \mathsf{r}\_{N-1} \right) \\ &= \lambda p\_{\mathbb{D}} \int q\_{\mathbb{D}} (t, \ t + \mathsf{r}\_{1}, \dots, t + \mathsf{r}\_{N-1}) dt \\ &= \lambda p\_{\mathbb{D}} \cdot \begin{cases} \prod\_{i=1}^{N-1} \alpha\_{i+1} & \mathsf{r}\_{i} \geq \mathsf{r}\_{i-1} \\ \cdot \exp \left[ -\alpha\_{i+1} (\mathsf{r}\_{i} - \mathsf{r}\_{i-1}) \right] & i = 1, \dots, N-1, \\ 0 & \text{otherwise} \end{cases} \end{aligned} \tag{5}$$

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

$$
\omega\_{13}^{N}(\mathfrak{r}) = \lambda p\_{\mathbb{D}} \Theta(\mathfrak{r}) \cdot \begin{cases}
\frac{\alpha\_{2}\alpha\_{3}}{\alpha\_{3} - \alpha\_{2}} \left\{ \exp\left[ -\alpha\_{2}\mathfrak{r} \right] - \exp\left[ -\alpha\_{3}\mathfrak{r} \right] \right\} & \alpha\_{2} \neq \alpha\_{3} \\
\alpha\_{2}\alpha\_{3}\mathfrak{r} \exp\left[ -\alpha\_{2}\mathfrak{r} \right] & \alpha\_{2} = \alpha\_{3} \\
\end{cases}.
$$

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,

$$\begin{aligned} \mathsf{x}\_{124}^{\mathsf{X}}(\mathsf{\tau}\_{1},\mathsf{\tau}\_{2}) &= \lambda \mathsf{p}\_{\mathsf{D}} \\ \cdot \begin{cases} \frac{\alpha\_{2}\alpha\_{3}\alpha\_{4}}{\alpha\_{4} - \alpha\_{3}} \exp\left[-\alpha\_{2}\mathsf{\tau}\_{1}\right] \left\{\exp\left[-\alpha\_{3}(\mathsf{\tau}\_{2} - \mathsf{\tau}\_{1})\right] \right. \\\ -\; \exp\left[-\alpha\_{4}(\mathsf{\tau}\_{2} - \mathsf{\tau}\_{1})\right] \left\{ \\\ \alpha\_{2}\alpha\_{3}\alpha\_{4}(\mathsf{\tau}\_{2} - \mathsf{\tau}\_{1})\exp\left[-\alpha\_{2}\mathsf{\tau}\_{1} - \alpha\_{3}(\mathsf{\tau}\_{2} - \mathsf{\tau}\_{1})\right] \end{cases} & \quad \alpha\_{3} \neq \alpha\_{4}, \end{aligned} (6)$$

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 (Grün and Rotter, 2010).

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)

$$\kappa\_{i,\text{pop}}^{\mathbf{X}}(\mathfrak{r}) = \sum\_{j \neq i} \kappa\_{ij}^{\mathbf{X}}(\mathfrak{r}).$$

This function is linearly related to the spike-triggered average of the population activity conditioned on that of cell *i*. In **Figure 4** 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

$$
\kappa\_{ij,\text{pop}}^{\mathbf{X}}(\mathfrak{r}\_1,\mathfrak{r}\_2) = \sum\_{k \neq i,j} \kappa\_{ijk}^{\mathbf{X}}(\mathfrak{r}\_1,\mathfrak{r}\_2).
$$

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

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.

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., 2010, 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.

#### *2.2.3. 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 *Vi* obeying

$$\frac{dV\_i}{dt} = -V\_i + \sum\_{j=1}^{N} \omega\_{ij} (\mathbf{F} \ast \mathbf{z}\_j)(t) + \mathbf{w}^{\text{in}} \mathbf{x}\_i(t), \quad i = 1, \dots, N. (7)$$

When the membrane potential of cell *i* reaches a threshold *V*th, an output spike is recorded and the membrane potential is reset to zero, after which evolution of *Vi* resumes the dynamics in Equation (7). Here *wij* 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 *F* = τsyn−1*e*−(*t* − τd)/τsyn(*t* − τd) is a delayed, unit-area exponential synaptic kernel with time-constant τsyn and delay τd. The output of the *i*th neuron is

$$z\_i(t) = \sum\_j \delta(t - t\_i^j),$$

where *t j i* is the time of the *j*th spike of neuron *i*. In addition, the input {*xi*}*N i* = 1 is

$$\mathbf{x}\_i(t) = \sum\_j \mathbb{S}(t - s\_i^j),$$

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

$$w\_{ij} = \begin{cases} \nu^{\text{syn}} & i - j \mod N \equiv 1 \text{ or } N - 1 \\ 0 & \text{otherwise} \end{cases}$$

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.

#### **2.3. 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,..., *XN*) on R*N*, given a collection of Borel subsets *Ai* ∈ *B*(R),*i* = 1,..., *N*, then **X**(*A*1 ×···× *AN*) = (*X*1(*A*1), . . . , *XN* (*AN*)) ∈ N*N* is a random vector where the value of each coordinate *i* indicates the (random) number of points which fall inside the set *Ai*. Note that the GTaS model defines processes that are marginally Poisson. All GTaS model parameters and related quantities are defined in **Table 1**.

**Table 1 | Common notation used in the text.**

For each *D* ⊂ D = {1,..., *N*}, define the tail probability *p*¯*D* by

$$\bar{p}\_{\mathcal{D}} = \sum\_{D \subset D' \subset \mathcal{D}} p\_{D'}.\tag{8}$$

Since *pD* 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 *Xi* with probability *p*¯{*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 *QD* onto dimension *i*. Thus, the events in the marginal process *Xi* are shifted in an independent and identically distributed (IID) manner according to the mixture distribution *Qi* given by

$$Q\_i = \frac{\sum\_{D \ni i} p\_D Q\_D^{\{i\}}}{\sum\_{D \ni i} p\_D}.$$

Note that IID shifting of the event times of a Poisson process generates another Poisson process of identical rate. Thus, the process *Xi* 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

$$
\begin{pmatrix} X\_1(A\_1) \\ \vdots \\ X\_N(A\_N) \end{pmatrix} =\_{\text{distr}} \begin{pmatrix} \sum\_{D \ni 1} \xi(D; A\_1, \dots, A\_N) \\ \vdots \\ \sum\_{D \ni N} \xi(D; A\_1, \dots, A\_N) \end{pmatrix} . \tag{9}
$$

Here, each ξ(*D*; *A*1,..., *AN* ) 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 *Ai*. 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 *Ai*.

Suppose we are interested in calculating dependencies among a subset of daughter processes, {*Xij* }*ij*∈*D*¯ for some set *D*¯ ⊂ D, consisting of |*D*¯ | = *k* distinct members of the collection of counting processes **X**. Then the following alternative representation will be useful:

$$\begin{pmatrix} X\_{i\_1}(A\_{i\_1}) \\ \vdots \\ X\_{i\_k}(A\_{i\_k}) \end{pmatrix} =\_{\text{distr}} \begin{pmatrix} \sum\_{i\_1 \in D \subset \bar{D}} \xi\_D(A\_1, \dots, A\_N) \\ \vdots \\ \sum\_{i\_k \in D \subset \bar{D}} \xi\_D(A\_1, \dots, A\_N) \end{pmatrix} \tag{10}$$

where

$$\xi\_D(A\_1, \ldots, A\_N) = \sum\_{\substack{D' \supset D \\ (\vec{D} \backslash D) \cap D' = \emptyset}} \xi(D'; A\_1, \ldots, A\_N).$$

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 *ij* and contained within *D*¯ . 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 subset *D*¯ 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.

#### *2.3.1. 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)

$$
\kappa\_{12}^{\mathbf{X}}(\mathfrak{r}) = \lambda p \int f(t) f(t+\mathfrak{r}) dt = \lambda p (f \times f)(\mathfrak{r}).
$$

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**)

$$
\begin{pmatrix} X\_i(A\_i) \\ X\_j(A\_j) \end{pmatrix} =\_{\text{distr}} \begin{pmatrix} \xi\_{\{i,j\}}(A\_i, A\_j) + \xi\_{\{i\}}(A\_i) \\ \xi\_{\{i,j\}}(A\_i, A\_j) + \xi\_{\{j\}}(A\_j) \end{pmatrix}. \tag{11}
$$

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. First, we have

$$\begin{split} \mathsf{cov}[X\_{i}(A\_{i}), X\_{j}(A\_{j})] &= \mathsf{\kappa}[X\_{i}(A\_{i}), X\_{j}(A\_{j})] \\ &= \mathsf{\kappa}[\xi\_{\{i,j\}}(A\_{i}, A\_{j}), \xi\_{\{i,j\}}(A\_{i}, A\_{j})] \\ &+ \mathsf{\kappa}[\xi\_{\{i\}}(A\_{i}), \xi\_{\{i,j\}}(A\_{i}, A\_{j})] \\ &+ \mathsf{\kappa}[\xi\_{\{i,j\}}(A\_{i}, A\_{j}), \xi\_{\{j\}}(A\_{j})] \\ &+ \mathsf{\kappa}[\xi\_{\{i\}}(A\_{i}), \xi\_{\{j\}}(A\_{j})] \\ &= \mathsf{\kappa}\_{2}[\xi\_{\{i,j\}}(A\_{i}, A\_{j})] + 0 \\ &= \mathsf{E}[\xi\_{\{i,j\}}(A\_{i}, A\_{j})]. \end{split}$$

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):

$$\begin{aligned} \text{cov}[X\_i(A\_i), X\_j(A\_j)] \\ = \lambda \sum\_{D' \supset \{i, j\}} p\_{D'} \int P\left(t + Y\_i \in A\_i, t + Y\_j \in A\_j \mid \mathbf{Y} \sim \mathbf{Q}\_{D'}\right) dt. \end{aligned} \tag{12}$$

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 *QD* ≡ *Q*, Equation (12) reduces to the solution given in Bäuerle and Grübel (2005). Using the tail probabilities defined in Equation (8), if *QD* ≡ *Q* for all *D*, the integral in Equation (12) no longer depends on the subset *D* , and the equation may be written as

$$\begin{aligned} &\mathbb{E}\mathbb{cov}[X\_i(A\_i), X\_j(A\_j)] \\ &= \lambda \bar{p}\_{\{i,j\}} \int P\left(t + Y\_i \in A\_i, t + Y\_j \in A\_j \mid \mathbf{Y} \sim Q\right) dt. \end{aligned}$$

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

$$\begin{aligned} \mathsf{x}\_{ij}^{\mathbf{X}}(\mathbf{r}) &= \lim\_{\Delta t \to 0} \frac{\mathsf{cov}[X\_i(\mathbf{0}, \Delta t)), X\_j(\mathbf{r}, \tau + \Delta t))}{\Delta t^2} \\ &= \lambda \sum\_{D' \supset \{i, j\}} p\_{D'} \\ &\quad \int \lim\_{\Delta t \to 0} \frac{P\left(t + Y\_i \in [0, \Delta t), t + Y\_j \in [\tau, \tau + \Delta t) \mid \mathbf{Y} \sim Q\_D\right)}{\Delta t^2} dt, \end{aligned} \tag{13}$$

Before continuing, we note that given a random vector **Y** = (*Y*1,..., *YN*) ∼ *Q*, where *Q* has density *q*(*y*1,..., *yN*), the vector **Z** = (*Y*2 − *Y*1,..., *YN* − *Y*1) has density *qZ* given by

$$q\_Z(\mathbf{r}\_1, \dots, \mathbf{r}\_{N-1}) = \int q(t, t + \mathbf{r}\_1, \dots, t + \mathbf{r}\_{N-1}) dt. \quad \text{(14)}$$

Assuming that the distributions *QD* have densities *qD* , and denoting by *q* {*i*, *j*} *D* the bivariate marginal density of the variables *Yi*, *Yj* under *QD* , we have that

$$\kappa\_{ij}^{X}(\mathbf{r}) = \lambda \sum\_{D' \supset \{i, j\}} p\_{D'} \int q\_{D'}^{\{i, j\}}(t, t + \mathbf{r}) dt. \tag{15}$$

According to Equation (14), the integrals present in Equation (15) are simply the densities of the variables *Yj* − *Yi*, where **Y** ∼ *QD* .

Thus κ**X** *ij*(τ), which captures the additional probability for events in the marginal processes *Xi* and *Xj* 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" λ*pD* 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 arrangement specifically, the shift applied to the event as attributed to *Xi* precedes that applied to the event mapped to *Xj* 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.

#### *2.3.2. 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**)

$$
\begin{pmatrix} X\_i(A\_i) \\ X\_j(A\_j) \\ X\_k(A\_k) \end{pmatrix} =\_{\text{distr}} \begin{pmatrix} \xi\_{\{i,j,k\}} + \xi\_{\{i,j\}} + \xi\_{\{i,k\}} + \xi\_{\{i\}} \\ \xi\_{\{i,j,k\}} + \xi\_{\{i,j\}} + \xi\_{\{j,k\}} + \xi\_{\{j\}} \\ \xi\_{\{i,j,k\}} + \xi\_{\{i,k\}} + \xi\_{\{j,k\}} + \xi\_{\{k\}} \end{pmatrix}, \tag{16}
$$

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

$$\kappa(X\_i(A\_i), X\_j(A\_j), X\_k(A\_k)) = \kappa\_\Im(\zeta\_{\{i,j,k\}}) = \mathbf{E}[\zeta\_{\{i,j,k\}}] \dots$$

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

$$\begin{aligned} \mathbb{K}(X\_i(A\_i), X\_j(A\_j), X\_k(A\_k)) &= \lambda \sum\_{D' \supset \{i, j, k\}} p\_{D'} \int P(t + Y\_i \in A\_i, \\\\ t + \ \ Y\_j \in A\_j, t + Y\_k \in A\_k \mid \mathbf{Y} \sim Q\_{D'}) dt. \end{aligned}$$

The third cross-cumulant density is then given similarly to the second order function by

$$\kappa\_{ijk}^{\mathbf{X}}(\mathfrak{r}\_1, \mathfrak{r}\_2) = \lambda \sum\_{D' \supset \{i, j, k\}} p\_{D'} \int q\_{D'}^{\{i, j, k\}}(t, t + \mathfrak{r}\_1, t + \mathfrak{r}\_2) dt.$$

Here, we have again assumed the existence of densities *qD* , and denoted by *q* {*i*, *j*, *k*} *D* the joint marginal density of the variables *Yi*, *Yj*, *Yk* under *qD* . The integrals appearing in the expression for the third order cross-cumulant density are the probability densities of the vectors (*Yj* − *Yi*, *Yk* − *Yi*), where **Y** ∼ *QD* .

## *2.3.3. 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.

**Theorem 1.1.** *Let* **X** *be a joint counting process of GTaS type with total intensity* λ*, marking distribution* (*pD*)*D* ⊂ D*, and family of shift distributions* (*QD*)*D* ⊂ D*. Let A*1,..., *Ak be arbitrary sets in B*(R)*, and D*¯ = {*i*1,...,*ik*} ⊂ D *with* |*D*¯ | = *k. The cross-cumulant of the counting processes may be written*

$$\begin{aligned} \mathbb{K}(X\_{\dot{i}\_1}(A\_1), \dots, X\_{\dot{i}\_k}(A\_k)) \\ = & \lambda \sum\_{D' \supset \bar{D}} p\_{D'} \int P(t\mathbf{1} + \mathbf{Y}^{\bar{D}} \in A\_1 \times \dots \times A\_k | \mathbf{Y} \sim Q\_{D'}) dt \end{aligned} \tag{17}$$

*where* **Y***D*¯ *represents the projection of the random vector* **Y** *onto the dimensions indicated by the members of the set D. Furthermore,* ¯ *assuming that the shift distributions possess densities* (*qD*)*D* ⊂ D*, the cross-cumulant density is given by*

$$\begin{aligned} \mathsf{w}\_{i\_1\ldots i\_k}^X(\mathfrak{r}\_1,\ldots,\mathfrak{r}\_{k-1}) \\ = \lambda \sum\_{D'\supset\bar{D}} p\_{D'} \int q\_{D'}^{\bar{D}}(t, t+\mathfrak{r}\_1, \ldots, t+\mathfrak{r}\_{k-1}) dt,\end{aligned} \tag{18}$$

*where qD*¯ *D indicates the kth order joint marginal density of qD in the dimensions of D.* ¯

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

$$\gamma\_{i\_1\dots i\_k}^{\mathbf{X}}(\infty) = \int \cdots \int \kappa\_{i\_1\dots i\_k}^{\mathbf{X}}(\mathbf{r}\_1,\dots,\mathbf{r}\_{k-1})d\mathbf{r}\_{k-1}\cdots d\mathbf{r}\_1$$

$$= \lambda \sum\_{D'\supset D} p\_{D'} \cdot 1 = \lambda \bar{p}\_D.\tag{19}$$

This shows that the infinite time window cumulants for a GTaS process are non-increasing with respect to the ordering of sets, i.e.,

$$
\boldsymbol{\gamma}\_{i\_1\cdots i\_k}^{\mathbf{X}}(\infty) \ge \boldsymbol{\gamma}\_{i\_1\cdots i\_k i\_{k+1}}^{\mathbf{X}}(\infty).
$$

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)

$$\begin{aligned} &\mathbb{K}(X\_1(A\_1), X\_1(A\_2), X\_3(A\_3)) \\ &= \mathbb{K}\_{(2,1)}(X\_1(A\_1 \cap A\_2), X\_3(A\_3)) \quad \text{if} \quad \mathbf{X} \sim \text{GTaS.} \end{aligned} \tag{20}$$

For a general counting process **X**, it may be shown that

$$
\kappa\_{113}^{\mathbf{X}}(\mathbf{r}\_1, \mathbf{r}\_2) = \delta(\mathbf{r}\_1)\kappa\_{13}^{\mathbf{X}}(\mathbf{r}\_2) + \text{ \text{\textquotedbl{}non-singular\textquotedbl{}continuous\textquotedbl{}} }\text{\textquotedbl{}}\tag{21}
$$

In addition, the second order auto-cumulant density may be written (Cox and Isham, 1980)

$$
\kappa\_{ii}^{\mathbf{X}}(\mathbf{r}) = r\_i \delta(\mathbf{r}) + \text{ "non-singular contributions"},
$$

where *ri* 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.

## **3. 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. 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.

## **4. METHODS**

#### **4.1. 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,..., *XN*) by

$$\log(t\_1, \dots, t\_N) = \log\left(\mathbf{E}\left[\exp\left(\sum\_{j=1}^N t\_j X\_j\right)\right]\right).$$

The **r**-cross-cumulant of the vector **X** is given by

$$\kappa\_{\mathbf{r}}(\mathbf{X}) = \left. \frac{\partial^{|\mathbf{r}|}}{\partial t\_1^{r\_1} \cdots \partial t\_N^{r\_N}} g(t\_1, \dots, t\_N) \right|\_{t\_1 = \cdots = t\_N = 0}.$$

where **r** = (*r*1,...,*rN*) is a *N*-vector of positive integers, and |**r**| = -*N i* = 1 *ri*. 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,..., *XN*) as

$$\kappa(\mathbf{X}) := \kappa\_1(\mathbf{X}) = \left. \frac{\partial^N}{\partial t\_1 \cdots \partial t\_N} \mathbf{g}(t\_1, \dots, t\_N) \right|\_{t\_1 = \cdots = t\_N = 0}$$

where **1** = (1,..., 1).

This relationship may be expressed in combinatorial form:

$$\kappa(X\_1, \dots, X\_N) = \sum\_{\pi} (|\pi| - 1)! (-1)^{|\pi| - 1} \prod\_{B \in \pi} \mathbf{E} \left[ \prod\_{i \in B} X\_i \right] \tag{22}$$

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

$$g(t\_1, \ldots, t\_N) = \sum\_{\mathbf{r}} \frac{\kappa\_{\mathbf{r}}(X\_1, \ldots, X\_N)}{\mathbf{r}!} x\_1^{r\_1} \cdots x\_d^{r\_N} \quad \text{with} \quad \mathbf{r}!$$

$$\mathbf{r}! = \prod\_{i=1}^N r\_i!,$$

similarly expanding the moment generating function *M*(*t*) = *eg*(*t*) , and matching the polynomial coefficients. Note that the *n*th cumulant κ*n* of a random variable *X* may be expressed as a joint cumulant via

$$\kappa\_n(X) = \kappa(\underbrace{X, \dots, X}\_{n \text{ copies of } X})\,\,\,^\vee$$

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*,{*Zi*} *N i* = 2, we have

$$\begin{aligned} \kappa(aX + bY, Z\_2, \dots, Z\_N) &= a\kappa(X, Z\_2, \dots, Z\_N) \\ &+ b\kappa(Y, Z\_2, \dots, Z\_N). \end{aligned}$$

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,..., *XN*1 } and {*Y*1,..., *YN*2 } are sets of random variables such that each *Xi* is independent from each *Yj*, then

$$\mathbf{x}\_{\left(\mathbf{r}\_{X},\mathbf{r}\_{Y}\right)}\left(X\_{1},\ldots,X\_{N\_{1}},Y\_{1},\ldots,Y\_{N\_{2}}\right) = \mathbf{0}$$

$$\text{for all } \mathbf{r}\_{X} \in \mathbb{N}\_{+}^{N\_{1}}, \mathbf{r}\_{Y} \in \mathbb{N}\_{+}^{N\_{2}}.$$

To exhibit another key property of cumulants, consider a 4 vector **X** = (*X*1, *X*2, *X*3, *X*4) with non-zero fourth cumulant and a random variable *Z* independent of each *Xi*. Define **Y** = (*X*1 + *Z*, *X*2 + *Z*, *X*3 + *Z*, *X*4). Using properties (C1), (C2) above, it follows that

$$
\kappa(Y\_1, Y\_2, Y\_3) = \kappa(X\_1, X\_2, X\_3) + \kappa\_3(Z).
$$

On the other hand, it is also true that

$$
\kappa(\mathbf{Y}) = \kappa(\mathbf{X}),
$$

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 *r*th 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 κ(*Xi*) = **E**[*Xi*], and κ(*Xi*, *Xj*) = **cov** *Xi*, *Xj* . In addition, the third cumulant, like the second, is equal to the corresponding central moment:

$$\kappa(X\_i, X\_j, X\_k) = \mathbf{E}[(X\_i - \mathbf{E}[X\_i])(X\_j - \mathbf{E}[X\_j])(X\_k - \mathbf{E}[X\_k])\,\big]\,.$$

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:

$$\begin{split} \mathbb{K}(X\_i, X\_j, X\_k, X\_l) \\ &= E[(X\_i - \mathbb{E}[X\_i])(X\_j - \mathbb{E}[X\_j])(X\_k - \mathbb{E}[X\_k])(X\_l - \mathbb{E}[X\_l])] \\ &- \mathbf{cov}[X\_i, X\_j] \left[ \mathbf{cov}[X\_k, X\_l] - \mathbf{cov}[X\_i, X\_k] \right. \\ &\left. - \mathbf{cov}[X\_i, X\_l] \left[ \mathbf{cov}[X\_j, X\_k] \right. \right] . \end{split} \tag{23}$$

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).

#### **4.2. 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,..., *XN*) where each *Xi* is marginally Poisson.

Formally, a counting process **X** is an integer-valued random measure on *B*(R*N*). Evaluated on subset *A*1 ×···× *AN* of *B*(R*N*), the random vector (*X*1(*A*1), . . . , *XN*(*AN*)) counts events in *d* distinct categories whose times of occurrence fall in to the sets *Ai*. 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 *Ai* ∈ *B*(R), the random variable *Xi*(*Ai*) follows a Poisson distribution with mean λ*i* (*Ai*), where is the Lebesgue measure on R, and λ*i* is the (constant) rate for the *i*th process. The processes under consideration will further satisfy a joint stationarity condition, namely that the distribution of the vector (*X*1(*A*1 + *t*), . . . , *XN* (*AN* + *t*)) does not depend on *t*, where *Ai* + *t* denotes the translated set {*a* + *t* : *a* ∈ *Ai*}.

We now consider some common measures of temporal dependence for jointly stationary vector counting processes. We will refer to the quantity *Xi*[0, *T*] as the *spike count* of process *i* over [0, *T*]. The quantity γ**X** *i*1···*ik* (*T*) (which we will refer to as a *spike count cumulant*) is given by