Cross-Correlations and Joint Gaussianity in Multivariate Level Crossing Models

A variety of phenomena in physical and biological sciences can be mathematically understood by considering the statistical properties of level crossings of random Gaussian processes. Notably, a growing number of these phenomena demand a consideration of correlated level crossings emerging from multiple correlated processes. While many theoretical results have been obtained in the last decades for individual Gaussian level-crossing processes, few results are available for multivariate, jointly correlated threshold crossings. Here, we address bivariate upward crossing processes and derive the corresponding bivariate Central Limit Theorem as well as provide closed-form expressions for their joint level-crossing correlations.


Introduction
Various phenomena in the biological or physical sciences are amenable to the description by level crossings of random Gaussian processes [1,2]. Examples of these phenomena are spike coordination of neurons in the brain [3], insurance risk assessment [4] and stress levels generated by ocean waves [5]. Therefore a number of mathematical studies in recent decades have focused on the statistical properties of level E. Di Bernardino Laboratoire Cédric EA4629, Conservatoire National des Arts et Métiers, 292 rue Saint-Martin, Paris Cedex 03, France e-mail: elena.di_bernardino@cnam.fr crossings arising from stationary Gaussian processes [2]. However, largely this literature addresses the properties of one level-crossing process and rarely deals with the coordinated level crossings of multivariate Gaussian processes. A prominent application where correlated level crossings are of particular importance is neuroscience. Recent work has shown that the spikes of a cortical neuron can be approximated by a Gaussian level-crossing process [3,6]. The assumption of Gaussianity is prompted by the experimental observation that cortical neurons are on average connected to ∼10000 neurons and therefore receive a barrage of inputs that together lead to a near-Gaussian fluctuation at the cell body of any given cortical neuron [7]. The spikes of two neurons are then modeled as upward level-crossing times of two cross-correlated fluctuating Gaussian potentials.
In this article we aim to address two features of level crossings of multiple correlated Gaussian processes. First, we want to clarify whether level-crossing counts derived from multiple correlated processes are jointly Gaussian. Second, we want to understand how many more coincident level crossings in a given time instance are expected if the underlying Gaussian random processes are correlated. Let us provide an intuitive reason for these questions. Starting with the first question, we recognize that if level-crossing counts of two neurons were jointly Gaussian, then a simple measure of dependence is the covariance or the Pearson correlation coefficient. Measuring a vanishing correlation coefficient or vanishing covariance between two neuronal spike counts would in this case imply true statistical independence, because only in the case of multivariate Gaussian distribution is it permissible to conclude independence from vanishing count correlations. This implication is not permissible if the marginal distributions are not Gaussian or are Gaussian but the joint distribution is not a multivariate Gaussian distribution. While marginal Gaussianity has been shown for level-crossing counts in [2] for large bin sizes, joint Gaussianity is still an open question. It might seem natural to imply joint Gaussianity from marginal Gaussianity for multivariate level-crossing processes, however, numerous counter examples exist to prove this intuition wrong, see Sect. 5 in [8]. Here, we use a modified Breuer-Major Theorem to prove joint Gaussianity and show that any linear combination of level-crossing counts of the two processes is also Gaussian.
The second question we address in this article deals with the conditional probabilities of two level-crossing processes. We are interested in how the level crossings of one Gaussian process can be used to predict the level-crossing probability of the partner process in a specific time interval relative to the observed level crossing in one process. In neuroscience, coordinated neuronal firing drives changes in synaptic connectivity and calculating the spike count dependencies across neurons is therefore a topic of current research efforts (e.g. Chap. 8 in [9]). The available mathematical results for conditional upward crossings in Gaussian processes currently comprise mostly variance and moments for one level-crossing process (see Chaps. 3-5 in [2]) as well as the low and high correlation limit in pairs of processes [3,10]. As yet, a comprehensive closed-form solution covering the complete level-crossing crosscorrelation function is currently lacking. Here, we use a regression approach to derive, for all correlation strengths, the conditional level-crossing correlation functions in two continuous Gaussian processes. We hypothesize that the level-crossing correlations we provide in this article could also be valuable in other fields outside of neuroscience for example in risk assessment calculations to predict the risk of joint default for insurance purposes.
The article is structured as follows. In Sect. 2 we define the mathematical model setting and introduce the concept of level crossings and specifically the upward crossings. In Sect. 3 we use a regression approach to obtain a general closed-form solution for cross-correlations of level crossings in two correlated Gaussian processes. In Sect. 4 we prove the joint Gaussianity (Central Limit Theorem) for the correlated joint upward crossings for two correlated Gaussian processes. In the section on materials and methods (Sect. 6) we provide detailed derivations of the reported results. We assume throughout this article that both level-crossing processes arise from crossings of the same threshold level by two Gaussian processes with different variances. This is permissible because the number of level crossings, the Rice rate [11], depends only on the variance-to-threshold ratio, but not on these quantities individually. We therefore work with a pair of level-crossing processes where each process has a unique voltage variance and therefore the rate of crossings in the two neurons being considered are, unless stated otherwise, not the same. Let us note that this assumption is prompted by the observation that in a living brain typically no two neurons are identical in all their properties and differ at least in their firing rate.

