Spike Pattern Structure Influences Synaptic Efficacy Variability under STDP and Synaptic Homeostasis. II: Spike Shuffling Methods on LIF Networks

Synapses may undergo variable changes during plasticity because of the variability of spike patterns such as temporal stochasticity and spatial randomness. Here, we call the variability of synaptic weight changes during plasticity to be efficacy variability. In this paper, we investigate how four aspects of spike pattern statistics (i.e., synchronous firing, burstiness/regularity, heterogeneity of rates and heterogeneity of cross-correlations) influence the efficacy variability under pair-wise additive spike-timing dependent plasticity (STDP) and synaptic homeostasis (the mean strength of plastic synapses into a neuron is bounded), by implementing spike shuffling methods onto spike patterns self-organized by a network of excitatory and inhibitory leaky integrate-and-fire (LIF) neurons. With the increase of the decay time scale of the inhibitory synaptic currents, the LIF network undergoes a transition from asynchronous state to weak synchronous state and then to synchronous bursting state. We first shuffle these spike patterns using a variety of methods, each designed to evidently change a specific pattern statistics; and then investigate the change of efficacy variability of the synapses under STDP and synaptic homeostasis, when the neurons in the network fire according to the spike patterns before and after being treated by a shuffling method. In this way, we can understand how the change of pattern statistics may cause the change of efficacy variability. Our results are consistent with those of our previous study which implements spike-generating models on converging motifs. We also find that burstiness/regularity is important to determine the efficacy variability under asynchronous states, while heterogeneity of cross-correlations is the main factor to cause efficacy variability when the network moves into synchronous bursting states (the states observed in epilepsy).


Basic properties of efficacy variability under STDP and synaptic homeostasis
Here we summarize some basic properties of efficacy variability under STDP and synaptic homeostasis, which is important for understanding this paper.

Comparing E T (Var S (∆W)) with Var T (E S (∆W))
In Section 2.1 of the main text, we mentioned that the efficacy variability can be written as Var S,T (∆W) = Var T (E S (∆W)) + E T (Var S (∆W)), and that in our theoretical context, Var T (E S (∆W)) ≈ 0, so that Var S,T (∆W) ≈ E T (Var S (∆W)).
Here, we provide numerical validation to this statement.
We implement synaptic homeostasis during the synaptic evolution under STDP, so that the mean weight of all the synapses into a neuron is kept constant (eq.6 in the main text). Therefore, under synaptic homeostasis, E S (∆W) in different trials are all the same, so Var T (E S (∆W)) = 0, and eq.2 is strictly correct. To understand the situation that synaptic homeostasis is absent, in Supplementary Figure 1, we plot the ratio Var T (E S (∆W)) E T (Var S (∆W)) under the original spike patterns, in which ∆W representing the synaptic changes only caused by STDP. We see that this ratio is very small, which suggests that we still have eq.2 even without synaptic homeostasis. under the original spike patterns, in which ∆W representing the synaptic changes only caused by STDP. When τ I d = 3ms, 4ms or 5ms, this ratio is so small that we cannot clearly see the bar representing its value in the plot. Spike patterns lasted for 41s, with STDP starting after the initial 2s of transient period. We used the synaptic strengths at the end of the STDP evolution to construct the matrix ∆W. Simulations were run for 32 trials.

Efficacy variability in converging motifs and in a general network under synaptic homeostasis
In our previous paper (Bi and Zhou, 2016), we studied the efficacy variability in converging motifs. Here we explain the relation between the efficacy variability in converging motifs and that in a general network. Using the law of total variance (Weiss, 2005), the variance of all the synaptic changes in a general plastic network can be decomposited as Var ab (∆w ab ) = Var a (E b∈∂a (∆w ab )) + E a (Var b∈∂a (∆w ab )), with ∆w ab being the synaptic changes from neuron b to neuron a, and ∂a being all the neurons that input to neuron a in the plastic network. Under synaptic homeostasis eq.6 in the main text, the mean synaptic strength input to every neuron is the same, so Var a (E b∈∂a (∆w ab )) = 0. Therefore, Var ab (∆w ab ) = E a (Var b∈∂a (∆w ab )), which means that the efficacy variability totally comes from the efficacy variability within the converging motifs in the network (i.e., Var b∈∂a (∆w ab ) in the second term). Therefore, converging motifs are basic units to understand the efficacy variability of a general network under synaptic homeostasis, and our previous study paved the way to understand the efficacy variability of a general network.
1.3 Some properties of DriftV in converging motifs

