Effects of Firing Variability on Network Structures with Spike-Timing-Dependent Plasticity

Synaptic plasticity is believed to be the biological substrate underlying learning and memory. One of the most widespread forms of synaptic plasticity, spike-timing-dependent plasticity (STDP), uses the spike timing information of presynaptic and postsynaptic neurons to induce synaptic potentiation or depression. An open question is how STDP organizes the connectivity patterns in neuronal circuits. Previous studies have placed much emphasis on the role of firing rate in shaping connectivity patterns. Here, we go beyond the firing rate description to develop a self-consistent linear response theory that incorporates the information of both firing rate and firing variability. By decomposing the pairwise spike correlation into one component associated with local direct connections and the other associated with indirect connections, we identify two distinct regimes regarding the network structures learned through STDP. In one regime, the contribution of the direct-connection correlations dominates over that of the indirect-connection correlations in the learning dynamics; this gives rise to a network structure consistent with the firing rate description. In the other regime, the contribution of the indirect-connection correlations dominates in the learning dynamics, leading to a network structure different from the firing rate description. We demonstrate that the heterogeneity of firing variability across neuronal populations induces a temporally asymmetric structure of indirect-connection correlations. This temporally asymmetric structure underlies the emergence of the second regime. Our study provides a new perspective that emphasizes the role of high-order statistics of spiking activity in the spike-correlation-sensitive learning dynamics.


INTRODUCTION
Spike-timing-dependent plasticity (STDP) (Markram et al., 1997;Bi and Poo, 1998;Caporale and Dan, 2008) is one of the major forms of Hebbian learning. In canonical STDP, a synapse is potentiated (depressed) when the presynaptic spikes precede (follow) the postsynaptic ones. With this temporally asymmetric learning window, STDP can regulate both the mean rate and variability of postsynaptic firing (Song et al., 2000), enforce the competition between convergent synaptic inputs (Song et al., 2000), and support temporal sequence learning (Rao and Sejnowski, 2001).
While STDP of a single neuron with many synapses has been well characterized, issues related to the STDP learning dynamics in recurrent networks are yet to be fully addressed, in spite of recent progress (Morrison et al., 2007;Clopath et al., 2010;Gilson, 2010;Kozloski and Cecchi, 2010;Abbott, 2013, 2016;Litwin-Kumar and Doiron, 2014;Ocker et al., 2015;Zenke et al., 2015;Bi and Zhou, 2016;Ravid Tannenbaum and Burak, 2016). The role of the loworder statistical structure (i.e., firing rate) of spiking activity has been intensively investigated in the learning dynamics in recurrent networks (Morrison et al., 2007;Kozloski and Cecchi, 2010;Babadi and Abbott, 2013). Though it is known that, as a correlation-sensitive learning rule, STDP inherently incorporates the high-order statistics of spiking activity (Morrison et al., 2008), the role of the high-order statistics of spiking activity in STDP learning dynamics for recurrent networks remains largely unknown. In this work, we explore the issue of how the high-order firing statistics (for example, firing variability) affect the network connectivity structures learned through STDP in recurrent networks. We demonstrate that the information of firing variability is encoded in a linear response kernel that quantifies how a neuron responds to a weak incoming signal. By developing a self-consistent linear response theory pioneered in Pernice et al. (2012), Trousdale et al. (2012), and Helias et al. (2013), we show that this response kernel, together with the firing rate, determines the pairwise correlation of spiking activity. By decomposing the pairwise correlation into one component associated with local direct connections and the other associated with indirect connections, we can identify two distinct dynamical regimes. In the first regime, the local direct connections dominate over the indirect connections in determining the learned network structure. This regime corresponds to those studied in Babadi and Abbott (2013), Kozloski and Cecchi (2010), and Morrison et al. (2007), where the connections from neurons of higher firing rate to those of lower firing rate are strengthened whereas the connections with the opposite situation are weakened. The other is a regime in which the indirect connections dominate over the direct connections in learning dynamics. In this second regime, the connections from neurons of higher firing rate and higher firing variability to those of lower firing rate and lower firing variability are weakened whereas the connections with the opposite situation (i.e., connections from neurons of lower firing rate and lower firing variability to those of higher firing rate and higher firing variability) are strengthened. We find that there is a temporally asymmetric structure of commoninput-induced spike correlations between different neuronal populations. This structure is induced by the heterogeneity of firing variability across neuronal populations, which underlies this counterintuitive dynamical regime. Our study provides a theoretical framework for addressing the question of how the high-order statistics of spiking activity affect network structures learned through STDP. Using this framework, we demonstrate how the introduction of heterogeneity of firing variability across different neuronal populations can give rise to counterintuitive network connectivity structures in a plastic neuronal network.

Model
The model we used here is the current-based leaky integrate-andfire (LIF) model with finite synaptic time constant. It is the same as that in Babadi and Abbott (2013).
where τ m = 20 ms is the membrane time constant, V r = −60 mV is the resting potential, τ s = 5 ms is the synaptic time constant, w ij is the synaptic strength from neuron j to neuron i and s j (t) is the spike train of neuron j. The parameters µ ext,i and σ 2 ext,i are the mean and variance of the external input, respectively, and ξ i is the white noise satisfying ξ i = 0 and ξ i (t)ξ j (t ′ ) = δ ij δ(t − t ′ ). Every time when the membrane potential V i crosses the threshold V th = −40 mV, neuron i would emit a spike and V i would be reset to V r .
The network consists of N E excitatory neurons and N I inhibitory neurons. We use N E = 250, N I = 250. The connections between these neurons are of the all-to-all type. We keep the inhibitory-to-excitatory (EI), excitatory-to-inhibitory (IE) and inhibitory-to-inhibitory (II) synaptic strengths constant during the entire simulation while the excitatory-to-excitatory (EE) connections are subject to the STDP learning rule to be introduced in section "The effect of firing rate and firing variability on learning dynamics-numerical evidence". At the beginning of simulations, the strengths for the EE and IE are drawn at random from a uniform distribution from 0 mV to w max EE = 1 mV and w max IE = 2 mV, respectively, and the strengths for the EI and II connections are drawn at random from a uniform distribution from w max EI = −4 mV and w max II = −4 to 0 mV, respectively.
Based on the mean and variance of external inputs, we divide the excitatory neurons into three populations in which the αth (α = 1, 2, 3) population receives the external input of the mean µ p ext,α and standard deviation σ p ext,α . The number of neurons in the αth (α = 1, 2, 3) population is denoted as N α .