Mathematical Definitions of Multivariate Level Crossings
Here we address the statistics of coincident level crossings arising from two Gaussian processes that share a common latent source. This situation is illustrated in Fig. 1(a). We choose to illustrate the situation using a neuroscience perspective. Neurons in a mammalian brain receive synaptic inputs, both excitatory and inhibitory, from thousands of other neurons. Particularly in the visual cortex, the excitatory and inhibitory inputs largely cancel each other and lead to a net fluctuating residual current at the cell body of each neuron. These residual fluctuations drive the spikes of individual neurons. These voltage fluctuations arise from largely independent inputs so they are well approximated by a random Gaussian process with temporal correlations determined by the temporal structure of synaptic currents [7].

Definitions of Multivariate Voltage Distributions
We begin by defining the random, temporally correlated Gaussian zero mean process V j (t) which represents the voltage of a neuron j where f V is a combination of filters f V (λ) = γ (λ)g(λ), where γ represents the membrane filter and g the synaptic filter. Both of these filters can be chosen freely, but their combination should guarantee a continuously differentiable voltage trajectory.
where c(τ ) = ∞ −∞ e iτ λ f V (λ)dλ, and τ is the considered delay. The vector (V 1 (0), V 1 (0), V 2 (τ ), V 2 (τ )) comprising the voltages and their derivatives is Gaussian and has the covariance matrix where the variances are σ 2 and covariance functions are given by We use the correlation time τ s to quantify the width of the correlation function c(τ ): If the filters γ (λ) and g(λ) are classic low-pass filters, then the correlated voltage processes of Eq. (1) can be written in a differential form for each neuron j : where I j (t) is the residual Gaussian current fluctuation with variance σ , τ M the membrane time constant of the neuron, e.g. [12][13][14][15]. The synaptic drive I j (t) can be separated into two parts: a common and an individual noise component where C ξ (λ) is the synaptic noise spectral density. Using Eq. (6) we obtain the following spectral representation for the stationary solutions: In this form the spectral density of each V j is given by

