Home | About | Journals | Submit | Contact Us | Français |

**|**Front Comput Neurosci**|**v.4; 2010**|**PMC2995522

Formats

Article sections

Authors

Related links

Front Comput Neurosci. 2010; 4: 147.

Published online 2010 November 19. doi: 10.3389/fncom.2010.00147

PMCID: PMC2995522

Edited by: Jakob H. Macke, University College London, UK

Reviewed by: Taro Toyoizumi, RIKEN Brain Science Institute, Japan; Kresimir Josic, University of Houston, USA

*Correspondence: Shy Shoham, Faculty of Biomedical Engineering, Technion – Israel Institute of Technology, Haifa 32000, Israel. e-mail: li.ca.noinhcet.mb@mahohss

Received 2009 December 8; Accepted 2010 October 25.

Copyright © 2010 Krumin, Reutsky and Shoham.

This is an open-access article subject to an exclusive license agreement between the authors and the Frontiers Research Foundation, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are credited.

This article has been cited by other articles in PMC.

The correlation structure of neural activity is believed to play a major role in the encoding and possibly the decoding of information in neural populations. Recently, several methods were developed for exactly controlling the correlation structure of multi-channel synthetic spike trains (Brette, 2009; Krumin and Shoham, 2009; Macke et al., 2009; Gutnisky and Josic, 2010; Tchumatchenko et al., 2010) and, in a related work, correlation-based *analysis* of spike trains was used for blind identification of single-neuron models (Krumin et al., 2010), for identifying compact auto-regressive models for multi-channel spike trains, and for facilitating their causal network analysis (Krumin and Shoham, 2010). However, the diversity of correlation structures that can be explained by the feed-forward, non-recurrent, generative models used in these studies is limited. Hence, methods based on such models occasionally fail when analyzing correlation structures that are observed in neural activity. Here, we extend this framework by deriving closed-form expressions for the correlation structure of a more powerful multivariate self- and mutually exciting Hawkes model class that is driven by exogenous non-negative inputs. We demonstrate that the resulting Linear–Non-linear-Hawkes (LNH) framework is capable of capturing the dynamics of spike trains with a generally richer and more biologically relevant multi-correlation structure, and can be used to accurately estimate the Hawkes kernels or the correlation structure of external inputs in both simulated and real spike trains (recorded from visually stimulated mouse retinal ganglion cells). We conclude by discussing the method's limitations and the broader significance of strengthening the links between neural spike train analysis and classical system identification.

Linear system models enjoy a fundamental role in the analysis of a wide range of natural and engineered signals and processes (Kailath et al., 2000). Hawkes (Hawkes, 1971a,b; cf. Johnson, 1996) introduced the basic point processes equivalent of the linear auto-regressive and multi-channel auto-regressive process models, and derived expressions for their output correlations and spectral densities. The Hawkes model was later used as a model for neural activity in small networks of neurons (Brillinger, 1975, 1988; Brillinger et al., 1976; Chornoboy et al., 1988), where maximum likelihood (ML) parameter estimation procedures can be used to estimate the synaptic strengths between connected neurons, but where no external modulating processes were considered. Interestingly, the recent renaissance of interest in explicit modeling and model-based analysis of neural spike trains (e.g., Brown et al., 2004; Paninski et al., 2007; Stevenson et al., 2008), has largely disregarded the Hawkes-type models, focusing instead on their non-linear generalizations: the generalized linear models (GLMs), and related multiplicative models (Cardanobile and Rotter, 2010). GLMs are clearly powerful and flexible models of spiking processes, and are also related to the popular Linear–Non-linear encoding models (Chichilnisky, 2001; Paninski et al., 2004; Shoham et al., 2005). However, they do not enjoy the same level of mathematical simplicity as their Hawkes counterparts – only approximate analytical expressions for the correlation and the spectral properties of a GLM model were derived (Nykamp, 2007; Toyoizumi et al., 2009) under fairly restrictive conditions, while exact parameters for detailed, heterogeneous GLM models can only be evaluated numerically (Pillow et al., 2008).

The significance and applications of spike train models with closed-form expressions for the output correlation/spectral structure have begun to emerge in a number of recent studies. These include: (1) the ability to generate synthetic spike trains with a given auto- and cross-correlation structure (Brette, 2009; Krumin and Shoham, 2009; Macke et al., 2009; Gutnisky and Josic, 2010); (2) the ability to identify neural input-output encoding models “blindly” by analyzing the spectral and correlation distortions they induce (Krumin et al., 2010); (3) the ability to fit compact multivariate auto-regressive (MVAR) models to multi-channel neural spike trains (Krumin and Shoham, 2010); and (4) the ability to apply the associated powerful framework of Granger causality analysis (Granger, 1969; Krumin and Shoham, 2010). These early studies relied on the analysis of tractable non-linear spiking models such as threshold models (Macke et al., 2009; Gutnisky and Josic, 2010; Tchumatchenko et al., 2010) or the Linear–Non-linear-Poisson (LNP) models (Krumin and Shoham, 2009) driven by Gaussian input processes.

In this paper we revisit the Hawkes model within this new emerging framework for correlation-based, closed-form identification and analysis of spike trains models. The framework is thereby extended from the exclusive treatment of feed-forward models to treating more general and neuro-realistic (yet analytically tractable) models that also include feedback terms. In Section “Methods” we begin by reviewing some basic results for the correlation structure of the classical, homogenous (constant input) single and multivariate Hawkes model, derive new integral equations for the correlation structure of a Hawkes model driven by a *time-varying* (*inhomogeneous*) stationary random non-negative process input (see Figure Figure1),1), and propose a numerical method for solving them. In Section “Results,” we present the results of applying these methods to real neural recordings from isolated mouse retina, and the required methodological adaptations. We conclude with a discussion in Section “Discussion.”

In this section we begin by defining the Hawkes model, recalling its auto-correlation structure and then generalizing to multivariate (mutually exciting) non-homogeneous Hawkes model of point processes. Next, we propose a method for the solution of the resulting equations, and for the estimation of the different parameters of the model. In the final subsection the experimental methods of stimulation and data acquisition are presented.

Let us consider the intensity of a self-exciting point process to be defined by the following expression:

$$\mu (t)=\lambda +{\displaystyle \sum _{k}g\left(t-{t}_{k}\right)}$$

(1)

Here, the instantaneous firing intensity μ(*t*) is the exogenous input λ summed together with multiple shifted replicas of the self-excitation kernel *g*(*t*). The kernels are causal (*g*(*t*)=0, *t*<0), and *t _{k}* represents all the past spike-times. For technical reasons we will write the expression using the Stieltjes integral:

$$\mu (t)=\lambda +{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)dN(u)}$$

(2)

where *N*(*t*) is the counting process (number of spikes up to time *t*). The sum term in Eq. 1 is now replaced by a convolution of the spiking history with a linear kernel. The mean firing rate (denoted throughout the paper by *dN*) of this point process is given by:

$$\begin{array}{l}\langle dN\rangle \triangleq \mathbb{E}\left\{\frac{dN(t)}{dt}\right\}=\mathbb{E}\left\{\lambda +{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)dN(u)}\right\}=\\ =\lambda +{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)\mathbb{E}\left\{\frac{dN(u)}{du}\right\}du}=\lambda +\langle dN\rangle \cdot {\displaystyle \underset{0}{\overset{\infty}{\int}}g(u)du}\end{array}$$

(3)

Resulting in:

$$\langle dN\rangle =\frac{\lambda}{1-{\displaystyle \underset{0}{\overset{\infty}{\int}}g(u)du}}$$

(4)

The stability (and stationarity) condition for this model $({\int}_{0}^{\infty}g(u)du<1)$ can easily be inferred from this equation. An expression for the auto-covariance function of such a point process was derived in Hawkes (1971a), and we will briefly review here the main results (adapted from his auto-covariance notation into auto-correlation function notation used here for simplicity). We will distinguish between two different auto-correlation functions, the first:

$${\tilde{R}}_{dN}(\tau )\triangleq \frac{\mathbb{E}\left\{dN(t+\tau )dN(t)\right\}}{d{t}^{2}},$$

(5)

which has a delta function singularity *dN*·δ(τ) at τ=0 due to the nature of point processes, and the second:

$${R}_{dN}(\tau )\triangleq {\tilde{R}}_{dN}(\tau )-\langle dN\rangle \cdot \delta (\tau ),$$

(6)

from which this singularity was subtracted.

Using these definitions we get the following integral equation for the auto correlation of the output point process of the Hawkes model:

$${R}_{dN}(\tau )=\lambda \cdot \langle dN\rangle +g(\tau )\cdot \langle dN\rangle +{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}g(\tau -u){R}_{dN}(u)du}$$

(7)

This equation can be solved numerically (Mayers, 1962) or by using Wiener–Hopf related techniques (Noble, 1958; Hawkes, 1971b).

Similarly, Hawkes (1971a) generalized this solution (Eqs 4 and 7) to *multivariate* mutually exciting point processes by using matrix notation. The intensity of mutually exciting process becomes:

$$\mu (t)=\mathbf{\lambda}+{\displaystyle \underset{-\infty}{\overset{t}{\int}}\mathbf{\text{G}}(t-u)dN(u)}$$

(8)

with mean firing rates:

$$\langle d\mathbf{\text{N}}\rangle ={\left(I-{\displaystyle \underset{0}{\overset{\infty}{\int}}\mathbf{\text{G}}(u)du}\right)}^{-1}\cdot \mathbf{\lambda}$$

(9)

and the cross-correlation matrix as a solution of:

$${\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )=\mathbf{\lambda}\cdot {\langle d\mathbf{\text{N}}\rangle}^{T}+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}\mathbf{\text{G}}(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}$$

(10)

Let us now consider a more general case of a non-homogeneous Hawkes model, where the exogenous input λ(*t*) can be a time-varying (stationary) process:

$$\mu (t)=\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)dN(u)}$$

(11)

For example, this class of models includes the important special case (Figure (Figure1)1) where λ(*t*) is itself a non-negative stationary random process generated by a Linear–Non-linear cascade acting on a Gaussian process input (possibly a stimulus). Note the difference between the proposed linear–non-linear-Hawkes (LNH) model and the GLM-type models, in which the feedback term is summed with the *x* (*t*) and not with the **λ**(*t*) (according to the notation in Figure Figure1).1). This effectively changes the locus of the non-linearity present in the model and affects the model's properties and analytical tractability.

The mean firing rate of this point process can, in general, be found in a similar way as in Eqs 3 and 4:

$$\langle dN\rangle =\frac{\mathbb{E}\left\{\lambda (t)\right\}}{1-{\displaystyle \underset{0}{\overset{\infty}{\int}}g(u)du}}$$

(12)

Next, the auto-correlation function *R _{dN}*(τ) of this process can be derived using a similar procedure to the derivation of Eq. 7 (the detailed derivation can be found in Section “Correlation Structure of the LNH Model” of Appendix). This time, the auto-correlation function is governed by two coupled integral equations:

$$\begin{array}{l}{R}_{dN}(\tau )={R}_{\lambda dN}(\tau )+g(\tau )\cdot \langle dN\rangle +{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}g(\tau -u){R}_{dN}(u)du}\\ {R}_{\lambda dN}(\tau )={R}_{\lambda}(\tau )+{\displaystyle \underset{\tau}{\overset{\infty}{\int}}g(u-\tau ){R}_{\lambda dN}(u)du}\end{array}$$

(13)

These two equations provide the solution for the output auto-correlation function *R _{dN}*(τ) and for the cross-correlation ${\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )\triangleq \mathbb{E}\{\lambda (t+\tau )(d\mathbf{\text{N}}(t)/dt)\}$ between the exogenous input λ(

Equations 12 and 13 can be further generalized to a multivariate case (mutually exciting point processes), and be written using the matrix notation:

$$\begin{array}{l}\langle d\mathbf{\text{N}}\rangle ={\left(I-{\displaystyle \underset{0}{\overset{\infty}{\int}}\mathbf{\text{G}}(u)du}\right)}^{-1}\cdot \mathbb{E}\left\{\mathbf{\lambda}(t)\right\}\\ {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}\mathbf{\text{G}}(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}\\ {\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda}(\tau )+{\displaystyle \underset{\tau}{\overset{\infty}{\int}}{\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(u){\mathbf{\text{G}}}^{T}(u-\tau )du}\end{array}$$

(14)

Note that for constant **λ** these equations are reduced to Eqs 9 and 10.

The equations for the correlation structure of a single self-exciting point process and multivariate mutually exciting point processes (Eqs 13 and 14 respectively) can be solved numerically by switching from continuous time integral notation to discrete time matrix notation, and consequently performing matrix calculations. The integration operations in the Eqs 13 and 14 are thus converted to matrix multiplication operations. This allows a simple and straightforward way to solve the equations for the output correlation structure. Here, we only briefly present the main results. All the detailed explanations on the notation used, on how the appropriate matrices and vectors are built, and how the equations are solved in both single- and multi-channel cases can be found in Section “Solution of the Integral Equations” of Appendix. Using the new notation the output correlation is estimated by:

$$\begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}={\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}\right)}^{-1}\cdot \left({\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}+\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)\right),\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\right)}^{-1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}\end{array}$$

(15)

where ${\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}},\text{\hspace{0.17em}}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda},\text{\hspace{0.17em}}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T},{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\text{\hspace{0.17em}and\hspace{0.17em}}\underset{\xaf}{\mathbf{\text{G}}}$ and G are block column vectors that represent the sampled versions of the correlations **R**_{dN}(τ), **R**_{λ}(τ), **R**_{λdN}(τ), and the feedback kernel **G**(τ). Block matrices ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}$ and ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}$ are built from **G**(τ), and $\underset{\xaf}{\underset{\xaf}{I}}$ is the unity matrix of appropriate dimensions (see also Solution of the Integral Equations of Appendix). The generalized Hawkes model has three different sets of parameters – the input correlation structure **R**_{λ}(τ), the output correlation structure **R*** _{dN}*(τ), and the Hawkes feedback kernel