P-D imbalance and its contribution to DriftV
In a converging motif, the trial-averaged change of the ath synapse per unit time can be written as in which r 0 and r a respectively being the time-averaged firing rates of the central and non-central neurons, H(τ ) being the STDP time window, τ delay being the axonal delay, and C a (τ ) ≡ r0(t)ra(t+τ ) r0(t) ra(t) being the unit crosscorrelation between the central neuron and the ath non-central neuron. Then the DriftV of a converging motif can be quantified by If C a is the same for different as, then it can be shown that with Var a (r a ) quantifying the heterogeneity of firing rates, and |E a (v a )| quantifying the imbalance of the strengths of the potentiation and depression process under STDP (P-D imbalance). This equation means that P-D imbalance is able to induce DriftV through its interaction with heterogeneity of rates.
Because of the additive nature of our synaptic homeostasis (eq.6 in the main text), efficacy variability of converging motifs does not change after we consider the synaptic changes caused by synaptic homeostasis. However, we should notice that E a (v a ) in eq.7 means the synaptic changes only caused by STDP (without considering synaptic homeostasis). We have also pointed this out in the discussion below eq.11 in the main text.
1.3.2 Manifesting the contributions to DriftV by P-D imbalance and heterogeneity of cross-correlations Eq. 6 suggests that the heterogeneity of cross-correlation can influence DriftV, and our discussion above suggests that the interaction of P-D imbalance and heterogeneity of rates is another source of DriftV. These two contributions to DriftV can be explicitly split under the condition that r a and C a in eq.6 are independent.
Suppose X and Y are independent random variables, then it can be proved that Using this formula, assuming r a and C a are independent, eq.6 becomes From eq.5, E a [ ∞ −∞ dτ H(τ )C a (−τ delay − τ )] = Ea(va) Ea(ra)r0 , which makes the equation above becomes In this equation, the first term represents the contribution of the P-D imbalance (i.e. |E a (v a )|) under heterogeneity of firing rates (i.e. Var a (r a )), while the second term represents the contribution of heterogeneity of cross-correlations . When potentiation and depression balance each other (i.e. E a (v a ) = 0), we see that 2 Adjusting the firing rates of the excitatory population for different As we mentioned in Section 2.3 of the main text, in our simulations, the decay time constant τ I d of all the inhibitory synapses are the same in one simulation trial, but may be different for different trials. We used different c parameters for different τ I d s, so that the firing rates of the excitatory population are kept around 20Hz regardless of τ I d . In this section, we introduce the method we used to choose the c parameter under a given τ I d . During our simulations, we found that the firing rate of the excitatory population tends to monotonically increase with c in a relative large parameter range. Based on this fact, given the value of τ I d , we used the following pseudocode to find the suitable c parameter: After iterating the 2nd and 3rd steps above for a sufficient number of times, r new starts to randomly walk around 20Hz because of the trial-to-trial variability of the calculated firing rate. From the last 15 times of the iteration, we chose the c new value that corresponding to the r new most close to 20Hz by eye looking, and regarded it as the c parameter for the given τ I d . We plot the firing rate as a function of τ I d under the chosen c parameters in Supplementary Figure 2A. The firing rates in this figure are all larger than 20Hz, which is probably due to the human bias during the eye-looking choosing. However, this does not matter, because what we care is that the firing rates of the excitatory population should be kept around the same value for different τ I d s, so that the change of efficacy variability with τ I d shown in Figure 10AB in the main text is totally due to the high-order statistics of the spike patterns. From Supplementary Figure 2A, we see that when τ I d changes, the firing rate fluctuates within 0.5% of its mean value, which is significantly smaller than the relative changes of efficacy variability with τ I d ( Figure 10AB in the main text).  Figure  3 Technical details of fitting c dif f v and c drif tv As we explained in Section 2.4 of the main text, we estimated the strengths of DiffV and DriftV by fitting the formula using the recorded evolution of efficacy variance v(t) with time. In this section, we explain the technical details of this fitting.