Upward Crossing Definitions
Neurons communicate using brief pulses, the so called spikes, which are emitted whenever a voltage threshold is crossed [9]. The integrate-and-fire-type neuron models that are frequently used in computational neuroscience [12][13][14][15] generate a spike in neuron j any time a voltage V j (t) crosses a fixed threshold ψ and subsequently reset the voltage to a reset potential. Recently, it has been shown that in many physiologically relevant cases the leaky integrate-and-fire model can be equivalent to a level-crossing model without reset, where spikes are modeled as positive threshold crossings and are not followed by a reset [3,6,16]. Here, we therefore identify the spikes of a neuron j with the positive level crossings of its voltage V j (t) and quantify the cross-correlation between level crossings in neurons 1 and 2 by the following level functional: where Q(ε) = I 1 × I 2 is a bounded and finite rectangle in R 2 , ψ denotes the voltage threshold in both neurons (see Fig. 1). Here δ is introduced to quantify the infinitesimal interval around the threshold ψ where a spike takes place. We choose the same threshold for both neurons and two different variances (σ V 1 = σ V 2 ) and keep all other parameters the same. σ V 1 = σ V 2 represents the biological situation in which two neurons of the same neuronal type could have differences in the strength of their synaptic input and threshold-to-variance ratio but are exposed to the same temporal background statistics. We will consider the following random field Z : R 2 × Ω → R 2 , defined as (s 1 , s 2 ) → Z(s 1 , s 2 ) = (V 1 (s 1 ), V 2 (s 1 + s 2 )). Ω denotes the probability space; here Ω is the Gaussian probability space. The field Z is Gaussian and Z(s 1 , s 2 ) and Z(0, s 2 − s 1 ) are equal in distribution. We denoted by p s 2 −s 1 (·, ·) the bivariate Gaussian density of vector and the prerequisites of Theorem 6.2 in [2] are fulfilled we can write where the expectation value is denoted by E, and det(Z(s 1 , s 2 )) is the determinant of the correlation matrix for the vector field Z(s 1 , s 2 ). Now, we are left to prove the conditions of Theorem 6.2 in [2]. First, we find that conditions (i) and (ii) of Theorem 6.2 are satisfied by definition. Condition (iii) holds because p s 2 (ψ, ψ) is not degenerate. If we let I 1 and I 2 be two finite and bounded intervals in R, condition (iv) is satisfied because Here P denotes the probability measure. We can define the correlation of two spike trains as where The conditional firing rate ν cond (τ ) then is where is the firing rate of a neuron j , for j = 1, 2. In the next sections we provide closed-form expressions for s 1 (t)s 2 (t + τ ) and ν cond (τ ).

Cross-Correlations of Two Upward Level Crossings
Here, we address s 1 (t)s 2 (t + τ ) and provide a closed-form solution that is valid for any cross-correlation strength r between two level-crossing processes and any time delay τ .
where p τ (·, ·) is a joint Gaussian distribution for voltages V j as defined in Sect.

and C (a,b) (τ ) is the series given by
where and φ and Φ are the standard Gaussian density and distribution, respectively, and 2 )e z 2 /2 are the Hermite polynomials. Note that the first two terms in Eq. (18) correspond to truncation orders n = 0 and n = 1, respectively. To aid the explicit evaluation of C (a,b) (τ ) in Eq. (18) we provide code for a computer algebra system. 1 Figure 2(a), (b) demonstrate ν cond (τ ) obtained using Eq. (17) for progressively large truncation orders n. Notably, we find close correspondence between the first truncation order n = 1 and the large n limit (n = 10). Figure 2 Let us now briefly discuss the result we obtained in Eq. (17) within the context of previous level-crossing literature. One of the closely related results is the variance of level crossings and maxima provided in Proposition 4.4 in [2]. However, this result is derived for one level-crossing process, while we addressed a pair of levelcrossing processes. For multiple cross-correlated processes orthant probabilities that describe expressions of the form P(V 1 (t) > ψ 1 , V 2 (t) > ψ 2 ) have been obtained, e.g. Lemma 4.3 on p. 78 in [2]. However, specific results for the cross-correlations of upward crossings are sparse. For two correlated upward level-crossing processes previous studies have addressed the limiting cases of weak (r ≈ 0) or strong (r ≈ 1, r < 1) cross-correlations [3,10]. However, to address upcrossing correlations in the intermediate regimes where neither r ≈ 0 nor r ≈ 1 no analytical methods are available. Therefore, it was previously necessary to numerically evaluate the associated Gaussian probability densities in Eq. (15). The direct numerical evaluation of multidimensional Gaussian integrals can be computationally and algorithmically demanding, requires adaptive grid procedures and its accuracy can be hard to evaluate [17]. For the specific case of τ = 0 we show in the materials and methods section on 'Zero time lag correlations', Sect. 6.2, that a direct evaluation of the Gaussian double integral is possible via a variable substitution. The key to this variable substitution method was a manageable unity correlation matrix. For any other finite τ > 0 and a given finite correlation strength 0 r < 1 we could not identify a transformation that leads to a tractable integral and we therefore derived the series expansion in Eq. (17). This solution is explicit such that each series term of order n can be evaluated and studied separately. Furthermore, Eq. (17) consists of analytical functions with a wellstudied parameter dependence. This makes it possible to predict the influence of a specific parameter, such as time scale τ s , firing rate ν i or correlation strength r on the upward level-cross correlations. As an example, we evaluate Eq. (18) for an identical pair of neurons up to the third order in r via a Taylor expansion. We obtain We recognize that the linear and quadratic expressions are equivalent to the first and second order r-expansion reported in [3,10]. The cubic term has not been reported before, to the best of our knowledge. This demonstrates consistency with previous results and illustrates that expansions of any order can be obtained via Eq. (18). In this context, it is desirable to have an exact reference point for deciding how many n-orders are necessary for Eq. (17) to be sufficiently accurate. Such a reference point can be the zero lag value which we calculate exactly. Deviation of Eq. (17) for a specific n from this reference point can serve as an accuracy benchmark. Following the calculations in the methods Sect. 6.2, we derive Figure 2(a), (b) demonstrates ν cond (τ ) obtained using Eq. (17) for different truncation orders n alongside the zero lag correlation ν cond (0). Figure 2(c), (d) demonstrates ν cond (τ ) obtained using Eq. (17) as a function of the correlation strength r alongside the zero lag correlation ν cond (0). As previously, we chose c(τ ) = cosh(τ/τ s ) −1 and r ∈ [0, 1). We note that for two identical neurons (σ V 1 = σ V 2 ) ν cond (τ ) is a symmetric function. Yet, for a pair of neurons with different rates (σ V 1 = σ V 2 ) the spike correlation function ν cond (τ ) is asymmetric, indicating that the lower rate neuron spikes on average after the higher rate neuron.

