A Formalism for Evaluating Analytically the Cross-Correlation Structure of a Firing-Rate Network Model

We introduce a new formalism for evaluating analytically the cross-correlation structure of a finite-size firing-rate network with recurrent connections. The analysis performs a first-order perturbative expansion of neural activity equations that include three different sources of randomness: the background noise of the membrane potentials, their initial conditions, and the distribution of the recurrent synaptic weights. This allows the analytical quantification of the relationship between anatomical and functional connectivity, i.e. of how the synaptic connections determine the statistical dependencies at any order among different neurons. The technique we develop is general, but for simplicity and clarity we demonstrate its efficacy by applying it to the case of synaptic connections described by regular graphs. The analytical equations so obtained reveal previously unknown behaviors of recurrent firing-rate networks, especially on how correlations are modified by the external input, by the finite size of the network, by the density of the anatomical connections and by correlation in sources of randomness. In particular, we show that a strong input can make the neurons almost independent, suggesting that functional connectivity does not depend only on the static anatomical connectivity, but also on the external inputs. Moreover we prove that in general it is not possible to find a mean-field description à la Sznitman of the network, if the anatomical connections are too sparse or our three sources of variability are correlated. To conclude, we show a very counterintuitive phenomenon, which we call stochastic synchronization, through which neurons become almost perfectly correlated even if the sources of randomness are independent. Due to its ability to quantify how activity of individual neurons and the correlation among them depends upon external inputs, the formalism introduced here can serve as a basis for exploring analytically the computational capability of population codes expressed by recurrent neural networks. Electronic Supplementary Material The online version of this article (doi:10.1186/s13408-015-0020-y) contains supplementary material 1.


Introduction
The brain is a complex system whose information processing capabilities critically rely on the interactions between neurons. One key factor that determines interaction among neurons is the pattern of their anatomical or structural connectivity, namely the specification of all the synaptic wirings that are physically present between neurons. However, communication among neurons appears to change dynamically [1], suggesting the presence of not-yet understood network mechanisms that modulate the effective strength of a given connection. Understanding how the functional connectivity of a neural network (i.e. the set of statistical dependencies among different neurons or neural populations [2]) depends upon the anatomical connectivity and is further modulated by other network parameters has thus become a central problem in systems neuroscience [3][4][5][6][7][8].
In this article we introduce a new formalism for evaluating analytically the structure of dependencies among neurons in the finite-size firing-rate network with recurrent connections introduced in [9]. Although these dependencies are computed from neural activity in a number of ways [10], in most cases functional connectivity is inferred from computing the correlation among neurons or populations of neurons [2]. In this article, we therefore concentrate on computing the correlations among neurons in a firing-rate network, although we also discuss how to compute, with the same formalism, also other measures of functional connectivity (Sect. 5).
To our knowledge, the problem of determining analytically the correlation structure of a neural network has been begun to be investigated systematically only recently. This is in part due to the new experimental insights into functional connectivity among cortical neurons [3][4][5][6][7][8], and in part due to the focus on many previous mathematical studies of neural networks on the mean-field approximation. This approximation exploits the fact that (under certain hypotheses) neurons become independent in the thermodynamic limit when the number of neurons N in the network goes to infinity. This kind of mean-field approximation has been developed by Sznitman [11][12][13], Tanaka [14][15][16], McKean [17,18] and others. According to it, if the neurons are independent at time t = 0 (initial chaos), then in the thermodynamic limit this independence propagates to every t > 0. 1 This phenomenon of propagation of chaos has been studied in different kinds of neural network models [19][20][21][22]. However, recent studies have begun to investigate the more interesting and rich structure of correlations arising in descriptions of networks dynamics beyond the thermodynamic limit. For example, new studies considered finite-size networks with excitatory and inhibitory populations, where the firing rates are determined by a linear response theory [23][24][25]. These studies included in the network sources of Poisson randomness in the spike times [23,24], as well as randomness originating from normal white noise in the background for the membrane potentials [25]. Pioneer approaches [26] relied on estimating correlation by using a perturbative expansion around the thermodynamic limit in the inverse number of neurons in the network. The method was developed for binary neurons, where the sources of randomness were the transitions between the two states of each neuron and the topology of the synaptic connections, and a similar model was reintroduced recently in [27] for large networks. In [28] the author considered an alternative way to calculate correlations as a function of the inverse number of neurons (which is known as the linear noise approximation) and applied it to homogeneous populations of identical neurons with random fluctuations in the firingrates. In [29] the authors introduced a density functional approach adapted from plasma physics to study correlations in large systems, and applied it to a heterogeneous network of phase neurons with random initial conditions. Another effective approach is represented by large deviations theory. In [30][31][32] the authors considered a discrete-time network of rate neurons, whose sources of randomness were background Brownian motions for the membrane potentials and normally distributed synaptic weights.
Building on these previous attempts to study network correlations including finitesize effects that go beyond the mean-field approximation, here we develop an approach based upon a first-order perturbative expansion of the neural equations. We introduce randomness through three different sources: the background noise of the membrane potentials, their initial conditions and the distribution of the recurrent synaptic weights. These sources of variability are normally distributed and can be correlated, and their standard deviations are used as perturbative parameters. Using μ (N ) t is said to be ν t -chaotic if, for all integers k ≥ 1 and for all continuous and bounded test functions ϕ 1 , . . . , ϕ k , we have lim N →∞ R d×N ϕ 1 (x 1 ) · · · ϕ k (x k ) dμ (N ) t (x 1 , . . . , x N ) If the sequence μ (N ) t is ν t -chaotic, the neurons are said to be chaotic at time t . Intuitively, we can think of μ (N ) t as the joint probability of N (exchangeable) neurons. The exchangeability condition is equivalent to the above symmetry condition. We see that the case of independent and identically distributed neurons is a special one. Indeed, if μ (N ) t = ν ⊗N t (where ⊗ denotes the tensor product of measures), then it is obvious that the above condition is satisfied ∀N . In turn, the tensor product of measures can be interpreted as the factorization of the joint probability that defines independence in probability theory. The initial conditions are said to be ν 0 -chaotic if the previous condition holds at time t = 0. The propagation of chaos refers to the fact that if the initial conditions are ν 0 -chaotic, then, if the neurons are exchangeable, their joint law μ (N ) t is ν t -chaotic for some probability measure ν t on R d for all times t ∈ [0, T ].
this formalism and this model, we quantify analytically how synaptic connections determine statistical dependencies at any order (not only at the pairwise level, as in previous studies) among different neurons. The technique developed in this article is general, but for simplicity we demonstrate its efficacy by applying it to the case of synaptic connections described by regular graphs. A regular graph is a graph in which each vertex has the same number of neighbors, so this means that we consider networks where each neuron receives and makes the same number of connections. While this assumption is of course biologically implausible, it is sufficient to show interesting and non-trivial behaviors and will be relaxed to study more plausible connections in our future studies. We use this formalism to investigate in detail how the correlation structure depends on the strength of the external input to the network. We find that external input exerts profound and sometimes counterintuitive changes in the correlation among neurons: for example, a strong input can make the neurons almost independent. Moreover, we prove that in general it is not possible to find a mean-field description à la Sznitman of the neural network, due to the absence of chaos, if the anatomical connections are too sparse or our three sources of variability are correlated. This demonstrates the fairly limited range of applicability of the mean-field approximation. Finally, we also show a very counterintuitive phenomenon, which we call stochastic synchronization, through which neurons become almost perfectly correlated even if the sources of randomness are independent.
This article is organized as follows. In Sect. 2 we describe the details of the firingrate network we use. We then develop a first-order perturbative expansion (Sect. 3) that allows the approximate analytical calculation, for a generic anatomical connectivity matrix, of the membrane potentials and the firing rates of the network. (In this section we assume the reader to be familiar with stochastic calculus [33,34].) Then we use this formula for the membrane potentials and the firing rates in Sect. 4 to calculate analytically the pairwise and higher-order correlation structure of the network and the joint probability distribution for both the membrane potentials and the firing rates. In Sect. 5 we briefly discuss how other measures of functional connectivity can be evaluated analytically through our theory. In Sect. 6 we specialize to the case of regular graphs and we investigate network dynamics using some explicit examples of anatomical connectivity. We start by considering relatively simple cases, in particular a block-circulant graph with circulant blocks (Sect. 6.1) and a more general case of symmetric undirected graphs (Sect. 6.2). Then in Sect. 6.3 we conclude by showing how to extend the theory to highly complex regular graphs and by discussing also some possible extensions to irregular networks. In Sect. 7 we investigate the goodness of our perturbative approach by comparing it to the numerical simulation of the network's equations. In Sect. 8 we show that the correlation structure depends dynamically on the external input of the network. In Sect. 9 we demonstrate with counterexamples that in general Sznitman's mean-field approximation cannot be applied to the network in the case when the sources of randomness are correlated (Sect. 9.1) or when the anatomical connectivity matrix is too sparse (Sect. 9.2). In Sect. 10 we introduce the phenomenon of stochastic synchronization. Finally, in Sect. 11 we discuss the implications of our results as well as the strengths and limitations of our mathematical approach.