The starting time for fitting
Under STDP, a synapse is depressed when a pre-synaptic spike gets the axonal terminal, and is potentiated when a post-synaptic spike is back-propagated to the dendritic end. Therefore, the weight of the synapse undergoes a jump process driven by the pre-and post-synaptic spike trains. When the evolution time t is large, the number of jumps  is also large. In this case, DiffV ∝ t because of the uncorrelated part of these jumps (and central limit theorem), while DriftV ∝ t 2 because of the correlated part of these jumps. However, if t is small, the number of jumps is also small, then such relationships may not emerge. Because of this, it may not be suitable to fit the formula using the efficacy variances recorded at the very beginning of the simulations.
To estimate the suitable time t to start with, we plot the evolution of E T (Var (ab) (∆w ab ))/t with time t under several spike patterns (Supplementary Figure 3), with Var (ab) (∆w ab ) being the variance of synaptic changes within the excitatory population, and E T representing trial average. According to eq.9, we should expect that E T (Var (ab) (∆w ab ))/t is linear with time t. However, we found that this linearity is only significant for sufficiently large t, especially for spike pattens under which DriftV is small (see Supplementary Figure 3A2,A3,B2). After considering this fact, we chose t = 15s as the time when we started to use the recorded efficacy variances for our fitting.

The goodness of fit
In our simulations, we recorded the efficacy variance at a given time for 32 trials. Because of this, we use lack-of-fit F-statistic (Walpole et al, 2012) to assess the goodness of fit.
F-statistic is defined by with v i (t) being the variance of synaptic changes at time t in the ith simulation trial, v(t) being the average variance over all the trials, v(t) being the variance predicted by the model, N trial = 32 being the number of trials, and N time = 25 being the number of time points at which we recorded the variances. F-statistic quantifies the fit error per degree of freedom caused by the lack of fit of the model relative to the experimental error (here, the trial-to-trial variability). Small F-statistic means that the lack of fit of the model is weaker than the experimental error, which suggests that the proposed model well explains the experimental data; while large F-statistic has the opposite meaning.
In Supplementary Figure 4A1, B1, C1, we exemplify the F-statistics under the original spike patterns, P SS in asynchronous states and P T SiE in synchronous states; we see that most of the F-statistics are very small. For the highest F-statistics under each type of patterns, we also calculate the p-values for the null hypothesis that "there is no lack of linear fit" (Walpole et al, 2012), and get high p-values in all cases. These facts suggest that our model fits the data quite well. We also draw the scatter plot and the fitting line for the lowest p-value cases (i.e., highest F-statistic cases) under each type of patterns in Supplementary Figure 4A2-A3, B2, C2, exemplying the worst fittings.