$$\begin{array}{ll}\left(\text{I}\right)\hfill & {\text{R}}_{d\mathbf{\text{N}}}(\tau ),\mathbf{\text{G}}(\tau )\Rightarrow \hfill \\ \left(\text{II}\right)\hfill & {\text{R}}_{d\mathbf{\text{N}}}(\tau ),{\text{R}}_{\lambda}(\tau )\text{\hspace{0.17em}}\Rightarrow \hfill \\ \left(\text{III}\right)\hfill & {\text{R}}_{d\mathbf{\text{N}}}(\tau )\Rightarrow \hfill \end{array}$$

In the first scenario we are interested in the estimation of the input correlation structure, given the output correlation structure **R*** _{dN}*(τ) and the Hawkes kernel

$$\begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}=\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}\right)\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}-\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}=\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\right)\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\end{array}$$

(16)

After **R**_{λ}(τ)is estimated one can proceed, if interested, with the estimation of an LN cascade model for this correlation structure by applying the correlation pre-distortion procedures developed and detailed in (Krumin and Shoham, 2009) and (Krumin and Shoham, 2010). Estimation of the Linear–Non-linear cascade model, in addition to the connectivity kernels **G**(τ), can provide additional insights about the stimulus-driven neural activity.

The second possible scenario is to estimate the Hawkes kernels when the output and the input correlation structures are known (see, e.g., Figure Figure3B).3B). Here, once again, we can use the advantage of the same matrix notation (block column vector ${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}$ and block matrix ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{R}}}}}_{d\mathbf{\text{N}}}$ represent the **R**_{λdN}(τ) and **R*** _{dN}*(τ) correlation functions, respectively) and solve the following equations in an iterative manner to estimate

$$\begin{array}{l}{\underset{\xaf}{\mathbf{\text{G}}}}^{T}={\left(\underset{\xaf}{\underset{\xaf}{\langle d\mathbf{\text{N}}\rangle}}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{R}}}}}_{d\mathbf{\text{N}}}^{T}\right)}^{-1}\left({\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}^{T}-{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\right)\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\right)}^{-1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}\end{array}$$

(17)

where $\underset{\xaf}{\underset{\xaf}{\langle d\mathbf{\text{N}}\rangle}}$ stands for the block diagonal matrix with *diag* (*d N*) as its block elements on the main diagonal.

The iterative solution of this set of equations is explained in detail in Section “Solution of the Integral Equations” of Appendix, Eq. A23.

The third possible scenario is to estimate both the kernels **G**(τ) and the input correlation structure **R**_{λ}(τ), given only the output correlation structure **R*** _{dN}*(τ). In general, this problem is not well-posed and does not have a unique solution, and additional application-driven constraints on the structure of