Description of the Model
We suppose that the neural network is described by the following firing-rate model [9]: is the membrane potential of the ith neuron, I i (t) is its external time-varying input current and τ is the membrane time constant describing the speed of convergence of the membrane potential to its resting state. Note that the external input can assume both positive and negative values, modeling the effect of prevailingly depolarizing or hyperpolarizing external influences, respectively. Moreover, J ij (t) is the synaptic weight from the j th to the ith neuron, while M i is the number of incoming connections of the ith neuron. In graph theory this quantity is called incoming vertex degree and its role will be explained later in this section. A (·) represents a generic activation function, which converts the membrane potential V of a neuron into its corresponding firing rate ν = A (V ). A typical choice is to consider S-shaped (or sigmoidal) activation functions, because they are biologically plausible and their saturation for |V | → ∞ ensures the boundedness of the solutions of Eq. (2.1). Some classic examples are shown below: Tangent), where in the above ν max is the maximum firing rate, Λ determines the speed with which the neuron switches from a low (ν ≈ 0) to a high (ν ≈ ν max ) firing rate, and V T is the threshold between low and high firing rates, namely the value of the membrane potential such that ν = ν max 2 . An example of the functions (2.2) is shown in Fig. 1 for some particular values of ν max , Λ, and V T .
The functions B i (t), the first of the three sources of randomness introduced in the network, are non-fractional Brownian motions (or in other terms, Wiener processes with independent increments). They can be equivalently interpreted as a background noise for the membrane potentials V i (t) or as the stochastic component of the external input I i (t). σ 0 is the standard deviation (or intensity) of the noise, that for simplicity is supposed to be the same for all the neurons and constant in time. This is the first perturbative parameter that will be used in Sect. 3 to develop a first-order perturbative expansion of the neural equations. In general the Brownian motions are correlated Cov Here δ ij is the Kronecker delta, δ(·) is the Dirac delta function and C (0) represents the correlation between two different Brownian motions (the derivative of the Brownian motion with respect to time here is meant in the weak sense of distributions and is interpreted as white noise). The covariance matrix must be positive-semidefinite. Since it is symmetric, then it is positive-semidefinite if and only if its eigenvalues are non-negative. Moreover, with our choice, the covariance matrix is circulant, therefore its eigenvalues are 1 + C (0) (N − 1) and 1 − C (0) , with algebraic multiplicity 1 and N − 1, respectively. Therefore the matrix is positive-semidefinite if and only if 1 1−N ≤ C (0) ≤ 1. Note that there are no technical obstructions to increasing the complexity of this correlation structure, if desired.
The initial conditions V i (0) are normally distributed around their mean μ i with standard deviation σ 1 , the second of our perturbative parameters. The stochastic variables N i (see Eq. (2.1)) are normally distributed with zero mean and covariance matrix: (1) otherwise, (2.4) namely Cov(N i , N j ) = δ ij + C (1) (1 − δ ij ). The coefficient C (1) is the correlation between pairs of membrane potentials at time t = 0. As before, the covariance matrix must be positive-semidefinite, and this is true if and only if 1 1−N ≤ C (1) ≤ 1. Again, it is possible to increase the complexity of this correlation structure, if desired.
The third and last source of randomness in the network is represented by the synaptic connectivity J (t). We assume that each entry J ij (t) is normally distributed around its mean J ij (t) with standard deviation σ 2 (the third perturbative parameter used in Sect. 3), or in other terms:

5)
W ij are zero mean normal stochastic variables (their covariance structure is shown below, see Eq. (2.6)), while the matrix T represents the topology of the connectivity matrix, namely the mere absence (T ij = 0) or presence (T ij = 1) of the synaptic connection from the j th neuron to the ith neuron, for all the pairs of neurons (i, j ). So if the connection is present, its strength is given by J ij (t) + σ 2 W ij , otherwise it is equal to zero. Below we show an example of connectivity matrix in a network of four neurons and its corresponding topology: In graph theory, T is known as the adjacency matrix of the unweighted graph of the network, and in this article is supposed to be deterministic and time-independent. Therefore the only source of randomness in the synaptic matrix is represented by W ij , whose covariance structure is chosen as follows: ifT ij = 0 and/or T kl = 0, 1 ifi = k and j = l and T ij = 1, C (2) otherwise, (2.6) or in other terms Cov(W ij , W kl ) = T ij T kl [δ ik δ jl + C (2) (1 − δ ik δ jl )]. This simply means that W ij has zero (unit) variance if the connection i ← j is absent (present), while the covariance between W ij and W kl is zero (C (2) ) if at least one of the connections i ← j and k ← l is absent (they are both present). As for the covariance structures (2.3) and (2.4), also (2.6) can be made more complicated, if desired. With our choice, the range of allowed values of C (2) depends on the topology of the connectivity matrix. In order to find this range, we start by vectorizing the matrices W and T as follows: This allows us to reinterpret the N × N matrix-variate random variable W as a N 2dimensional multivariate vector W with a N 2 × N 2 covariance matrix Cov(W i , W j ). Now we call Z the number of absent connections (i.e. the number of zeros in the matrix T ), and we suppose that the vectorization is such that T i = 0 for i = 0, . . . , Z − 1.
According to (2.6), if we call Θ the covariance matrix of W, we obtain where Id K is the K × K identity matrix. Therefore Θ has eigenvalues 0, 1 + C (2) (N 2 − Z − 1) and 1 − C (2) , with algebraic multiplicity Z, 1 and N 2 − Z − 1, respectively. This means that Θ is a true covariance matrix if and only if 1 1+Z−N 2 ≤ C (2) ≤ 1, where Z depends on the topology of the network.
In order to avoid biologically unrealistic sign changes of the synaptic weights, J ij (t) should not change sign during time evolution. Moreover, |J ij (t)| must be much larger than σ 2 for every i, j , and t, because in this way the probability that J ij (t) + σ 2 W ij changes sign from trial to trial is small: having used an asymptotic expansion of the error function for |J ij (t)| σ 2 1. Since in Sect. 3 we will solve perturbatively the system of equations (2.1), it is clear that this cannot be accomplished with the current formulation of the synaptic weights. Actually, even if our differential equations were linear, in general it would not be possible to solve them exactly, since their coefficients are time dependent, due to J ij (t). Linear differential equations with variable coefficients can be solved perturbatively in terms of a Magnus expansion [35], but this is not the approach followed in our work. In order to unify the perturbative expansion introduced in this article with the problem of the variation of the coefficients, we rewrite the matrix J ij (t) as follows: where J c ij is constant in time, while J v ij (t) is variable. σ 3 is assumed to be small and will be used as a perturbative parameter in Sect. 3. In this way the time variability of the synaptic matrix is treated perturbatively, as for the three sources of randomness in the network, and this will allow us to find analytical solutions for the perturbative expansion developed in the next section. Moreover, the variable part of the synaptic weights should satisfy the constraint |J v ij (t)| ∈ [0, 1] for all t, because according to (2.7) this implies max t |J ij (t) − J c ij | ≤ σ 3 . In this way we ensure that at every time instant J ij (t) is not too different from J c ij , and therefore that a first-order perturbative expansion of the neural equations provides a good approximation of the real dynamics of the network.
It is also important to observe that when a neuron receives more and more connections from the other neurons (i.e. when M i → ∞, see Eq. (2.1)), the sum ) in (2.1) is divergent, therefore the limit of large and densely connected networks would not be well defined. This explains the need to introduce the factor 1 M i to normalize the divergent sum. For later purpose, it is useful to express M i in terms of the topology of the network: (2.8) Finally, as we did for J ij (t), we suppose that the external input current is deterministic (if we interpret B i (t) as the background noise of the membrane potentials) and given by where I c i is constant in time, while I v i (t) is variable. σ 4 is our last perturbative parameter, and together with σ 3 quantifies the time variability of the network. As for the synaptic weights, we have the constraint |I v i (t)| ∈ [0, 1], because according to (2.9) this implies max t |I i (t) − I c i | ≤ σ 4 . For simplicity, the three sources of randomness are supposed to be independent from each other, namely: (2.10) This assumption reduces considerably the complexity of the formula for the correlation structure that we will calculate in Sect. 4, but can be relaxed if desired. This concludes our description of the neural equations, so now we have all the ingredients to develop a perturbative expansion of the system. This method is introduced in the next section, and will allow us to find a series of new results for the behavior of our stochastic neural network.