Statistical Properties of a Single LIF Neuron
We first describe the statistical properties of a single LIF neuron. The dynamics of a single LIF neuron with finite synaptic time constant is governed by Equation (1) and the following equation: where µ and σ are the mean and standard deviation of the external input respectively, and ξ (t) is the Gaussian white noise satisfying ξ = 0 and ξ (t)ξ (t ′ ) = δ(t − t ′ ).
When τ s is much smaller than τ m , the mean firing rate r has the following approximation: where β = −ς(1/2) τ s /(2τ m ), erf(·) is the error function and ζ (·) is the Riemann-Zeta function (Fourcaud and Brunel, 2002). We use a response kernel to characterize how the LIF neuron ensemble responds to a perturbation. Specifically, we consider Equation (1) and the following equation: where u(t) is a perturbation. When u(t) is weak, the difference between the instantaneous firing rate r(t) of neuronal ensemble that is governed by Equations (1, 5) and the mean firing rate r of neuronal ensemble that is governed by Equations (1, 3) is small. Therefore, we can approximate the instantaneous firing rate r(t) by the mean firing rate r plus a linear response term with respect to u(t), i.e., where h(τ ) is the linear response kernel of Equations (1, 3) and the symbol * stands for temporal convolution. The analytical results for the linear response kernel are available for the case of zero synaptic time constant (Fourcaud and Brunel, 2002). Recently, for the case of nonzero synaptic time constant, the analytical result was derived (Schuecker et al., 2015). Here, we use the reverse correlation method to compute it (Dayan and Abbott, 2001;Ostojic et al., 2009). Specifically, we inject a weak Gaussian white noise signal ξ s (t) with correlation ξ s (t + τ )ξ s (t) = σ 2 s δ(τ ) to the LIF neuron, i.e., adding ξ s (t) on the right hand side of Equation (3). Note that this testing white noise signal ξ s (t) is different from the external noise input ξ (t) in Equation (3). By using the spike-triggered average technique, we calculate the response kernel in the following way: where is the spike timing of the LIF neuron in the duration T and the firing rate r is given by Equation (4) (Dayan and Abbott, 2001). In our simulation, σ s = 2 mV and T = 10 6 ms. We use one realization of ξ s (t) and perform the average over 1,000 realizations of external noise input ξ (t).

Linear Response Theory
Here, we describe how we obtain r i and h i (τ ) for i = 1, · · · , N in a self-consistent manner for the network described in Equations (1, 2). By splitting the spike train s j (t) in Equation (2) into the mean firing rate and its fluctuations, we can rewrite Equation (2) in the following form: where the term N j=1 w ij (s j − r j ) can be regarded as a perturbation. By enforcing the first-order Volterra approximation of Equations (1, 8), we arrive at Equation (31), where r i and h i (τ ) are the mean firing rate and response kernel, respectively, of Equation (1) and the following equation: Here, we can observe that, even with the same mean and variance value of external inputs, different neurons may still have different mean firing rates due to the randomness of synaptic strengths.
To simplify the self-consistent derivation of mean firing rate, we make a further assumption that the mean firing rates of single neurons within the same population are approximated by the population-averaged mean firing rate. Denote the populationaveraged mean firing rates of the different populations as r p α , α = 1, 2, 3, I, where I stands for the inhibitory population. Then, r p α is the mean firing rate of the following equations: where w p αβ stands for the average synaptic strength from Population β to Population α. Therefore, r p α can be determined in the following self-consistent way: The Power Spectrum for the Hawkes Process The first derivation of the power spectrum of Hawkes process (i.e., Equation 31) was carried out by Hawkes (1971a), in which the Wiener-Hopf theory was used. Here, we follow the method in Grytskyy et al. (2013). First, for the spike correlation function C(τ ), the following equations hold: for τ ≥ 0, where s(t) and r(t) are the spike train and the instantaneous firing rate of Hawkes process, respectively, and D = diag{r 1 , · · · , r N }. In Equation (14b), we use the fact that s(t + τ )s T (t) = r(t + τ )s T (t) + Dδ(τ ). In Equation (14c), we insert the Hawkes process (Equation 31). In Equation (14d), we take advantage of the definition of spike correlation. Next, define and x(t) is the Gaussian white noise with correlation x(t + τ )x(t) = D x δ(τ ). For y(t), for τ ≥ 0. In Equation (16b), we have made use of the definition of y(t). In Equation (16c), we take advantage of the definition of C y (t). In Equation (18), we use the fact that x(t + τ )r T y (t) = x(t + τ ) r T y (t) = 0. By setting D x = D, we obtain C(τ ) = C y (τ ). Therefore, to understand the correlation structure of the Hawkes process, it suffices to study the correlation structure of y(t). Definer y (ω) as the Fourier transform of r y (t). From Equation (16b), we havê r y =ĥ(ω)w(r y +x).