In general, the connectivity between the different units (**G**(τ) feedback terms in the Hawkes model) is not limited to non-negative values. Hence, the firing intensity μ(*t*) defined in Eqs 1 or 8 can occasionally become negative. However, the analytical derivations for the output mean rate and correlation structure are based on the assumption that μ(*t*) is non-negative for all *t*. The violation of this assumption results in a discrepancy between the actual and the analytical results. Simulation of the estimated LNH model [while using the effective firing intensity μ_{eff}(*t*)= max{μ(*t*), 0}] yields output spike trains with a correlation structure ${\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{\text{sim}}(\tau )$ that is different from the desired output correlation structure *R _{dN}*(τ) (used for the estimation of the model parameters). To address this issue an additional procedure was developed for the estimation of the actual feedback kernel

- Simulate a Hawkes point process using the original input correlation
*R*_{λ}(τ) and the estimated kernel*g*(τ). Use μ_{eff}(*t*)= max{μ(*t*), 0}. - Estimate the output correlation ${R}_{dN}^{\text{sim}}(\tau )$ of the simulated spike train. The violation of the μ(
*t*)≥ 0 assumption will result in a difference between the desired (*R*(τ)) and the estimated (${R}_{dN}^{\text{sim}}(\tau )$) output correlation structures._{dN} - Use the estimated ${R}_{dN}^{\text{sim}}(\tau )$ instead of the input correlation
*R*_{λ}(τ) in the Eq. 13 to estimate the kernel δ*g*(τ). The output correlation that should be used is the desired*R*(τ) throughout the iterative solution, only the input correlation_{dN}*R*_{λ}(τ) changes from iteration to iteration. - Update
*g*(τ)←*g*(τ)+ α·δ*g*(τ). The scalar α≤ 1 is used for controlling the speed and/or smoothness of the convergence. In Section “Application to Neural Spike Trains – Single Cells” we have used a relatively small α=0.1 to ensure smooth convergence to the solution. - Loop through steps 2–5 until the actual ${R}_{dN}^{\text{sim}}(\tau )$ of the simulated spike train converges to the desired
*R*(τ)._{dN}

The above procedure uses the difference between the model-based (simulated) and the desired (data-estimated) correlation structures of the output spike trains to systematically update the feedback kernel *g*(τ) until the difference between these two correlation structures becomes small enough. The resulting model allows to relax the assumption of μ(*t*)≥ 0 and to use μ_{eff}(*t*)= max{μ(*t*), 0} instead.

Animal experiments and procedures were approved by the Institutional Animal Care Committee at the Technion – Israel Institute of Technology and were in accordance with the NIH Guide for the Care and Use of Laboratory Animals. Six-week-old wild type mice (C57/BL) were euthanized using CO_{2} and then decapitated. Eyes were enucleated and immersed in Ringer's solution containing (in mM): NaCl, 124; KCl, 2.5; CaCl_{2}, 2; MgCl_{2}, 2; NaHCO_{3}, 26; NaH_{2}PO_{4}, 1.25; and Glucose, 22 (pH 7.35–7.4 with 95% O_{2} and 5% CO_{2} at RT). An incision was made at the ora serrata using a scalpel and the anterior chamber of the eye was separated from the posterior chamber cutting along the ora serrata with fine scissors. The lens was removed and the retina was gently cleaned of the remaining vitreous. Retinal tissue was isolated from the retinal pigmented epithelium. Three radial cuts were made and the isolated retina was flattened with the retinal ganglion cells facing the multi electrode array (MEA). During the experiment the retina was continuously perfused with oxygenated Ringer's solution.

The retina was stimulated by wide-field intensity-modulated light flashes using a DLP-based projector. The stimulus intensities were normally distributed and updated at the rate of 60Hz. Resulting activity was recorded using 60-channel MEA with 10 μm diameter, planar electrodes spaced at 100μm. The data was acquired with custom written data acquisition software using Matlab 7.5.0 data acquisition toolbox.

We performed a number of simulation studies to validate the methods proposed for the solution of the integral Eqs 13 and 14.

In Figures Figures2A–C2A–C the forward model solution by the Eq. 15 is compared to the auto-correlation function estimated from single simulated point processes with different self-excitation kernels, *g*(τ), under two different conditions – constant input λ (pure Hawkes model), or time-varying input λ(*t*) with an exponentially shaped auto-correlation function (LNH model). In Figure Figure2D2D an example of a bivariate case is presented with a more complex correlation structure of the input **λ**(*t*) and a set of self- and mutually exciting kernels **G**(τ).

As can be seen in all of these examples, the analytically predicted correlation functions had a near-perfect match with the mean correlation functions of the simulated spike trains (correlation coefficient ≥0.99). Individual correlation functions calculated from 10-min traces were more noisy, thus the forward analytical prediction vs. simulation correlation coefficients for single traces were significantly lower: 0.83 ±0.06.

Figure Figure3A3A shows the result of applying the “scenario I” solution (Eq. 16) to spike trains generated by the model presented in Figure Figure2D;2D; the mean identified input correlations have an excellent match with the ones used for generating the data (correlation coefficients: 0.99 and 0.92 respectively for the auto- and cross-correlations).

Figure Figure3B3B shows the result of applying the “scenario II” solution (Eq. 17) to spike trains generated by the model presented in Figure Figure2D;2D; the mean identified kernels greatly match the ones used in generating the data (correlation coefficients >0.99 for all kernels.

Next, we applied the method on the data recorded from the retina (see Methods for the experimental protocol). We started by analyzing the spike trains using reverse-correlation techniques (Ringach and Shapley, 2004) based on a feed-forward Linear–Non-linear–Poisson (LNP) model. The LNP-based estimates of the linear filter, and the static non-linearity (Figure (Figure4A)4A) were further used for the calculation of the expected output auto-correlation function of the estimated LNP model. This LNP-based output auto-correlation function was found to be noticeably different from the actual auto-correlation function of the measured spike trains (Figure (Figure44B).

This LN cascade was then used for generating the input (λ(*t*) in Figure Figure1)1) to the Hawkes feedback stage of the LNH model. The auto-correlation function of λ(*t*) is exactly that of the LNP model's output estimated previously and found inconsistent with the real recordings. Now, the input auto-correlation function *R*_{λ}(τ) was used together with the measured output auto-correlation function *R _{dN}*(τ) to estimate the Hawkes feedback kernel

Finally, we validated that the improved fit of the LNH model to the data compared with the LNP model, does not result from a model overfitting due to the larger number of parameters in the LNH model. For each unit, we computed an LN-Hawkes for a different data set from the same unit (Gaussian distribution, different mean intensity). Next, we simulated an output spike train using a “hybrid” LNH model (“original” LN model +“new” feedback kernel *g*(τ)), and estimated its correlation function. This output correlation function was compared to the correlation function of the original data by calculating the correlation coefficient between the two functions ρ_{LNH}. This procedure was applied to the nine units in our data set where the mean firing rates were >2 Hz. In eight out of these nine units the hybrid LNH model provided considerably better fits to the output correlation function than the corresponding LNP model, providing in those cases an average improvement of Δρ= ρ_{LNH}− ρ_{LNP}=0.19 with Δρ/ρ_{LNP}=30%. Note that this procedure is over-conservative, since there is no guarantee that kernels calculated for different input stimulus ensembles will be the same or conversely, that neural models will generalize across different stimulus ensembles.

In this paper, we extended previous work on the correlation-based simulation, identification and analysis of multi-channel spike train models with a feed-forward Linear–Non-linear (LN) stage driven by Gaussian process inputs (Krumin and Shoham, 2009; Krumin et al., 2010), by allowing the non-negative process to drive a feedback stage in the form of a multi-channel Hawkes process. The move from doubly stochastic Poisson (Cox) models in our previous work to doubly stochastic Hawkes models employed here vastly expands the range of realizable correlation structures, thus relaxing the main limitation of the previous results, and allowing for a superior, excellent fit (ρ 0.98) of the auto-correlation structures of spike trains recorded from real visually driven retinal ganglion cells. At the same time, it preserves the analytical tractability and closed-form correspondence between model parameters and the second-order statistical properties of the output spike trains, and thus, essentially, all of the advantages and potential applications of the general model-based correlation framework, which was limited, thus far, to feed-forward models. These currently include the synthetic generation of spike trains with a pre-defined correlation structure (Brette, 2009; Krumin and Shoham, 2009; Macke et al., 2009; Gutnisky and Josic, 2010; Tchumatchenko et al., 2010), “blind” correlation-based identification of single-neuron encoding models (Krumin et al., 2010), the compact representation of multi-channel spike trains in terms of multivariate auto-regressive processes and the framework of causality (Granger) analysis (Nykamp, 2007; Krumin and Shoham, 2010). As noted above, the LNH model is related to the commonly used GLM model, with the LNH feedback kernels paralleling the GLM history terms. Both ways of altering the underlying feed-forward LNP model lead to more flexible models capable of fitting more complex correlation structures, but the preferred fitting procedures for the two models differ: the GLM model is typically fit using a maximum likelihood approach, but this does not suit the LNH model (due to possible zero firing rates), where a method of moments (like the one introduced here) is more appropriate for the estimation of the linear kernels. A systematic study on the differences between the statistical properties of the two approaches falls beyond the scope of the current manuscript.

The model and analysis presented here also provide a new context and results to a significant body of related previous work on the second-order statistics of Hawkes models, which we will now review very briefly. The basic properties of the output correlation structure and the spectrum of a univariate self-exciting and a multivariate mutually exciting linear point process model without an exogenous drive were derived in the original works of Hawkes (1971a,b) using the linear representation of this process (Eq. 2). Brillinger (1975) also analyzes linear point process models and uses spectral estimators for the kernels, which he applies to the analysis of synaptic connections (Brillinger et al., 1976). Bremaud and Massoulie (2002) and Daley and Vere-Jones (2003) (exercise 8.3.4) present expressions for the output spectrum of a univariate Hawkes model excited by an exogenous correlated point process derived using an alternative, cluster process representation of the Hawkes process:

$$d\tilde{\mathbf{\text{N}}}(\omega )=\frac{\Gamma \cdot \mathbb{E}\left\{\lambda (t)\right\}/\left(1-\Gamma \right)+\tilde{\lambda}(\omega )}{{\left|1-\tilde{g}(\omega )\right|}^{2}},$$

where $\Gamma \triangleq {\int}_{0}^{\infty}g(u)du$ and $d\tilde{\mathbf{\text{N}}}(\omega ),\text{\hspace{0.17em}}\tilde{\lambda}(\omega ),\text{\hspace{0.17em}}\tilde{g}(\omega )$ represent the respective spectra of *dN*(*t*), λ(*t*), *g*(*t*). Our derivation in the Section “Methods” and “Correlation Structure of the LNH Model” of Appendix focused on expressions for the correlation structure of exogenously driven Hawkes process and was based on the linear representation, similar to Hawkes (1971a). Adding the exogenous input introduces a new term into the Hawkes integral Eq. 10, and a second integral equation for the cross-covariance term between the exogenous input and the output spike trains **R**_{λdN}(τ). The parameters of these generalized models, i.e., the kernels **G**(τ) and/or the input correlation structure **R _{λ}**(τ), can be directly estimated from the output process correlation structure using an iterative application of this set of equations, as illustrated in Section “Results,” or they could, alternatively, be estimated from the spectral expressions.

We next turn to discuss certain limitations of the proposed framework. First, the analytical equations for the auto-correlation structure of the point processes (Eqs 7, 10, 13, and 14) are *exactly* true under the assumption μ(*t*)≥ 0 (Eqs 2, 8, and 11) or when the stochastic intensity is always non-negative. These exact results could also provide an excellent agreement to many practical cases wherein the self-exciting Hawkes kernel *g*(τ) is only weakly negative (e.g., Figure Figure2),2), leading in such cases to slight systematic deviations at “negative” peaks. In cases of strong refractoriness or other inhibitory interactions, *g*(τ) becomes strongly negative, and the rectification of the stochastic intensity around zero leads to strong deviations from the assumptions underlying Eqs 7 and 13. For such cases we introduced an intuitive iterative procedure for computing *g*(τ) (see Refractoriness and Strong Inhibitory Connections), and it is likely that related alternatives are also possible. Although the convergence of this procedure is not proven, in practice, it was capable of estimating kernels for real neural spike trains that not only dramatically improved the auto-correlation fits relative to LNP cascades, but also generalized across different stimulus ensembles (a very conservative cross-validation test). Second, we have not addressed the important but complex issue of uniqueness of the different identification problems encountered here. Interestingly, in the examples we have examined, an excellent match was found, in practice, between the Hawkes kernels and their estimates (Figure (Figure3B),3B), although we are not aware of any guarantees of uniqueness here (these may perhaps be related to the nature of point processes). In the more general problem where both $\widehat{\mathbf{\text{G}}}(\tau ),\text{\hspace{0.17em}}{\widehat{\mathbf{\text{R}}}}_{\lambda}(\tau )$ are simultaneously estimated, it seems obvious that unique solutions can only be obtained by imposing additional constraints on the solutions (i.e., degree of smoothness and/or sparseness). In section “Application to Neural Spike Trains – Single Cells” we presented an example of the “scenario III”-type problem, where only the output correlation structure is actually observable. In this example we used additional application-driven constraints on the input correlation structure **R _{λ}**(τ) to infer the feedback kernels

When considering the broader relevance of this work, and the directions to which it may develop in the future, it is worth noting that some of the most fundamental and widely applied tools for the identification of systems rely on the use of second-order statistical properties (Ljung, 1999) (correlation or spectral). The increasing arsenal of tools for identifying spike train models from their correlations, rather than from their full observed realizations could form a welcome bridge between “classical” signal processing ideas and tools and the field of neural spike train analysis.

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.

This work was supported by Israeli Science Foundation grant #1248/06 and European Research Council starting grant #211055. We thank A. Tankus and the two anonymous reviewers for their comments on the manuscript.

We consider the Hawkes point process driven by a time-varying exogenous input, with the intensity defined in Eq. 11:

$$\mu (t)=\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)d\mathbf{\text{N}}(u)}$$