Perturbative Expansion
As we said in the previous section, we want to develop a first-order perturbative expansion of the neural network in terms of the perturbative parameters σ 0 -σ 4 . To this purpose we define the following first-order expansion of the membrane potentials: can be determined by substituting the perturbative expansion (3.1) and the expressions (2.5) + (2.7) and (2.9) for, respectively, the synaptic weights and the external input current, into the system (2.1). If all the parameters σ m are small enough, we can expand the activation function A (V i ) in a Taylor series about μ i . In order to be rigorous, we have to determine the radius of convergence r(μ i ) of the Taylor expansion for every value of μ i and to check if the radius is big enough compared to σ m , because otherwise the series does not converge. Actually, the various σ m determine the order of magnitude of the fluctuations of V i around μ i , therefore it is important to check if V i lies inside the interval of convergence of the Taylor expansion (this will be quantified more rigorously at the end of Sect. 4). In Appendix A we evaluate r(μ i ) for two examples of A (V i ), namely the logistic and the inverse tangent activation functions (see (2.2)). In both cases we obtain that the radius decreases with the slope parameter Λ, and since all the sigmoidal functions are qualitatively similar, it is reasonable to assume that this result holds for all of them. Therefore, supposing that Λ is small enough, the Taylor expansion of A (V i ) truncated at the first order is Now we substitute this expansion inside the neural equations (2.1) and we equate the terms with the same σ coefficients, obtaining

4)
dY (1) dY (4) J is the Jacobian matrix of the network, while J eff can be interpreted as the real anatomical connectivity matrix that the system would have if it were linear. For this reason we call it the effective connectivity matrix of the network, and it should not be confused with the effective connectivity discussed in [2]. Now we observe that equations (3.3) are algebraic and non-linear, therefore in general they must be solved numerically. Eventually, it is possible to obtain analytical solutions when the activation function is approximated by a piecewise linear function. Moreover, (3.4) ((3.5)-(3.8)) are linear stochastic (ordinary) differential equations with constant coefficients, therefore can be solved analytically as a function of μ i , which are supposed to be known from the solution of (3.3). The fundamental matrix Φ(t) of the system is where J is given by (3.9). In this article we consider only cases when J is diagonalizable, so we can calculate the matrix exponential as follows: where D(t) = diag(e λ 0 t , . . . , e λ N−1 t ) and λ i are the eigenvalues of J , while P is an N ×N matrix whose columns are composed of the eigenvectors of J . The differential equations (3.4)-(3.8) are linear with constant coefficients since, as explained in the previous section, we have used the perturbative approach to fix the problem of the time variability of the coefficients. So their solutions are obtained straightforwardly as follows: By substituting the solutions (3.12)-(3.16) inside (3.1), we obtain an approximate formula for the membrane potentials of the neural network. Moreover, since ν = A (V ), (3.2) provides a perturbative expression for the firing rates. Now with these results we can determine analytically the behavior of the neural network, starting from its correlation structure and probability density, which are discussed in the next section.

Cross-Correlation and Probability Density
The aim of this section is to compute the statistical dependencies among the activity of different neurons.
We first calculate the Pearson cross-correlation among pairs of neurons, which is the simplest and most commonly used measure of functional connectivity [10]. This is defined as follows: The subscript "2" has been introduced to stress the fact that it is a pairwise correlation between two neurons. We note that the above expression quantifies time-dependent instantaneous correlations at any given time t. This equation can easily be extended to higher-order correlations between triplets, quadruplets, etc. of neurons. The most straightforward generalization of the pairwise covariance to the case of n neurons would be where the bar represents the statistical mean over trials computed at time t. This is the multivariate moment M 1,...,1 of the functions . However, this measure is not yet normalized to lie in the range [−1, 1]. To achieve this purpose, we observe that having used the fact that |x + y| ≤ |x| + |y| at the first step and a special case of the Hölder inequality at the second. This means that the function Therefore, we will use Eq. (4.1) to quantify correlations at any order. Note that Eq. (4.1) is equivalent for n = 2 to the pairwise Pearson coefficient, and thus Eq. (4.1) includes also the pairwise correlation as a special case. From (3.12)-(3.14) and remembering that B is a zero mean multivariate normal process, whose covariance matrix is given by where (4.7) Using the Isserlis theorem [36], and noting that we assumed that our sources of randomness are normally distributed, we obtain that the numerator of (4.1) is equal to zero when n is odd (in general this is false for non-normal processes), otherwise where means summing over all the distinct n! 2 n/2 (n/2)! ways of partitioning N i 0 (t), . . . , N i n−1 (t) into pairs. This completes the calculation of the numerator of (4.1).
For the denominator we use the formula of the absolute moments of a normal process, therefore for n even we obtain where N 2 i (t) is given by (4.3) and (4.4)-(4.6) for i = j . Finally, by substituting all these results into the definition (4.1), we obtain the final formula for the higher-order correlation of the network.
We observe that computation of Eq. (4.8) leads to a combinatorial problem, whose complexity is related to n and to the structure of the effective connectivity matrix J eff . To simplify matters, in Appendix B we will calculate the higher-order correlation for a generic n in the case of a complete graph (i.e. a fully connected network).
We also observe that in the special case n = 2 our formula reduces simply to where again the terms Y (m) 2 are given by (4.4)-(4.6), so in this case the combinatorial problem is absent. Due to its simplicity and in order to keep the article as short as possible, in the next sections we will evaluate explicitly only the pairwise correlation structure through (4.10) (therefore the subscript "2" in the notation Corr 2 (·, ·) will be omitted). The interested reader could apply the general Eq. (4.1) for the calculation of the correlation structure when n > 2.
Neuroscientists make use of measures of correlation between firing rates, rather than between membrane potentials, to study cross-neuron communication. This is because only spiking events (and not subthreshold membrane fluctuations) are transmitted to other neurons. For this reason we also derive a formula for the correlation of the firing rates ν. Since in this model having used the fact that A (μ i j ) is always positive. This proves that However, it is important to observe that the correlation structures of the firing rates and the membrane potentials are equivalent only at the first perturbative order, namely when all the parameters σ m are relatively small. At higher orders this equivalence does not hold anymore. Now we have all the ingredients required to evaluate the joint probability distribution of the neural network. Since we have linearized the differential equations (2.1), at the first perturbative order the joint probability density of the system is a multivariate normal distribution. Denoting by the matrix transposition operator and defining In a similar way, if we define ν def = (ν 0 , . . . , ν N −1 ) and we remember that This completes the description of the system at the first perturbative order. Now, if we suppose that, for given values of σ 0 -σ 4 , the perturbative corrections of order higher than one are negligible, Eq. (4.11) can be used to evaluate the probability P(t) that, at the time instant t, all the activation functions in (2.1) can be expanded in a Taylor series according to (3.2) is the radius of convergence of the activation function evaluated at the point V i = μ i (see Appendix A), then P(t) is defined as follows: where × represents the Cartesian product of subsets. For a multivariate normal distribution, an analytical expression of P(t) is not known, therefore it must be evaluated numerically (see Sect. 7). So if P(t) ≈ 1, we can safely expand the activation function by using Eq. (3.2), therefore under this constraint all the results found in this article are justified. In other terms, this can be considered as a self-consistency check of the theory, which can be further refined if higher-order corrections are taken into account.

Other Measures of Functional Connectivity
In order to illustrate the generality of our approach, here we briefly describe how it can be extended to compute two other symmetric quantities commonly used to mea-sure the functional connectivity, namely the mutual information and the magnitudesquared coherence [10,37,38]. The mutual information between the membrane potentials (a similar formula holds for the firing rates) of two neurons i, j is defined as follows: where the last identity holds only for normal probability distributions, which is indeed our case. This shows that the mutual information is a simple function that depends trivially on the pairwise correlation between the neurons, which in turn implies that it can be evaluated directly from the results obtained in the previous section. A similar result holds for the magnitude-squared coherence between the membrane potentials (or the firing rates) of two neurons i, j . If we call C ij (t, ω) the Fourier transform of the time-shifted cross-correlation: (the imaginary unit is denoted by ι, to avoid confusion with the neural index i), then the magnitude-squared coherence is defined as follows: It is straightforward to extend Eqs. (4.4)-(4.6) to include the temporal shift s, which allows us to calculate Cov(V i (t), V j (t + s)). This means that the functional connectivity inferred from the correlation, the mutual information and the coherence is qualitatively the same. This further justifies our decision to focus this article only on cross-correlations. We note that our formalism lends itself in principle also to the calculation of directed asymmetric measures of functional connectivity, such as those based upon transfer entropy [39,40] or the Granger causality [41][42][43]. However, an analytical calculation of these directed quantities is beyond the scope of this article.