Population-Averaged Spike Correlation
Assume that there are N α (α = 1, 2, 3, I) neurons in the αth population. The spike correlation between population α and population β is defined as follows: The corresponding spike cross correlation is defined as follows: Note that C p,cr αβ (τ ) is used in the caption of Figure 4. Inserting Equation (16a) into Equation (19), we can obtain where h p α (τ ) stands for the linear response kernel of Population α. Similar to the result (Equation 18), the population-averaged power spectrumĈ p (·) = Ĉ p αβ (·) has the following expression: where

Decomposition of Cross Correlation Structures
Formally, the power spectrum (Equation 18) can be expanded in the following way: (23) Therefore, the cross correlation C(τ ), which is the inverse Fourier transform ofĈ(ω) by the Wiener-Khinchin theorem, can be approximated by Equation (33). The same decomposition can also be carried out for the population-averaged cross correlation as discussed previously.

RESULTS
We first present the numerical evidence of the complexity of mutual interactions between plastic network structures and spiking activity, as it relates to the effect of high-order statistical structures implied in the STDP learning dynamics. In particular, we emphasize the importance of firing variability of spiking activity for the network structures learned through STDP. This fact is particularly emphasized by the sharp contrast of learned network structures between two scenarios (Case I in Figure 1 and Case II Figure 2); they share the same trend for firing rate change across different neuronal populations, but have the opposite trend for their firing variability change. We then develop a selfconsistent linear response theory to account for the effects of both firing rate and firing variability on learning dynamics and explain the role observed in our numerical study of high-order statistics in learning dynamics. The initial values of excitatory-to-excitatory connection strengths are drawn from a uniform distribution ranging from 0 to 1 mV. (A,C) Raster plots before learning and after learning, respectively. (B,D) Firing rates and coefficients of variation (CVs) of the different populations before learning and after 10 6 s of learning, respectively. Error bars indicate the standard deviation of corresponding values within each population. (E) Firing rate dynamics of the different populations during learning. (F) Weight matrix after 10 6 s of learning. The largest (i.e., brightest) and the smallest (i.e., darkest) value of the gray bar is 1 and 0 mV, respectively.

The Effect of Firing Rate and Firing Variability on Learning Dynamics-Numerical Evidence
For our numerical experiments, we use the current-based leaky integrate-and-fire (LIF) neuron model with finite synaptic time constant (Babadi and Abbott, 2013). The network consists of 250 excitatory neurons and 250 inhibitory neurons (See section Methods for details). Based on the values of mean and variance of external inputs, we divide the excitatory neurons into three populations-Population 1 (denoted as P 1 , including neurons indexed from 1 to 50), Population 2 (denoted as P 2 , including neurons indexed from 51 to 200) and Population 3 (denoted as P 3 , including neurons indexed from 201 to 250). The initial values of excitatory-to-excitatory (EE) connection strengths are drawn at random from a uniform distribution ranging from 0 to w max EE = 1 mV. During the simulation, the EE connections are subject to the canonical all-to-all pairwise additive STDP learning rule (Song et al., 2000;Babadi and Abbott, 2013). Specifically, let τ = t i − t j , where t i and t j are a pair of spike times of excitatory neuron i and excitatory neuron j, respectively. This pair of spike times will affect a change of the synaptic strength w ij from neuron j to neuron i as follows: w ij → w ij + L(τ ), where the STDP learning window L(τ ) is given by Similarly, for the synaptic strength w ji from neuron i to neuron j, the update is w ji → w ji + L(−τ ). As in Song et al. (2000) FIGURE 2 | The effect of firing rate and variability on STDP learning dynamics-Case II. Different populations receive external inputs with different means as well as different variances. Population 1 (including neurons labeled by 1-50, purple) receives the smallest mean (µ p ext,1 = 27.5 mV) but the largest variance (σ p ext,1 = 31.6 mV), Population 2 (including neurons labeled by 51-200, blue) receives an intermediate mean (µ p ext,2 = 30 mV) and variance (σ p ext,2 = 22.4 mV), and Population 3 (including neurons labeled by 201-250, black) receives the largest mean (µ p ext,3 = 32.5 mV) but the smallest variance (σ p ext,3 = 11.2 mV). The initial values of excitatory-to-excitatory connection strengths are drawn from a uniform distribution ranging from 0 to 1 mV. (A,C) Raster plots before learning and after learning, respectively. (B,D) Firing rates and coefficients of variation (CVs) of the different populations before learning and after 10 6 s of learning, respectively. Error bars indicate the standard deviation of corresponding values within each population. (E) Firing rate dynamics of the different populations during learning. (F) Weight matrix after 10 6 s of learning. The largest (i.e., brightest) and the smallest (i.e., darkest) value of the gray bar is 1 and 0 mV, respectively. and Babadi and Abbott (2013), we enforce the hard bound to the value of w ij and w ji , i.e., set In this work, we will restrict ourselves to the balanced case in which A + = A − = 0.005 w max EE and τ + = τ − = 20 ms. We characterize spiking activity in the network by its firing rate and coefficient of variation (CV), where CV is the ratio of the standard deviation of interspike interval (ISI) to the mean of ISI. As is well known, for a given neuron under fixed input conditions, CV is often used as a measure of firing variability -the larger CV, the more irregular the spiking activity.
If these populations receive external inputs of different means but an identical variance, we observe that, as the mean external input increases, the firing rate increases whereas the CV decreases. Figure 1B illustrates the corresponding firing statistics of these populations. For this case, as shown in Figure 1F, the connections from the populations of higher firing rate to those of lower firing rate are strengthened after learning whereas the connections with the opposite situation (i.e., connections from those of lower firing rate to those of higher firing rate) are weakened. This result is consistent with results of the previous studies (Morrison et al., 2007;Kozloski and Cecchi, 2010;Babadi and Abbott, 2013). Here, STDP is found to act as a firing rate equalizer, i.e., decreasing the rate of the neuronal populations of higher firing rate and increasing the rate of the neuronal populations of lower firing rate. Figure 1E displays an example of this behavior. Furthermore, by comparing the firing statistics before learning ( Figure 1B) with that after learning (Figure 1D), we observe that STDP also acts as a firing variability equalizer, i.e., decreasing the CV of higher-CV neuronal populations and increasing the CV of lower-CV neuronal populations. For ease of discussion below, this case will be termed as Case I.
Next, we consider the case in which both the mean and variance of external inputs vary across these populations. We adjust both the mean and variance of external inputs (see the caption of Figure 2 for the specific values) such that Population 1 has the highest firing rate and the highest CV, Population 2 has an intermediate firing rate and CV, and Population 3 has the lowest firing rate and the lowest CV. The corresponding firing statistics of these populations before learning is displayed in Figure 2B. From Figure 2F, it can be seen clearly that after learning the connections from Population 1, which has the highest firing rate, to Population 3, which has the lowest firing rate, are weakened. We note that this phenomenon still exists even when the projections from the excitatory populations to the inhibitory populations become plastic (Supplementary Material). Furthermore, by comparing Figure 2B with Figure 2D, we observe that STDP tends to increase the rate of the population of the highest firing rate and decrease the rate of the population of the lowest firing rate. That is, STDP is no longer acting as a firing rate equalizer. However, it still acts as a weak equalizer for firing variability. This phenomenon cannot be explained solely by the firing rate theory (Morrison et al., 2007;Kozloski and Cecchi, 2010;Babadi and Abbott, 2013), in which the connections from the population of higher firing rate to the population of lower firing rate should be strengthened and STDP should act as a firing rate equalizer. In the following, we develop a linear response theory to characterize the importance of high-order firing statistics in learning dynamics as illustrated in Figure 2. For ease of discussion below, the case in Figure 2 will be termed as Case II.