For the mean firing rate we receive:

$$\begin{array}{l}\langle d\mathbf{\text{N}}\rangle \triangleq \mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}\right\}=\mathbb{E}\left\{\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)d\mathbf{\text{N}}(u)}\right\}\\ =\mathbb{E}\left\{\lambda (t)\right\}+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)\mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}\right\}du}=\mathbb{E}\left\{\lambda (t)\right\}+\langle d\mathbf{\text{N}}\rangle \cdot {\displaystyle \underset{0}{\overset{\infty}{\int}}g(u)du}\end{array}$$

(A1)

resulting in:

$$\langle d\mathbf{\text{N}}\rangle =\frac{\mathbb{E}\left\{\lambda (t)\right\}}{1-{\displaystyle \underset{0}{\overset{\infty}{\int}}g(u)du}}$$

(A2)

Next, we expand the expressions for the correlation structure of the output spike trains, following a similar formalism to the derivation found in Hawkes (1971a) for the correlations of homogeneous Hawkes processes:

$$\begin{array}{l}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )\triangleq {\tilde{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}(\tau )-\langle d\mathbf{\text{N}}\rangle \cdot \delta (\tau )\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}=\mathbb{E}\left\{d\mathbf{\text{N}}(t+\tau )d\mathbf{\text{N}}(t)\right\}/d{t}^{2}-\langle d\mathbf{\text{N}}\rangle \cdot \delta (\tau )\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}=\mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}\cdot \left[\lambda (t+\tau )+{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}g(t+\tau -u)d\mathbf{\text{N}}(u)}\right]\right\}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}=\mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}\cdot \lambda (t+\tau )\right\}+\mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}g(t+\tau -u)d\mathbf{\text{N}}(u)}\right\}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}g(t+\tau -u){\tilde{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}(t-u)du}\end{array}$$

(A3)

Now, substituting ${\tilde{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )+\langle d\mathbf{\text{N}}\rangle \cdot \delta (\tau )$ we get:

$${\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+g(\tau )\cdot \langle d\mathbf{\text{N}}\rangle +{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}g(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}$$

(A4)

We have arrived to a solution similar to Eq. 7 with one additional term *R*_{λdN}(τ) that will be derived in Part II.

The derivation of *R*_{λdN}(τ) has much in common with the derivations in Part I above.

$$\begin{array}{l}{\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )\triangleq \mathbb{E}\left\{\lambda (t+\tau )\cdot \frac{d\mathbf{\text{N}}(t)}{dt}\right\}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}}=\mathbb{E}\left\{\lambda (t+\tau )\cdot \left[\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)d\mathbf{\text{N}}(u)}\right]\right\}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}}=\mathbb{E}\left\{\lambda (t+\tau )\cdot \lambda (t)\right\}+\mathbb{E}\left\{\lambda (t+\tau )\cdot {\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)d\mathbf{\text{N}}(u)}\right\}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}}={\mathbf{\text{R}}}_{\lambda}(\tau )+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u)\mathbb{E}\left\{\lambda (t+\tau )\frac{d\mathbf{\text{N}}(u)}{du}\right\}du}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}}={\mathbf{\text{R}}}_{\lambda}(\tau )+{\displaystyle \underset{-\infty}{\overset{t}{\int}}g(t-u){\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(t+\tau -u)du}\end{array}$$

(A5)

To summarize, the derivations in Part I and Part II of the current Appendix result in two coupled integral equations:

$$\begin{array}{l}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+g(\tau )\cdot \langle d\mathbf{\text{N}}\rangle +{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}g(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}\\ {\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda}(\tau )+{\displaystyle \underset{\tau}{\overset{\infty}{\int}}g(u-\tau ){\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(u)du}\end{array}$$

(A6)

Let us now consider a multivariate inhomogeneous Hawkes process:

$$\mu (t)=\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}\mathbf{\text{G}}(t-u)d\mathbf{\text{N}}(u)},$$

(A7)