Examples
Now we consider some explicit examples of calculation of the correlation structure. First of all, it is important to observe that in this article we consider only cases when the Jacobian matrix J (see (3.9)) is diagonalizable, so we can calculate the fundamental matrix Φ as shown by Eq. (3.11). For this reason we need to know the eigenquantities of J . However, due to the eventual inhomogeneities of the static synaptic weights J c ij and of the incoming vertex degree M i , and to the non-linearity of the network introduced by the activation function A , in general it is not possible to find a simple relation between the eigenquantities of J and those of the underlying topology T . This means that, even if the matrix T has some special structure like circularity or symmetry, in general this cannot be exploited to calculate the eigenquantities of J , because the same structure is not preserved in J due to the term 1 M i J c ij A (μ j ) (see (3.9) + (3.10)). However, it is important to observe that this is not a problem per se. Actually the method introduced in this article has been applied to the study of relatively generic connectivity matrices, but these results will be published in other papers. For the sake of clarity, here we want to avoid complicated algebraic calculations of the eigenquantities, therefore we will consider the simplest case possible, namely neural networks where the term 1 where Γ is a free parameter that describes the strength of the synaptic connection (if present), while M is the number of incoming connections per neuron. Under these assumptions, the condition μ i = μ ∀i can be satisfied for symmetry reasons if we also suppose that , it is easy to verify that the parameter μ is given by does not depend on i, j anymore, the eigenvalues and eigenvectors of T , which we call, respectively, λ i and v i , and those of J , respectively, λ i and v i , are trivially related to each other as follows: therefore now the fundamental matrix Φ can be calculated in terms of the eigenquantities of T . It is important to observe that (directed) regular graphs with uniform input satisfy the assumptions above, and for this reason they will be considered from now on, even if the hypothesis of regularity is not strictly required, since we do not need to have also the same number of outgoing connections for each neuron. We also observe that even if under our assumptions the neurons have the same J c ij , I c i (and, as a consequence, also the same μ i ), this does not mean that they all behave identically. For example, from (4.11) we see that the mean of the membrane potentials is are not uniform. This proves that V i (t) depends on the index i, or in other terms that the neurons are not identical in this network.
To conclude, it is interesting to observe that if we choose A (μ) to be the algebraic activation function (see (2.2)), then Eq. (6.1) can be solved analytically, since it can easily be reduced to a fourth-order polynomial equation. Notwithstanding, in every numerical simulation of this article we will use the logistic activation function, since its properties are ideal for studying the phenomenon of stochastic synchronization introduced in Sect. 10. Now we are ready to start with the first example.

Block-Circulant Matrices with Circulant Blocks
Given two positive integers F and G, with 1 ≤ F, G ≤ N , here we suppose that the topology of the network is given by an All the entries b According to Eqs. (4.4)-(4.6), the correlation structure depends on Φ ij (t) and ij , therefore now we want to calculate the matrices Φ(t) and Φ(t)Φ (t) in terms of the eigenquantities of T . Since T is blockcirculant, its eigenvalues are those of the following matrices [44]: In a similar way, since the matrices B (i) are circulant, we can compute their eigenvalues λ (i) j as follows: Furthermore the matrix of the eigenvectors of T is where ⊗ is the Kronecker product of matrices. Since in this article we suppose that the matrix exponential that defines Φ(t) could be calculated according to (3.11), we obtain * is the element-by-element complex conjugation, and D(t) = diag(e λ 0 t , . . . , e λ N−1 t ), where λ k for k = 0, . . . , N − 1 are the eigenvalues of J (namely the collection of all the λ (i) j , as given by (6.2) in terms of the λ (i) j , for k = iG + j ). Here we have used the fact that D(t) and P are symmetric matrices (see (6.4)) and also the identity: due to the mixed-product property of the Kronecker product and to the elementary We also observe that and finally through (4.10) we obtain an explicit expression for the pairwise correlation structure of the network. It is interesting to observe that Eq. (6.6) is a consequence of the regularity of the graph. Actually, it is well known that v 0 = (1, . . . , 1) is an eigenvector of any regular graph, and that the other eigenvectors are orthogonal to v 0 , Since P is the matrix whose columns are composed of the eigenvectors of J , this means that of which Eq. (6.6) is a particular case. Now we show an explicit example of this technique, namely the case when the blocks of the matrix T have the following symmetric circulant band structure: where, supposing for simplicity that G ≥ 3, the first row of B (i) (excluding the term b (i) 0 , which is 0 for i = 0 and 1 for i > 0) can be written explicitly as The graph with this special block-circulant adjacency matrix will be represented by the notation BC F,G (M 0 , . . . , M F −1 ), and some examples are shown in Fig. 2 for different values of F and ξ . This can be considered as a toy model for describing a network of F cortical columns containing G neurons each. The parameters ξ i can be adjusted in order to generate M 0 local and M i long-range connections compatible with recent neuroanatomical studies [45], providing a rough description of a wide area of neural tissue. This idea will be extended to the case of irregular graphs in Sect. 6.3.2. Moreover, it is important to observe that even if in this case all the matrices B (i) are symmetric, the matrix T is not, since the number of connections in every block is different (the case of symmetric connectivity matrices is studied in Sect. 6.2). Now, by using Eqs. (6.3) + (6.9), we obtain   (1, 2, . . . , ξ) Ci N (1, 2), and finally K N = Ci N (1, 2, . . . , N 2 ) (complete graph) However, many different special cases can be studied. The simplest one is obtained for ξ 0 = · · · = ξ F −1 def = ξ (see the example BC 3,10 (4, 5, 5) in Fig. 2, obtained for F = 3, G = 10, and ξ 0 = ξ 1 = ξ 2 = 2), and in this case Eq. (6.10) gives: . Therefore in this case all the eigenvalues are real, as it must be, since with this special choice of the parameters the matrix T is symmetric. For F = 1 and ξ < N 2 we have M = 2ξ and Eq. (6.11) gives the eigenvalues of the circulant graph (see the example Ci N (1, 2, . . . , ξ) in Fig. 2): Now, the cyclic graph C N is obtained in the special case ξ = 1, and due to the Dirichlet kernel identity, (6.12) reduces to: Instead for ξ = G 2 (full band) and ∀F, G we have M = N − 1 and Eq. (6.11) gives the eigenvalues of the complete graph K N : By replacing Eqs. (6.10)-(6.14) in (6.7), we obtain the pairwise correlation structure of the corresponding network topology. In general there is no closed form for the finite sums in (6.7), with only the exception of the complete graph, for which we obtain Some interesting consequences of these formulas, for the complete graph and other kinds of topologies, will be analyzed in Sects. 8,9,10. However, before that, in the next section we want to show the effectiveness of this perturbative method by applying it to another class of topologies, that of symmetric connectivity matrices.

Symmetric Matrices
Another case where the matrices Φ(t) and Φ(t)Φ (t) can be computed easily is when T is a general symmetric matrix. Since its entries are real, it can be diagonalized by an orthogonal matrix P (namely a matrix such that P −1 = P ), therefore we have Φ(t) = P D(t)P . Since in this case the matrix J is symmetric, we also obtain Fig. 3 Three examples of the hypercube Q n and so Now we show an explicit example, by applying equations in (6.16) to the case when the neurons are connected according to a hypercube graph Q n . The hypercube can be defined in terms of the Cartesian product of graphs [46] (see also Sect. 6.3.1): where n is an integer and K 2 is the complete graph with 2 vertices. Some examples of Q n for different values of n are shown in Fig. 3. Clearly in this case M = n, and from (6.17) and by definition of the Cartesian product, the topology of the hypercube can be expressed as follows: It is easy to check that the eigenvalues of the matrix T Q n are n − 2m, for m = 0, . . . , n and with algebraic multiplicity n m . If we rewrite these eigenvalues with the following order: then the corresponding eigenvectors are the columns of the matrix: where H N is an N × N Hadamard matrix, defined as follows: , The reader can check that N −1 j =0 (−1) S ij k = Nδ 0k , as it must be according to (6.8), so we get We observe that Eqs. (6.7) and (6.21) are very similar. This is clearly a consequence of the regularity of the corresponding graphs.

Examples with More Complex Connectivity Matrices
Now we briefly discuss some more complex examples of connectivity. In particular, in Sect. 6.3.1 we focus on examples of more complex regular graphs, while in Sect. 6.3.2 we relax the hypothesis of regularity.

Product of Regular Graphs
In Sects. 6.1 and 6.2 we showed some relatively simple examples of regular graphs. It is possible to build more complicated topologies by means of graph operations that transform a graph into another while allowing us to calculate easily the spectrum of the new graph from that of the old one. There are two main kinds of graph operations: unary and binary. An example of unary operation is the graph complement that transforms a graph G into its complement G, namely in the graph with the same vertices of G and such that two distinct vertices of G are connected if and only if they are disconnected in G. For example, the complement of C 4 is the disjoint union of two graphs K 2 . On the other side, binary operations create a new graph from two initial graphs G, H. In this section we discuss only graph products, namely a particular kind of binary operations that prove very useful for studying networks made of different interconnected populations.

