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.

1 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 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 cross-correlation 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.

2 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].

Fig. 1

Cross-correlations in the Gaussian variables lead to correlations of coincident level crossings. a Spike correlations can arise from common input in a neuronal network. b We consider coincident level crossings arising from two Gaussian processes that share a common latent source. Whenever the voltage crosses a threshold ψ from below a spike is emitted. Spikes are indicated by vertical solid lines. The vertical dotted lines indicate the width of a time bin T used to compute spike counts ,

2.1 Definitions of Multivariate Voltage Distributions

We begin by defining the random, temporally correlated Gaussian zero mean process which represents the voltage of a neuron j

where is a combination of filters , 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. with and are complex random measures with independent increments, such that for all Borelian sets we have , the Lebesgue measure of A. By we denote the common noise component. Moreover if then and are independent Gaussian random variables. The correlation strength r, , denotes the presynaptic overlap of neurons 1 and 2 generated in a neuronal network, and it is illustrated in Fig. 1. If the voltages and are independent if the voltages and are identical. The auto- and cross-correlation functions between and are, respectively

where , and τ is the considered delay. The vector comprising the voltages and their derivatives is Gaussian and has the covariance matrix

where the variances are , and covariance functions are given by

We use the correlation time to quantify the width of the correlation function :

If the filters and 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 is the residual Gaussian current fluctuation with variance σ, the membrane time constant of the neuron, e.g. [12-15]. The synaptic drive can be separated into two parts: a common and an individual noise component

where 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 is given by . Analogously to Eq. (3), we obtain

2.2 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-15] generate a spike in neuron j any time a voltage 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 and quantify the cross-correlation between level crossings in neurons 1 and 2 by the following level functional:

where is a bounded and finite rectangle in , ψ 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 () and keep all other parameters the same. 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 , defined as . Ω denotes the probability space; here Ω is the Gaussian probability space. The field Z is Gaussian and and are equal in distribution. We denoted by the bivariate Gaussian density of vector . If and the prerequisites of Theorem 6.2 in [2] are fulfilled we can write

where the expectation value is denoted by , and is the determinant of the correlation matrix for the vector field . 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 is not degenerate. If we let and be two finite and bounded intervals in ℝ, condition (iv) is satisfied because

Here ℙ denotes the probability measure. We can define the correlation of two spike trains as

where is the joint Gaussian density of the vector . The conditional firing rate then is

where is the firing rate of a neuron j, for . In the next sections we provide closed-form expressions for and .

3 Cross-Correlations of Two Upward Level Crossings

Here, we address and provide a closed-form solution that is valid for any cross-correlation strength r between two level-crossing processes and any time delay τ.

Proposition 3.1Following the steps outlined in the methods section, Sect. 6.1, we can apply a regression model and Mehler’s Formula (see Lemma 10.7 in[2]) to prove that

whereis a joint Gaussian distribution for voltagesas defined in Sect. 2.2 andis the series given by

where

and

with

ϕandΦare the standard Gaussian density and distribution, respectively, and. are the Hermite polynomials. Note that the first two terms in Eq. (18) correspond to truncation ordersand, respectively.

To aid the explicit evaluation of in Eq. (18) we provide code for a computer algebra system.a Figure 2(a), (b) demonstrate obtained using Eq. (17) for progressively large truncation orders n. Notably, we find close correspondence between the first truncation order and the large n limit (). Figure 2(c), (d) show vs. τ as in Eq. (17) for varying correlation strength r. For simplicity, we chose and . We note that for two identical neurons () is a symmetric function while for a pair of neurons with different rates (), is asymmetric.

Fig. 2

Convergence of pairwise level crossing correlations. a Spike correlation function vs. time lag τ for different truncation orders n in Eq. (17); , , , , . b vs. τ for a pair of rate heterogeneous neurons with , , , , , , . c vs. time lag for varying correlation strengths r in a pair of neurons with , , , . d vs. time lag for varying correlation strengths r, in a pair of rate heterogeneous neurons with , , , , , . The truncation order of -series in Eq. (17) in c–d is . In a–dthe filled circles at indicate the predicted (as in Eq. (23)) and colored squares denote the corresponding numerical simulations obtained with independent realizations of 20 s length

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 level-crossing processes. For multiple cross-correlated processes orthant probabilities that describe expressions of the form 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 () or strong (, ) cross-correlations [3,10]. However, to address upcrossing correlations in the intermediate regimes where neither nor 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 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 and a given finite correlation strength 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 well-studied parameter dependence. This makes it possible to predict the influence of a specific parameter, such as time scale , firing rate 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 obtained using Eq. (17) for different truncation orders n alongside the zero lag correlation . Figure 2(c), (d) demonstrates obtained using Eq. (17) as a function of the correlation strength r alongside the zero lag correlation . As previously, we chose and . We note that for two identical neurons () is a symmetric function. Yet, for a pair of neurons with different rates () the spike correlation function is asymmetric, indicating that the lower rate neuron spikes on average after the higher rate neuron.