AT WithinEvent in P T SiE , P T SiE+SSiE and P T SiE+N RC
In the main text, we use CV W ithinEvent to quantify AT WithinEvent , which is the mean CV value of the pieces of spike trains that contain more than 2 spikes in synchronous events. However, if a spike train contains only a small number of spikes, its CV may also depends on the number of spikes it contains. For example, the CV of a long Poisson spike train with stationary firing rate is typically 1, but if we calculate the CV of a short piece of this spike train which contains, say, 3 or 4 or 5 spikes, we will find that its CV is typically smaller than 1. To get more understanding on the AT WithinEvent in P T SiE , P T SiE+SSiE and P T SiE+N RC , we will discuss CV as a function of the spike number n in a synchronous event under these three patterns in this section.
We found that in all these three patterns, CV (n) tends to increase with n (Supplementary Figure 5, upper panels). CV in P T SiE is smaller than that in P T SiE+SSiE for all n values, which reflects the increase of the burstiness of AT WithinEvent caused by SSiE. We also found that distribution of the spike number n is narrow under P T SiE and P T SiE+SSiE , and gets broadened in P T SiE+N RC (Supplementary Figure 5, middle and lower panels). However, for the n values whose probability is nonzero under P T SiE+SSiE , its CV (n) is almost the same as the CV (n) under P T SiE+N RC (Supplementary Figure 5,   In A1, we show the F-statistics for different τ I d s, and also the p-values for the null hypothesis that "there is no lack of linear fit" for the two τ I d s with highest F-statistics. In A2 and A3, we draw the scatter plots and the fitting lines for these two τ I d s to exemplify the worse fittings. (B1-B2) The F-statitics, p-value, scatter plot and fitting line under P SS in asynchronous states. (C1-C2) The F-statistics, p-value, scatter plot and fitting line under P T SiE in synchronous states. The CV is only defined when n ≥ 3, and the probability that n ≥ 5 in P T SiE and P T SiE+SSiE is zero, so CV (n) is only meaningful when n = 3, 4 in P T SiE and P T SiE+SSiE . Note that for the meaningful n values, the CV in P T SiE+SSiE is very close to that in P T SiE+N RC . Middle: The probability distribution of n under P T SiE+SSiE . We only shows the values when n ≥ 3. And this distribution under P T SiE+SSiE is exactly the same with that under P T SiE . Lower: The distribution of n under P T SiE+N RC . (B) The same as A, except that τ I d = 11ms. (C) The same as A, except that τ I d = 14ms. Simulation details are the same as Figure 9 in the main text. Error bars are now shown for simplicity. 5 Comparing the results of this paper with the results of our previous paper In this paper, we investigate the influence of spike pattern structure onto efficacy variability by implementing spike shuffling methods onto spike patterns self-organized by a LIF networks; in our previous paper (Bi and Zhou, 2016), we investigated a similar problem using spike generating models on converging motifs. In this section, we briefly review the results of our previous paper, and compare them with the results of this paper. We find that the results of this paper provide further validation to some of the results of our previous paper.
For the convenience of the following discussions, we call the neuron that receives from many neurons in a converging motif to be the central neuron, and call the other neurons non-central neurons.

Results of the previous paper:
If spike trains of all the neurons (including the non-central and central neurons) in a converging motif are stationary (i.e., firing rates do not change with time) with equal firing rate, then DiffV is smallest when CV is within the range 0.3 ∼ 0.7, and gets large either when the spike trains get burstier or more regular. If different non-central neurons have different firing rates, then DiffV tends to monotonically increases with CV .

Results of this paper:
For asynchronous spike patterns, we compared P T M with P T M +SS (Term 3 in Section 3.3 of the main text).
We found that CV > 1 in P T M , but SS reduces CV to almost 1. At the same time, c dif f v is larger under P T M than under P T M +SS , which suggests that the burstiness of stationary spike trains increases DiffV. This is consistent with our previous result in the CV > 1 case.

Results of the previous paper:
For stationary spike trains, introducing heterogeneity of rates to non-central neurons does not strongly influence DiffV when the spike trains are bursty, but discounts the increase of DiffV with the decrease of CV under strong regularity.

Results of this paper:
For asynchronous spike patterns, we compared P T M +SS with P T M +N RC (Term 2 in Section 3.3 of the main text). We found that CV ≈ 1 in both these two patterns, and that c dif f v is almost the same under these two patterns, which suggests that heterogeneity of rates does not influence DiffV when CV is around 1. This is consistent with our previous result in the burstiness case.

Synchronous firing
Results of the previous paper: Synchronous firing usually increases DiffV, except for spike-to-spike synchrony and good temporal precision.

Results of this paper:
For asynchronous spike patterns, we compared P SS with P SS+T M to investigate the effect of synchronous firing (Term 1 in Section 3.3 of the main text). Here, for asynchronous spike patterns, synchronous firing means the fluctuation of population firing rate with time. We found that c dif f v is larger under P SS than under P SS+T M , which suggests that synchronous firing facilitates DiffV. This is consistent with our previous result.
For synchronous spike patterns, we compared P SSiE with P SSiE+ST R , and found that c dif f v is larger under P SSiE than under P SSiE+ST R (Term 1 in Section 3.4 of the main text). In our spike patterns, there is no spike-to-spike synchrony, so this result is also consistent with our previous result.

Auto-correlation structure under synchronous firing
Results of the previous paper: Under synchronous firing, auto-correlation structure may have a number of types, such as AT WithinEvent , AT SpikeNum and AT events (see Figure 4 in the main text). Burstier spike trains within synchronous events (for AT WithinEvent ), broader distribution of the spike numbers a neuron fires in different synchronous events (for AT SpikeNum ), and burstier occurrences of synchronous events (for AT events ) tend to increase DiffV. Results of this paper: 1) To investigate the effect of AT WithinEvent , we compared P T SiE with P T SiE+SSiE (Term 3 in Section 3.4 of the main text). We found that pieces of spike trains within synchrous events are burstier under P T SiE+SSiE than under P T SiE , and c dif f v is also larger under P T SiE+SSiE than under P T SiE . This is consistent with our previous result.
2) To investigate the effect of AT SpikeNum , we compared P T SiE+SSiE with P T SiE+N RC (Term 4 in Section 3.4 of the main text). We found that the spike number distribution is broader under P T SiE+N RC than under P T SiE+SSiE , and c dif f v is also larger under P T SiE+N RC than under P T SiE+SSiE . This is consistent with our previous result.
3) To investigate the effect of AT events , we compared P SSiE with P SSiE+ET S (Term 5 in Section 3.4 of the main text). We found that the occurrence of synchronous events is burstier under P SSiE+ET S than under P SSiE , but c dif f v is almost the same under these two patterns. However, we argue here that this may be because that the contribution to DiffV from the interactions between spikes in different synchronous events are small due to the far separation between synchronous events and the decaying STDP time window.
To study the effect of AT events onto DiffV in more details, we compare c dif f v under P T SiE+SSiE with that under P T SiE+SSiE+ET S . Here we first implement TSiE to homogenize firing rate, so that DriftV = 0 in both these patterns. In this case, the efficacy variability is totally caused by DiffV, so that the estimated value of c dif f v will become more precise. We find that c dif f v under P T SiE+SSiE+ET S is larger than that under P T SiE+SSiE when τ I d ≤ 10ms, but there is hardly any difference when τ I d ≥ 11ms (Supplementary Figure 6A). This is possibly because that when τ I d is large, the frequency of the occurrence of synchronous events becomes lower, which separates two adjacent synchronous events farther away in time (see Figure 3A in the main text). This weakens the STDP interactions between the spikes in adjacent synchronous events, so that changing AT events does not significantly influences the synaptic changes. To validate this idea, we enlarge the STDP time scale τ ST DP (eq. 5 in the main text) from the default value 20ms to 50ms, to increase the interactions between spikes in adjacent synchronous events. Indeed, we found that in this case, c dif f v under P T SiE+SSiE+ET S becomes significantly larger than that under P T SiE+SSiE (Supplementary Figure 6A, inset). This validates the result in our previous paper (Bi and Zhou, 2016) that the burstiness of the occurrence of synchronous events increases DiffV.  Figure 9 in the main text.