In all the examples that follow, the new graph resulting from the product of G and H has a vertex set V(G) × V(H), where × is the Cartesian product of sets and V(G), V(H) represent the collection of vertices of G, H, respectively. A well-known example that has already been introduced in Sect. 6.2 is the Cartesian product G H. This represents a new graph, where any two vertices (g, h) and (g , h ) in V(G) × V(H)
are connected if and only if either g = g and h is connected to h in H, or h = h and g is connected to g in G. Due to this rule, G H has the following topology:  (i, j ).
. This result is true for every pair of graphs that are combined through the Cartesian product. However, if G, H are regular with vertex degrees M G , M H , respectively, then G H is also regular, with degree M G + M H . This is a consequence of the fact that a graph is regular if and only if (1, . . . , 1) is an eigenvector (with the vertex degree as corresponding eigenvalue), and the fact that the tensor product v G i ⊗ v H j between two all-ones vectors v G i , v H j is itself an all-ones vector with λ G i + λ H j = M G + M H as corresponding eigenvalue. Therefore we conclude that, given graphs with known spectra, it is possible to build more complex graphs whose spectra are easily determined through the rules shown above. This proves that the theory introduced in this article can easily be used to calculate analytically the correlation structure of neural networks with highly complex connectivity matrices. Typical examples of the Cartesian product are the hypercube (see Eq. (6.17)), the circular ladder CL N = C N K 2 (also known as prism graph), the d-dimensional torus T (N 0 , . . . , N d−1 ) = C N 0 · · · C N d−1 , and so on.
Another case is the tensor product G ⊗ H, where any two vertices (g, h) and (g , h ) are connected if and only if g is connected to g in G and h is connected to h in H. From this rule we get the following topology: are their corresponding eigenvectors. Again, this result is true for any G, H, but if the graphs are both regular, then G ⊗ H is also regular, with vertex degree M G M H . An example of tensor product is the crown graph S 0 N = K N ⊗ K 2 . Now we consider the strong product G H, where any two vertices (g, h) and (g , h ) are connected whenever g and g are equal or connected in G, and h and h are equal or connected in H. So we get

From this formula it follows that the eigenvalues of G H are (λ
are their corresponding eigenvectors. Again, this result is true for any G, H, but if the graphs are both regular, then G H is also regular, with vertex degree (M G + 1)(M H + 1) − 1. A trivial example is K N G +N H = K N G K N H , from which it is possible to prove in an alternative way Eq. (6.14) by iteration.
Finally, we show the lexicographic product G • H, where any two vertices (g, h) and (g , h ) are connected if and only if either g is connected to g in G, or g = g and h is connected to h in H. Therefore the topology matrix is In general there is no simple expression for the spectrum of G • H. However, if H is regular, from (6.8) it is easy to prove that . . . , 1). If also G is regular, then G • H is regular with vertex degree M G N H + M H . An example of lexicographic product is the so called double graph of G, namely D[G] = G • K 2 [47], where K 2 is the complement of K 2 , or in other terms the graph on 2 vertices without edges.
A more complex example of the graph products introduced so far is shown in Fig. 4. This example clearly shows that the products can be used to generate easily networks with sub-populations connected in different ways, increasing the biological plausibility of the connectivity matrix. In other terms, this can be interpreted as a way to build more complex connections between the neural populations compared to the case BC F,G (M 0 , . . . , M F −1 ) studied in Sect. 6.1. We conclude by observing that it is also possible to define many other kinds of products, which are not considered here. The interested reader is referred to the literature. In general, a product between two graphs G, H can be interpreted as a system of N G neural populations with N H neurons each, interconnected in different ways according to the graph product that has been chosen

Irregular Graphs
Up to now we have studied only regular graphs, because for this class it is possible to calculate easily the eigenquantities of J from those of T by means of Eq. (6.2). In this section we show that this is not a strict requirement of our theory and that it can be applied also to irregular graphs. Regularity can be broken in two different ways: either by introducing non-uniform weights (since, by definition, regular graphs are unweighted), or by considering vertices with different (incoming or outgoing) degrees. In both cases we show how to calculate the eigenquantities of the Jacobian matrix in a relatively simple way.
First of all, in Sect. 6 we observed that Eq. (6.2) could be applied more widely also to irregular graphs with uniform weights and the same number of incoming connections, but with different outgoing connections for each neuron. In this section we generalize that result. Indeed, if for a given collection of input currents, we consider those graphs with a generally irregular topology T and a non-uniform weight ma- does not depend on the indices i, j , we can easily see that J ij = RT ij J c ij . Therefore for this class of graphs the eigenquantities of the Jacobian matrix depend trivially on those of the (unperturbed) weight matrix T • J c (here • represents the Hadamard product of matrices), which are supposed to be known. An example of such connectivity is represented by the ring model of Hansel and Sompolinsky [48], which is a well-known model for feature selectivity in primary visual cortex. In this model each cortical hypercolumn is modeled as a collection of F minicolumns with G neurons each that respond to a particular feature of the stimulus, namely the orientation of bars in the visual scene. If we introduce the function p(·), which maps each neuron i to the minicolumn it belongs to, then we call θ p(i) the preferred orientation of that neuron. In this way all the neurons in the same minicolumn have the same preferred orientation. According to experimental evidence, Hansel and Sompolinsky proposed the following connectivity matrix for the hypercolumn, where the strength of the synaptic interaction between two neurons depends on the difference between their preferred orientations: Here Γ , Δ are free parameters that define the strength of the synaptic connections within and among the minicolumns. We also observe that this formula defines a nonuniform weight matrix, therefore the corresponding graph is irregular. Now, in the primary visual cortex the preferred orientations are organized in a circular scheme around special points of the orientation map, known as pinwheels, in order to represent all the possible bar orientations in the range [0, π). For this reason, we can choose θ p(i) = ϑ + π F i G where ϑ is an arbitrary angle, so the connectivity matrix of the system can be rewritten as follows: for k = 0, . . . , F − 1. In [48] the authors also considered an external current of the form where C is the maximum amplitude of the external input, 0 ≤ ε ≤ 0.5 measures the degree of modulation of I i , and θ is the orientation for which the external input is maximum. We also observe that, if ε is small enough, ε[−1 + cos(2(θ p(i) − θ))] in the formula of I i can be interpreted as the term σ 4 I v i in Eq. (2.9), so that we can identify I c i = C ∀i. Clearly this is an extension to the irregular case of the topology studied in Sect. 6.1, which can be re-obtained for Δ = 0. It is easy to verify that M i = N − 1 ∀i (so in this case the topology T is regular, but the graph is not, due to the non-uniform weight matrix J c ), and, moreover, Finally, this neural network can be extended to the case of multiple populations with different sizes and vertex degrees (of which a very special example is the complete k-partite graph, whose topology is generally irregular). The analysis is beyond the purpose of this work and will be developed in upcoming articles.

Numerical Comparison
In this section we show that our first-order perturbative expansion is in good agreement with the real behavior of the neural network obtained from the simulation of the system (2.1). These stochastic differential equations have been solved numerically 10,000 times with the Euler-Maruyama scheme, and this collection of trials has been used to calculate the correlation by a Monte Carlo method (the code, running under Python 2.6, is available in the Supplementary Material). This result is then compared to the perturbative formula of the correlation obtained in the previous sections. The topologies that have been chosen for this comparison are C 10 , K 10 , BC 3,10 (4, 5, 5) and Q 4 (see Figs. 2 and 3), while the values of the parameters used in the numerical simulations are shown in Table 1. Moreover, the variable part of the synaptic weights and the external input currents have been chosen as follows:    Table 1 Parameters used for the numerical simulations of Figs. 5, 6, 7 and the right-hand side of Fig. 8.
For the left-hand side of Fig. 8 and for Fig. 9 the parameters are the same, with only the exception of C (0) , C (1) and C (2) , which have been set to zero Neuron Initial conditions Synaptic weights External input Logistic function We plot this comparison as a function of time (Fig. 5) and also the percentage-relative error ε% = 100 × numerical Corr − first order perturbative Corr numerical Corr as a function of the perturbative parameters (left-hand side of Fig. 6). In order to avoid high dimensional plots, we assume that σ 0 = · · · = σ 4 def = σ . Figure 5 has been obtained for σ = 0.1 and it clearly shows that the membrane potential follows very closely its numerical counterpart, while for the correlation the difference between the numerical simulation and the perturbative formula is of order 10 −2 . This is compatible with the law of large numbers, according to which the statistical error introduced by a Monte Carlo method with T trials is of order √ T . The error ε% has been calculated as a function of the perturbative parameter, for σ = 10 −3 − 1. Since we want to take into account also the error introduced by the perturbative expansion with respect to the initial conditions, whose effect quickly vanishes due to the time constant τ , the error ε% has been calculated at a small time instant, namely t = 1. The result is shown in the left-hand side of Fig. 6, which confirms the goodness of the perturbative approximation, since the error is always smaller than 3.5% if calculated over 10,000 trials. ε% could be even smaller if T is increased.
The right-hand side of Fig. 6 shows the numerical evaluation of the probability P(t) for t = 1 (see (4.13)) according to the algorithm introduced in [49]. From the figure it is easy to check that for σ = 10 −3 − 1 we obtain P ≈ 1, which further confirms the validity of our results.
To conclude, in Fig. 7 we show a comparison between the numerical and analytical probability density for both the membrane potential and the firing rate, in networks with topologies K 8 and Q 3 . Again, the parameters used in the simulations are t = 1, σ = 0.1, and those of Table 1 and Eq. (7.1). For the sake of clarity we have considered only the single-neuron marginal probability, since it facilitates the comparison. The numerical probability has been calculated by solving the system (2.1) 1,000,000 times and by applying a Monte Carlo method, while the analytical density has been evaluated by integrating Eqs. (4.11) + (4.12) over all but one dimension. The figure confirms that at the first order the neural network can be described by a normal process, even if small deviations from the normal distribution, due to the non-linearity introduced by A (V ), can be observed. Comparison between the first-order perturbative expansion and the real behavior of the network, for the topologies C 10 , K 10 , BC 3,10 (4, 5, 5) and Q 4 . The parameters used for the simulation are σ = 0.1 and those shown in Table 1 and by Eq. (7.1). Correlation has been calculated by simulating equations in (2.1) 10,000 times with the Euler-Maruyama scheme and then by applying a Monte Carlo method. Finally, this result is compared to the first-order analytical formula of the correlation. The figure shows good agreement, which validates the use of the perturbative approach