where **μ**(*t*), **λ**(*t*), and *d***N**(*t*)are now column vectors, and **G**(τ) is a square matrix. The values in the row *#r* and column *#s* of the matrix **G**(τ) correspond to the mutual-excitation kernel that explains the effect of the firing history of the process *#s* on the stochastic intensity of the process *#r*.

The expression for the mean firing rate *d***N** of the process is derived in the following way:

$$\begin{array}{l}\langle d\mathbf{\text{N}}\rangle \triangleq \mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}\right\}=\mathbb{E}\left\{\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}\mathbf{\text{G}}(t-u)d\mathbf{\text{N}}(u)}\right\}\\ =\mathbb{E}\left\{\lambda (t)\right\}+{\displaystyle \underset{-\infty}{\overset{t}{\int}}\mathbf{\text{G}}(t-u)\mathbb{E}\left\{\frac{d\mathbf{\text{N}}(t)}{dt}\right\}du}=\mathbb{E}\left\{\lambda (t)\right\}+\langle d\mathbf{\text{N}}\rangle \cdot {\displaystyle \underset{0}{\overset{\infty}{\int}}\mathbf{\text{G}}(u)du},\end{array}$$

(A8)

resulting in

$$\langle d\mathbf{\text{N}}\rangle ={\left(I-{\displaystyle \underset{0}{\overset{\infty}{\int}}\mathbf{\text{G}}(u)du}\right)}^{-1}\cdot \mathbb{E}\left\{\lambda (t)\right\}$$

(A9)

The output correlation structure is now defined by:

$$\begin{array}{l}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )\triangleq {\tilde{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}(\tau )-diag\left(\langle d\mathbf{\text{N}}\rangle \right)\cdot \delta (\tau )\\ =\mathbb{E}\left\{d\mathbf{\text{N}}(t+\tau )d{\mathbf{\text{N}}}^{T}(t)\right\}/d{t}^{2}-diag\left(\langle d\mathbf{\text{N}}\rangle \right)\cdot \delta (\tau )\\ =\mathbb{E}\left\{\left[\lambda (t+\tau )+{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}\mathbf{\text{G}}(t+\tau -u)d\mathbf{\text{N}}(u)}\right]\cdot \frac{d{\mathbf{\text{N}}}^{T}(t)}{dt}\right\}\\ =\mathbb{E}\left\{\lambda (t+\tau )\cdot \frac{d{\mathbf{\text{N}}}^{T}(t)}{dt}\right\}+\mathbb{E}\left\{{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}\mathbf{\text{G}}(t+\tau -u)d\mathbf{\text{N}}(u)}\cdot \frac{d{\mathbf{\text{N}}}^{T}(t)}{dt}\right\}\\ ={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )+{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}\mathbf{\text{G}}(t+\tau -u){\tilde{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}(t-u)du}\\ ={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{t+\tau}{\int}}\mathbf{\text{G}}(t+\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(t-u)du}\\ ={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}\mathbf{\text{G}}\left(\tau -u\right)\text{\hspace{0.05em}}}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du\end{array}$$

(A10)

Similarly to the Eq. A5 we can also derive:

$$\begin{array}{l}{\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )\triangleq \mathbb{E}\left\{\lambda (t+\tau )\cdot \frac{d{\mathbf{\text{N}}}^{T}(t)}{dt}\right\}\\ \text{}=\mathbb{E}\left\{\lambda (t+\tau )\cdot {\left[\lambda (t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}\mathbf{\text{G}}(t-u)d\mathbf{\text{N}}(u)}\right]}^{T}\right\}\\ \text{}=\mathbb{E}\left\{\lambda (t+\tau )\cdot {\lambda}^{T}(t)\right\}+\mathbb{E}\left\{\lambda (t+\tau )\cdot {\displaystyle \underset{-\infty}{\overset{t}{\int}}d{\mathbf{\text{N}}}^{T}(u){\mathbf{\text{G}}}^{T}(t-u)}\right\}\\ \text{}={\mathbf{\text{R}}}_{\lambda}(t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}\mathbb{E}\left\{\lambda (t+\tau )\frac{d{\mathbf{\text{N}}}^{T}(u)}{du}\right\}{\mathbf{\text{G}}}^{T}(t-u)du}\\ \text{}={\mathbf{\text{R}}}_{\lambda}(t)+{\displaystyle \underset{-\infty}{\overset{t}{\int}}{\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(t+\tau -u){\mathbf{\text{G}}}^{T}(t-u)du}\\ \text{}={\mathbf{\text{R}}}_{\lambda}(\tau )+{\displaystyle \underset{\tau}{\overset{\infty}{\int}}{\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(u){\mathbf{\text{G}}}^{T}(u-\tau )du}\end{array}$$

(A11)

The following coupled equations govern the relationship between the input correlation structure **R _{λ}**(τ), the output correlation structure

$$\begin{array}{l}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}\mathbf{\text{G}}(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}\\ {\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda}(\tau )+{\displaystyle \underset{\tau}{\overset{\infty}{\int}}{\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(u){\mathbf{\text{G}}}^{T}(u-\tau )du}\end{array}$$

(A12)

We can rewrite these equations in the following manner:

$$\begin{array}{l}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}\mathbf{\text{G}}(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}={\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}+{\displaystyle \underset{0}{\overset{\tau}{\int}}\mathbf{\text{G}}(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}+{\displaystyle \underset{0}{\overset{\infty}{\int}}\mathbf{\text{G}}(\tau +u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}\\ {\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}^{T}(\tau )={\mathbf{\text{R}}}_{\lambda}^{T}(\tau )+{\displaystyle \underset{\tau}{\overset{\infty}{\int}}\mathbf{\text{G}}(u-\tau ){\mathbf{\text{R}}}_{\lambda d\mathbf{\text{N}}}^{T}(u)du}\end{array}$$

(A13)

To solve these equations numerically we use the following discretized representation:

$$\begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}+\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{T}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{H}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}+\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T},\end{array}$$

(A14)

where *d N* – is a block column vector representing the mean firing rates of the output spike trains

$\underset{\xaf}{\mathbf{\text{G}}},\text{\hspace{0.17em}}{\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}},\text{\hspace{0.17em}}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda},\text{\hspace{0.17em}}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}$ – block column vectors of *N* block elements with the first block element representing τ =0, and the last block element representing τ =τ_{max}. The choice of the discretization time-step *d*τ depends on the desired time resolution of the solution.

${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T},\text{\hspace{0.17em}}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}$ – also block column vectors, but with their block elements transposed (in the univariate case ${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda ,\lambda d\mathbf{\text{N}}}^{T}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda ,\lambda d\mathbf{\text{N}}}$)

${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1},\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{T},\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{H},\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}$ – square block matrices of size *N* ×*N* blocks that match the dimensions of the block column vectors.

To convert the integration operations into matrix multiplication operations we define the matrices ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}$ and ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}=\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{T}+\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{H}$ (*d*τ – time resolution) in the following way:

$${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}=d\tau \cdot \left[\begin{array}{ccccc}{\mathbf{\text{G}}}_{0}& {\mathbf{\text{G}}}_{1}& {\mathbf{\text{G}}}_{2}& \cdots & {\mathbf{\text{G}}}_{\mathbf{\text{N}}-1}\\ 0& {\mathbf{\text{G}}}_{0}& {\mathbf{\text{G}}}_{1}& \vdots \\ \vdots & \ddots & \ddots & \ddots \\ {\mathbf{\text{G}}}_{1}\\ 0& \cdots & 0& {\mathbf{\text{G}}}_{0}\end{array}\right]$$

(A15)

is a block Toeplitz matrix with the elements of the block vector G in the first row, and zeros in the first block column (excluding the main diagonal). The block elements of the matrix are;

$${\mathbf{\text{G}}}_{k}\triangleq \mathbf{\text{G}}\left(\tau =k\cdot d\tau \right)$$

(A16)

${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}$ is a sum of two other matrices: ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}=\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{T}+\text{\hspace{0.17em}}{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{H}$, where

$${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{T}=d\tau \cdot \left[\begin{array}{cccc}{\mathbf{\text{G}}}_{0}& 0& \cdots & 0\\ {\mathbf{\text{G}}}_{1}& {\mathbf{\text{G}}}_{0}& \vdots \\ {\mathbf{\text{G}}}_{2}& \ddots & \ddots \\ \vdots & \ddots & 0\\ {\mathbf{\text{G}}}_{\mathbf{\text{N}}-1}& \cdots & {\mathbf{\text{G}}}_{1}& {\mathbf{\text{G}}}_{0}\end{array}\right]$$

(A17)

is a block Toeplitz matrix with the elements of the block vector $\underset{\xaf}{\mathbf{\text{G}}}$ in the first block column, and zeros in the first block row (excluding the main diagonal).

$${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{H}=d\tau \cdot \left[\begin{array}{ccccc}{\mathbf{\text{G}}}_{0}& {\mathbf{\text{G}}}_{1}& {\mathbf{\text{G}}}_{2}& \cdots & {\mathbf{\text{G}}}_{\mathbf{\text{N}}-1}\\ {\mathbf{\text{G}}}_{1}& {\mathbf{\text{G}}}_{2}& \u22f0& 0\\ {\mathbf{\text{G}}}_{2}& \u22f0& \vdots \\ \vdots & \u22f0\\ {\mathbf{\text{G}}}_{\mathbf{\text{N}}-1}& 0& \cdots & 0\end{array}\right]$$

(A18)

is a block Henkel matrix with the elements of the block vector $\underset{\xaf}{\mathbf{\text{G}}}$ in the first block column, and zeros in the last block row (excluding the secondary diagonal).

The solution of the Eq. A14 for the output correlation structure ${\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}$ (the forward model) is straightforward:

$$\begin{array}{l}\begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}+\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\end{array}\}\Rightarrow \\ \begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}={(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2})}^{-1}\cdot \left({\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}+\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)\right)\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1})}^{-1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}\end{array}\},\end{array}$$

(A19)

where the second equation is solved in the beginning and then substituted into the first (after the appropriate rearrangement of ${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}$).

For scenario $(I)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau ),\hspace{0.17em}\mathbf{\text{G}}(\tau )\Rightarrow {\stackrel{\wedge}{\mathbf{\text{R}}}}_{\lambda}(\tau )$ the solution is also straightforward:

$$\begin{array}{l}\begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}+\underset{\xaf}{\mathbf{\text{G}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\end{array}\}\Rightarrow \\ \begin{array}{l}{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}=\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{2}\right)\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}-\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}=\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\right)\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\end{array}\}\end{array}$$

(A20)

For scenario $(II)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau ),{\mathbf{\text{R}}}_{\lambda}(\tau )\text{\hspace{0.17em}}\Rightarrow \stackrel{\wedge}{\mathbf{\text{G}}}(\tau )$ we will reorganize the equations and the matrix notation. Let us rewrite the first equation of Eq. A12 in the following way:

$$\begin{array}{l}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{-\infty}{\overset{\tau}{\int}}\mathbf{\text{G}}(\tau -u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u)du}\\ \text{}={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(\tau )+\mathbf{\text{G}}(\tau )\cdot diag\left(\langle d\mathbf{\text{N}}\rangle \right)+{\displaystyle \underset{0}{\overset{\infty}{\int}}\mathbf{\text{G}}(u){\mathbf{\text{R}}}_{d\mathbf{\text{N}}}(u-\tau )du}\\ {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(\tau )={\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(\tau )+diag\left(\langle d\mathbf{\text{N}}\rangle \right)\cdot {\mathbf{\text{G}}}^{T}(\tau )+{\displaystyle \underset{0}{\overset{\infty}{\int}}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(u-\tau ){\mathbf{\text{G}}}^{T}(u)du}\end{array}$$

(A21)

This equation, written in the matrix form is:

$${\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}^{T}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}+\underset{\xaf}{\underset{\xaf}{\langle d\mathbf{\text{N}}\rangle}}\cdot {\underset{\xaf}{\mathbf{\text{G}}}}^{T}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{R}}}}}_{d\mathbf{\text{N}}}^{T}\cdot {\underset{\xaf}{\mathbf{\text{G}}}}^{T},$$

where the matrix $\underset{\xaf}{\underset{\xaf}{\langle d\mathbf{\text{N}}\rangle}}$ is a block diagonal matrix with blocks of *diag* (*d N*) replicated

