Frontiers in Computational Neuroscience Correlation-based Analysis and Generation of Multiple Spike Trains Using Hawkes Models with an Exogenous Input

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 (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. 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 IntroductIon Linear system models enjoy a fundamental role in the analysis of a wide range of natural and engineered signals and processes 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, where maximum likelihood (ML) parameter estimation procedures can be used to estimate the syn-aptic 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 The correlation structure of neural activity is believed to play a major role in the …

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 autoand 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 ; (3) the ability to fit compact multivariate auto-regressive (MVAR) models to multi-channel neural spike trains ; and (4) the ability to apply the associated powerful framework of Granger causality analysis (Granger, 1969;. 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

IntroductIon
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 autoregressive 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(Brillinger, , 1988Brillinger 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 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 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."

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

theoretIcal Background
Let us consider the intensity of a self-exciting point process to be defined by the following expression: Here, the instantaneous firing intensity μ(t) is the exogenous input λ summed together with multiple shifted replicas of the selfexcitation 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: 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: Resulting in: The stability (and stationarity) condition for this model 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: which has a delta function singularity 〈dN〉·δ(τ) at τ = 0 due to the nature of point processes, and the second: 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: These two equations provide the solution for the output auto-correlation function R dN (τ) and for the cross-correlation between the exogenous input λ(t) and the point process whose intensity is defined by Eq. 11. Here, the input auto-correlation function R λ (τ) and the self-exciting kernel g(τ) serve as given parameters (see also Identification of the LNH Model).
Equations 12 and 13 can be further generalized to a multivariate case (mutually exciting point processes), and be written using the matrix notation: Note that for constant λ these equations are reduced to Eqs 9 and 10.

IdentIfIcatIon of the Lnh modeL
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: where and are block column vectors that represent the sampled versions of the correlations R dN (τ), R λ (τ), R λdN (τ), and the feedback kernel G(τ). Block matrices G 1 and G 2 are built from G(τ), and 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 G(τ). Thus, in addition to the forward problem solution presented in Eq. 15, there are three other possible basic scenarios for the identification of the different parts of the proposed generalized Hawkes model from the correlation structure of the observed spike train(s).
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: with mean firing rates: and the cross-correlation matrix as a solution of: the LInear-non-LInear-hawkes modeL and Its correLatIons Let us now consider a more general case of a non-homogeneous Hawkes model, where the exogenous input λ(t) can be a timevarying (stationary) process: For example, this class of models includes the important special case (Figure 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 1). This effectively changes the locus of the nonlinearity 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: 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: Correlation-based Hawkes model analysis feedback kernel g(τ) from the input and the output correlations (R λ (τ) and R dN (λ), respectively). The procedure is summarized in the following algorithm: 1. Estimate initial g(τ) from R λ (τ) and R dN (λ) by solving Eq. 13 (in its matrix form of Eq. 17). 2. Simulate a Hawkes point process using the original input correlation R λ (τ) and the estimated kernel g(τ). Use μ eff (t) = max{μ(t), 0}. 3. Estimate the output correlation R dN sim ( ) τ of the simulated spike train. The violation of the μ(t) ≥ 0 assumption will result in a difference between the desired (R dN (τ)) and the estimated ( ( )) R dN sim τ output correlation structures. 4. Use the estimated R dN sim ( ) τ 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 dN (τ) throughout the iterative solution, only the input correlation R λ (τ) changes from iteration to iteration. 5. 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. 6. Loop through steps 2-5 until the actual R dN sim ( ) τ 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.

Retina preparation
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.

Electrophysiology
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 60 Hz. Resulting activity was recorded using 60-channel MEA with 10 μm diameter, 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 G(τ). By using the aforementioned matrix notation the solution can be achieved in a straightforward manner, akin to the forward problem: 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 . 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 3B). Here, once again, we can use the advantage of the same matrix notation (block column vector R N d and block matrix R N d represent the R λdN (τ) and R dN (τ) correlation functions, respectively) and solve the following equations in an iterative manner to estimate G(τ): where dN stands for the block diagonal matrix with diag (〈dN〉) 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 applicationdriven constraints on the structure of G(τ) and/or R λ (τ) should be considered. We will leave additional discussion on the uniqueness of the solution to the results (see Application to Neural Spike Trains -Single Cells) and in Sections "Discussion."

Refractoriness and strong inhibitory connections
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 R dN sim ( ) τ 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 In Figures 2A-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 2D 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(τ).
planar electrodes spaced at 100 μm. The data was acquired with custom written data acquisition software using Matlab 7.5.0 data acquisition toolbox.

results sIMulatIon studIes
We performed a number of simulation studies to validate the methods proposed for the solution of the integral Eqs 13 and 14. This LN cascade was then used for generating the input (λ(t) in Figure 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 g(τ) (Figure 4C) from Eq. 13 (including the procedure described in the Refractoriness and Strong Inhibitory Connections). Interestingly, the output auto-correlation function of the newly estimated LNH model (as measured from the simulated spike trains) was in excellent agreement with the auto-correlation function of the actual neural data (Figure 4D). The addition of the linear Hawkes feedback stage to the classical feed-forward LNP model proved beneficial to the model's capability of explaining more complex spike train correlation structures of real neural recordings ( Figure 4E).
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 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 3A shows the result of applying the "scenario I" solution (Eq. 16) to spike trains generated by the model presented in Figure 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 3B shows the result of applying the "scenario II" solution (Eq. 17) to spike trains generated by the model presented in Figure 2D; the mean identified kernels greatly match the ones used in generating the data (correlation coefficients >0.99 for all kernels.

applIcatIon to neural spIke traIns -sIngle cells
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 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 autocorrelation function of the measured spike trains (Figure 4B). 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 , the compact representation of multi-channel spike trains in terms of multivariate auto-regressive processes and the framework of causality (Granger) analysis (Nykamp, 2007;. 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 feedforward 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 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.

dIscussIon
In this paper, we extended previous work on the correlationbased 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;, 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 autocorrelation 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 references Bremaud, P., and Massoulie, L. (2002).
Power spectra of general shot noises and Hawkes point processes with a random excitation. Adv. Appl. Probab. 34, 205-222. Brette, R. (2009). Generation of correlated spike trains. Neural. Comput. 21, 188-215. Brillinger, D. (1988). Maximum likelihood analysis of spike trains of interacting nerve cells. Biol. Cybern. 59,[189][190][191][192][193][194][195][196][197][198][199][200] 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 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 ˆ,ˆ( ) G( ) R τ τ  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 G(τ). Interestingly, the exact same "scenario III"-type framework can be used for generating synthetic spike trains with a controlled correlation structure. This application will benefit from using the LNH feedback model by harnessing the capability of generating spike trains with a much richer ensemble of possible correlation structures in comparison with the feed-forward-only models like LNP. Additionally, once ˆ( ) R  τ is determined there is an additional level of non-uniqueness in the determination of the underlying LN structure, which can also be overcome by imposing constraints (e.g., a minimum phase constraint ).
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.

acknowledgMents
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. 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: ω 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 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 appendIx correlatIon structure of the lnh Model

Part I -Derivation of the output correlation of the inhomogeneous Hawkes point process
We consider the Hawkes point process driven by a time-varying exogenous input, with the intensity defined in Eq. 11: For the mean firing rate we receive: 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: We have arrived to a solution similar to Eq. 7 with one additional term R λdN (τ) that will be derived in Part II.

Part II -Derivation of the cross-correlation between the exogenous input λ (t) and the output point process
The derivation of R λdN (τ) has much in common with the derivations in Part I above.
To summarize, the derivations in Part I and Part II of the current Appendix result in two coupled integral equations:

Part III -Derivation of the output correlation structure for the multidimensional LNH model
Let us now consider a multivariate inhomogeneous Hawkes process: where m(t), λ(t), and dN(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 〈dN〉 of the process is derived in the following way: We can rewrite these equations in the following manner: To solve these equations numerically we use the following discretized representation: where 〈dN〉 -is a block column vector representing the mean firing rates of the output spike trains -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.
-also block column vectors, but with their block elements transposed (in the univariate case R R T H -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 G 1 and G G G 2 = + T H (dτ -time resolution) in the following way: 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; G 2 is a sum of two other matrices: Similarly to the Eq. A5 we can also derive:

Part I -Developing the discrete time matrix notation formalism for the integral equations
The following coupled equations govern the relationship between the input correlation structure R l (τ), the output correlation structure R dN (τ), and the feedback linear kernel G(τ) of the generalized Hawkes model:   These can be solved iteratively: (i) Start with a random R N d (ii) Find G from (*) (iii) Build matrix G 1 (iv) Find R N d from (**) (v) Goto ii) We can alternatively set the initial condition to R R N   d = , which corresponds to G = 0.
This iterative solution converges very rapidly and, in practice, a single iteration brings us very close to the final solution.
is a block Toeplitz matrix with the elements of the block vector G in the first block column, and zeros in the first block row ( excluding the main diagonal). is a block Henkel matrix with the elements of the block vector G in the first block column, and zeros in the last block row (excluding the secondary diagonal).

Part II -Solution of the equations for different scenarios
The solution of the Eq. A14 for the output correlation structure