Relation to the Leaky Integrate-and-Fire Model
Here, we address the relation between spike statistics and spike correlations in the level-crossing setting in our article and previous results in the leaky integrate-andfire framework [12][13][14][15]. In a current-based leaky integrate-and-fire model driven by fluctuating noise the voltage of a neuron V j (t) is given by where I j (t) is the input current of a neuron, ξ(t) a white noise, unit variance drive. The voltage power spectrum for this model is a combination of low-pass filters f V (λ) ∼ [(1 + τ 2 M λ 2 )(1 + τ 2 I λ 2 )] −1 and its correlation function can be determined according to Eqs. (8). If the voltage V j reaches the threshold φ the neuron j emits a spike and the voltage is subsequently reset to a reset value V r . The integrate-andfire model differs only in one important detail from the level-crossing approach-the presence of a reset after a spike. A recent article by Laurent Badel systematically compared the validity of upward level-crossing approximation for the firing rate, spike correlations and frequency response of a leaky integrate-and-fire neuron [16]. This study reached the conclusion that the upward level-crossing approach accurately represents the leaky integrate-and-fire model if two conditions are fulfilled: (1) the firing rate is much lower than the typical relaxation time of the voltage, (2) the synaptic filtering time constants remain of the same order of magnitude as the membrane time constant (τ I /τ M ≈ 1). Numerically, the validity of the approximation remained highly accurate even for synaptic time constants 0.4 τ I /τ M 2.6.
A number of spike correlation results have been derived in the leaky integrate-andfire model for the limit of weak correlations [18][19][20]. They include the observation that the spike correlation coefficient increases with firing rate [18,19]. The equivalent firing rate dependent increase in spike correlations and correlation coefficients for low correlation strengths has been reported for the level-crossing model, see [3] and Fig. 3(A) (right) and Fig. 2(A) (top) in [10]. Furthermore, leaky integrate-and-fire model exhibits a sublinear dependence of correlation coefficients on input strength r [18,21], which we see confirmed in Fig. 4.