3.1 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-and-fire framework [12-15]. In a current-based leaky integrate-and-fire model driven by fluctuating noise the voltage of a neuron is given by

where is the input current of a neuron, a white noise, unit variance drive. The voltage power spectrum for this model is a combination of low-pass filters and its correlation function can be determined according to Eqs. (8). If the voltage reaches the threshold ϕ the neuron j emits a spike and the voltage is subsequently reset to a reset value . The integrate-and-fire 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 (). Numerically, the validity of the approximation remained highly accurate even for synaptic time constants .

A number of spike correlation results have been derived in the leaky integrate-and-fire model for the limit of weak correlations [18-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.

Fig. 4

Finite asymptotic limit of covariancesand correlation coefficient. a Count covariance vs. time bin T for varying correlation strengths r in the case of two identical neurons. The covariances converge for increasing time bin T towards the value predicted in Theorem 4.1 (indicated by small, thick lines). Here , and . b Limit value of the covariance vs. r for (from a) vs. r. The case of (black) corresponds to the variance , indicated by the dashed horizontal line. c Correlation coefficient vs. r for large time bins (). The dashed line indicates the equality line

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-28]. Spike counts and their cross correlations in neuroscience are often computed for a variety of bin sizes varying from [22] to [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 converges for to a one-dimensional 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 and are uncorrelated, then only if the vector 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 correlation between the variables X and Y does not imply independence: and . We obtain , 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 one-dimensional 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.

Fig. 3

Univariate and multivariate count Gaussianity. a Univariate Gaussianity of counts for large bin size . A solid black line represents the corresponding zero mean Gaussian fit. b Shapiro test p-value of the projections , for all . c Probability density of (left). QQ-plot between the theoretically predicted -quantiles and the empirical quantiles of in Eq. (35). Figures b and c are both validations of the multivariate Gaussianity of the bivariate vector X in Eq. (34). In all panels, ,

Theorem 4.1Let, be two processes satisfying Eq. (1), with covariancewhereand a common spiking thresholdψ. To take advantage of the available Gaussianity proofs that are typically derived for unit variance processes, we rescale the voltagesand the thresholdψto obtain processeswith unit variance and unit variance of the derivatives. then has the correlation functionand the spiking thresholds are. The number of-level upcrossings in a time intervalTfor processis given by. We assume that the following necessary and sufficient conditions hold. First, . This holds if and only ifsatisfies Geman’s Condition (see, e.g., Theorem 3 in[30]). Second, whereis the order of derivation. Then

where the count covarianceswithare three convergent series. Each is then given byand, both of which are finite. The first two terms in these series for, are

and the first and second order covariances are

We note that only forwe obtain. The asymptotic Pearson correlation coefficient, defined by

will also converge to the respective ratio of the asymptotic covariances and variances.

4.1 Numerical Confirmation of Joint Gaussianity and Limit Covariances

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 , and , 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 the consecutive sample number. We project each two-dimensional element in all directions by calculating the scalar product , with . 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 , where d is the dimensionality, μ is the mean and Σ the standard deviation, then the Mahalanobis distance with entries

is distributed according to a -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 and compare it with a -distribution, using the QQ-plot method (see Fig. 3(c)).

Figure 3 demonstrates the results of the joint Gaussianity tests for a bin size , where . Figure 3(a) shows the empirical univariate distribution of spike counts in one level-crossing process derived from independent count realizations. Figure 3(b) demonstrates that in independent count realizations of Xp-values for all θ are above the 10 % significance level. Figure 3(c) (left) illustrates that the Mahalanobis distance of a two-dimensional spike count variable X are well approximated by the -distribution (solid line). Figure 3(c) (right) demonstrates in a QQ-plot of the empirically measured -quantiles and the theoretical -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 introduced in Theorem 4.1 in Eq. (26). We choose , , for independent count realizations. To numerically compute the covariances from spike count simulations we used the covariance-matrix estimator proposed by [31]. Figure 4(a) demonstrates the convergence of covariance 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 is close to zero for and it is close to the variance for . Figure 4(c) demonstrates that the interplay between covariances and variances and leads to a sublinear dependence of the asymptotic correlation coefficient (Eq. (33)) in the large time bin limit, .

5 Conclusions

Level-crossing phenomena occur in a variety of physical and biological sciences. In many of these situations coordination between level crossings of multiple cross-correlated Gaussian processes is of interest. Here, we focused on neuroscience and modeled the spikes of two cross-correlated neurons by two cross-correlated level-crossing 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 cross-correlations and (2) the joint Gaussian limit of level-crossing counts. Our first result provides an explicit solution to 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 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 () 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.

6 Materials and Methods

6.1 Proof of Proposition 3.1

For simplicity of notation we adopt the following convention: denotes the vector . In order to calculate , defined in (14), we use the following regression model:

where and are independent. We use the covariance matrix in (4) and obtain

The conditional distribution is a bivariate Gaussian distribution

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 and are the Hermite coefficients associated with the expectation . Then these Hermite coefficients are given by

In particular,

– for , , ,

– for , , ,

– for , , .

Analogously we obtain

with

In particular,

– for , , ,

– for , , ,

– for , , .

We also have

with

In particular,

– for , , ,

– for , , ,

– for , , .

Finally,

with

In particular,

– for , , ,

– for , , ,

– for , , .

Then, combining the four expectations (38)–(41), the associated Hermite coefficients, and Eq. (37), we obtain

Note that the first two terms in Eq. (42) correspond to orders and , respectively. Denoting we find . Here, is a uniformly convergent series.  □

6.2 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 we can write

Now, we substitute the variables:

The correlation matrix for is a four-dimensional identity matrix.

Solving this integral we obtain Eq. (23).

6.3 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 and . Defining the level count deviation for neuron i by 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 are independent and have finite variance. From Mehler’s Formula we recognize that the contributions for are independent. The finite variance follows from the observation that for all q the variance of 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 is asymptotically Gaussian and satisfies a Central Limit Theorem, for all . We start from Eq. (45) and consider a truncated series for consisting of the first Q terms and denote by the norm in . First, we show that the remainder is bounded

As Q grows, this difference diminishes such that if is Gaussian for large Q then this will imply that is Gaussian, too. Now, we only need to show Gaussianity of . Using a modified version of the Breuer–Major Theorem (Sect. 6.4), we know that for each q

where is given in Eq. (51). The same theorem implies that and are asymptotically independent if . Thus we obtain for any

To calculate the count covariances in Eq. (26) we use Lemma 10.7 in [2] and obtain

Applying the Mehler Formula to a four-dimensional Gaussian random vector (Lemma 10.7 in [2]) we find

Here, is the Kronecker delta, Z is the set of ’s satisfying: , , , , and . We thus can express the covariances of the spike counts as

This is the result reported in Theorem 4.1.  □

6.4 Modified Breuer–Major Theorem

Here, we adapt the Breuer–Major Theorem [34] to show that the bivariate vector is Gaussian.

Theorem 6.1We consider two zero mean and unit variance Gaussian processes, which are described by the properties in Sect. 2 and a correlation function, where. For all functionsthat fulfillandand two real constantsthe following integral is Gaussian and we have

where

Proof We start by considering the Hermite expansion of the functions and write

where denotes the sum over q in Eq. (52) truncated at Q, and the convergence is in . We write

Now, to prove the Gaussianity of Eq. (54) it is sufficient to prove the asymptotic Gaussianity of . This is sufficient because

where is the norm in . Each of the terms is bounded,

where and a is a positive real number such that whenever . Since Eq. (60) is a vanishing series, . Now, we address the Gaussianity of . A result of Kratz and León, Lemma 3, p. 653 in [35], implies that

can be approached in by the sequence of ϵ-dependent processes . The -dependent processes for , are defined as

where ∗ denotes a convolution, is , β being a positive function with , and such that its Fourier transform has support in , where these processes are -dependent. Then

This follows from the Central Limit Theorem for -dependent random vectors and concludes the proof. □

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

The main idea of this paper was proposed jointly by EB, JL and TT. All authors contributed equally to the writing of this paper performed all the steps of the proofs in this research. All authors read and approved the final manuscript.