STDP Learning Dynamics
To characterize the learning dynamics, we consider two mutually-connected neurons in a recurrent network, say neuron i and neuron j. The STDP learning rule introduced in the last section can be simply written in the following weight dynamics: where s i (t) and s j (t) respectively stand for the spike trains of neuron i and neuron j. Expressing x + j and x − i as the integral of s j (t) and s i (t), respectively, we can obtain By assuming that the change of w ij (t) is sufficiently slow (Ocker et al., 2015), s j (t − τ )s i (t) can be approximated by C ij (τ ) + s i (t) s j (t) where C ij (τ ), defined as s i (t + τ )s j (t) − s i (t) s j (t) , is the temporal spike cross correlation between neuron i and neuron j in the recurrent network with a fixed synaptic strength matrix. Combining the dynamics of w ji in Equation (25), we arrive at the following equations: Note that unlike the scenario studied in Babadi and Abbott (2013) and Ocker et al. (2015), there is no offset term in Equations (29, 30). This is due to the balanced condition A + = A − and τ + = τ − in our STDP learning rule. Therefore, to understand the learning dynamics (Equations 29, 30), it is facilitative to study the correlation structure in recurrent networks with static couplings. In the following sections, we will follow the linear response theory proposed in Pernice et al. (2012), Trousdale et al. (2012), and Helias et al. (2013) and develop a self-consistent framework for the learning dynamics (Equations 29,30). To this end, we first characterize how the firing variability information is encoded in the linear response kernel for the single LIF neuron case in section Incorporation of Firing Variability in Linear Response Kernel. Then, in section Hawkes Process as an Approximation of an LIF Network Dynamics, we approximate the LIF network dynamics with a stochastic point process called Hawkes process in which the instantaneous firing rate depends on the spiking history of interacting neurons (Hawkes, 1971b). With this Hawkes process approximation, in section Spike Correlation Structure of Hawkes Process, we discuss the issue of how firing statistics of spiking activity and network connectivity structure shape the spike correlation structure of the LIF network dynamics. In the end, we combine the correlation structure analysis with the STDP learning dynamics (Equations 29,30) to arrive at the final learning dynamics (Equations 35, 36), which explain the observed numerical results, in section Learning Dynamics with Spike Correlation Decomposition.

Incorporation of Firing Variability in Linear Response Kernel
For a single LIF neuron described by Equations (1, 3), the dependence of mean firing rate and CV on input parameters (mean input value µ and standard deviation value σ ) is illustrated in Figures 3A,B, respectively. By using the reverse correlation method (Dayan and Abbott, 2001;Ostojic et al., 2009) (see section Methods for details), We also computed the response kernel h(τ ) that characterizes how the LIF neuron ensemble responds to a perturbation. We found that response kernel h(τ ) exhibits an exponential decay profile, A exp{−τ /τ eff } (τ ), where (τ ) is the Heaviside function, provided that fluctuations of the external input are sufficiently high, i.e., comparable to the difference between firing threshold and reset voltage V th − V r (for example, σ is about 10 ∼ 30 mV). We can use two quantities -the integral of response kernel h(τ )dτ , which is τ eff A after the substitution of the exponential form of h(τ ) = A exp{−τ /τ eff } (τ ), and the decay time constant τ effto characterize h(τ ). Figures 3C,D illustrate the dependence of h(τ )dτ and τ eff on µ and σ , respectively. In Figure 3, we can observe that, for any fixed mean firing rate, when the standard deviation σ increases, both CV of ISI and τ eff increase whereas h(τ )dτ decreases. Therefore, if there are two neurons with the same mean firing rate but different firing CVs, both the decay time constant and the integral of response kernel of these two neurons will differ from each other. In this sense, the information of firing variability is incorporated in the linear response kernel.