Joint Gaussianity of Upcrossing Counts
Spike count cross correlations and correlation coefficient measurements in pairs of neurons are ubiquitous in neuroscience and are often used to measure the strength of interdependencies in a pair of neurons, e.g. in cortical neurons [18,19,22], in model neurons [23] and in theoretical and experimental studies of net correlations emerging in recurrent networks [23][24][25][26][27][28]. Spike counts and their cross correlations in neuroscience are often computed for a variety of bin sizes varying from T = 0.1-1 ms [22] to T = 2 s [29]. Here, we are interested in the question when spike count correlations of two neurons computed in a bin size T are jointly Gaussian such that their cross correlations are unbiased measures of statistical dependence or independence.
In this section we will prove that the spike counts of two neurons, approximated by up crossings of a Gaussian process, approach a joint multivariate Gaussian distribution for large bin sizes T . We start by considering the marginal distributions of upcrossing counts. From the one-dimensional Central Limit Theorem proven in [2] we know that 1 converges for T → ∞ to a onedimensional centered normal variable with finite variance. We provide a direct illustration of this classical univariate result in Fig. 3(a). Now, it is tempting to conclude that because the distribution of counts in each neuron is Gaussian, the joint distribution for the vector is also a multivariate Gaussian distribution. Yet, this conclusion is mathematically forbidden. While a joint Gaussian distribution implies marginal Gaussianity it is general not possible to inverse this relation [8]. The joint Gaussianity of spike counts is a highly desired property. If two counts are jointly Gaussian zero count correlation directly implies statistical independence. If count correlations between neuron 1 and neuron 2 are zero such that U V 1 [0,T ] (ψ) and U V 2 [0,T ] (ψ) are uncorrelated, then only if the vector (U V 1 [0,T ] (ψ), U V 2 [0,T ] (ψ)) is from a multivariate Gaussian distribution is it possible to infer independence of neuron 1 and neuron 2. Let us consider a teaching counter example where vanishing correla-tion between the variables X and Y does not imply independence: X ∈ N (0, 1) and Y = X 2 . We obtain Cov(X, Y ) = 0, but the two random variables are strongly linked. Contrasting examples of where X and Y variables are both marginally but not jointly Gaussian, have a zero correlation but are not independent can be found in Sect. 5 in [8]. To benefit from joint Gaussianity and be able to infer true statistical independence from vanishing correlations, we prove here the joint Gaussianity of upward level crossings/spike counts for large T . In the following we derive the Central Limit Theorem for the spike counts of two neurons, using two steps. First, we use the onedimensional Central Limit Theorem proven in [2]. Second, we apply a version of the Breuer-Major Theorem adapted to our upward crossing setting which we present in Sect. 6.4.
where i, j = 1, 2 and a common spiking threshold ψ. To take advantage of the available Gaussianity proofs that are typically derived for unit variance processes, we rescale the voltages V i (t) and the threshold ψ to obtain processes X i (t) with unit variance and unit variance of the derivatives.
then has the correlation function c(τ ) =