Auto-correlation structure in stationary spike trains
Results of the previous paper: Auto-correlation structure in stationary spike trains does not influence DriftV.

Results of this paper:
For asynchronous spike patterns, we compared P T M with P T M +SS (Term 3 in Section 3.3 of the main text).
We found that CV > 1 in P T M , but SS reduces CV to almost 1. However, we found that c drif tv is almost the same under these two patterns, which suggests that the burstiness of stationary spike trains does not change DriftV. This is consistent with our previous result.

Results of the previous paper:
Under the absence of heterogeneity of cross-correlations, the interaction of heterogeneity of rates and P-D imbalance induces DriftV.

Results of this paper:
For asynchronous spike patterns, we compared P T M +SS with P T M +N RC (Term 4 in Section 3.3 of the main text). We found that c drif tv is always near to zero under P T M +N RC (whose firing rates are homogeneous), and is positively correlated with the strength of P-D imbalance under P T M +SS (whose firing rates are heterogeneous).
This shows that DriftV can be induced by the interaction of heterogeneity of rates and induces DriftV, which is consistent with our previous result.
For synchronous spike patterns, we compared P SSiE with P SSiE+T SiE (Term 2 in Section 3.4 of the main text). We found that c drif tv is always near to zero under P SSiE+T SiE (whose firing rates are homogeneous), and is positively correlated with the strength of P-D imbalance under P SSiE (whose firing rates are heterogeneous). This also shows that DriftV can be induced by the interaction of heterogeneity of rates and induces DriftV, which is again consistent with our previous results.

Synchronous firing
Results of the previous paper: Synchronous firing can influence P-D imbalance |E a (v a )|, thereby influencing DriftV under heterogeneity of rates.