Correlation as a Function of the Strength of the Network's Input
In this section we consider how the cross-correlation among neurons depends upon a crucial network parameter, namely the strength of the external input current I c . As explained above, I c represents the external input to the network (for example, Fig. 6 Percentage-relative error of the correlation calculated between the first-order perturbative expansion and the numerical simulation of the neural network (left) and the probability P defined by (4.13) (right), for σ = 10 −3 − 1. The error is small (<3.5%) even for relatively large values of the perturbative parameter (σ ∼ 1), which proves the goodness of the perturbative approach. ε% increases considerably for σ 1, but this result has not been shown, since such values correspond to biologically unrealistic levels of randomness for a neural network. On the other hand, the figure shows that P ≈ 1, which further confirms the legitimacy of the Taylor expansion (3.2) and therefore the validity of our results. Clearly P decreases with σ because a larger variance brings the membrane potential closer to the borders defined by the radius of convergence a feed-forward input from the sensory periphery, or a top-down modulatory input) that drives or inhibits the activity of our network. Studying how the network properties depend on the parameter I c is important for many reasons. From the mathematical and theoretical point of view, this is important because this parameter may profoundly affect network dynamics. For example, the input can change the dynamical behavior of the system from a stationary to an oscillatory activity, because the eigenvalues of the Jacobian matrix (3.9) depend on μ, which in turn is determined by I c through Eq. (6.1). So changing I c can transform real eigenvalues into imaginary ones (in non-symmetric connectivity matrices) and therefore generate oscillations, or change the sign of the real part of an eigenvalue from negative to positive, giving rise to an instability. From the neural coding point of view, characterizing the dependence of different aspects of network activity upon the external input is necessary to understand and quantify how different aspects of network activity take part in the encoding of external stimuli [50][51][52][53]. Here we investigate specifically how the correlations among neurons depend on I c .
The dependence of correlation on I c is shown in Fig. 8. In this figure, the top panels show correlations for any pair of neurons in a network with a complete connectivity graph (in which case, the correlation has the same value for all pairs of neurons and so is independent of the neural indices i, j ). The bottom panels show the correlation values for a pair of directly connected neurons in a hypercube graph (in this network, the correlation value depends only on the distance between two vertices, i.e. the number of edges in a shortest path connecting them, which can range between the value of 1 which corresponds to directly connected vertices, and the maximal value of log 2 N ).
We first examined the case when the sources of variability are independent (left panels of Fig. 8), i.e. when C (0) , C (1) , and C (2) are equal to zero. Considering (3.10), it is apparent that this behavior originates from the sigmoidal shape of the activation Fig. 7 Single-neuron marginal-probability density for the membrane potential (left) and the firing rate (right) in a network with topology K 8 (top) and Q 3 (bottom). The parameters used for the simulation are t = 1, σ = 0.1, and those of Table 1 and Eq. (7.1). The numerical probability density has been calculated by simulating equations in (2.1) 1,000,000 times with the Euler-Maruyama scheme and then by applying a Monte Carlo method, while the analytical density has been evaluated by integrating Eqs. (4.11) + (4.12) over all but one dimension. From the comparison it is easy to observe that the mean and the variance of the numerical simulations are in good agreement with the corresponding analytical quantities, even if the numerical probability density is not perfectly normal, due to relatively small higher-order corrections that have been neglected in our first-order perturbative approach function: when |I c | is large, then |μ| is large as well, therefore A (μ) and the entries of the effective connectivity matrix are small. In other words, the neurons become effectively disconnected, due to the saturation of the sigmoidal activation function. An important consequence of this phenomenon is that the neurons become independent, even if the size of the network is finite. This result holds for both the complete (top-left panel) and the hypercube graph (bottom-left panel of Fig. 8). An important implication of this result is that, taking into account that ν increases with I c , in general Corr(ν i (t), ν j (t)) is not a monotonic function of the firing rate.
When the sources of variability are correlated, we found (for both network topologies; see right panels of Fig. 8) that the dependence of the correlation upon the parameter I c was very different from the case of uncorrelated sources of variability. In this case, for both considered topologies, Corr(ν i (t), ν j (t)) increases with the firing rate provided that the sources of randomness were sufficiently correlated and the network is large enough (see the case N = 32 in the right panels of Fig. 8).  Table 1, while for the independent case the parameters are the same, with only the exception of C (0) , C (1) , and C (2) , which have been set to zero. This figure shows that the correlation is strongly modulated by I c , confirming its relation with the effective connectivity matrix J eff

Failure of Sznitman's Mean-Field Theory
In this section we take advantage of our ability to study generic networks to investigate the ranges of applicability of Sznitman's mean-field theory for the mathematical analysis of a neural network. A neural network is generally described by a large set of stochastic differential equations, which makes it hard to understand the underlying behavior of the system. However, if the neurons become independent, their dynamics can be described with the mean-field theory using a highly reduced set of equations that are much simpler to analyze. For this reason the mean-field theory is a powerful tool that can be used to understand the network. One of the mechanisms through which the independence of the neurons can be obtained is the phenomenon known as propagation of chaos [19][20][21][22]. Propagation of chaos refers to the fact that, if we choose chaotic initial conditions for the membrane potentials, then any fixed number of neurons are independent ∀t > 0 in the so called thermodynamic limit, namely when the number of neurons in the system grows to infinity. Therefore the term propagation refers to the "transfer" of the chaotic distribution of the membrane potentials from t = 0 to t > 0. Under simplified assumptions as regards the nature of the network (namely that the other sources of randomness in the system, in our case the Brownian motions and the synaptic weights, are independent), propagation of chaos does occur. However, in Sects. 9.1, 9.2 and 10 we show that in many cases of practical interest, e.g. for a system with either correlated Brownian motions, initial conditions and synaptic weights, or with a sufficiently sparse connectivity matrix, or with an arbitrarily large (but still finite) size, the correlation between pairs of neurons can be high. Therefore in general any fixed number of neurons are not independent, which invalidates the use of Sznitman's mean-field theory for analyzing such networks.

Chaos Does not Occur if the Sources of Randomness Are not Independent
Here the proof is provided through a simple counterexample, namely the complete graph. From (6.15) we obtain, in the limit N → ∞: From this formula it is easy to see that if at least one of the parameters C (0) , C (1) , and C (2) is not equal to zero, then Corr(V i (t), V j (t)) = 0 (absence of chaos), even if we are in the thermodynamic limit. In particular, this means that: • if C (0) , C (2) = 0, then C (1) = 0 does not imply Corr(V i (t), V j (t)) = 0 (i.e. there is no propagation of initial chaos); • at every finite t, if C (1) = 0, then C (0) , C (2) = 0 does not imply Corr(V i (t), V j (t)) = 0 (i.e. absence of initial chaos does not lead to chaos).
Therefore Corr(V i (t), V j (t)) = 0 can be obtained only for C (0) = C (1) = C (2) = 0, which is compatible with Sznitman's mean-field theory. However, in the next section we will see that even under the last condition, namely even if all the sources of randomness are independent, propagation of chaos may not occur if the neurons are not densely connected. Clearly the fully connected network has the largest number of connections possible, for this reason it does show propagation of chaos in the thermodynamic limit. Other topologies may not satisfy this requirement.