Hawkes Process as an Approximation of an LIF Network Dynamics
In general, the response of a neuron to recurrent spike inputs is nonlinear. Here, we investigate the asynchronous state, as shown in Figures 1A, 2A, in which spikes from different neurons arrive at the synapses of the target neuron at different times. By considering that the magnitude of a postsynaptic potential elicited by the firing event of a single presynaptic neuron is in general small, i.e., the synaptic strength is weak, it is reasonable to treat the asynchronous recurrent inputs in a linear manner (Pernice et al., 2012;Trousdale et al., 2012;Helias et al., 2013). Specifically, by regarding N j=1 w ij (s j − r j ) in the LIF neuronal network as a perturbation, we arrive at the following effective description: where r i and h i (τ ) are the mean firing rate and response kernel of the ith neuron that is driven by both the external input and the mean-driven recurrent input (see section Methods for details), respectively, and N is the number of neurons in the network. Note that (i) the mean firing rate r i is determined by the selfconsistent Equations (12, 13); and (ii) the response kernel h i (τ ) is computed based on the single LIF neuron Equations (1, 3) with the total mean input value µ p i in Equation (13) and the standard deviation value σ p ext,i in Equation (12). Equation (31) can be understood in the following way. On the one hand, given spike trains s j (t) for j = 1, · · · , N, Equation (31) provides a way to evolve the instantaneous firing rate of the ith neuron r i (t). On the other hand, from Equation (31), the spike train s i (t) can be regarded as being generated by a jump process with the instantaneous firing rate r i (t). Therefore, we have obtained a self-consistent stochastic jump process that produces both the instantaneous firing rate r i (t) and the stochastic spike train s i (t). This stochastic process Equation (31) is a Hawkes process (Hawkes, 1971b), which has been extensively used in the neuroscience literature (Pernice et al., 2012;Trousdale et al., 2012;Helias et al., 2013). From the discussion above, such a Hawkes process can be viewed as the first-order (i.e., linear) Volterra approximation of the LIF network dynamics.
To examine the validity of our linear response theory, we compute the population-averaged cross correlation (see section Methods for details). We find that our theoretical prediction of the population-averaged cross correlation is in very good agreement with the cross correlation measured from our LIF neuronal network simulation. This agreement is clearly seen in Figure 4, demonstrating that our linear response theory can well capture the spike correlation structure of the LIF network dynamics.
We note that Equation (32) also provides a starting point to study how network connectivity affects the spike correlation structure since we can decompose the connections into direct connections, common inputs, and other higher-order connections. In particular, the spike correlation has the following decomposition (see section Methods for a derivation): The first term on the right hand side (RHS) of Equation (33) can be obtained by setting the connection strength w in the Hawkes process (Equation 31) to zero. There is no cross correlation in this term.
For the second term on the RHS of Equation (33), the corresponding component with respect to the correlation between neuron i and neuron j is h i (τ )w ij r j . This involves firing rate r j of presynaptic neuron j, the connection strength w ij from neuron j to neuron i and the linear response kernel h i (τ ) of postsynaptic neuron i. Therefore, this component corresponds to the correlation induced by the direct connection from neuron j to neuron i. Similarly, the component h j (−τ )w ji r i (the third term on the RHS of Equation 33), corresponds to the correlation induced by the direct connection from neuron i to neuron j. Combining these two terms, we obtain the spike correlation induced by direct connections between neuron i and neuron j. For ease of discussion below, we denote these terms as C di ij (t).
We next discuss the terms associated with the indirect connections. For the fourth term on the RHS of Equation (33), the corresponding component with respect to the correlation between neuron i and neuron j is k w ik r k w jk [h i * h −j ](τ ). This involves the firing rate r k of neuron k, the connection strengths w αk and the response kernels h α (·) for α = i, j. Therefore, this term corresponds to the correlation induced by the common inputs from neuron k to both neuron i and neuron j. As we have demonstrated above, h α (·), α = i, j, exhibits an exponential decay profile, which will be expressed as B α exp(−t/τ α ) (t). Then, this component becomes For ease of discussion below, we denote this term as C co ij (t), the remaining higher-order terms in Equation (33) as C hi ij (t), and the summation of C co ij (t) and C hi ij (t) as the indirect term C in ij (t). We can observe from Equation (34) that, if there is a difference between response time constants τ i and τ j , C co ij (t) will show an asymmetric property with respect to the axis t = 0. More precisely, if we assume τ i > τ j , C co ij (t) > C co ij (−t) for any t > 0. This asymmetry can be understood intuitively as follows. When there is a common input from neuron k into both neurons i and j, the neuron with the smaller response time constant (neuron j) will respond more rapidly than the neuron with the larger response time constant (neuron i) on average. Therefore, there is a higher likelihood that neuron j will fire before neuron i. In terms of the spike correlation, this implies that C co ij (t) > C co ji (t) = C co ij (−t) for any t > 0. We note that, as shown in Equation (34), the common inputs from both excitatory and inhibitory neurons into neuron i and neuron j induce a positive spike correlation between neuron i and neuron j. Therefore, if there is a significant difference between τ i and τ j , there will be a substantial asymmetry between the positive-t component of C co ij (t) and the corresponding negative-t component. This will have a profound impact on the structure of C ij (t) as discussed below.
Now we study how different connectivity patterns affect the cross correlation structure between different neuronal populations for Cases I and II in Figures 1, 2. As an example, we pick neuron i from Population 3 and neuron j from Population 1. Note that, for Case I in Figure 1, neuron i (j) has a lower (higher) firing rate and a higher (lower) firing variability. Whereas, for Case II in Figure 2, neuron i (j) has a lower (higher) firing rate and a lower (higher) firing variability.
For Case I in Figure 1, we find that the net mean inputs to neuron j in Population 1 and neuron i in Population 3 are 31.45 and 11.55 mV, respectively. Note that the external mean inputs µ p ext,j and µ p ext,i are 40 and 20 mV, respectively. Since the net mean input is the summation of external mean input and recurrent input, the net recurrent inputs to neuron j and neuron i are −8.55 mV and −8.45 mV, respectively. This means the overall recurrent inputs to both neurons are inhibitory. In Case I, the standard deviation values for both neurons are 15.81 mV. With these net FIGURE 4 | Population-averaged cross correlation-comparison between theoretical prediction and numerical results. (A-F) Normalized population-averaged cross correlation between the different populations. Note that C p,cr αβ for α, β = 1, 2, 3 is defined in Equation (20) and r p α for α = 1, 2, 3 is computed with Equations (12, 13) in a self-consistent way. The black dots are the cross correlation measured from a simulation of the LIF neuronal network with a static weight matrix. The red dashed line is the theoretical prediction of our linear response theory.
mean input values and standard deviation values, we compute the response kernel by using the reverse correlation method (see section Methods for more details). We find that h j (t)dt > h i (t)dt and τ j < τ i , which is evident in Figures 3C,D. From the left panels of Figure 5A, we can observe that (i) the positive-t component of C di ij (t) is larger than the corresponding negativet component; and (ii) the positive-t component of C in ij (t) is slightly smaller than the corresponding negative-t component. To further understand how different connectivity patterns shape C in ij (t), we split C in ij (t) into the common-input correlation C co ij (t) and the higher-order-connection correlation C hi ij (t). It can be seen clearly from the right panels of Figure 5A that the positive-t component of C hi ij (t) is smaller than the corresponding negative-t component whereas the positive-t component of C co ij (t) is larger than the corresponding negative-t component. A combination of these two facts yields a small difference between the positive-t component of C in ij (t) and the corresponding negative-t component. We note that the asymmetry of C co ij (t) with respect to the axis t = 0 exhibited in the top right panel of Figure 5A is a consequence of τ j < τ i .
For Case II in Figure 2, we find that the net mean inputs to neuron j in Population 1 and neuron i in Population 3 are 15.58 and 20.50 mV, respectively. Note that the external mean inputs µ p ext,j and µ p ext,i are 27.5 and 32.5 mV, respectively. Since the net mean input is the summation of external mean input and recurrent input, the net recurrent inputs to neuron j and neuron i are −11.92 and −12.00 mV, respectively. This means the overall recurrent inputs to both neurons also are inhibitory. In Case II, the standard deviation values for neuron j and neuron i are 31.6 and 11.1 mV, respectively. With these net mean input values and standard deviation values, we compute the response kernel by using the reverse correlation method (see section Methods for FIGURE 5 | Decomposition of the spike cross correlation between neurons from the different populations before learning. Here, we choose neuron i from Population 3 and neuron j from Population 1. The connection strengths w ij and w ji are fixed at 0.5 mV. The top left, top right, bottom left and bottom right panels of both (A,B) correspond to the correlations associated with direct connections, common inputs, indirect connections, and higher-order connections, respectively. (A) Compared to neuron i, neuron j receives a larger mean and an identical variance of external inputs. This corresponds to Case I in Figure 1. (B) Compared to neuron i, neuron j receives a smaller mean but a larger variance of external inputs. This corresponds to Case II in Figure 2. more details). We find that h j (t)dt < h i (t)dt and τ j > τ i , which is also evident in Figures 3C,D. From the left panels of Figure 5B, we can observe that (i) the positive-t component of C di ij (t) is again larger than the corresponding negative-t component; and (ii) the positive-t component of C in ij (t) is much smaller than the corresponding negative-t component. From the right panels of Figure 5B, we find that the positive-t components of both C co ij (t) and C hi ij (t) are smaller than the corresponding negative-t components. A combination of these two effects yields a large difference between the positive-t component of C in ij (t) and the corresponding negative-t component. Note that the fact that C co ij (t) < C co ij (−t) for t > 0 in Figure 5B is a consequence of the inequality τ j > τ i . In the next section, we will demonstrate the extent of asymmetry between the positive-t component of C in ij (t) and the corresponding negative-t component underlies the difference of network structures learned through STDP presented in Figures 1, 2.