Results of this paper:
For asynchronous spike patterns, we compare P SS with P SS+T M (Term 1 in Section 3.3 of the main text). We found that when A p = A d , both c drif tv and |E a (v a )| are larger than zero in P SS but almost zero in P SS+T M , which suggests that synchronous firing influences DriftV under heterogeneity of rates by changing P-D imbalance. This is consistent with our previous result.
For asynchronous spike patterns, we compare P SSiE with P SSiE+ST R (Term 1 in Section 3.4 of the main text). We found that when A p = A d , both c drif tv and |E a (v a )| are larger than zero in P SSiE but almost zero in P SSiE+ST R , which also suggests that synchronous firing influences DriftV under heterogeneity of rates by changing P-D imbalance. This is again consistent with our previous result.

Auto-correlation structure under synchronous firing
Results of the previous paper: As for the three types of auto-correlation structure under synchronous firing ( Figure 4 in the main text), AT WithinEvent and AT SpikeNum do not change P-D imbalance, thereby nor DriftV under heterogeneity of rates; AT events changes P-D imbalance, thereby changing DriftV under heterogeneity of rates. The influence of AT events onto P-D imbalance is complicated, but if τ delay τ ST DP and synchronous events do not overlap with each other, E a (v a ) will increase (or decrease) with the burstiness of the occurrence of synchronous events if A p > A d (or Results of this paper: 1) To investigate the effect of AT WithinEvent , we compared P T SiE with P T SiE+SSiE (Term 3 in Section 3.4 of the main text). We found that P-D imbalance is almost the same in these two patterns, which suggests that AT WithinEvent does not change DriftV under heterogeneity of rates. This is consistent with our previous result.
2) To investigate the effect of AT SpikeNum , we compared P T SiE+SSiE with P T SiE+N RC (Term 4 in Section 3.4 of the main text). We found that P-D imbalance is almost the same in these two patterns, which suggests that AT SpikeNum does not change DriftV under heterogeneity of rates. This is consistent with our previous result.
3) To investigate the effect of AT events , we compared P SSiE with P SSiE+ET S (Term 5 in Section 3.4 of the main text). We carefully designed the method of ETS so that synchronous events do not overlap with each other in the spike patterns treated by ETS. We found that P SSiE+ET S has larger CV events value than P SSiE , and E a (v a ) in P SSiE+ET S is larger (or smaller) than that in P SSiE . This suggests that E a (v a ) will increase (or decrease) with CV events if A p > A d (or A p < A d ), which is consistent with our previous result.

Results of the previous paper:
Heterogeneity of cross-correlations induces DriftV together with heterogeneity of rates (eq.6).

Results of this paper:
For asynchronous spike patterns, we compared P ST R and P ST R+T M (Term 2 in Section 3.3 of the main text). In this two patterns, potentiation and depression are almost balanced with each other, and P ST R has larger heterogeneity of cross-correlations than P ST R+T M . We found that c drif tv is larger in P ST R than in P ST R+T M . This is consistent with the result derived from eq.8.
For synchronous spike patterns, we compared P W T S and P SSiE (Section 3.5 of the main text). We found that DriftV in P SSiE decreases with τ I d , while both DriftV and heterogeneity of cross-correlations increase with τ I d , which suggests that heterogeneity of cross-correlations is the main reason of DriftV when the network works in synchronous bursting states.
6 Several other topics 6.1 The change of efficacy variability with τ I d in the original spike patterns The change of efficacy variability with τ I d in the original spike patterns is shown in Figure 7. Τ d I ms Figure 7: Efficacy variability as a function of τ I d in the original spike patterns. The dashed vertical line indicates the transition between asynchronous state and synchronous state. We used the synaptic strengths at the end of the STDP evolution to calculate the efficacy variability. Simulation details are the same as those in Figure  10 in the main text.

The influence of STR onto CV and heterogeneity of cross-correlations in asynchronous states
We compare the CV and index of heterogeneity of cross-correlations (HCC) under P W T S and P ST R in Supplementary Figure 8. We see that STR slightly reduces both of these two quantities.  Figure 9: The efficacy variability caused by the spike patterns after WTS (blue) and after EOS (red) for synchronous states.

The influence of EOS onto the efficacy variability
EOS hardly changes the efficacy variability in synchronous spike patterns (Supplementary Figure 9), which suggests that the dependences (on, say, spike times or spike numbers) between the pieces of spike trains in adjacent synchronous events do not have significant influence onto the efficacy variability our model.