Propagation of Chaos Does not Occur in Sufficiently Sparse Networks
Again, we show this through a counterexample. Since in this section we are interested in sparse systems, we study propagation of chaos in the thermodynamic limit as a function of the number of connections in a circulant and block-circulant network. To this purpose, we set C (0) = C (1) = C (2) = 0 (see previous section). For N → ∞ and finite M, the right-hand sides of equations in (6.7) do not converge to zero, therefore for every finite value of M propagation of chaos does not occur.
These results have been obtained by using Eq. (6.7) with C (0) = C (1) = C (2) = 0 (while all the remaining parameters are those of Table 1). The figure shows that correlation does not go to zero in the thermodynamic limit (absence of propagation of chaos) if lim N →∞ M is finite, namely if the network is sufficiently sparse However, from Fig. 9 we see that correlation decreases with M, therefore propagation of chaos occurs only in the thermodynamic limit and if M is an increasing function of N , namely if lim N →∞ M = ∞. For example, in the complete graph M = N − 1, so it explains why in this case correlation goes to zero in the thermodynamic limit. Instead in a network with a cyclic topology, propagation of chaos is never possible, also for N → ∞, since M = 2. In other words, having infinitely many neurons is not a sufficient condition for getting independence, because also infinite connections per neuron are required.

Stochastic Synchronization
Finally, we use our formalism to demonstrate a theoretically interesting regime of network dynamics. In particular, we show that for every finite and arbitrarily large number of neurons in the network, it is possible to choose special values of the parameters of the system such that, at some finite and arbitrarily large time instant, correlation is (approximately) equal to one. In other terms, the stochastic components of the membrane potentials become perfectly synchronized, therefore from now on we refer to this phenomenon as stochastic synchronization. This is a very counterintuitive behavior of the network, since it does occur even when all the sources of randomness are independent (namely C (0) = C (1) = C (2) = 0). It is important to observe that this phenomenon requires a precise tuning of the parameters of the network, which is really hard to find by chance through numerical simulations. For this reason we need a rigorous theory that tells us how to set the parameters: such a theory is developed in the next section.

The General Theory
More precisely, here we show that even when C (0) = C (1) = C (2) = 0, if the Jacobian matrix (3.9) has an eigenvalue of algebraic multiplicity one with non-negative real part, while all the other eigenvalues have negative real parts, then correlation goes to one for t → ∞, for every finite N . This is proved for a generic anatomical connectivity, therefore the assumption of regularity is relaxed. To prove this result, we suppose that J has an eigenvalue λ max with non-negative real part and with a generic algebraic multiplicity m > 0, while all the other eigenvalues have negative real parts. Now from (3.11) we recall that D(t) is the diagonal matrix of the eigenvalues of e J t , and P is the matrix of its eigenvectors. If λ max s are the first m eigenvalues of J , for t → ∞ we have D(t) ≈ diag e λ max t , . . . , e λ max t m-times , 0, . . . , 0 because all the eigenvalues have negative real part but λ max . Therefore and, moreover, According to Eqs. (4.4)-(4.6) for C (0) = C (1) = C (2) = 0, this means that where for λ max = 0 we mean e γ λmaxt −1 γ λ max = t, given γ ∈ {1, 2}. Therefore . Now, in the special case m = 1 we obtain so we conclude that lim t→∞ Corr(V i (t), V j (t)) = 1. In other terms, the neurons become perfectly correlated even if the sources of randomness are independent, which is what we wanted to prove. It is interesting to observe that, due to the Perron-Frobenius theorem [54], if the matrix with entries 1 M i J eff ij (see Eq. (3.9)) is non-negative and irreducible (namely if its corresponding directed graph is strongly connected, which means that it is possible to reach each vertex in the graph from any other vertex, by moving on the edges according to their connectivity directions), then it has a unique largest positive eigenvalue, which can be used to generate stochastic synchronization.
To conclude, it is important to observe that we must be careful when we use the perturbative expansion to describe stochastic synchronization. Actually the divergence of the term e γ λ max t implies a fast growth of the variance of the membrane potential, therefore the first-order approximation may not be good enough due to a possibly larger magnitude of the higher-order perturbative corrections. However, this problem can easily be fixed by choosing sufficiently small values of σ m that ensure the variance is still small when the correlation is close to one. Another possibility is to choose the parameters of the network in order to have λ max negative but very close to zero. For continuity, in this case correlation will be very close to one, and the variance cannot diverge since λ max < 0. Now we are ready to see an explicit example of stochastic synchronization, which will be developed in the next section for the complete and the hypercube graphs.

Examples: The Complete and the Hypercube Graphs
For both these topologies, the largest eigenvalue is λ max = − 1 τ + Γ A (μ) with algebraic multiplicity one. According to Sect. 10.1, we have to set λ max ≥ 0 in order to obtain stochastic synchronization. In particular, we consider the case λ max = 0 and we use the logistic function A (V ) = X(V ), since we can take advantage of the following property: Now, the condition λ max = 0 can be rewritten as Γ X (μ) = 1 τ , namely The solutions of this algebraic equation are where μ 1,2 are two possible stationary solutions of the membrane potential. Moreover, from Eq. (6.1) we know that Putting together Eqs. (10.1) and (10.2) we obtain Replace this value of μ 1,2 in (10.2) to obtain the final result: This non-linear algebraic equation is the constraint that must be satisfied by all the parameters of the system in order to have correlation equal to 1 in the limit t → ∞.

An example of solution of this equation is
In this case μ 1,2 = 0 and it should be used as initial condition in order to ensure the stationarity of the system. In Fig. 10 we show the phenomenon of stochastic synchronization only in the case of the complete graph (for the hypercube the results are qualitatively similar). As we can see, correlation goes to one more and more slowly if we increase the number of neurons N in the network or if we decrease the current I c . Therefore in the limit N → ∞ and/or I c → 0 the system has correlation 0 at every finite time instant. Actually, from (6.15) it is possible to prove that, given t 1, the time instant t * such that Corr(V i (t * ), V j (t * )) = C is having used the fact that Γ A (μ) = 1 τ and τ = − 2 I c . From this result we see that, for C fixed, t * increases linearly with N for large networks and is inversely proportional to I c , as obtained numerically in Fig. 10. In particular, this proves that in the thermodynamic limit there is still propagation of chaos at every finite time instant. This is in agreement with Sznitman's mean-field theory and the results on propagation of chaos proved in [20][21][22].
Moreover, from (B.1), it is interesting to observe that if there is a perfect stochastic synchronization between pairs of neurons, then it is "transmitted" to all the higher-order correlations with even order, at least for the complete graph. In other terms, if the neurons are all-to-all connected, then Corr 2 (V i (t), V j (t)) = 1 implies Corr n (V i 0 (t), . . . , V i n−1 (t)) = 1, ∀n even.

Discussion
In this article we developed a novel formalism for evaluating analytically the crosscorrelation structure of a finite-size firing-rate network with recurrent connections, using a first-order perturbative expansion of the neural equations. Importantly, the network we considered is stochastic and includes three distinct sources of randomness, namely the background noise of the membrane potentials, their initial conditions and the distribution of the recurrent synaptic weights. With this approach we succeeded in calculating analytically correlations at any order among all groups of neurons in the network. This formalism is general and in principle can be applied to networks with any kind of topology of the anatomical connections, but here we applied it to the case of regular graphs. In upcoming articles this technique will be employed to study more general kinds of anatomical connections. In other terms, the present article represents a proof of concept of the ability of our theory to relate analytically the anatomical and functional connectivity.
The cases we have decided to study are networks with block-circulant and hypercube topologies. Clearly some of the results we have obtained could be specific for these special graphs. Nevertheless, our formalism applied to these cases has shown a series of (to our knowledge) new results, whose generality or specificity can be later determined by comparison with other kinds of anatomical connections.