Learning Dynamics with Spike Correlation Decomposition
Combining the STDP learning dynamics (Equations 29, 30) with the direct-and indirect-connection decomposition of the correlation structure in the previous section, we arrive at the following learning dynamics: where A α = ∞ 0 A + exp{−τ/τ + }h α (τ )dτ for α = i, j and B = L(τ )C in ij (τ )dτ . In Equations (35, 36), the first and second terms on the RHS correspond to the effect of direct connections between neuron i and neuron j, whereas the third term corresponds to the effect of indirect connections between neuron i and neuron j. Note that here we have used the balanced STDP assumption, i.e., A + = A − and τ + = τ − .
As argued in Babadi and Abbott (2013), while a linear system of differential equations such as Equations (35, 36) cannot have more than two attractors, the hard lower and upper bounds enforced by the STDP learning rule, i.e., the condition 0 ≤ w ij ≤ w max EE , can lead to the following two attractors: (w ij , w ji ) = (w max EE , 0) and (w ij , w ji ) = (0, w max EE ). For the sake of discussion, we consider A i r j > A j r i below. Then, the line w ji = 1 A j r i (A i r j w ij + B) has the slope greater than 1, as illustrated by the black dash line in Figure 6. If neglecting the indirect connection effect, that is, setting B = 0, we can observe from Equations (35, 36) that, for any (w ij , w ji ) pair satisfying A i w ij r j − A j w ji r i > 0, the pair will converge to (w max EE , 0). Whereas, for any (w ij , w ji ) pair satisfying A i w ij r j − A j w ji r i < 0, the pair will converge to (0, w max EE ). This is illustrated in Figure 6A with the black dash line separating the two basins. In this case, the size of basin of the attractor (w max EE , 0) is larger than that of the attractor (0, w max EE ). Therefore, w ij would be more likely to be potentiated, whereas w ji would be more likely to be depressed.
When the indirect connection effect is included, there exist two distinct regimes depending on the value of B. For the sake of discussion, we denote the size of basin of (w max EE , 0) (i.e., the area below the line A i w ij r j −A j w ji r i +B = 0) as S ij and the size of basin of (0, w max EE ) (i.e., the area above the line A i w ij r j − A j w ji r i + B = 0) as S ji . One is a regime in which S ij > S ji . This occurs when B ≥ 0 or B is negative but small. As illustrated in Figure 6A, the red line separates the two basins. In this regime, w ij would be more likely to be potentiated whereas w ji would be more likely to be depressed. For neuronal populations, the synapses from the neuron j's population to neuron i's population will have a higher likelihood to be strengthened. In contrast, the synapses with the opposite direction will be more likely to be weakened. The other is a regime in which S ij < S ji . This occurs when B is sufficiently negative, as indicated by the red line in Figure 6B. In this regime, w ij would be more likely to be depressed whereas w ji would be more likely to be potentiated. For neuronal populations, the synapses from neuron i's population to neuron j's population will have a higher likelihood to be strengthened, whereas the synapses with the opposite direction will be weakened with a larger probability.
We now examine the cases in Figures 1, 2. For Case I in Figure 1 and Table 1 displays the corresponding values of A i r j , A j r i , B, S ij , S ji as well as the average connection strengths w ij and w ji after 10 6 s of learning. It can be seen clearly that for all the three (i, j) pairs, namely (i ∈ P 2 , j ∈ P 1 ), (i ∈ P 3 , j ∈ P 1 ) and (i ∈ P 3 , j ∈ P 2 ), S ij > 0.5mV 2 > S ji , indicating that w ij for all these three pairs would be strengthened. Indeed, this strengthening is confirmed by the values of w ij presented in Table 1 (i.e., w ij > 0.5mV > w ji for all three pairs). To evaluate the contribution of common inputs and higher-order connections separately, we also include B co (defined as L(τ )C co ij (τ )dτ ) and B hi (defined as L(τ )C hi ij (τ )dτ ) in Table 1. We can observe that for all three pairs B co > 0 whereas B hi < 0. Since B = B co + B hi , this gives rise to a negative but small value of B, yielding a regime in which S ij > S ji . Since the connections from the populations of higher firing rate to those of lower firing rate are strengthened whereas the connections with the opposite direction are weakened, the rate of the populations of lower firing rate will increase whereas that of the populations of higher firing rate will decrease. The strengthening of connections from the populations of higher firing rate to those of lower firing rate leads to the increase of the mean of total inputs into the populations of lower firing rate. As illustrated in Figure 3B, when the variance of total inputs is fixed, the larger the mean total inputs, the smaller the CV. Therefore, this explains why in Figure 1 the CV of the populations of lower firing rate, which is larger before learning, decreases after learning. Therefore, STDP here acts as a weak equalizer for both firing rate and firing variability, as demonstrated in Figure 1.
For Case II in Figure 2 and Table 2 shows that for all three (i, j) pairs, namely (i ∈ P 2 , j ∈ P 1 ), (i ∈ P 3 , j ∈ P 1 ) and (i ∈ P 3 , j ∈ P 2 ), S ij < 0.5 mV 2 < S ji and w ij < 0.5 mV < w ji . Since the connections from the populations of higher firing rate to the populations of lower firing rate are weakened whereas the connections with the opposite direction are strengthened, the rate of the populations of higher firing rate will increase further whereas that of the populations of lower firing rate will decrease further. Meanwhile, the weakening of connections from the populations of higher firing rate to those of lower firing rate leads to the decrease of the mean of total inputs into the populations of lower firing rate. As illustrated in Figure 3B, the smaller the mean total inputs, the larger the CV. Therefore, the CV of these populations of lower firing rate, which is smaller before learning, increases after learning. Note that, in Case II, those populations of lower firing rate also have the lower firing FIGURE 6 | Phase plane analysis of learning dynamics. This learning dynamics has two attractors -(w ij , w ji ) = (1, 0) and (w ij , w ji ) = (0, 1). The black dashed line satisfies A i w ij r j − A j w ji r i = 0 while the red line satisfies A i w ij r j − A j w ij r i − B = 0. We assume that A i r j > A j r i . Therefore, the slope of both the black dashed line and the red line is greater than 1. Any (w ij , w ji ) pair in the area above the red line will converge to (0, 1) whereas any pair in the area below the red line will converge to (1, 0), as indicated by the blue arrows. Therefore, the area above the red line is the basin of the attractor (0, 1) whereas the area below the red line is the basin of the attractor (1, 0). Note that w ij and w ji evolve along the line w ij + w ji = Const as indicated by the arrows. (A) The size of basin of the attractor (0, 1) is smaller than that of the attractor (1, 0). In this case, synapses from the neuronal population to which neuron j belongs, to the neuronal population to which neuron i belongs, would be more likely to be strengthened, whereas the synapses with the opposite direction will be weakened. (B) The size of basin of the attractor (0, 1) is larger than that of the attractor (1, 0). In this case, synapses from the neuronal population to which neuron j belongs, to the neuronal population to which neuron i belongs, would be more likely to be weakened, whereas the synapses with the opposite direction will be strengthened. A i r j × 10 6 A i r j × 10 6 B × 10 6 B co × 10 6 B hi × 10 6 S ij w ij w ji (ms −1 ) (ms −1 ) (ms −1 · mV) (ms −1 · mV) (ms −1 · mV) (mV 2 ) (mV) (mV) (i ∈ P 2 , j ∈ P 1 ) 0 variability. This explains why STDP acts as a weak equalizer of firing variability, instead of firing rate, in Figure 2. Furthermore, for Case II in Figure 2, by examining the contribution of common inputs and higher-order connections, we can find from Table 2 that, compared to B co , B hi is negligible for the pairs (i ∈ P 2 , j ∈ P 1 ) and (i ∈ P 3 , j ∈ P 1 ). Therefore, for these two pairs, it is the contribution of commoninput correlations that yields a large negative value of B and the corresponding network structures in Figure 2. Note that the asymmetric property of common-input correlations arises from the difference of response time constants between different neurons. As demonstrated in Figure 3, the response time constant encodes the information of the firing variability of spiking activity. In this sense, we have established a direct link between firing variability and learning dynamics by analyzing the spike correlation structures, and have demonstrated explicitly how firing variability shapes the network structures learned through STDP.