and the spiking thresholds are ψ i = ψ/ C V ii (0). The number of ψ i -level upcrossings in a time interval T for process X i is given by U X i [0,T ] (ψ i ).
We assume that the following necessary and sufficient conditions hold. First, E{(U where the count covariances a ij with i, j ∈ {1, 2} are three convergent series. Each is then given by 0 < a ii = ∞ q=1 σ 2 X i (q) and 0 < a 12 = ∞ q=1 σ X 1 ,X 2 (q), both of which are finite. The first two terms in these series for ψ = ψ i , and the first and second order covariances are We note that only for r = 1 we obtain σ X 1 ,X 2 (j ) = σ 2 X 1 (j ). The asymptotic Pearson correlation coefficient ρ T , defined by will also converge to the respective ratio of the asymptotic covariances and variances.

Numerical Confirmation of Joint Gaussianity and Limit Covariances a ij
In the last section we showed that spike counts of two neurons in large bins approach a bivariate Gaussian distribution with finite variances. Here, we illustrate this theoretical result in simulations of level-crossing processes. We choose two methods based on linear combinations and the Mahalanobis distance to numerically confirm joint Gaussianity. First, we numerically confirm the joint Gaussianity by showing that all linear combinations of two simulated spike counts for a large bin are Gaussian. We consider a vector where are the spike counts and their average in neuron 1 and 2, respectively. X consists of N × two-dimensional samples. N denotes the number of sample realizations and i ∈ [1, N] the consecutive sample number. We project each two-dimensional element X i in all directions (cos(θ ), sin(θ )) by calculating the scalar product A = X T i · (cos(θ ), sin(θ )), with θ ∈ [0, 2π). Subsequently, we apply a Shapiro test on A to verify Gaussianity for all univariate projections. The p-value of this Shapiro test conveys the certainty with which the Gaussian hypothesis cannot be rejected. As a second numerical test of joint Gaussianity we use the Mahalanobis distance. This test is based on the fact that if X ∼ N d (μ, Σ), where d is the dimensionality, μ is the mean and Σ the standard deviation, then the Mahalanobis distance D 2 with entries D 2 is distributed according to a χ 2 d -distribution with d degrees of freedom (see, e.g., Sect. 3.1.4 and Eq. (3.16) in [4]). By numerically estimating the count sample average μ and Σ we calculate in our case D 2 i and compare it with a χ 2 2 -distribution, using the QQ-plot method (see Fig. 3(c)).  Figure 3(a) shows the empirical univariate distribution of spike counts in one level-crossing process derived from N = 10000 independent count realizations. Figure 3(b) demonstrates that in N = 10000 independent count realizations of X p-values for all θ are above the 10 % significance level. Figure 3(c) (left) illustrates that the Mahalanobis distance D 2 of a two-dimensional spike count variable X are well approximated by the χ 2 2 -distribution (solid line). Figure 3(c) (right) demonstrates in a QQ-plot of the empirically measured D 2 -quantiles and the theoretical χ 2 2 -quantiles that they are linearly related. This is an indication that both distributions are equal. Figure 4 addresses the numerical confirmation of the constant asymptotic covariances a ij introduced in Theorem 4.1 in Eq. (26). We choose ψ = 0.3, τ s = 1, for N = 5000 independent count realizations. To numerically compute the covariances a ij from spike count simulations we used the covariance-matrix estimator proposed by [31]. Figure 4(a) demonstrates the convergence of covariance a 12 to a finite value that is predicted by Theorem 4.1. The asymptotic large T limit is denoted by a brief colored horizontal line. Figure 4(b) demonstrates the dependence of this asymptotic limit on the correlation strength r. As expected, we find that the asymptotic estimated covariance a 12 is close to zero for r = 0 and it is close to the variance for r = 1. Figure 4(c) demonstrates that the interplay between covariances a 12 and variances a 11 and a 22 leads to a sublinear dependence of the asymptotic correlation coefficient ρ T = a 12 /( √ a 11 a 22 ) (Eq. (33)) in the large time bin limit, T = 25τ s .

Conclusions
Level-crossing phenomena occur in a variety of physical and biological sciences. In many of these situations coordination between level crossings of multiple crosscorrelated Gaussian processes is of interest. Here, we focused on neuroscience and modeled the spikes of two cross-correlated neurons by two cross-correlated levelcrossing processes. While crossings and extrema of one level-crossing process have been the focus of mathematical research, results describing the coordination of multiple level-crossing processes are sparse and typically available only in specific and limited cases. Limits where level-crossing cross-correlations have been previously calculated are the weak and strong input correlation limit [3]. Here, we studied the case of two cross-correlated upward crossing processes and derived closed-form expressions for their joint level-crossing coordination as well their joint count Gaussianity. Importantly, the results we present in this article are consistent with previously reported limits but we now extended and generalized them. The two main results of our article are (1) closed-form explicit solution of the level-crossing crosscorrelations and (2) the joint Gaussian limit of level-crossing counts. Our first result provides an explicit solution to ν cond (τ ) = s 1 (t)s 2 (t + τ ) / √ ν 1 ν 2 that is valid for all correlation strengths and which comprises previously obtained limits, see discussion in Sect. 3. The rate of level crossings by a one-dimensional Gaussian process is given by the prominent Rice's equation derived by Rice in the 1950s [11]. The solution we obtained for the level-crossing cross-correlation ν cond (τ ) extends the Rice rate to the joint rate of two correlated processes. Our second result proves the joint Gaussianity of level crossings for large bin sizes. The joint Gaussianity of spike counts is a highly desired property because if and only if two level-crossing counts are jointly Gaussian can zero count cross-correlation imply statistical independence. Notably, marginal Gaussianity of spike counts in each neuron combined with zero count cross-correlation is not sufficient to imply independence. Contrasting examples of where X and Y variables are both marginally but not jointly Gaussian, have a zero cross-correlation but are not independent can be found in Sect. 5 in [8]. Count covariance and measures derived from it, such as the Pearson correlation coefficient, are computationally inexpensive and widely used as measures of statistical interdependencies [8]. Therefore, it is highly desirable to investigate the joint Gaussianity of level counts and thereby delimit the parameter space and mathematical conditions ensuring that independence can be implied from zero correlation coefficient. Notably, the joint Gaussianity of spike counts in bins of size T where T is much larger than the intrinsic time constant τ s (T τ s ) also implies that models of multi-neuronal dynamics only need to consider the mean and variance of spike counts because all higher cumulants are zero.

Materials and Methods
6.1 Proof of Proposition 3.1 For simplicity of notation we adopt the following convention: (X 1 , Y 1 , X 2 , Y 2 ) denotes the vector (V 1 (0), V 1 (0), V 2 (τ ), V 2 (τ )). In order to calculate s 1 (t)s 2 (t + τ ) , defined in (14), we use the following regression model: where ( 1 , 2 ) and (X 1 , X 2 ) are independent. We use the covariance matrix in (4) and obtain Using the regression system above, we write the conditional expectation where . Now, we calculate the four different expectations in Eq. (37) using Mehler's Formula (Lemma 10.7 in [2]). First, we write where c n (a) and c n (b) are the Hermite coefficients associated with the expectation Then these Hermite coefficients are given by In particular,

Zero Time Lag Correlations
In this section we derive the zero lag spike correlation using the Gaussian probability integrals in Eq. (15). Following the previously introduced notation the spiking threshold level in a neurons is ψ and variance σ V i we can write Now, we substitute the variables: .

Proof of Theorem 4.1
To prove joint Gaussianity, it is instructive to briefly recapitulate the one-dimensional Central Limit Theorem. We start by writing with d Defining the level count deviation for neuron i by S i (T ) we obtain where . A Gaussian distribution is a stable limit distribution for a sum of independent finite variance variables. Therefore, all that is left to prove is that contributions q = q are independent and have finite variance. From Mehler's Formula we recognize that the contributions for q = q are independent. The finite variance follows from the observation that for all q the variance of G i q (X i (s), X i (s)) is proportional to the expectation of a product of four Hermite polynomials, which has been proven to be finite (Theorem 10.10 in [2]) if the conditions of Theorem 4.1 are satisfied. Now, we address the joint Gaussianity. A random vector is jointly Gaussian if and only if any linear combination of its components has a univariate normal distribution (see [4,32,33]). Thus, we need to prove that the sequence α 1 S 1 (T ) + α 2 S 2 (T ) is asymptotically Gaussian and satisfies a Central Limit Theorem, for all α 1 , α 2 ∈ R. We start from Eq. (45) and consider a truncated series for S i Q (T ) = Q q=1 J 1 q (T , X 1 , X 1 ) consisting of the first Q terms and denote by · the norm in L 2 (Ω). First, we show that the remainder is bounded As Q grows, this difference diminishes such that if α 1 S 1 Q (T ) + α 2 S 2 Q (T ) is Gaussian for large Q then this will imply that α 1 S 1 (T ) + α 2 S 2 (T ) is Gaussian, too. Now, we only need to show Gaussianity of α 1 S 1 Q (T ) + α 2 S 2 Q (T ). Using a modified version of the Breuer-Major Theorem (Sect. 6.4), we know that for each q where σ 2 G q is given in Eq. (51). The same theorem implies that α 1 J 1 q (T , X 1 , X 1 ) + α 2 J 2 q (T , X 2 , X 2 ) and α 1 J 1 q (T , X 1 , X 1 ) + α 2 J 2 q (T , X 2 , X 2 ) are asymptotically independent if q = q . Thus we obtain for any Q ≥ 1 To calculate the count covariances in Eq. (26) we use Lemma 10.7 in [2] and obtain a ij = 2 ∞ q=1 q k=0 q k =0 d (1) q−k (ψ i )d (2) q−k (ψ j )a k a k  (1) q−k (ψ)d (2) q−k (ψ)a k a k This is the result reported in Theorem 4.1.
Theorem 6. 1 We consider two zero mean and unit variance Gaussian processes X i (t), which are described by the properties in Sect. 2 and a correlation function E[X i (0)X j (t)] = c ij (t), where i, j ∈ {1, 2}. For all functions G i (·, ·) that fulfill E[G i (X i (0), X i (0))] = 0 and E[G 2 i (X i (0), X i (0))] < ∞ and two real constants α i the following integral is Gaussian and we have Proof We start by considering the Hermite expansion of the functions G i (·, ·) and write G i X 1 (t), X 2 (t) = lim Q→∞ Q q=1 k 1 +k 2 =q c k 1 k 2 ,G i H k 1 X 1 (t) H k 2 X 2 (t) can be approached in L 2 (Ω) by the sequence of -dependent processes X ε