{'fraction_non_alphanumeric': 0.051708664135312395, 'fraction_numerical': 0.01318605453917846, 'mean_word_length': 4.216240547353259, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytic Investigation for Spatio-temporal Patterns Propagation in Spiking Neural Networks

Based upon the moment closure approach, a Gaussian random field is constructed to quantitatively and analytically characterize the dynamics of a random point field. The approach provides us with a theoretical tool to investigate synchronized spike propagation in a feedforward or recurrent spiking neural network. We show that the balance between the excitation and inhibition postsynaptic potentials is required for the occurrence of synfire chains. In particular, with a balanced network, the critical packet size of invasion and annihilation is observed. We also derive a sufficient analytic condition for the synchronization propagation in an asynchronous environment, which further allows us to disclose the possibility of spatial synaptic structure to sustain a stable synfire chain. Our findings are in good agreement with simulations and help us understand the propagation of spatio-temporal patterns in a random point flied.


Introduction
Towards fully understanding an evolutionary random point field, the (joint) probability distribution density is often very hard, if it is not impossible, to be calculated and estimated. However, in many cases, the first few moments are sufficient to present a holistic picture of the evolution of this random point field [1,2,3,4]. To this end, the moment closure approach was proposed to approximate random fields by employing their first few moments [5,6,7,8], which transforms the spatio-temporal random point field into a few non-random dynamical systems of the moments respectively. Then, this idea enables us to investigate these moment dynamical systems for the depiction of the asymptotical properties of the random point field. For instance, the spatio-temporal pattern of random point field can be regarded as attracting dynamics of the corresponding moment dynamical system.
Mathematically, neural activities can be formulated as random point processes (spike trains), and a large ensemble of neural spiking trains with an underlying geometric structure (manifold) composes of a random point field [9]. In particular, the propagation of spatio-temporal patterns of neural spiking activities, as the essence of cortical function, has been widely studied in recent decades [10,11,12,13,14]. A typical example of spatio-temporal pattern, the so-called pulse-packet, to encode a piece of information and reliably transmit it from one layer of nervous system to the other is a synfire chain: at each layer, spikes are synchronized inside a packet of neurons but asynchronous out of the packet [15,16,17,18] (see Fig. 1 for illustration). There are large numbers of studies showing that synfire chain can transmit information between neuronal populations experimentally and theoretically [19,20,21]. As an information transmitter, it has been successfully applied in many computational tasks in recent decades [22,23,24,25]. This sort of coexistence of synchronous and asynchronous states of neural activities can also be categorised as "chimera" phenomenon [26,27,28].
The synfire chain can be naturally described in the framework of random point fields. Traditionally, dynamical models of mean firing rate [29,30] were used to depict the evolution of spike spatio-temporal patterns [31,32,33,34,35] in cortical circuits. However, such models are only true under the condition of independent or weak correlation of spiking activities between neurons [36,37] and they have limitations to fully account for the evolution of the spiking patterns [38,39,40,41]. On the other hand, temporal dispersion in terms of pulse-packets [42] and correlation map [43,44] represented the degree of synchrony well but failed to realistically describe the pattern in a random point field. How the survival of a pulse packet in a multi-layered neural network depends on the structural and physiological characteristics of the network has not been answered analytically in the literature.
Under the framework of moment closure, Pearson correlation coefficients were utilized to describe synchrony [44,45]. In the present paper, we developed a novel and general theoretical framework to investigate the synchrony dynamics of random point field of multilayer feedforward neural network (FNN) with spiking leaky integrate-and-fire neurons. Besides, we also find similar dynamical features in discrete-time recurrent spiking neural network using our model. By constructing the moment dynamical system that includes the first (mean firing rate) and second-order statistics (variance and correlation) of spiking random point fields, the synfire chain dynamics can be regarded as the existence of certain sort of attracting set of the correlation map. Combining with the mean field approach, we discover the necessity of the balance of network for the stable synchronization propagation. In addition, we analytically derive a sufficient condition for the existence of this attracting set that enables to obtain the appropriate size of this synfire chain and the proper synpatic density of the neural network. These results are in a good agreement with the numerical results and help understand the evolution of synchronisation pattern in spiking neural network.
For each ξ k (x, t), the number of events occurring in an infinitesimal interval of time, [t, t + dt], denoted by dξ k (x, t), follows a probability distribution, P k (ξ; x, dt). Inspired by the moment closure approach, with assuming that ξ(x, t) is a renewal process, which disregards the correlation with time, we are in the stage to approximate P k (ξ; x, dt) by the first moment, mean, the second statistics, variation, and correlation. Thus, we approximate dξ k (x, t) by a Gaussian field as follows: where µ k (x, t) and σ k (x, t) are the functions of mean and variance, η k (x, t) is a Gaussian field with zero means and unit variances: E(η k (x, t)) = 0, var(dη k (x, t)) = dt, but correlated. Here, E(·) stands for the expectation. Equivalently, we can use the following coefficient of variation (CV) to stand for the second moment: The Pearson correlation coefficient (CC) between locations x and y at t: Assume that dξ k (x, t) is asymptotically stationary. Then, with a sufficiently large ∆, the long-term duration behaviours of µ k (x, t), σ k (x, t) and ρ k (x, y, t) can be represented by the asymptotics, namely, µ k (x), σ k (x) and ρ k (x, y). Thus, (1) can be rewritten as: Here, the correlation between dη k (x, t) and dη k (y, t) is ρ k (x, y)dt. In this stationary case, CC can be equivalently defined in the following way. The shift-correlation between ξ k (x, t) and ξ k (y, t) in a sliding window with length T is defined as: where n k (x, τ ) is the number of events occurring at location x in the time interval [τ, τ + T ], cov(·, ·) stands for the covariance and the covariance (cov) and variance (var) are calculated with respect to τ . Thus, the correlation coefficient between this spike train pair is defined as its limit ρ k (x, y) = lim T →∞ ρ k T (x, y), which equals formula (3) in stationary state [44].

Moment clourse method
The above analysis constructs a random point field system (4)-(5) to study the dynamical characteristics of spiking neural network. Discribing spike activities with probability distribution function directly seems intractable, so it is more efficient to focus on their moments of distribution. Specifically, suppose Λ is a finite index set and M m i : i ∈ Λ represents the moments we need. Then, we can transfer a stochastic system into moment system, m i : i ∈ Λ = F ( m i : i ∈ Λ ), using moment closure methods [8], where F (·) represents corresponding moment map. Different from the famous model proposed by Wilson and Cowan [46] which only includes the first order moment (firing rate), we hope to find a more complex model based on higher order, such as the second order moments (variance, correlation). Gaussian part in the right-hand side of (4) can naturally provide the informatin of the second order moments and yield a moment closure.
We regard the long-term dynamics in the time bin [(k−1)∆, k∆], ξ k−1 (x, t), as the input to that of the successive time bin, [k∆, (k + 1)∆], ξ k (x, t). We aim to formulate the iteration of these first and second moments/statistic in the form of where cov k−1 (x, y) = ρ k−1 (x, y)σ k−1 (x)σ k−1 (y) is the covariance of ξ k−1 (·, t) between the location x and y. Here, S x 1,2 stand for the functions of mean and variance with respect to the last moment mean and covariance, named mean and variance maps respectively. Ψ x,y stands for the correlation function, named correlation map. The details of S x 1,2 and Ψ x,y will be introduced in Section 2.3. Thus, the set of mean, variance and correlation maps compose of the moment dynamic system (6) that represents the random point field ξ(x, τ ) in the moment closure fashion.
The problem of spatio-temporal pattern of ξ(x, τ ), which can be described by the first and second moments/statistics, can be transformed into the dynamics of (6). For instance, synchronization pattern of a random point process can be described by the CC (3). ρ k (x, y) = 1 implies the processes at x and y completely synchronise; ρ k (x, y) = 0 means uncorrelated and ρ k (x, y) = −1 means complete desynchronization. The asymptotic dynamics, namely, asymptotic attractor, can depict the pattern related to synchronization. If ξ(x, t) is asymptotically steady, we can equivalently study the equilibrium functions: µ k−1 (x) = µ k (x), σ k−1 (x) = σ k (x) and ρ k−1 (x, y) = ρ k (x, y), towards understanding the asymptotic pattern of the random point process. We highlight that this approach can be extended to include the high-order moments and correlations of the random point field.

Multilayer feedforward spiking network
As an application, let us specify a random point field described by multilayer FNN of spiking neurons with a sparse random coupling structure [18,47]. At each layer, there are exactly N = 5n neurons. Among them, there are N E = 4n excitatory (E-) neurons and N I = n inhibitory (I-) neurons. For the E-neuron group at the kth layer, a synfire chain is carried by a packet of E-neurons, denoted by W k with the identical size #{W k } = W for all k, as the packet size, where #{·} is the number of elements in a finite set. Between any two successive layers, for instance, from the (k − 1)th to the kth layer, each neuron in W k receives inputs from all neurons in W k−1 . This constructs a completely connected FNN in the packet, which is believed as the basic model for synfire chain [47]. This completely connected FNN is embedded into a sparse network [18] and the couplings for other neurons outside W k are randomly picked with an equal probability so that the total number of excitation and inhibition synaptic links are K E = λN E and K I = λN I respectively (see Fig. 1). Here, λ ∈ (0, 1) is the indicator of the sparsity of the synaptic density. To ensure that every neuron receives an equal number of synapses, W K E is necessary.
We use a discretized random point field N k i (t) to stand for the spike counts from the neuron i at the kth layer, i.e., N k is the Kronecker-delta function and {t i,k l } are the time points of the presynaptic spikes from neuron i (at the kth layer). Here, the neuron label i stands for the spatial variable, t for the continuous-time variable at kth layer.
The evolution of the random point field follows the integration-and-fire model. The potential activity of each neuron i at layer k is described as: for i = 1, · · · , N . Here, τ m is the capacitance constant, I k ext,i is the external current stimulus at layer k, and I k syn,i is the synaptic stimulus from the neurons at the (k − 1)th layer: with w k ij standing for the strength of the excitatory post-synaptic potential (EPSP) or inhibitory post-synaptic potential (IPSP) from neuron j at the (k − 1)th layer to neuron i at the kth layer. In this paper, we take values of EPSP/IPSP as follows. If there is a synaptic link from neuron j at the (k − 1)th layer to neuron i at the kth layer, then we set w k ij = w 0 if j is an E-neuron and w k ij = −rgw 0 if j is an I-neuron, for some constant w 0 > 0, where g = 4 equalizes the ratio between the numbers of E-synapses over I-synapses and thus r serves as the ratio of ISPS over ESPS. If there is no link from j to i, w ij = 0.
Once V k i (t) reaches a threshold (V th ), neuron i is depolarized and emits a spike. Then, V k i (t) is reset to V r after a period of refractory time T ref . By this way, neuron i at layer k emits a spike train, which is the input to the neurons at the (k + 1)th layer that are linked with neuron i. Therefore, the random point field of spiking trains, N k i (t), can be generated. By the approach mentioned above, we are to approximate the spiking train of neuron j (a random point process) by a Gaussian process: , j = 1, · · · , N are correlated Brownian motions [45,48]. µ k−1 j and σ k−1 j are the mean and variance of neuron j (at the (k − 1)th layer), which can be derived by the renew process theory [49]. Here, τ m stands for the time-scale constant. According to (8), (7) becomes correlated Ornstein-Uhlenbeck (OU) processes: are the mean and variance of the sum of the post-synapses received by neuron i respectively, where ρ k−1 jl is the correlation coefficient between the j-th and l-th synapses of neuron i.
The output spike trains derived from (9) (with potential threshold V th ) are also approximated as Gaussian processes. This establishes maps of the first (mean map), second-order statistics (variance and correlation maps) of the random point field N k i (t) of the successive layers. By the Siegert's expression [50], these three maps can be formulated: where S 1 and S 2 are the mean and variance maps. Specifically, they can be written as [45,48]: Here T ref represents the refractory period, V r is the rest potential and V th is the threshold of membrane potential to emit a spike. Besides, where D − (u) is exactly the Dawson's integral. For the details, we refer the readers to Ref. [45,48]. As pointed out in Ref. [44], the evolution of correlation coefficient from the (k − 1)th layer to the kth is formulated as: whereρ k i,j is the correlation coefficient between the collection of synpatic inputs of neuron i and j at layer k: Φ is the correlation map, which is assumed to be the identity map, namely, Φ(ρ) = ρ in this paper, according to the arguments in Ref. [45]. However, the following results still hold when Φ(µ) is monotonous increasing with respect to µ as shown in Ref. [44,48] with some minor modifications.
To sum up, the spike trains of neurons in the FNN are modeled by the discrete random point field N k i (t). Their evolution through layers is formulated by the iteration of the first and second moments/statistics, which is completely described by the mean, variance, correlation maps and called the moment neural network (MNN) [45].
In particular, evolution equations of correlation map through a multilayer FNN can be utilized to analytically and quantitatively study the synchronization propagation. The (stable) synfire chain is defined by an attractor of the correlation map (13) together with the mean and variance maps (10) with two properties: (P 1 ) the correlation coefficients between the neurons in the packet of synchrony are large; (P 2 ) the correlation coefficients between the other pairs of neurons are low.
Hence to characterize the synfire chain, we focus on the correlation map, which turns out to be a useful way to understand the complex dynamics. Let us use the subscript + to denote all variables in the packet W k and − out the packet at each layer. Then, for given two small constants 0 < < 1 and 0 < δ < 1, define Λ ,δ = {(ρ + , ρ − ) : ρ + > 1 − , |ρ − | < δ}. If Λ ,δ is an attracting set of the feedforward network of LIF neurons or the theoretical model (13), then we call it synfire attractor and ρ + − ρ − the synfire gap. Thus, a synfire chain can be defined as: The multilayer FNN is said to possess a synfire chain (with respect to and δ) if Λ ,δ is an attracting set of the iteration map (13) accompanied with (10).

Spatial mean-field approximation
With the mean field approximation, the correlation map allows to provide analytic inference and insights on the stability of synfire chain by some algebras, and so greatly simplifies the analysis of the transmission of spiking point process in multilayer FNN. According to whether the link or neuron belongs to the packet or not, we substitute the specific correlation between neurons p and q, ρ k pq , and the variance of neuron p, σ k p , by ρ k ± and σ k ± . In detail, the mean in-packet can be calculated through (13)- (14) as: Here · S represents the average calculation over the set S, noting that the correlation is assumed to be the identity map in this paper. Similarly, the out-of-packet correlation coefficients are approximately written as: where ≈ 0, when N >> W , are two weights corresponding to the proportions of neuron pairs with all outof-packet and neuron pairs with one in-packet and the other out-of-packet respectively. Here, Through a squared decomposition in Fig. 2, we have the following expression of A k+1 ± and B k+1 ± in mean field form. See Appendix A for the detailed derivation.   Figure 2: The graphical computation of p,q w k+1 ip σ k p w k+1 jq σ k q ρ k pq in (14). We lined up the subscripts p = 1, 2, . . . , W, . . . , K E , . . . , K E +K I and q = 1, 2, . . . , W, . . . , K E , . . . , K E +K I in two mutually perpendicular directions to form two squares for (A) (i, j) ∈ W k+1 ×W k+1 and (B) (i, j) ∈ W k+1 × W k+1 . Both of the squares can be decomposed into 3 diagonal intervals and 7 blocks according to the value of ρ k pq , the location of p and q(in the packet or not) and the attribute of p and q(excitatory or inhibitory). The value of number of neuron pairs, variance weight, correlation coefficient (CC) on each region and the expression of A k+1 ± and B k+1 ± are listed in the right tables.

Non-constant coefficient of variation
To our best knowledge, in most existing literature, the spike trains of LIF neurons were simulated by Poission or sub-/supra-Poisson processes through FNNs, and the coefficient of variation (CV) of the spike trains were assumed constant [18,51]. Thus, the first order moment (mean) can be utilized to depict the second-order dynamics of synfire propagation. However, this assumption could be away from the facts. To justify it, we simulate 2000 LIF neurons with initial Poisson input in synfire propagation to see the evolution of CV, which was defined in (2). As shown in Fig. 3(A), CV is non-stationary for layers. Moreover, even if the CVs go steady after many layers, they vary for parameters. For instance, there is a moderate rise of CVs with increasing cluster size W or decreasing the proportion of ISPS and ESPS, denoted by r, as shown in Fig. 3(B). This motivates us to employ Gaussian process to approximate the spike trains, instead of assuming fixed CV.

Balanced network
The parameter r measures the ratio of ISPS over ESPS. r = 1 means that the excitatory and inhibitory are exactly balanced. Consider the correlation map in the mean field sense, (15)- (16) and their expressions in detail. Noting K E = λN E with N E = 4n = 4 5 N , the terms of the highest order in the denominator, B k+1 − , and the numerator, So, if r = 1, the identical term of the highest order in both denominator and numerator implies that ρ k+1 as k goes to infinity when N is sufficiently large. Therefore, with a large size, r = 1 is the critical value of the I-ESPS ratio for the existence of synfire chain in terms of property P 2 . As illustrated in Fig. 4(A), one can disclose that Λ ,δ exists with some small values of and δ only in a balanced network, namely, r ≈ 1, which implies the largest synfire gap (ρ + − ρ − > 0.7). This finding can be verified by simulating LIF neuron in a feedforward network, as shown in Fig.4(A) as well. However, when r > 1, the out-of-packet mean firing rate decreases quickly so that the out-of-packet spike frequency almost disappears (lower than 1 Hz) when r > 1.2. Therefore, to maintain a stable pattern of synchronization, the balanced network, namely, r = 1, is necessary. . Simulation results of the propagation of µ ± , σ ± and ρ ± (the embedded subfigure) through layers by LIF network (7). At the initial layer, N Possionian spiking trains are generated in a time interval of 20 sec with the firing rate of 20Hz, among which W trains in the packet are strongly correlated (ρ ≈ 1) that represent and the rest out-of-packet N − W trains are weakly correlated (ρ < 0.2).

Synfire chain condition and packet sizes
We discuss the packet size in the synfire chain: Λ ,δ with = δ = 0.3 as an attracting set of the correlation map (14) in a balanced network [52,53]. As illustrated in Fig. 5(A), to maintain Λ ,δ , the packet of synfire chain should have an appropriate size. A large W will enhance synchronization between neurons not only in-packet but also out-of-packet, and thus the in-packet synchrony will invade over the whole network, which leads the synfire gap to disappear due to the large packet size. That is, property P 2 fails to hold (invasion). There exists a corresponding threshold value of W by which synchronization invading occurs, named invasion packet size and denoted by W i . In comparison, a very small W may depress synchronization propagation in the packet and at the same time destruct the maintaining of the mean firing rates of the whole network, and thus the in-packet synfire chain will annihilate, which also leads the synfire gap disappears due to the small packet size. That is, property P 1 cannot hold (annihilation). There also exists a corresponding threshold value of W , by which synchronization annihilation occurs, named annihilation packet size and denoted by W a . Therefore, an ideal W to maintain a stable synfire chain is medium, belonging to the interval (W a , W i ), which leads both synchronous synfire packet and an asynchronous background with regular firing rates. As illustrated by Fig. 5(B), the simulation result by LIF neurons in the feedforward network has a good qualitative agreement with the theoretical results.
With the methods introduced in Section 2, we can derive an analytic result (a sufficient condition) with respect to , δ and λ for the existence of attractor Λ ,δ in a balanced network (r = 1). We consider the asymptotic stationary state of system (15)- (16), in which all the mean-field variables such as σ t ± , ρ t ± , A t ± and B t ± have reached their steady states as shown in Fig. 4(B), whose values are denoted by σ ± , ρ ± , A ± and B ± respectively. For simplicity, let With r = 1, according to (15)-(21), by some algebras, we have and Thus, it can be seen that if and λP + λ 2 QW 2 + λRW P + λ 2 QW 2 + λRW < δ, hold for some W , then the existence of Λ ,δ can be guaranteed in the meanfield sense. By some algebras, one can derive that if holds, named synfire condition, then (24)-(25) hold. We are to identify the values of W that satisfied (24)-(25) under the synfire condition. In fact, this synfire condition (26) is equivalent to 1−λ −1 < δ−λ λ(1−δ) . Thus, for any τ ∈ ( 1−λ − 1, δ−λ λ(1−δ) ), picking W as the solution of equation λ 2 QS 2 +λRS λP = τ with respect to S, one can easily verify that (24)- (25) hold. Hence, we can derive a region of synfire size: For each W ∈ Reg, Λ ,δ is an attracting set of the model (13) in the meanfield sense. In addition, max Reg and min Reg give the lower bound of the invasion packet size and the upper bound of the annihilation packet size respectively.

Sparse synaptic density
Furthermore, the synfire condition (26) reveals the dependence of synfire propagation on the synaptic density λ. The synfire condition (26) and the region of synfire chain packet (27) imply that a smaller λ leads a larger interval of synfire region in terms of the existence of Λ ,δ , which may cause a larger value of the maximum synfire gap max W (ρ + − ρ − ). However, since the mean map in (10) can be written as µ k+1 in a balanced network (r = 1), the mean out-of-packet firing rate will converge to 0 through layers, becauseμ k+1 − → 0 as λ → 0. That is to say, very small λ results in the mean out-of-packet firing rate going extremely low, even near zero. As shown in Fig. 6, the maximum synfire gap max W ρ + − ρ − decreases with λ but an extremely small λ (less than 0.075) makes the spiking activities out of the packet vanish (< 0.3 Hz). Therefore, to sustain a synfire chain and the asynchronous out-of-packet spiking activity , λ should be taken a "modestly" small value. For instance, λ = 0.1 was suggested in Refs. [47,54].  Figure 6: Barplots of the maximum synfire gap max W ρ + −ρ − (blue) and the corresponding mean out-of-packet firing rate µ − (red) with respect to linking density λ in balanced LIF network.

Synfire chain in recurrent neural network
Besides FNN, we claim that above analysis still works in recurrent neural network (RNN) [48]. If we consider the dynamics of a single layer of neurons in discrete time (see Fig. 7(A)), then (7) can be written as: where s denotes time index. Similar to the Section 2.2, we iterate moment maps and get the moment closure: where the mapping functions are still defined as (11) and (13). Therefore FNN and RNN share almost the same evolution equations in our framework. For simplicity to compare, we denote time index t in RNN as k in the following paragraphs.
As shown in Fig. 7(A), to analyze RNN, firstly, we should unfold the network structure in discrete time. Then, RNN can be treated as FNN with constrained weight [55], which means packets W k (k = 1, · · · ) share the same position and other random sparse connects are also fixed in different layers (weight sharing). If other parameters (e.g., I-ESPS proportion r, synaptic density λ, packet size W ) take the same value as above, the dynamics of RNN and FNN are almost identical since we only focus on the forward propagation and do not consider the backpropagation like supervised learning. We denote the symbols using in FNN (Section 2.3) as A k ± , B k ± , ρ k ± in RNN so that (17) turns to: Assuming that the input spikes in RNN and FNN are identical, with the mean field approximation, it is easy to get: which leads to ρ k ± = ρ k ± . From Fig. 7(B), we can find that both networks have almost the same attractors. Furthermore, other dynamical characteristics referred in Section 3.1-3.4. can also be varified in the same way. The mean in-out packet correlation of recurrent ( ρ ± ) and forward (ρ ± ) LIF networks with the same Poission input (20Hz, 20 seconds) and other parameters (W = 50, r = 1) in ten runs.

Conclusion
We have developed a theoretical framework of Gaussian random field to study how synchronization pattern propagates in feedforward or recurrent spiking neural network. Combined with the mean field approach, we analytically proved that the balance network is necessary for the stability of synfire chain in terms of suppressing the outside correlation and maintaining the in-packet spiking activities. And we derived a sufficient condition for a stable synfire chain by providing an estimation of the appropriate packet size region which avoids both invasion and annihilation of synchronization, and revealing the role of spatial synaptic structure to sustain this spiking spatio-temporal pattern. Our analytic results show good agreements with the simulations of LIF network. We highlight that this approach based on moment closure is powerful and general to investigate propagation and stability of spatio-temporal patterns in random point field. One step further, we have not included the complex interactions like brain manifold structure in the random field here, as developed in Ref. [56], which is our future work and such a theoretical framework should be a valuable tool for investigating the varied dynamics of the spike patterns observed in experiments.
through the same index set of length W . In region b, p and q go through different index set of length K E − W , thus p connects q for K E − W times if i = j otherwise λ(K E − W ) times, in which λ describes the connection probability out of the cluster. In region c, p and q go through different index set of length K I , thus p connects q for K I times if i = j otherwise λK I times. As for the variance weight, in region a, both p and q are in the cluster and excitatory, thus the variance weight is (σ k + ) 2 ; in region b, both p and q are out of the cluster and excitatory, thus the variance weight is (σ k − ) 2 ; in region c, both p and q are out of the cluster and inhibitory, thus the variance weight is (−4rσ k − ) 2 . Therefore, the contribution of these region to A k+1 To sum up the components in each region, we have