DISCUSSION
We have investigated how the high-order statistics of spiking activity affect the learning dynamics in a plastic LIF neuronal network. We have developed a self-consistent linear response theory, which incorporates the information of both firing rate and firing variability. By decomposing the spike correlation structure into the component associated with local direct connections and that associated with indirect connections, we have identified two distinct regimes regarding the network structures learned through STDP. In one regime, the contribution of the correlation associated with local direct connections dominates in STDP learning dynamics, leading to network structures consistent with the previous work (Morrison et al., 2007;Kozloski and Cecchi, 2010;Babadi and Abbott, 2013). In the other regime, the contribution of the correlation associated with indirect connections, instead of local direct connections, dominates, yielding network structures at variance with the prediction by the simple firing rate theory. The dominance of the indirect-connection correlation mainly arises from the asymmetric common-input correlations, which is induced by the heterogeneity of firing variability across different neuronal populations. Therefore, our study highlights the importance of firing variability of spiking activity for the evolution of plastic neural networks.
In a series of work (Gilson et al., 2009a(Gilson et al., ,b,c,d, 2010, the STDP learning dynamics in recurrent networks has been investigated with Hawkes processes. However, the response kernel of the Hawkes process used there is identical across neuronal populations, thus the effect of firing variability was not addressed. Here, we study the STDP learning dynamics with the LIF neuron, in which the effect of firing variability is addressed by varying the mean and variance of external inputs across different neuronal populations. In Babadi and Abbott (2013), the authors developed a theory to explain the formation of in-and out-hubs. However, the effect of firing variability has not been discussed. Our theory can be viewed as an extension of the work (Babadi and Abbott, 2013). More recently, motif dynamics have been studied with linear response theory for one homogeneous population (Ocker et al., 2015). Instead, our work focuses on the evolution of connection strengths across different neuronal populations.
We have focused here on the balanced STDP rule. However, it is straightforward to extend our analysis to the unbalanced case by including additional terms in learning dynamics (Equations 35,36). An unbalanced STDP rule can be expected to act as a loop-generating or loop-eliminating mechanism, depending on whether the rule is potentiation-dominated or depressiondominated (Babadi and Abbott, 2013). The neuron model we used here is the current-based LIF neuron. As we have demonstrated, the key property of a single neuron model to have an indirect-connection-dominated regime is that the firing variability of spiking activity can be encoded in the time constant of linear response kernel. Therefore, it can be expected that the indirect-connection-dominated regime can exist with other neuron models, such as conductance-based models, exponential integrate-and-fire models or more realistic Hodgkin-Huxley neurons. The connectivity here is of the allto-all type. We can readily extend our self-consistent linear response theory to the sparsely connected network case. For larger network size, such as the network of 500 excitatory neurons and 500 inhibitory neurons, we have also found that there is an indirect-connection-dominated regime, provided that there is a significant difference of firing variability between different neuronal populations. In this work, the input statistics (i.e., input means and input variances) are kept constant, indicating that the circuit we studied receives the same sensory stimuli during the whole simulation. Under this assumption, the neural population activity can be simply quantified by mean firing rate and CV of ISI. It would be very interesting to study the learning dynamics when the stimulation periods are interleaved by spontaneous states. Another interesting issue to be addressed in the future work is how the high-order firing statistics of spiking activity affect the connectivity structures in a circuit with more complex learning rules, such as triplet STDP learning rule (Gjorgjieva et al., 2011).
In conclusion, we have demonstrated that the introduction of heterogeneity of firing variability across different neuronal populations can give rise to counterintuitive network structures in a plastic neural circuit. Our work suggests that the high-order statistics of spiking activity can play an important role in shaping the connectivity patterns in our brain.

AUTHOR CONTRIBUTIONS
Conceived and designed the research, performed experiments, analyzed data, and wrote the paper: BM, DZ, and DC.