Dependence of the Correlation Structure on the Parameters of the System
First of all we quantified analytically how the correlation depends dynamically on the external input of the network. This has revealed a number of new and partly counterintuitive insights. We have shown that a strong input can make the neurons almost independent, and this reveals a simple mechanism to achieve network decorrelation that adds to those, such as the balance of excitation and inhibition (e.g. [27,55]) or the use of purely inhibitory feedback (e.g. [56]), that were recently proposed. Moreover, we have shown that it is not possible to obtain a mean-field description à la Sznitman of the neural network, if the anatomical connections are too sparse or our three sources of variability are correlated. We have also proved that correlation depends not only on the input, but also on the topology of the network and on the correlation structure of the sources of randomness. To conclude, we have shown that for very special values of the parameters, the neurons become almost perfectly correlated even if the sources of randomness are independent. We have called this phenomenon stochastic synchronization, and we stress the fact that the formalism developed in this article is able to prove its existence for a completely generic anatomical connectivity whose eigenvalues satisfy a bland condition.
The dependence of network correlations on the neuron's firing rates has been the subject of extensive investigations in recent years [57][58][59]. Our study of the dependence of the correlation on the strength of the external input allowed us to consider analytically this problem in our network. It is interesting to compare our results to those obtained in [57] for in-vitro real networks and for model integrate-and-fire networks. They reported that Corr(ν i (t), ν j (t)) increases with the geometric mean of the firing rates. However, in our model, this is not always the case. This happened in our case for strongly correlated inputs and relatively large networks (a scenario compatible with the cases studied in [57]). However, in our model the network showed a non-monotonic dependence of the correlation on the firing rates in other instances. A consequence of this non-monotonic dependence is that rates and correlations expressed by recurrent networks can indeed act as separate information channels for the encoding of the strength of the external stimuli. We would also like to underline the fact that, according to those authors, the correlation between the firing rates is bounded by the correlation between the inputs. According to our model, this is generally correct, but in some cases the neural network is able to generate almost perfectly correlated firing rates even if the inputs are independent. This is the phenomenon of stochastic synchronization discussed in Sect. 10.

Strengths and Weaknesses of the Presented Approach
As discussed in Sect. 1, our approach presents some advantages when compared to other methods based on linear response theory [23][24][25], networks of stochastic binary neurons [26,27], the linear noise approximation [28], the density functional approach [29], and large deviations theory [30][31][32]. These advantages consist in the possibility to use different sources of variability, to study synchronization and the effect of axonic delays, and to quantify finite-size effects also for small-size networks. This means that our formalism lends itself to the possibility of multiple generalizations and extensions. Additional sources of stochasticity, such as a random threshold V T in the activation function or a stochastic membrane time constant τ , can be introduced in the model even including correlations among different sources. As we stated above, delays in the transmission of the electric signal through the axons can be taken into account as well, following [60,61]. Another possibility of further extensions of this study is the introduction of Hebbian learning. In this article we assumed for simplicity that the dynamics of the synaptic weights is already known, through the functions (2.5). However, in the case of synaptic plasticity the time evolution of the matrix J (t) depends on the membrane potentials V (t), so the system of differential equations (2.1) should be extended to include the differential version of Hebb's learning rule. We also observe that in this article we have considered a deterministic topology T for the anatomical connectivity, which means that T is fixed from trial to trial. An interesting extension is the study of random topologies, in particular random regular graphs [62], but this problem will be tackled in another article.
A detailed analysis of the limits of our formalism for different values of all the parameters of the model and many graph topologies is beyond the purpose of the article. Nevertheless, being a perturbative approach, in general it is possible to assert that our method presents the same limits and advantages elucidated by (non-singular) perturbation theory, to which the interested readers are referred. Our formalism can be applied also to other neural equations, such as the Wilson-Cowan model [63]. However, it is important to observe that it requires the existence of a stable equilibrium point, around which the neural equations are linearized. Therefore this technique cannot be used to study the correlation structure of spiking neurons, like those described by FitzHugh-Nagumo [64,65] or the Hodgkin-Huxley [66] or integrate-and-fire [67] neurons, because in these systems spikes are generated by periodic orbits. For example, for FitzHugh-Nagumo and Hodgkin-Huxley neurons, stable periodic orbits occur around unstable equilibria, therefore our method predicts the divergence of the covariance matrix for t → ∞, which is clearly a consequence of the linearization of the neural equations. This also means that our formalism cannot be used to evaluate the correlation structure when equations (2.1) undergo neural oscillations generated through Hopf bifurcations, but can still describe damped oscillations around a stable focus in the phase space when the connectivity matrix has complex eigenvalues.
Another difficulty of our formalism is the need for an analytical expression of the eigenquantities of the Jacobian matrix J , of which we have shown a biologically relevant example in Sect. 6.3.2. Clearly spectra of brain areas that accomplish complex functions are difficult to evaluate analytically. For this reason we are forced to introduce some simplifications of the structural connectivity that we want to study. Another possibility is to determine the eigenquantities numerically, and then Eqs. (4.4)-(4.6) provide an algorithm for evaluating numerically the correlation structure of the network. Clearly even with this method the eigenquantities cannot be calculated for very large networks, since the matrix J is N × N and therefore grows quickly with the network size. However, the advantage of evaluating numerically Eqs. (4.4)-(4.6) is evident, compared to the Monte Carlo approach. Actually, if the randomness of the synaptic weights is taken into account (namely if σ 2 = 0), one needs to generate numerically by a random generator the N 2 entries of the matrix W , according to the covariance matrix (2.6), which has N 4 entries. This calculation must be repeated for a sufficiently high number of trials, according to the Monte Carlo method, so it is computationally much more expensive in terms of time and memory consumption.
It is important to observe that in this article we focused mainly on regular graphs for the sake of clarity, since for this class of connectivity matrices the eigenquantities of J can be evaluated easily from those of T • J c through Eq. (6.2). For a general connectivity this relation is harder to find, but we underline that this is in part due to our choice to use a biologically realistic activation function A (·) (see Eqs. (3.9) and (3.10)). Usually, in order to obtain analytical results, in the literature there is a wide use of piecewise linear activation functions (e.g. in [48,68,69]). Clearly in this case it is much easier to evaluate the eigenquantities of J from those of T • J c , taking some care at the connection points between the segments of A (·), where the piecewise linear function is not differentiable. Another useful feature of our approach is that it allowed the calculation of the dependence on the strength of the external input of correlations of arbitrary order (not only pairwise correlations). This feature will be useful for the evaluation of the ability of networks to encode genuinely additional information in the variations with inputs of higher-order correlations, a subject that has been under intense theoretical [70] and experimental debate in recent years [71,72].

Analyzing the Consequences of Structural Damage
Similarly to spectral graph theory, where the properties of a graph are studied in relationship to its characteristic polynomial and eigenquantities, in this article we have found the relation between the functional connectivity and the spectrum of the underlying structural connectivity. This, in principle, allows one to study the effect on the functional connectivity caused by lesions to the synaptic connections. These structural damages can be modeled as perturbations to the topology matrix. Thus, in principle they can be studied by perturbative techniques such as those described in [73][74][75][76][77]. This branch of graph theory deals with discrete perturbations (such as the removal of connections or vertices from a given graph), as opposed to the Rayleigh-Schrödinger theory from quantum mechanics, that studies the effect of continuous perturbations to the generalized eigenvalue problem. This approach would help to understand abnormal functional behavior, complementing other studies of the consequences of structural damage, e.g. [78].

Possible Extensions to Other Measures of Communication Among Neurons
It is also interesting to observe that the correlation structure can be used to estimate causal relations between neurons or neural populations. This can be achieved in many ways. However, in our view a promising direction is to take advantage of hierarchical clustering techniques already used in economics, whose potential application is briefly described as follows. According to [79], the correlation structure can be used to define a distance measure d ij (t) def = 2(1 − Corr(V i (t), V j (t))) between every pair of neurons. Clearly we are not interested in the hierarchical structure of single neurons, but rather in that of mesoscopic or macroscopic areas. For this reason, from d ij (t) we have to define an arbitrary distance between these areas of the brain (e.g. the mean distance between all the pairs of neurons). Then, from the distance matrix of the areas, we can determine the minimum spanning tree of the system, a concept introduced in the context of graph theory to find the most relevant (or more informative) connections in a network. Finally, on the minimum spanning tree it is possible to define an ultrametric distance, which in turn allows us to build a dendrogram (i.e. a hierarchical tree) in an unambiguous way, by using techniques such as UPGMA [80].

Concluding Statement
We have shown that the formalism introduced in this article can be effectively used to calculate the functional connectivity of neurons within a firing-rate network model. In this article we concentrated mostly on computing the Pearson correlation among all pairs of neurons in the network. However, the work reported in this paper also lays the basis for computing more refined measures of functional connectivity (such as those based on information theory). This in turn will allow in future studies the analytical quantification of the transmission of information among the elements of this recurrent network and of how information transmission is modulated by factors such as the strength and dynamics of external inputs. Fig. 11 Radius of convergence r of the Taylor series of the logistic function X(V ), as a function of the point V = μ about which the expansion is performed. For large μ the radius of converge increases linearly since the logistic function is asymptotically flat. Instead for Λ → ∞ we obtain r(0) → 0, because in this limit the logistic function becomes a Heaviside step function with a discontinuity in V = 0
Therefore the radius of convergence increases with μ, as it must be. Moreover, in the limit Λ → ∞ it gives r(μ) = |μ|, as with the logistic function. The same result can be proved for other sigmoidal functions and is left as an exercise for the reader.