$${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{R}}}}}_{d\mathbf{\text{N}}}^{T}\triangleq d\tau \cdot \left[\begin{array}{cccc}{\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(0)& {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(1)& {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(2)& \cdots \\ {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(1)& {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(0)& \ddots & \ddots \\ {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(2)& \ddots & \ddots \\ \vdots & \ddots & {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(1)\\ {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(1)& {\mathbf{\text{R}}}_{d\mathbf{\text{N}}}^{T}(0)\end{array}\right]$$

(A22)

This, together with the matrix form of the second equation of Eq. A12 brings us to a couple of equations:

$$\begin{array}{l}{\underset{\xaf}{\mathbf{\text{G}}}}^{T}={\left(\underset{\xaf}{\underset{\xaf}{\langle d\mathbf{\text{N}}\rangle}}+{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{R}}}}}_{d\mathbf{\text{N}}}^{T}\right)}^{-1}\left({\underset{\xaf}{\mathbf{\text{R}}}}_{d\mathbf{\text{N}}}^{T}-{\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}\right)\text{}(\ast )\\ {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}^{T}={\left(\underset{\xaf}{\underset{\xaf}{I}}-{\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}\right)}^{-1}\cdot {\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda}^{T}\text{}(\ast \ast )\end{array}$$

(A23)

These can be solved iteratively:

- (i) Start with a random ${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}$
- (ii) Find $\underset{\xaf}{\mathbf{\text{G}}}$ from (*)
- (iii) Build matrix ${\underset{\xaf}{\underset{\xaf}{\mathbf{\text{G}}}}}_{1}$
- (iv) Find ${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}$ from (**)
- (v) Goto ii)

We can alternatively set the initial condition to ${\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda d\mathbf{\text{N}}}={\underset{\xaf}{\mathbf{\text{R}}}}_{\lambda},$ which corresponds to $\underset{\xaf}{\mathbf{\text{G}}}=\underset{\xaf}{0}.$

This iterative solution converges very rapidly and, in practice, a single iteration brings us very close to the final solution.

- Bremaud P., Massoulie L. (2002). Power spectra of general shot noises and Hawkes point processes with a random excitation. Adv. Appl. Probab. 34, 205–22210.1239/aap/1019160957 [Cross Ref]
- Brette R. (2009). Generation of correlated spike trains. Neural. Comput. 21, 188–21510.1162/neco.2009.12-07-657 [PubMed] [Cross Ref]
- Brillinger D. (1988). Maximum likelihood analysis of spike trains of interacting nerve cells. Biol. Cybern. 59, 189–20010.1007/BF00318010 [PubMed] [Cross Ref]
- Brillinger D. R. (1975). The identification of point process systems. Ann. Probab. 3, 909–92410.1214/aop/1176996218 [Cross Ref]
- Brillinger D. R., Bryant H. L., Segundo J. P. (1976). Identification of synaptic interactions. Biol. Cybern. 22, 213–22810.1007/BF00365087 [PubMed] [Cross Ref]
- Brown E. N., Kass R. E., Mitra P. P. (2004). Multiple neural spike train data analysis: state-of-the-art and future challenges. Nat. Neurosci. 7, 456–46110.1038/nn1228 [PubMed] [Cross Ref]
- Cardanobile S., Rotter S. (2010). Multiplicatively interacting point processes and applications to neural modeling. J. Comput. Neurosci. 28, 267–28410.1007/s10827-009-0204-0 [PubMed] [Cross Ref]
- Chichilnisky E. J. (2001). A simple white noise analysis of neuronal light responses. Network 12, 199–213 [PubMed]
- Chornoboy E., Schramm L., Karr A. (1988). Maximum likelihood identification of neural point process systems. Biol. Cybern. 59, 265–27510.1007/BF00332915 [PubMed] [Cross Ref]
- Daley D. J., Vere-Jones D. (2003). An Introduction to the Theory of Point Processes, Vol. 1 New York: Springer
- Granger C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424–43810.2307/1912791 [Cross Ref]
- Gutnisky D. A., Josic K. (2010). Generation of spatio-temporally correlated spike-trains and local-field potentials using a multivariate autoregressive process. J. Neurophysiol. 103, 2912–293010.1152/jn.00518.2009 [PubMed] [Cross Ref]
- Hawkes A. G. (1971a). Spectra of some self-exciting and mutually exciting point processes. Biometrika 58, 83–9010.1093/biomet/58.1.83 [Cross Ref]
- Hawkes A. G. (1971b). Point spectra of some mutually exciting point processes. J. R. Stat. Soc. Series B Methodol. 33, 438–443
- Johnson D. H. (1996). Point process models of single-neuron discharges. J. Comput. Neurosci. 3, 275–29910.1007/BF00161089 [PubMed] [Cross Ref]
- Kailath T., Sayed A. H., Hassibi B. (2000). Linear Estimation. Upper Saddle River, NJ: Prentice Hall
- Krumin M., Shimron A., Shoham S. (2010). Correlation-distortion based identification of Linear-Nonlinear-Poisson models. J. Comput. Neurosci. 29, 301–30810.1007/s10827-009-0184-0 [PubMed] [Cross Ref]
- Krumin M., Shoham S. (2009). Generation of spike trains with controlled auto- and cross-correlation functions. Neural. Comput. 21, 1642–166410.1162/neco.2009.08-08-847 [PubMed] [Cross Ref]
- Krumin M., Shoham S. (2010). Multivariate auto-regressive modeling and granger causality analysis of multiple spike trains. Computat. Intell. Neurosci. 2010, Article ID 752428. [PMC free article] [PubMed]
- Ljung L. (1999). System Identification – Theory for the User, 2nd Edn Upper Saddle River, NJ: Prentice Hall PTR
- Macke J. H., Berens P., Ecker A. S., Tolias A. S., Bethge M. (2009). Generating spike trains with specified correlation coefficients. Neural. Comput. 21, 397–42310.1162/neco.2008.02-08-713 [PubMed] [Cross Ref]
- Mayers D. F. (1962). “Part II. Integral equations, Chapters 11–14,” in Numerical Solution of Ordinary and Partial Differential Equations, ed. Fox L., editor. (London: Pergamon; ), 145–183
- Noble B. (1958). Methods Based on the Wiener-Hopf Technique. London: Pergamon
- Nykamp D. Q. (2007). A mathematical framework for inferring connectivity in probabilistic neuronal networks. Math. Biosci. 205, 204–25110.1016/j.mbs.2006.08.020 [PubMed] [Cross Ref]
- Paninski L., Pillow J., Lewi J. (2007). Statistical models for neural encoding, decoding, and optimal stimulus design. Prog. Brain Res. 165, 493–50710.1016/S0079-6123(06)65031-0 [PubMed] [Cross Ref]
- Paninski L., Shoham S., Fellows M. R., Hatsopoulos N. G., Donoghue J. P. (2004). Superlinear population encoding of dynamic hand trajectory in primary motor cortex. J. Neurosci. 24, 8551–856110.1523/JNEUROSCI.0919-04.2004 [PubMed] [Cross Ref]
- Pillow J. W., Shlens J., Paninski L., Sher A., Litke A. M., Chichilnisky E. J., Simoncelli E. P. (2008). Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature 454, 995–99910.1038/nature07140 [PMC free article] [PubMed] [Cross Ref]
- Ringach D., Shapley R. (2004). Reverse correlation in neurophysiology. Cogn. Sci. 28, 147–16610.1207/s15516709cog2802_2 [Cross Ref]
- Shoham S., Paninski L. M., Fellows M. R., Hatsopoulos N. G., Donoghue J. P., Normann R. A. (2005). Statistical encoding model for a primary motor cortical brain-machine interface. IEEE Trans. Biomed. Eng. 52, 1312–132210.1109/TBME.2005.847542 [PubMed] [Cross Ref]
- Stevenson I. H., Rebesco J. M., Miller L. E., Koerding K. P. (2008). Inferring functional connections between neurons. Curr. Opin. Neurobiol. 18, 582–58810.1016/j.conb.2008.11.005 [PMC free article] [PubMed] [Cross Ref]
- Tchumatchenko T., Malyshev A., Geisel T., Volgushev M., Wolf F. (2010). Correlations and synchrony in threshold neuron models. Phys. Rev. Lett. 104, 058102.10.1103/PhysRevLett.104.058102 [PubMed] [Cross Ref]
- Toyoizumi T., Rad K. R., Paninski L. (2009). Mean-field approximations for coupled populations of generalized linear model spiking neurons with Markov refractoriness. Neural. Comput. 21, 1203–124310.1162/neco.2008.04-08-757 [PubMed] [Cross Ref]

Articles from Frontiers in Computational Neuroscience are provided here courtesy of **Frontiers Media SA**