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

^{1}State Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing, China^{2}Department of Physics, Hong Kong Baptist University, Kowloon Tong, Hong Kong^{3}Centre for Nonlinear Studies, Beijing-Hong Kong-Singapore Joint Centre for Nonlinear and Complex Systems, Institute of Computational and Theoretical Studies, Hong Kong Baptist University, Kowloon Tong, Hong Kong^{4}Beijing Computational Science Research Center, Beijing, China^{5}Research Centre, Hong Kong Baptist University Institute of Research and Continuing Education, Shenzhen, China

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).

## 1. Introduction

Variability is a prominent feature of the neuronal activities. The neurons in the same population may respond quite differently to the same stimulus (structural variability), and the responses of a neuron to the same stimulus can also differ in different trials (trial variability). Structural variability comes from the heterogeneity of neuronal responsive properties and the randomness of inter-neuronal connections. It is found that even the same type of neurons may have different responsive properties due to the difference in the gene expression of membrane ion channels (Schulz et al., 2006; Padmanabhan and Urban, 2010); and the strengths of synapses may span several magnitudes (Song et al., 2005; Buzsáki and Mizuseki, 2014), continuously changing with time (Zucker and Regehr, 2002; Keck et al., 2008). Trial variability partly comes from biomolecular noises such as the open and close of ion channels and the release of synaptic vesicles (see Faisal et al., 2008 for review). Such noises may enter any stage of information processing in the brain, from perception and decision making to motion generation (Faisal et al., 2008), influencing the reliability and timing of action potentials (Allen and Stevens, 1994; Zador, 1998; Dorval and White, 2005; Faisal and Laughlin, 2007), especially in neurons with thin axons (Faisal et al., 2005). Additionally, a neuronal network may have internal states such as slow synaptic currents, the strengths of the synapses under short-term plasticity, and the phase of the internal oscillation (Mongillo et al., 2008; Buonomano and Maass, 2009; VanRullen et al., 2011); and its response may vary depending on these internal states (Cohn et al., 2015; Daie et al., 2015), combined with being modulated or gated by inputs from the other brain areas (Masquelier, 2013; Lin et al., 2015). If the dynamics of the network exhibits deterministic chaos, such as theoretically suggested for the networks under excitatory-inhibitory balanced state (van Vreeswijk and Sompolinsky, 1998; Monteforte and Wolf, 2010; Ostojic, 2014), then the responsive variability will be exacerbated due to the high sensitivity to noises and initial conditions.

The nervous system is able to adapt its response to external stimuli by changing the strengths of its synapses. As the synaptic changes depend on the spike timings of the pre- and post-synaptic neurons (Dan and Poo, 2006; Caporale and Dan, 2008; Markram et al., 2012), the variability of neuronal activities must result in variability of synaptic changes. Besides, synaptic plasticity is also influenced by other factors such as neuromodulators (Bailey et al., 2000; Berke and Hyman, 2000) and series of pre- and post-synaptic biomolecular mechanisms (Graupner and Brunel, 2010; Yang and Calakos, 2013), and variability may come into any of these factors. Overall, synaptic plasticity is a “noisy” process, and we call the variability of the synaptic changes during plasticity induced by the stochasticity and randomness of neuronal networks to be the *efficacy variability*. See Section 2.1 for more discussions on the definition of efficacy variability.

In this paper, we focus on discussing the influence of variability of neuronal activities on efficacy variability, in the context of spike-timing dependent plasticity (STDP) (Gerstner et al., 1996; Dan and Poo, 2006; Caporale and Dan, 2008; Markram et al., 2012). Under STDP, synaptic plasticity is driven by spike trains, so the efficacy variability is caused by the temporal stochasticity and spatial randomness of spike trains. Additionally, spike trains may exhibit a variety of statistical features, which form rich spike pattern structures. Groups of neurons may spurt firing activity (*synchronous firing*) (Kamioka et al., 1996; Buzsáki and Draguhn, 2004; Bartos et al., 2007), the spike train of a single neuron can be bursty or regular (*auto-correlation structure*) (Softky and Koch, 1993; Schwindt and Crill, 1999; Jacob et al., 2012), firing rates of cortical neurons are typically heavily-skewed distributed *in vivo* (*heterogeneity of rates*) (Shafi et al., 2007; O'Connor et al., 2010; Buzsáki and Mizuseki, 2014), and the spike trains of different neurons also display rich inter-dependences (*heterogeneity of cross-correlations*) (Funahashi and Inoue, 2000; Schneidman et al., 2006; Ostojic et al., 2009; Trousdale et al., 2012). See Figure 1A for the concepts above. Under STDP, synaptic plasticity is induced by spike trains, so these spike pattern structures must have strong influence on efficacy variability, inducing neuronal networks with sharply different structures even under the same population rate. How these spike pattern structures may influence the efficacy variability is the topic we are interested in.

**Figure 1. Schematic of the key concepts in our modeling work**. **(A)** The four aspects of pattern structure studied in this paper. “Synchronous firing” typically means the spurt of firing activity of a population; in this paper, it also represents the time fluctuation of the population rate in asynchronous spike patterns. For asynchronous spike patterns, “auto-correlation structure” reflects the burstiness/regularity of the spike trains, which is quantified by coefficient of variance (*CV*) in this paper. Here, by “burstiness,” we typically mean the irregular structure of spike trains, instead of the regular burstiness in the spike patterns of, say, central pattern generator. For spike trains in synchronous states, we consider three types of “auto-correlation structure” to reflect the burstiness/regularity features of the spike patterns (see Figure 4). By “heterogeneity of rates,” we mean that the time-averaged firing rates are different for different neurons. By “heterogeneity of cross-correlations,” we mean that different pre-synaptic neurons of a neuron tends to fire spikes at different times relative to the spikes of the neuron. For example, in the right-bottom subplot, before a spike of neuron 2, neuron 1 tends to fire before neuron 3. **(B)** The STDP time window used in our work. Note that the axons in our work have time delay τ_{delay}, and the synapses are updated according to the spike time of the post-synaptic neuron and the time that the pre-synaptic spike arrives at the terminal. The STDP updatings of all spike pairs are summed together. **(C)** Synaptic homeostasis. The synapses input to a neuron are subject to a bound on their mean strength: when their mean strength is different from this bound, all the incoming synapses of that neuron will undergo an adjustment. **(D)** Converging motif, on which we conducted all the simulations in the previous paper (Bi and Zhou, 2016). Modeling details are presented in Section 2.

Theoretical and experimental results suggest that the dynamic pattern of a neuronal network during plasticity may strongly influence the functional performance of the resulting network after plasticity by influencing the efficacy variability. For example, large efficacy variability may blur the connection patterns that are required to successfully embed memory into the network (Figure 2A), thereby destroying memory. Experimentally, it is found that gamma oscillations are important for memory formation under normal physiological conditions (Sederberg et al., 2007; Jutras et al., 2009; Yamamoto et al., 2014); but too strong synchrony, such as that in epilepsy (Gulyás and Freund, 2015), is instead detrimental to memory (Butler and Zeman, 2008). These observations suggest that efficacy variability may get its smallest value under weak synchrony, and get larger at both asynchronous and strong synchronous states. As another example, the efficacy variability may influence the competition-and-elimination process of the synapses under development (Cancedda and Poo, 2009), resulting in networks with different sparsities and synaptic strengths (Figure 2B). Experimentally, it is found that if the spontaneous activity of the medial nucleus of the trapezoid body (MNTB) in the auditory pathway during early development is modified using genetic methods, then its feedforward projection to the lateral superior olive (LSO) becomes denser and weaker, which hinders the refinement of receptive fields of the LSO neurons, thereby destroying the hearing of the animal (Clause et al., 2014). This suggests that the genetically modified spontaneous activity in Clause et al. (2014) induces smaller efficacy variability than the normal one.

**Figure 2. Biological implications of efficacy variability**. **(A)** A network of excitatory neurons stores a memory using its attractor dynamics after the intra-connections within a sub-population (here, neurons 1–4) are strengthened (the inhibitory population that keeps the total activity of the network is not shown). When the efficacy variability is small (upper row), this subpopulation will exhibit persistently high activity if a sufficiently large number of neurons in the subpopulation have high activities initially, so the memory is retrieved. When the efficacy variability is large (lower row), this memory retrieval will fail even if the mean strength of the intra-connections (red) is stronger than that of the other ones (blue). The widths of arrows represent synaptic strengths. The subplots on the right represent the weight distributions of the blue and red synapses shown in the subplots on the left. **(B)** Efficacy variability causes different network structures by controlling the degree of synaptic competition. When the efficacy variability is small (upper), only a few synapses are weaker than the elimination threshold (black dashed vertical line) and get eliminated during neural development, so most synapses are left and their strengths tend to be uniform; when the efficacy variability is large (lower), more synapses are eliminated, and the left ones are more heterogeneous and also stronger than the upper case on average. Dashed arrows represent eliminated synapses. The subplots on the right represent the weight distributions of the synapses before the elimination process.

In our previous paper (Bi and Zhou, 2016), we studied how the four aspects of spike pattern structure shown in Figure 1A influence the efficacy variability under a conventional pair-wise additive STDP (Figure 1B) using converging motifs (Figure 1D). Additionally, we also added synaptic homeostasis which conserves the mean strength of the synapses input to a neuron (Figure 1C), which may be physiologically used to maintain the activity level in a plastic network (Turrigiano and Nelson, 2004; Turrigiano, 2011). In that paper, we first generated spike patterns using statistical models with tunable parameters, and then investigated how the efficacy variability of the converging motifs would change if the neurons fire according to the generated spike patterns with different statistics. We separated the efficacy variability (TotalV, short for “total variance”) into two parts:

with DriftV (short for “drift variance”) being the drift part induced by the heterogeneity of change rates of different synapses caused by the spatial heterogeneity of spike trains (mainly related to heterogeneity of rates and heterogeneity of cross-correlations shown in Figure 1A), and DiffV (short for “diffusion variance”) being the diffusion part induced by the weight diffusion caused by stochasticity of spike trains (mainly related to synchronous firing and auto-correlation structure shown in Figure 1A). Our main conclusions are that (1) synchronous firing generally increases DiffV, except for spike-to-spike synchrony with good temporal precision, (2) burstiness of auto-correlation structure tends to increase DiffV, (3) heterogeneity of rates induces DriftV when potentiation and depression in STDP are not balanced, and (4) heterogeneity of cross-correlations induces DriftV together with heterogeneity of rates.

However, the research strategy of our previous paper has its limitations. For example, the spike patterns were generated by statistical models, which prevents us from understanding the contributions of different pattern statistics to the efficacy variability under biologically more plausible spike patterns. Additionally, in practice, people may want to know how the spike trains experimentally observed may influence the efficacy variability of a local neural circuit, thereby understanding the functional meanings of the statistical features of the spike trains during learning or neural development. Therefore, it is desirable to develop an approach to manipulate the spike pattern statistics in the recorded spike patterns. In this paper, we propose to use spike shuffling methods to solve this problem.

Spike shuffling methods are commonly used experimental techniques to destroy inter-spike, inter-neuron or inter-trial dependencies of spike patterns, thereby establishing significance of dependencies. It has been used when, for example, studying functional interactions between neuronal population (Narayanan and Laubach, 2009), investigating replay of spike sequence (Nádasdy et al., 1999; Ji and Wilson, 2007), identifying spatio-temporal correlations on the background of noises (Amarasingham et al., 2012), and discovering information-containing spatio-temporal correlations in neural codes (Panzeri et al., 2002; Nirenberg and Latham, 2003; Ganmora et al., 2011). For example, in Ji and Wilson (2007), to validate the replay of hippocampal and cortical neuronal activities during sleep to the pre-sleep activities, the authors compared the pre-sleep neuronal firing sequences to those during sleep as well as to the shuffled during-sleep sequences in which the neuronal firing orders are randomized. They found that the pre-sleep sequences have a higher similarity with the original during-sleep sequences than with the shuffled ones. As another example, in practice, correlations between neurons may come either from the external stimuli or from the inter-neuronal connections. In Ganmora et al. (2011), the authors assessed these two contributions by comparing the inter-neuronal correlations in the original spike pattern and in the spike pattern after trial-shuffling: here, trial-shuffling means that the authors repeated the stimulus for several trials, and different neurons in the shuffled pattern fired according to the spike trains in different trials. They found that inter-neuronal connections contribute significantly to the high-order neuronal correlations.

Spike shuffling methods may provide a good opportunity for understanding the influence of spike pattern statistics to synaptic plasticity under experimental conditions. For example, people can treat the recorded spike patterns using a spike shuffling method, to evidently change a specific pattern statistical feature while keeping the others largely intact, and then get understanding on the impact of the statistical feature on plasticity by comparing the synaptic changes after optically stimulating the neuronal population according to the spike patterns before and after shuffling. For pre-synaptic neurons, people may straightforwardly stimulate their axons. For post-synaptic neurons, to observe the synaptic changes caused by the post-synaptic neuronal activities controlled by optical stimulations instead of evoked by synaptic couplings, people may have to inject perisomatic shunting inhibition while at the same time stimulate the dendritic arbors to mimick the back-propagated action potentials caused by the imaginary firing of the post-synaptic neurons. Recent progress on the spatio-temporal precision of optogenetics prospects the feasibility of this operation simultaneously onto a number of neurons at present time or in the near future (see e.g., Fenno et al., 2011; Peron and Svoboda, 2011; Hochbaum et al., 2014).

To manifest this idea and also get understanding on the efficacy variability in biologically plausible spike patterns, in this paper, we implement a variety of spike shuffling methods to the spike patterns self-organized by an excitatory-inhibitory LIF network model. It is known that the LIF network can generate synchronized oscillations through two mechanisms (Tiesinga and Sejnowski, 2009; Buzsáki and Wang, 2012), which may depend on the excitatory and inhibitory synaptic time scale τ_{E} and τ_{I} (Brunel and Wang, 2003). When τ_{E} ≪ τ_{I}, oscillations come from the interaction of the excitatory and inhibitory populations (E-I mechanism): neuronal firing is first driven up by fast excitation, and then dragged down by slow inhibition; when the excitation driving the interneurons wanes, the network recovers from inhibition and the next oscillation cycle starts. When τ_{E} ≫ τ_{I}, oscillations come from the interaction within the inhibitory population (I-I mechanism): in this case, the excitation driving the inhibitory interneurons can be regarded as constant in a relatively small time scale, and the interneurons may synchronously oscillate due to the interactions among them (Wang and Buzsáki, 1996; Brunel and Hakim, 1999), which in turn entrains the excitatory population into oscillation. In this paper, we focus on the synchronized oscillation induced by E-I mechanism, which is suggested to be the mechanism of the fast oscillation in the cortex (Salkoff et al., 2015), the hippocampal ripple oscillations (Stark et al., 2014), and the spike-and-wave electroencephalography (EEG) pattern observed in absence seizures (Destexhe, 2007, 2008). To do this, we fixed the time scale of the excitatory synapses, and increased the decay time scale ${\tau}_{d}^{I}$ (see Equation 8) of the inhibitory synaptic currents; and then found this network transits from asynchronous state to weak synchronous state and at last to synchronous bursting state (Figure 3A). In reality, the dynamics of a plastic network co-evolve with the synaptic weights. To only investigate the influence of the network dynamics onto the efficacy variability without worrying about the feedback to network dynamics from synaptic changes (Figure 3B), we take the following stategy in this paper (Figure 3C): We first record the spike patterns of the excitatory population of the LIF network, then shuffle the recorded spike patterns using different methods to change different statistical features, and at last evolve the E-E links under STDP (Figure 1B) and synaptic homeostasis (Figure 1C) when the excitatory population are supposed to fire according to the recorded or shuffled spike patterns. By comparing the statistics of the patterns as well as the efficacy variability under the patterns before and after implementing a spike shuffling method, we can gain understanding on how different aspects of the pattern structure may influence the efficacy variability.

**Figure 3. Overview of our research. (A)** Spike patterns of our LIF network (Section 2.3) at asynchronous (left, ${\tau}_{d}^{I}=3\mathrm{\text{ms}}$), weak synchronous (middle, ${\tau}_{d}^{I}=7\mathrm{\text{ms}}$) and synchronously bursting (right, ${\tau}_{d}^{I}=14\mathrm{\text{ms}}$) states. In weak synchronous states (middle), a neuron usually fires no more than one spike in a synchronous event; and if a neuron fires in a synchronous event, it may be silent in another one. In synchronously bursting states (right), a neuron typically fires more than one spikes burstly in a synchronous event. **(B)** In a plastic network, dynamics and synapses interact and co-evolve with each other. We would like to cut this loop to only investigate the influence of dynamics onto efficacy variability under STDP and synaptic homeostasis without worrying about the change of dynamics caused by synaptic changes. **(C)** To achieve the strategy shown in **(B)**, we first record the spike patterns of the LIF network when the synapses are fixed, then shuffle the spike patterns with a variety of methods to change specific pattern statistics, and at last let the neurons in the network fire according to the recorded or shuffled patterns with STDP and synaptic homeostasis being imposed onto the synapses.

We have three aims in this paper. Firstly, we would like to develop a systematic spike-shuffling approach to alter the statistical features of a given spike pattern. Secondly, we will apply this approach to investigate how different pattern statistics (Figure 1A) influence the efficacy variability under STDP and synaptic homeostasis (Figures 1B,C) in the spike patterns self-organized by the LIF network. Thirdly, we will estimate the contributions of different pattern statistics to the efficacy variability under the spike patterns generated by our LIF network, both in asynchronous and synchronous states. For example, we would like to know which spike pattern statistics influences significantly to the efficacy variability, and which statistics is the main reason for the change of the efficacy variability with ${\tau}_{d}^{I}$.

In the Results section, we will first illustrate the statistics of the spike patterns generated by the LIF network (Section 3.1), and then study the impact of the four statistical features (Figure 1A) onto the efficacy variability by implementing spike shuffling methods onto the spike patterns in both asynchronous and synchronous states (Sections 3.2–3.4). Finally, we will understand the contributions of different pattern statistics to the efficacy variability under these spike patterns (Section 3.5). In the Discussion section, we will discuss the connections of our results to known theoretical and experimental results. We compare the results of this paper with the results of our previous paper (Bi and Zhou, 2016) in Supplementary Materials Section S5, and find their consistency.

## 2. Materials and Methods

### 2.1. The Definition of Efficacy Variability

Suppose the weights of a set of synapses ${W}$ are to be changed by almost the same value during a plasticity process to make the network function normally. We run the plasticity process on the network for several trials, and construct a matrix Δ**W**, each column of which represents the synaptic changes of ${W}$ in one trial at a given time point *t*, and different columns represent different trials. To quantify the variability of synaptic changes during plasticity, we define *efficacy variability* of ${W}$ at time *t* to be the variance of the elements of the matrix Δ**W**, i.e., Var_{S, T}(Δ**W**). Here, the subscript *S* represents integrating over row index, i.e., structural index, and *T* represents integrating over column index, i.e., trial index. Using the law of total variance (Weiss, 2005), it can be shown (see Section 2.1 in our previous paper Bi and Zhou, 2016) that

and that

with E and Var representing mean and variance respectively.

In Equation (2), *E*_{T}(Δ**W**) represents the trial expectations of the changes of all the synapses in ${W}$; and Var_{S}(*E*_{T}(Δ**W**)) is the variance of these trial expectations, representing DriftV. Var_{T}(Δ**W**) represents the trial-to-trial variances caused by diffusion, and *E*_{S}(Var_{T}(Δ**W**)) is the average of these variance over all the synapses, representing DiffV. This equation is the formal writing of Equation (1) in the introduction.

In Equation (3), Var_{T}(*E*_{S}(Δ**W**)) represents the trial-to-trial variability of the mean synaptic change of the whole network. But a real biological process only allows a single trial, so this trial-to-trial variability cannot contribute to biological functions except for individual differences. Additionally, in this paper, we aim to understand the influence of spike pattern structures onto efficacy variability, so the firing rates and second-order statistics of spike patterns need to be controlled in our model. In this theoretical context, under STDP, Var_{T}(*E*_{S}(Δ**W**)) is presumably of the order of ${O}$(1∕|${W}$|), with |${W}$| being the number of synapses in ${W}$. So when |${W}$| is large enough, Equation (3) becomes

This tells that we can approximate the efficacy variability by the trial average of the variance of the synaptic changes in the network, which, as shown in Figure 2, may have strong biological implications. In this paper, we use *E*_{T}(Var_{S}(Δ**W**)) to quantify the efficacy variability in our simulations. As we show in Supplementary Figure 1, Var_{T}(*E*_{S}(Δ**W**)) is negligible comparing to *E*_{T}(Var_{S}(Δ**W**)), so *E*_{T}(Var_{S}(Δ**W**)) is indeed nearly the same with the full version Var_{S, T}(Δ**W**) in our simulations.

### 2.2. STDP and Synaptic Homeostasis

The STDP updating caused by a pair of pre- and post-synaptic spike at *t*_{pre} and *t*_{post} is

with τ_{delay} being the axonal delay. The contributions of all pairs of pre- and post-synaptic spikes are added together. τ_{STDP} = 20 ms, τ_{delay} = 1 ms throughout the paper, and *A*_{p} = *A*_{d} = 1 by default.

As explained in Figures 3B,C, we did not directly embed STDP and synaptic homeostasis into the self-organizing dynamics of the LIF network, but instead evolved the E-E links under STDP and synaptic homeostasis according to the recorded or shuffled spike patterns of the LIF network with fixed synaptic weights. We recorded the variance of synaptic weights during this evolution every 1 s of biological time, and before each recording, we implemented synaptic homeostasis onto the synapses within excitatory population as

with *w*_{ab} being the weight of the synapse from excitatory neuron *b* to *a*, *N*_{a} being the excitatory in-degree of the *a*th neuron, and *w*_{bound} = 0 being the ground line of synaptic homeostasis. In this way, the mean excitatory synaptic input to each excitatory neuron was fixed at *w*_{bound} = 0 before each recording.

### 2.3. The LIF Neuronal Network

The network consists of 2000 excitatory and 500 inhibitory conductance-based LIF neurons, with the links being randomly connected with probability 0.2. Each neuron in the network also receives excitatory external inputs.

The dynamics of the sub-threshold membrane voltage ${V}_{i}^{\alpha}$ of the *i*th neuron in the αth (α = *E, I* representing excitatory or inhibitory) population is (see e.g., Brunel and Wang, 2003)

with *C*^{α} being the membrane capacity, ${g}_{L}^{\alpha}$ the leakage conductance, *V*_{leak} the leakage voltage, *g*^{α, ext} the strength of the synapses from external inputs to the neurons in the αth neuronal population, ${s}_{i}^{\alpha ,ext}$ the synaptic conductance of unit synaptic strength from the external inputs, *g*^{αβ} (β = *E, I*) the strength of the synapses from the βth population to the αth population, ${s}_{ij}^{\alpha \beta}$ the synaptic conductance of unit synaptic strength from the *j*th neuron in the βth population connected to the *i*th neuron in the αth population, ${A}_{ij}^{\alpha \beta}=1,0$ indicating whether the *j*th neuron in the βth population connects to the *i*th neuron in the αth population or not, and *E*_{E} (*E*_{I}) being the inverse voltage for the excitatory (inhibitory) synaptic current. The membrane voltage ${V}_{i}^{\alpha}$ is reset to *V*_{r} as soon as it crosses threshold θ, and will stay at *V*_{r} for a refractory period ${\tau}_{ref}^{\alpha}$ after this reset.

The dynamics of the synaptic conductance of unit synaptic strength ${s}_{ij}^{\alpha \beta}$ follows (Brunel and Wang, 2003)

with {*t*_{1}, *t*_{2}, ⋯ } being the spike train that the synapse receives, ${\tau}_{m}^{\alpha}\equiv {C}^{\alpha}\u2215{g}_{L}^{\alpha}$ the membrane time constant of the αth population, ${\tau}_{r}^{\alpha}$ (${\tau}_{d}^{\alpha}$) the rising (decaying) time scale of the synaptic conductance in response to an incoming spike, and τ_{delay} = 1 ms the axonal delay. The normalization factor $\frac{{\tau}_{m}^{\alpha}}{{\tau}_{d}^{\alpha}-{\tau}_{r}^{\alpha}}$ is chosen so that varying the synaptic time constant does not affect the time integral of a postsynaptic current (Brunel and Wang, 2003). The dynamics of ${s}_{i}^{\alpha ,ext}(t)$ is determined by the external input spike trains, with the other parameters being that same as those for ${s}_{ij}^{\alpha E}(t)$.

In our simulations, ${g}_{L}^{E}={g}_{L}^{I}=10\mathrm{\text{nS}}$, ${C}^{E}=20\mathrm{\text{ms}}\xb7{g}_{L}^{E}$, ${C}^{I}=10\mathrm{\text{ms}}\xb7{g}_{L}^{I}$; *E*^{E} = 0 mV, ${V}_{leak}={E}^{I}=-70\mathrm{\text{mV}}$; *g*^{EE} = 0.4 nS, *g*^{EI} = 5.8 nS, *g*^{IE} = 0.74 nS, *g*^{II} = 9.6 nS; *V*_{r} = −60 mV, θ = −50 mV; ${\tau}_{ref}^{E}=2\mathrm{\text{ms}}$, ${\tau}_{ref}^{I}=1\mathrm{\text{ms}}$; ${\tau}_{r}^{E}={\tau}_{r}^{I}=0.5\mathrm{\text{ms}}$, ${\tau}_{d}^{E}=4\mathrm{\text{ms}}$. The decaying time constant ${\tau}_{d}^{I}$ of all the inhibitory synapses are the same in a simulation trial, but may differ in different trials, resulting in different network dynamics (see Figure 3A). In a single trial, ${\tau}_{d}^{I}$ may take one value in the 12 integers from 3 to 14 ms. Each neuron also receives 1000 Hz external Poisson input, with external conductance *g*^{E, ext} = *c* × 0.53 nS, *g*^{I, ext} = *c* × 0.75nS, with *c* being a coefficient whose value depends on ${\tau}_{d}^{I}$. To keep the excitatory population almost at 20 Hz for different trials, we choose *c* to be 3.1459, 3.28212, 3.38699, 3.46922, 1.60585, 1.46867, 1.27884, 1.04738, 0.82432, 0.616324, 0.410436, 0.323722 for ${\tau}_{d}^{I}$ as integer values from 3 to 14 ms (see Supplementary Materials Section S2 for more details on how we chose *c*). Simulations were performed using a second order Runge–Kutta scheme with fixed time step δ*t* = 0.05 ms; and an interpolation scheme was also used for the determination of the firing times of the neurons (Hansel et al., 1998).

The purpose of this work is to understand how the dynamic patterns influence the efficacy variability, instead of how the dynamic properties change with model parameters; so averaging configurations of the random LIF networks does not help to gain more insight to the problem being addressed, only increasing complexity. Therefore, our study focused on a single typical network configuration, except that we chose different initial states and seeds of random generators for different trials, which resulted in trial-to-trial variability. We did check our results using other network configurations, and found qualitatively the same results.

### 2.4. Fitting *c*_{diffv} and *c*_{driftv}

Because of the additive nature of our STDP model, DiffV ∝ *t* and DriftV ∝ *t*^{2}. Therefore, the time evolution of the efficacy variance *v*(*t*) can be written as $v(t)={c}_{diffv}t+{c}_{driftv}{t}^{2}$ in a long run, with *c*_{diffv} and *c*_{driftv} respectively quantifying the strength of DiffV and DriftV. To estimate the values of *c*_{diffv} and *c*_{driftv}, we let the excitatory population fire according to the recorded or shuffled spike pattern, with STDP and synaptic homeostasis starting after the initial 2 s of transient period. We then recorded the efficacy variance at time 15, 16,…, 39 s after the transient period, and did linear regression using the formula *v*(*t*)∕*t* = *c*_{diffv}+*c*_{driftv}*t*. The estimated values and standard errors of *c*_{diffv} and *c*_{driftv} were then used to plot the relevant panels in Figures 8, 9. See Supplementary Materials Section S3 for more details on the fitting procedure and the goodness of fit.

### 2.5. Spike Pattern Analysis

Here are the methods we used to analyze the statistics of the spike patterns of the LIF network, both originally recorded and after being shuffled. Each trial of our simulation lasted for 41 s biological time, with the first 2 s regarded as transient period and excluded from the following analysis.

#### 2.5.1. Analyzing the Statistics of Synchronous Firing

For spike patterns in synchronous states (${\tau}_{d}^{I}\ge 7\mathrm{\text{ms}}$), *synchronous events* are numerically defined as follows: the firing rate of the excitatory population is first calculated in bins of 0.1 ms, then filtered using Gaussian window of σ_{window} = 1 ms, and *synchronous events* are then defined as the sequential bins in which the filtered rates are above a threshold 1 Hz. The size of the Gaussian window and the value of the threshold are chosen so that the duration of a synchronous event is large enough to include nearly all the spikes that spurt synchronously (see Figure 3A, middle and right panels), and at the same time as small as possible. It is possible that some spikes are not included in any synchronous event, and these spikes will be excluded from the statistical analysis of synchronous events.

The firing profiles of the neurons in synchronous events shown in Figure 10E are calculated as follows. For a synchronous spike pattern, we first order the firing rates of all the neurons (2000 in total) from low to high, and then put the 290–299th neurons into “low rate” bracket, the 1190–1199th neurons into “middle rate” bracket, 1990–1999th into “high rate” bracket, and all the neurons into “whole” bracket. Suppose *r*_{m,i,s}(*t*) to be the firing rate (calculated by counting spike numbers within bins of 0.1 ms) of the *i*th neuron in the *m*th bracket in the *s*th synchronous event, then the firing profile of the *m*th bracket is defined as

with *N*_{i} being the number of neurons in the bracket, *N*_{s} being the number of synchronous events, and *t*_{s} being the middle time (i.e., the mean time of all the spikes) of the *s*th synchronous event. The idea of this equation is to first translationally move all the synchronous events so that the middle times of these synchronous events are all located at 0, and then average the firing rates in these synchronous events for all the neurons in the same bracket.

In Figure 7A, we use *p* and τ_{cross} to quantify the mean strength and duration of synchronous events. *p* is defined as the mean spike number per neuron per synchronous event. To calculate τ_{cross}, we use the firing profile ${\stackrel{\u0304}{r}}_{whole}(t)$ defined by Equation (9), and τ_{cross} is then defined as the duration between the two times at which ${\stackrel{\u0304}{r}}_{whole}(t)$ drops below 10% of its peak value from its peak time for the first time, along both the positive and negative directions.

#### 2.5.2. Analyzing the Statistics of Auto-Correlation Structure

For spike patterns in asynchronous states (τ_{d, I} ≤ 6 ms), auto-correlation structure is quantified by the mean coefficient of variance (*CV*, which is the ratio of the standard deviation and mean of the inter-spike intervals) over all the spike trains, representing the burstiness/regularity of these spike trains.

For spike patterns in synchronous states (τ_{d, I} ≥ 7 ms), auto-correlation structure is separated into three aspects (Figure 4): (1) the broadness of the distribution of the spike numbers a neuron fires in different synchronous events (AT_{SpikeNum}), (2) the burstiness/regularity of pieces of spike trains within synchronous events (AT_{WithinEvent}), and (3) the burstiness/regularity of the occurrence of synchronous events (AT_{events}). AT_{SpikeNum} is quantified by the variance of the distribution of the spike numbers a neuron fires in different synchronous events. AT_{WithinEvent} is quantified by *CV*_{WithinEvent}, which is calculated by averaging over the *CV* values of the pieces of spike trains that contain more than 2 spikes in synchronous events. AT_{events} is quantified by the *CV* of the middle times (i.e., the mean time of all the spikes) of the synchronous events.

**Figure 4. The three types of auto-correlation structure we consider under synchronous firing. (A)** The broadness of the distribution of the spike numbers a neuron fires in different synchronous events. Note that in the left panel, a neuron fires quite different number of spikes during different synchronous events; while in the right panel, the spike numbers of a neuron during different synchronous events are almost the same. **(B)** The burstiness/regularity of the pieces of spike trains within synchronous events. **(C)** The burstiness/regularity of the occurrence of synchronous events.

#### 2.5.3. Analyzing the Statistics of Heterogeneity of Rates

Heterogeneity of rates means the heterogeneity of time-averaged firing rates for different neurons in the spike patterns. It is quantified by the variance of the time-averaged firing rates.

#### 2.5.4. Analyzing the Statistics of Heterogeneity of Cross-Correlations

We define the *unit cross-correlation* between neuron *a* and neuron *b* to be *C*_{ab}(τ) = 〈*r*_{a}(*t*)*r*_{b}(*t* + τ)〉∕(〈*r*_{a}(*t*)〉〈*r*_{b}(*t*)〉), with *r*_{a}(*t*) being the firing rate of the *a*th neuron, and 〈·〉 representing averaging over time. In this paper, heterogeneity of cross-correlations typically means that the a post-synaptic neuron has different unit cross-correlations with its different pre-synaptic neurons. It reflects that different pre-synaptic neurons tend to fire at different times relative to a post-synaptic spike.

To quantify the strength of heterogeneity of cross-correlations, we define *Index of heterogeneity of cross-correlations (HCC)* to quantify ${\mathrm{\text{E}}}_{b}({\mathrm{\text{Var}}}_{a\in \partial b}({\int}_{-\infty}^{\infty}H(\tau ){C}_{ab}(\tau )\text{d}\tau ))$, with ∂*b* representing all the pre-synaptic neurons of *b*, and *H*(τ) being the STDP time window (see Equation 5). It is estimated as follows: for a synapse from neuron *a* to neuron *b*, we denote Δ*n*_{a→b} as the synaptic change per unit time under STDP alone (without considering synaptic homeostasis), and denote $\Delta {m}_{a\to b}=\frac{\Delta {n}_{a\to b}}{{r}_{a}{r}_{b}}$ (with *r*_{a} and *r*_{b} being the time-averaged firing rates). The index of HCC is then calculated as *E*_{b}(Var_{a∈∂b}(Δ*m*_{a→b})).

### 2.6. Spike Shuffling Methods

Here we list out the spike shuffling methods we used for spike patterns in asynchronous states, and discuss how each of them changes spike pattern statistics. We use different spike shuffling methods for the spike patterns in asynchronous and synchronous states, because of the sharp difference of their pattern structure. Here we first list out the shuffling methods that are used for both asynchronous and synchronous states, then list out those only used for asynchronous states, and finally those only for synchronous states. See Figure 5 and Supplementary Movie for illustration of these spike shuffling methods.

**Figure 5. The spike shuffling methods we use to study the spike patterns of our LIF network, and their influences onto spike pattern structures**. The most left panel explains the spike shuffling methods we use to treat spike patterns, and the rest columns in the right show the influence of these spike shuffling methods onto synchronous firing (SF), auto-correlation structure (AT), heterogeneity of rates (HR) and heterogeneity of cross-correlations (HCC). We consider three types of auto-correlation structure for synchronous states AT_{SpikeNum}, AT_{WithinEvent} and AT_{events} (see Figure 4), and the auto-correlation structure under asynchronous states are shortly represented by AT_{async}. ✓ means that a shuffling method keeps a pattern structure unchanged; × means that a shuffling method *destroys* a pattern structure. Here, “destroy” has a sense of “completely randomize.” For SF, “destroy” means that there is no time fluctuation of population firing rate. For AT_{async}, it means that the spike train of the *a*th neuron can be regarded as an inhomogeneous Poisson process with rate *r*_{a}(*t*) = *r*_{a}*x*(*t*), with *r*_{a} being the time-averaged firing rate, and *x*(*t*) being the same for all the neurons. Because of the weak fluctuation of the population firing rate in asynchronous states, this also makes *CV* ≈ 1. For AT_{SpikeNum}, it means that the spike numbers of the neurons within a synchronous event follows Poisson distribution of parameter *p*, with *p* being the mean spike number per neuron within the synchronous event. For AT_{WithinEvent}, it means that the spike train of the *a*th neuron can be regarded as an inhomogeneous Poisson process with rate *r*_{a}(*t*) = *r*_{a}*x*(*t*), with *r*_{a} being the time-averaged firing rate, and *x*(*t*) being the same for all the neurons. For AT_{events}, it means that the occurrence of synchronous events can be regarded as a Poisson process. For HR, it means that the time-averaged firing rates of all the neurons are the same. For HCC, it means that the unit cross-correlations ${C}_{ab}(\tau )=\frac{\langle {r}_{a}(t){r}_{b}(t+\tau )\rangle}{\langle {r}_{a}(t)\rangle \langle {r}_{b}(t)\rangle}$ (with *r*_{a}(*t*) representing the firing rate of the *a*th neuron) are the same for different neuronal pairs. ◯ means that a shuffling method may change, but does not “completely randomizes,” a pattern structure. Note that STR completely flattens the rate fluctuation with time, so AT_{SpikeNum}, AT_{WithinEvent} and AT_{events} are not applicable to the synchronous spike patterns after STR (indicated by the squares in the figure). However, note that AT_{SpikeNum} and AT_{WithinEvent} of the spike patterns before STR can influence the burstiness/regularity of the asynchronous spike patterns after STR. All these spike shuffling methods are further illustrated in Supplementary Movie.

#### 2.6.1. Spike Shuffling Methods for Both Asynchronous and Synchronous States

##### 2.6.1.1. Whole-train Swap (WTS)

This method randomly shuffles the neuronal indexes of spike trains in the pattern. For example, if we denote ${T}$_{a} to be the spike train of the *a*th neuron, then the whole spike pattern can be denoted as a set of neuron-train pairs {(*a*, ${T}$_{a}), (*b*, ${T}$_{b}), (*c*, ${T}$_{c}), ⋯ }; then after WTS, the spike pattern may become {(*a*, ${T}$_{c}), (*b*, ${T}$_{a}), (*c*, ${T}$_{b}), ⋯ }. By definition, WTS keeps all the statistics of a spike pattern, but destroys the possible correlation between the spike trains and the structure of the underlying neuronal network. To get rid of this pattern-network coupling thereby focusing on the influence to the efficacy variability by spike pattern statistics (see the discussions in Section 3.2), we treated all the recorded spike patterns by WTS before any other shuffling method.

##### 2.6.1.2. Spike-time Rescaling (STR)

The idea of this method is that the spike times are first projected to the rescaled time defined as the accumulative function of the population firing rate *r*(*t*)

and then are projected back to the normal time using ${\Lambda}_{0}^{-1}(s)$, where Λ_{0}(*t*) is the linear function connecting (0, 0) with (*T*, Λ(*T*)), with *T* being the duration of the spike pattern. In this way, inter-spike intervals are rescaled according to the population firing rate, so that the population firing rate is kept constant in the spike pattern after shuffling. Technically, STR is realized by first ordering all the *M* spikes in the pattern, then setting the time of the *i*th spike at *iT*∕*M*. By definition, this shuffling method flattens population firing rate, while conserving the time-averaged firing rate of each neuron. As it keeps the order of spikes, the burstiness/regularity of spike trains and cross-correlations between spike trains in the original pattern can be, to some extent, kept in the pattern after shuffling, especially if the rate fluctuation in the original spike pattern is weak.

##### 2.6.1.3. Neuron Re-choosing (NRC)

In this method, each spike in the pattern is assigned to a randomly selected neuron. When the population size is large, this method makes all the neurons to fire as Poisson processes with equal time-dependent firing rate (i.e., *r*_{a}(*t*) = *r*_{b}(*t*) for two different neurons *a* and *b*). NRC destroys auto-correlation structure, heterogeneity of cross-correlations and heterogeneity of firing rates, but keeps the time fluctuation of population firing rate.

Note that under synchronous states, NRC also changes the distribution of the spike numbers a neuron fires in different synchronous events (i.e., AT_{SpikeNum}). When the population size is large, this distribution after NRC is a Poisson distribution with parameter *p*, with *p* being mean spike number per neuron within a synchronous event (i.e., the strength of the synchronous event).

#### 2.6.2. Spike Shuffling Methods Only for Asynchronous States

##### 2.6.2.1. Translational Move (TM)

In this method, each spike train is translationally moved by a random displacement, and periodic boundary condition is used to deal with the spikes which are moved out of the boundaries of time. By definition, TM keeps the auto-correlation structure of spike trains, and the time-averaged firing rate of each neuron. It flattens the cross-correlations between any pair of spike trains, thereby destroying both synchronous firing and heterogeneity of cross-correlations.

##### 2.6.2.2. Spike Swap (SS)

The idea of this method is to swap pairs of randomly chosen spikes of different neurons many times. A spike pattern can be denoted as a set of number pairs {(*a, t*_{1}), (*b, t*_{2}), (*c, t*_{3}), ⋯ }, with *a, b, c*, ⋯ being neuronal indexes, and *t*_{1}, *t*_{2}, *t*_{3}, ⋯ being spike times. Technically, SS shuffles the order of the first fields of these number pairs, so that the spike pattern after SS may be {(*b, t*_{1}), (*c, t*_{2}), (*a, t*_{3}), ⋯ }. SS does not change the spike number of a neuron, but randomizes the occurrence of these spikes. When the size of the network is large, the probability that a neuron fires near time *t* is proportional to the global firing rate at *t*, and the times of different spikes are independent with each other. Because of this, SS makes the spike trains Poisson processes, and the rate fluctuations of all these spike trains are simultaneously time-modulated: i.e., the firing rate of the *a*th spike train can be written as *r*_{a}(*t*) = *r*_{a}*x*(*t*), with *r*_{a} being the time-averaged firing rate, and *x*(*t*) being the same for all the spike trains. By definition, SS destroys auto-correlation structure and heterogeneity of cross-correlations, but keeps the time fluctuation of population firing rate and the time-averaged firing rate of each neuron.

#### 2.6.3. Spike Shuffling Methods Only for Synchronous States

##### 2.6.3.1. Spike Swap in Events (SSiE)

The idea of SSiE is to swap pairs of randomly chosen spikes of different neurons many times, with each pair of chosen spikes being within the same synchronous event. It can be understood as implementing SS (see Section 2.6.2.2) onto each synchronous event in the spike pattern. Our numeric definition of “synchronous event” is presented in Section 2.5. By definition, SSiE keeps the spike number of a neuron within a synchronous event. And the same as SS, SSiE also makes the neurons to fire like Poisson processes when the size of the neuronal population is large, and the rate fluctuation of all the spike trains are simultaneously time-modulated. By definition, SSiE destroys heterogeneity of cross-correlations and AT_{WithinEvent}.

##### 2.6.3.2. Train Swap in Events (TSiE)

The idea of TSiE is to swap the pieces of spike trains of randomly selected neurons pairs in the same synchronous event many times, and the pieces of spike pattern in different synchronous events are shuffled independently. By definition, TSiE destroys heterogeneity of rates and heterogeneity of cross-correlations, but keeps AT_{WithinEvent}. The spike number distribution per neuron per synchronous event for the whole neuronal population is kept unchanged under TSiE, but the distribution for a single neuron in different synchronous events is changed dramatically.

##### 2.6.3.3. Event-time Shuffle (ETS)

The idea of ETS is that all the spikes within the same synchronous events are translationally moved by a random displacement, at the same time (1) avoiding the overlapping of different synchronous events, and (2) keeping the order of synchronous events unchanged. Here, by “avoiding synchronous-events overlapping”, we have two meanings. Firstly, we mean that if a synchronous event ${S}$ happens in the time interval [*t*_{1}, *t*_{2}], then no other synchronous events should happen in this interval. Secondly, because of the axonal delay τ_{delay}, a post-synaptic neuron receives its pre-synaptic spikes in ${S}$ during the interval [*t*_{1} + τ_{delay}, *t*_{2} + τ_{delay}], and we would like the post-synaptic neuron to receive no spikes from another synchronous event between the time that it fires in ${S}$ and the time that it finishes receiving all its pre-synaptic spikes in ${S}$. In a word, these two conditions mean that if a synchronous event happens in the time interval [*t*_{1}, *t*_{2}], then no other synchronous event should appear within [*t*_{1}−τ_{delay}, *t*_{2}], with τ_{delay} being the axonal delay. Technically, this is realized by first randomly selecting *N*_{event} points in the duration $\left[0,T-\sum _{i=1}^{{N}_{events}}({T}_{i}+{\tau}_{delay})\right]$ (with *N*_{event} being the number of the synchronous events, *T* being the time duration of the spike pattern, and *T*_{i} being the duration of the *i*th synchronous event, see Section 2.5 for the numeric definition of “synchronous event”), then set the beginning time of the *j*th synchronous event at ${x}_{j}+\sum _{i=1}^{j-1}({T}_{i}+{\tau}_{delay})$ (with *x*_{j} being the *j*th selected points). The reason why we avoid the synchronous-events overlapping is that firstly getting rid of the first constraint above may change AT_{WithinEvent} in some synchronous events, thereby changing the efficacy variability through this side effect; and secondly, our previous paper (Bi and Zhou, 2016) shows that under the second constraint above the change of DriftV caused by *CV*_{events} is simple to understand, and we would like to focus on this simple situation here. The reason why we keep the order of synchronous events is that there may be dependencies (on, say, spike times, or spike numbers) between the pieces of spike trains in adjacent synchronous events, and we would like to keep these possible dependencies to the most extent during ETS. Note that ETS not only change AT_{event}, but may also change heterogeneity of cross-correlations by changing the STDP interactions of spikes in adjacent synchronous events.

##### 2.6.3.4. Event-order Shuffle (EOS)

EOS shuffles the occurrence order of the synchronous events, so that the synchronous events appear at the same time points as in the original pattern, but in a shuffled order. For example, if the mean spike time of the *i*th synchronous event is at ${\stackrel{\u0304}{t}}_{i}$, then all the synchronous events in the whole spike pattern may be denoted as a set of number pairs $\left\{(i,{\stackrel{\u0304}{t}}_{i}),(j,{\stackrel{\u0304}{t}}_{j}),(k,{\stackrel{\u0304}{t}}_{k}),\cdots \phantom{\rule{0.3em}{0ex}}\right\}$; after EOS, this pattern may become $\left\{(i,{\stackrel{\u0304}{t}}_{k}),(j,{\stackrel{\u0304}{t}}_{i}),(k,{\stackrel{\u0304}{t}}_{j}),\cdots \phantom{\rule{0.3em}{0ex}}\right\}$. EOS destroys possible dependencies (on, say, spike times or spike numbers) between the pieces of spike trains in adjacent synchronous events. By comparing the efficacy variability under spike patterns before and after EOS, we found that such dependencies have little influences onto efficacy variability in our model (Supplementary Figure 9).

## 3. Results

### 3.1. The Dynamic Patterns of the LIF Network

When the decaying time scale ${\tau}_{d}^{I}$ of the inhibitory synaptic inputs (Equation 8) is within the range $3\mathrm{\text{ms}}\le {\tau}_{d}^{I}\le 6\mathrm{\text{ms}}$, the LIF network (Section 2.3) works in asynchronous states (Figure 3A). In this case, we found the population firing rate exhibits strong fluctuation with time (Figure 6A, upper panels). The firing rate of individual neurons also fluctuate strongly (Figure 6A, lower panels); and the coefficient of variance (*CV*) is larger than 1 (Figure 6B), reflecting the burstiness of the spike patterns. These features suggest that the network works in chaotic asynchronous states, which may be caused by the unstability of the network dynamics to heterogeneous perturbations (Ostojic, 2014). The firing rates are heavily skewed (Figure 6C), which may be caused by the heterogeneous input degrees and nonlinear current-rate relationship of neurons (Roxin et al., 2011). There also exists non-trivial cross-correlations between neurons in these asynchronous patterns (Figure 6D), which may be caused by connectivity details such as unidirectional connections and input sharing (Ostojic et al., 2009; Trousdale et al., 2012).

**Figure 6. Statistics of the spike patterns of the LIF network in asynchronous states. (A)** Upper panels: population firing rate as a function of time when ${\tau}_{d}^{I}=3$ ms (left) and 6 ms (right). Lower panels: firing rates for three example neurons, when ${\tau}_{d}^{I}=3$ ms (left) and 6 ms (right). We estimated the instantaneous population firing rates using bins of 0.1 ms, and estimated the instantaneous firing rates of individual neurons by convolving their spike trains with a 40-ms-wide Gaussian filter. **(B)** Mean coefficient of variance (*CV*) of all the spike trains in a spike pattern as a function of ${\tau}_{d}^{I}$. **(C)** Left panel: the standard deviation of population rate as a function of ${\tau}_{d}^{I}$. Right panel: the distributions of firing rates when ${\tau}_{d}^{I}=3$ ms and 6 ms. **(D)** The index of heterogeneity of cross-correlations (HCC) (see Section 2.5) as a function of ${\tau}_{d}^{I}$. **(A)** and the right panel of **(C)** are plotted from one simulation trial, while the other panels are plotted from averaging over 32 simulation trials. Simulation details are explained in Section 2.3. Error bars represent s.e.m.

When ${\tau}_{d}^{I}\ge 7\mathrm{\text{ms}}$, the LIF network works in synchronous states, and gradually goes from weak synchronous state to synchronously bursting state with the increase of ${\tau}_{d}^{I}$ (Figure 3A). The spike pattern structure in synchronous states is more complex than that in asynchronous states, and some key points are listed below:

(1) *Synchronous firing*. Both the strength *p* and the duration τ_{cross} of synchronous events increase with ${\tau}_{d}^{I}$ (Figure 7A).

(2) *Heterogeneity of rates*. The skewness and variance of the firing rates decrease with ${\tau}_{d}^{I}$ (Figure 7B).

(3) *Auto-correlation structure*. For AT_{SpikeNum}, because of heterogeneity of rates, different neurons fire quite different number of spikes in a synchronous event; but the spike number distribution for a neuron in different synchronous events is not broad (Figure 7C, left panel). For AT_{WithinEvent}, if a neuron is to fire more than one spikes in a synchronous event, then these spikes tends to appear regularly (Figure 7C, middle panel): this is because that the refractory periods of the excitatory neurons in our model is fixed at ${\tau}_{ref}^{E}=2\mathrm{\text{ms}}$, so the neurons will fire regularly at rate close to $1\u2215{\tau}_{ref}^{E}$ if they receive strong excitatory inputs during synchronous events. For AT_{events}, synchronous events in our model tends to appear regularly (Figure 7C, right panel), which reflects the oscillation nature of the E-I circuit (Brunel, 2000; Brunel and Wang, 2003).

(4) *Heterogeneity of cross-correlations*. We found that the heterogeneity of cross-correlations tends to increase with ${\tau}_{d}^{I}$ (Figure 7D), especially when ${\tau}_{d}^{I}\ge 12\mathrm{\text{ms}}$. The reason of this phenomenon and its influence onto the efficacy variability will be discussed in **Section 3.5**.

**Figure 7. Statistics of the spike patterns of the LIF network in synchronous states. (A)** The strength *p* and duration τ_{cross} of synchronous events as functions of ${\tau}_{d}^{I}$ (see Section 2.5 for details). **(B)** Left panel: Standard deviation of firing rates as a function of ${\tau}_{d}^{I}$. Middle and right panels: the distributions of firing rates when ${\tau}_{d}^{I}$ takes the indicated four values. **(C)** Left panel: σ_{SpikeNum} as a function of μ_{SpikeNum} when ${\tau}_{d}^{I}=7$ ms (blue), 10 ms (red) and 14 ms (black), with μ_{SpikeNum} and σ_{SpikeNum} being the mean and standard deviation of the spike number distribution for a neuron in different synchronous events. Middel panel: *CV*_{WithinEvent} (see Section 2.5) as a function of ${\tau}_{d}^{I}$, which quantifies the burstiness/regularity of pieces of spike trains within synchronous events. Right panel: *CV*_{events} (see Section 2.5) as a function of ${\tau}_{d}^{I}$, which quantifies the burstiness/regularity of the occurrence of synchronous events. **(D)** Index of HCC (see Section 2.5) as a function of ${\tau}_{d}^{I}$. The left panel of **(C)** is plotted from one simulation trial, while the other panels are plotted from averaging over 32 simulation trials. Error bars represent s.e.m.

### 3.2. Overview on the Strategy of Implementing Spike Shuffling Methods onto the Spike Patterns of the LIF Network

In the following two subsections, we will explain in details how we investigate the influence of spike pattern structure onto efficacy variability by implementing a variety of spike shuffling methods onto the spike patterns generated by the LIF network. But first of all, we give an overview on our strategy.

We use different spike shuffling methods for the spike patterns in asynchronous and synchronous states, because of the sharp difference of their pattern structure. These methods are explained in details in Section 2.6, Figure 5, and Supplementary Movie. For all the recorded spike patterns, we first treat them using a shuffling method called Whole-train Swap (WTS), before implementing any other method for further studies. The reason is that our spike patterns are generated by an underlying network with a specific structure, so the statistics of the spike patterns may be correlated with the structure of the underlying network. For example, if neuron *a* connects to neuron *b* in a network ${N}$, then the cross-correlation between these two neurons in the spike pattern ${P}$ generated by ${N}$ is more likely to be strong; so if the network ${N}$ undergoes STDP according to ${P}$, then the synapse *a*→*b* is more likely to be potentiated strongly. Now suppose that there is another network ${N}$′ with the exactly the same structure with ${N}$, except that the neuronal labels are different. In this case, ${N}$ and ${N}$′ will generate spike patterns with exactly the same statistics. However, neuron *a* may not connect to neuron *b* in ${N}$′, so if ${N}$′ also undergoes STDP according to ${P}$, then the strong cross-correlation between neuron *a* and neuron *b* in pattern ${P}$ will not contribute to the efficacy variability of ${N}$′. Therefore, the efficacy variability of a network not only depends on the spike pattern statistics, but also depends on the possible pattern-network coupling. By comparing the efficacy variability under the original spike pattern and that under the spike pattern treated by WTS, we can understand how strongly this pattern-network coupling may contribute to the efficacy variability. We found that in both asynchronous and synchronous spike patterns of the LIF network, this pattern-network coupling can only slightly increase the efficacy variability (see Figures 10A,B). This suggests that the efficacy variability under both asynchronous and synchronous spike patterns can be largely understood by the spike pattern statistics alone. For all the spike patterns used in the following subsections, we first implement WTS to remove this pattern-network coupling.

From Figure 5, we see that a shuffling method may simultaneously destroy more than one aspects of pattern structure. Therefore, we must carefully design the order to implement them when trying to understand how each aspect of pattern structure influences efficacy variability. For example, to understand the effect of population rate fluctuation with time in asynchronous patterns, we would like to treat the patterns using TM. However, from Figure 5, TM destroys not only synchronous firing but also heterogeneity of cross-correlations, so to investigate the effect of synchronous firing, we must implement TM onto the spike patterns whose heterogeneity of cross-correlations are already destroyed. Therefore, we first treat the spike patterns using SS, thereby destroying heterogeneity of cross-correlations, and then investigate the change of the efficacy variability after further implementing TM. This is the basic idea we used to design our research.

We outline all the spike patterns we compare and the main results we obtained from the following sections in Table 1. For the convenience of the following discussions, we denote *P*_{S1+S2+⋯} to be the spike patterns sequentially shuffled by methods *S*_{1}, *S*_{2}, … from the spike patterns firstly shuffled by WTS. We denote *P*_{0} to be the original spike patterns, and *P*_{WTS} to be the patterns shuffled by WTS.

**Table 1. The spike patterns that we compare to understand the influence onto efficacy variability by different spike pattern structures, and the outline of main results**.

### 3.3. The Influence of Spike Pattern Structure onto Efficacy Variability in Asynchronous Patterns

Here we explain in details how we investigate the influence of spike pattern structure onto efficacy variability by implementing spike shuffling methods onto the spike patterns of the LIF network working in asynchronous states (when ${\tau}_{d}^{I}\le 6\mathrm{\text{ms}}$).

(1) To examine the effect of synchronous firing, we observed the time evolution of the efficacy variability caused by *P*_{SS} and *P*_{SS+TM}. From Figure 5, TM not only destroys synchronous firing, but also influences heterogeneity of cross-correlations, so we first destroy heterogeneity of cross-correlations using SS before implementing TM to study the effect of synchronous firing. As DiffV ∝ *t* and DriftV ∝ *t*^{2}, we fit the time evolution of the efficacy variability with ${c}_{diffv}t+{c}_{driftv}{t}^{2}$ under *P*_{SS} and *P*_{SS+TM}, with *c*_{diffv} and *c*_{driftv} being two to-be-fitted coefficients that respectively manifest the strengths of DiffV and DriftV. See Section 2.4 for fitting details. We found that *c*_{diffv} is larger under *P*_{SS} than under *P*_{SS+TM} (Figure 8A1), which manifests the contribution of synchronous firing onto DiffV. We also found that *c*_{driftv} is approximately zero under *P*_{SS+TM}, but not (although weak) in *P*_{SS} (Figure 8A2, blue lower panel). To understand this, from Supplementary Materials Section S1.3.1, in the absence of heterogeneity of cross-correlations,

with $\Delta {\stackrel{~}{w}}_{ab}$ being the change of the synapse from neuron *b* to neuron *a* only caused by STDP (without considering synaptic homeostasis), *E*_{(ab)} representing averaging over all the synapses, and Var_{a}(*r*_{a}) being the variance of the firing rates. In this equation, $\left|{\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})\right|$ quantifies the imbalance of the strengths of the potentiation and depression process under STDP (*P-D imbalance*), and we see that DriftV can be induced by the interaction of P-D imbalance and heterogeneity of rates. Because of the absence of synchronous firing in *P*_{SS+TM} and *A*_{p} = *A*_{d} in our model (Equation 5), the potentiation and depression of STDP are balanced in *P*_{SS+TM} (Figure 8A2, red upper panel), but not in *P*_{SS} because of the population rate fluctuation with time (Figure 6A, upper panels). This weak P-D imbalance under *P*_{SS} makes DriftV non-zero under heterogeneity of firing rates, which makes *c*_{driftv} non-zero. These findings suggest that in asynchronous spike patterns, the fluctuation of population rate increases DiffV, and influences DriftV under heterogeneity of rates by changing P-D imbalance.

(2) To manifest the effect of heterogeneity of cross-correlations, we compared the efficacy variability under *P*_{STR} and *P*_{STR+TM}. Here we use STR to flatten the population rate fluctuation with time, before further using TM to study the effect of heterogeneity of cross-correlations. We found that *c*_{driftv} is significantly larger than zero under *P*_{STR} but close to zero under *P*_{STR+TM} (Figure 8B2, blue lower panel). From Supplementary Materials Section S1.3.2, DriftV can be largely considered to be contributed by two factors: heterogeneity of cross-correlations and the interaction between P-D imbalance and heterogeneity of rates. To understand the non-zero DriftV in *P*_{STR}, note that *P*_{STR} is free of population rate fluctuation, so because of *A*_{p} = *A*_{d} in our model (Equation 5), we should expect P-D *balance* in *P*_{STR}. Although P-D *imbalance* may slightly exist in *P*_{STR} (Figure 8B1), it is fairly weak. To understand how strong DriftV this weak P-D imbalance contributes, we can compare the P-D imbalance and DriftV under *P*_{STR} with those under *P*_{SS} (Figure 8A2, blue lower panel and Figure 8B2, blue lower panel). We can see that the P-D imbalance under *P*_{SS} is much stronger than that under *P*_{STR}, but the DriftV under *P*_{SS} is much weaker than that under *P*_{STR}. These suggest that the strong DriftV under *P*_{STR} is not caused by the weak P-D imbalance, but instead by the heterogeneity of cross-correlations remnant from the original patterns (Figure 6D). TM destroys this heterogeneity of cross-correlations (Figure 8B), and *c*_{driftv} accordingly becomes almost zero (Figure 8B2, blue lower panel).

(3) To examine the effect of auto-correlation structure, we compared the efficacy variability under *P*_{TM} with that under *P*_{TM+SS}. In these two patterns, synchronous firing and heterogeneity of cross-correlations are destroyed by TM; and SS then is used to reduce *CV* to almost 1 (Figure 8C1). We studied both the case *A*_{p} = *A*_{d} and the case *A*_{p} ≠ *A*_{d} here to investigate whether the influence of auto-correlation structure onto DriftV depends on P-D imbalance. We found that SS hardly influences *c*_{driftv} but significantly reduces *c*_{diffv} (Figures 8C2,C3). This suggests that the burstiness of stationary spike trains does not change DriftV, but increases DiffV.

(4) To study the effect of heterogeneity of rates, we compared the efficacy variability under *P*_{TM+SS} and *P*_{TM+NRC}. Both SS and NRC destroy heterogeneity of cross-correlations, and result in spike trains with *CV* ≈ 1 (Figure 8D1, blue lower panel), except that SS keeps the heterogeneity of rates, while NRC homogenizes the firing rates (Figure 8D1, red upper panel). We found that *c*_{diffv} is almost the same under *P*_{TM+SS} and *P*_{TM+NRC}, regardless whether *A*_{p} = *A*_{d} or *A*_{p} ≠ *A*_{d} (Figure 8D2), which suggests that heterogeneity of rates does not significantly change DiffV if the *CV* of the spike trains is kept around 1. To understand the influence of heterogeneity of rates onto DriftV, we kept *A*_{p} = 1, and observed how *c*_{driftv} and P-D imbalance change with *A*_{d} under these two patterns. We found that *c*_{driftv} is always near to zero under *P*_{TM+NRC}, but is positively correlated with the strength of P-D imbalance under *P*_{TM+SS} (Figure 8D3). This shows that P-D imbalance can induce DriftV under heterogeneity of rates.

**Figure 8. How spike shuffling methods influence the efficacy variability and spike pattern statistics under asynchronous states**. **(A1,A2)** Comparing *P*_{SS} (solid line) with *P*_{SS+TM} (dashed line). **A1**: *c*_{diffv} is larger under *P*_{SS} than under *P*_{SS+TM}. **A2:** *c*_{driftv} (blue) and ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ (red) are both nonzero under *P*_{SS}, and are both almost zero under *P*_{SS+TM}. Note that there are two panels marked in different colors in this subplot and some subplots below. Here, ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ represents the mean synaptic changes only caused by STDP (i.e., if the synaptic homeostasis is absent) (see Equation 11), and we use $\left|{\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})\right|$ to quantify P-D imbalance. The method for calculating *c*_{diffv} and *c*_{driftv} is explained in Section 2.4. **(B1,B2)** Comparing *P*_{STR} (solid line) with *P*_{STR+TM} (dashed line). **B1**: ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ is weak under *P*_{STR} and almost zero under *P*_{STR+TM}. **B2**: *c*_{driftv} (blue) and the index of HCC (red) are both larger under *P*_{STR} than under *P*_{STR+TM}. Here, the index of HCC quantifies the heterogeneity of cross-correlations (see Section 2.5), and the non-zero index of HCC under *P*_{STR+TM} manifests its background value caused by chance in the spike patterns. **(C1–C3)** Comparing *P*_{TM} (solid line) with *P*_{TM+SS} (dashed line). **C1**: The mean *CV* of spike trains under these two patterns. **C2**: *c*_{diffv} (blue) is larger under *P*_{TM} than under *P*_{TM+SS}, but *c*_{driftv} (red) is almost the same. *A*_{p} = *A*_{d} = 1. **C3**: The same as **C2**, except that *A*_{p} = 4∕3 and *A*_{d} = 1. **(D1–D3)** Comparing *P*_{TM+SS} (solid line) with *P*_{TM+NRC} (dashed line). **D1**: Mean *CV* of spike trains (blue) and the standard deviation of firing rates σ_{rate} (red) under these two patterns. **D2**: *c*_{diffv} is almost the same under *P*_{TM+SS} with under *P*_{TM+NRC}, both when *A*_{p} = *A*_{d} = 1 (blue), and when *A*_{p} = 4∕3 and *A*_{d} = 1 (red). **D3**: *c*_{driftv} (blue) and ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ (red) as functions of *A*_{d}, keeping *A*_{p} = 1, ${\tau}_{d}^{I}=6\mathrm{\text{ms}}$. Note that under *P*_{TM+NRC}, *c*_{driftv} ≈ 0; while under *P*_{TM+SS}, *c*_{driftv} increases with $\left|{\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})\right|$, and *c*_{driftv} ≈ 0 only when $\left|{\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})\right|\approx 0$ (indicated by the vertical dashed line). In **(A1–D3)**, *A*_{p} = *A*_{d} = 1 by default. Simulations lasted for 41 s of biological time, and STDP and synaptic homeostasis started after 2 s of transient period. Error bars represent s.e.m. over 32 trials.

### 3.4. The Influence of Spike Pattern Structure onto Efficacy Variability in Synchronous Patterns

Here we explain in details how we investigate the influence of spike pattern structure onto efficacy variability by implementing spike shuffling methods onto the spike patterns of the LIF network working in synchronous states (when ${\tau}_{d}^{I}\ge 7\mathrm{\text{ms}}$).

(1) To examine the effects of synchronous firing, we fitted the time evolution of efficacy variability with ${c}_{diffv}t+{c}_{driftv}{t}^{2}$ under *P*_{SSiE} and *P*_{SSiE+STR}. Here we would like to use STR to flatten the population firing rate. But as STR not only destroys synchronous firing, but also influences auto-correlation structure and heterogeneity of cross-correlations (Figure 5), we first implement SSiE to result in inhomogeneous Poisson processes before investigating the effect of synchronous firing using STR. We found that *c*_{diffv} under *P*_{SSiE} is larger than that under *P*_{SSiE+STR} (Figure 9A1), which manifests the contribution of synchronous firing onto DiffV. We also found that *c*_{driftv} is approximately zero under *P*_{SSiE+STR}, but not in *P*_{SSiE} (Figure 9A2, blue lower panel). The reason is that the potentiation and depression of STDP are balanced in *P*_{SSiE+STR} , because of the absence of synchronous firing in *P*_{SSiE+STR} and *A*_{p} = *A*_{d} in our model (Equation 8); but not in *P*_{SSiE}, because of the population rate fluctuation with time (Figure 9A2, red upper panel). And this P-D imbalance results in DriftV under heterogeneity of rates. These findings suggest that synchronous firing generally increases DiffV, and influences DriftV under heterogeneity of rates by changing P-D imbalance, which are consistent with our results for asynchronous patterns (see Term 1 in the previous section).

(2) To examine the effects of heterogeneity of rates, we compared the efficacy variability under *P*_{SSiE} with that under *P*_{SSiE+TSiE}. SSiE destroys heterogeneity of cross-correlations, so that the interaction of heterogeneity of rates and P-D imbalance becomes the only possible source of DriftV. TSiE further destroys heterogeneity of rates, which reduces the DriftV (i.e., *c*_{driftv}) under *P*_{SSiE+TSiE} almost to zero (Figure 9B, blue lower panel). When we keep *A*_{p} = 1 but assign *A*_{d} different values, the P-D imbalance is accordingly changed. We found that *c*_{driftv} under *P*_{SSiE} is positively correlated with the strength of P-D imbalance (Figure 9B). These observations again support the idea that P-D imbalance can induce DriftV under heterogeneity of rates (see Term 4 in the previous section).

(3) To examine the effects of AT_{WithinEvent}, we compared the efficacy variability under *P*_{TSiE} with that under *P*_{TSiE+SSiE}. By definition (Figure 5), SSiE not only randomizes the spike times, but also destroys the heterogeneity of cross-correlations. So before using SSiE to investigate the effect of AT_{WithinEvent}, we first implemented TSiE to destroy heterogeneity of cross-correlations in these two patterns. We found that after SSiE, the pieces of spike trains within synchronous events becomes burstier (Figure 9C1, red upper panel), and *c*_{diffv} also becomes larger (Figure 9C1, blue lower panel). This suggests that burstier spike trains within synchronous events increases DiffV.

In these two patterns, TSiE destroys both heterogeneity of cross-correlations and heterogeneity of rates (which are the two sources of DriftV), so *c*_{driftv} ≈ 0. To understand the possible influence of AT_{WithinEvent} onto DriftV under heterogeneity of rates, we compared the P-D imbalance under *P*_{TSiE} with that under *P*_{TSiE+SSiE}. We found that the P-D imbalance under these two patterns are almost the same (Figure 9C2), which suggests that AT_{WithinEvent} has no influence onto DriftV under heterogeneity of rates.

(4) To examine the effects of AT_{SpikeNum}, we compared the efficacy variability under *P*_{TSiE+SSiE} with that under *P*_{TSiE+NRC}. In these two patterns, TSiE homogenizes firing rates, so that different neurons have almost the same spike number distribution across synchronous events. Both SSiE and NRC randomize spike times (so that *P*_{TSiE+SSiE} and *P*_{TSiE+NRC} have similar AT_{WithinEvent}, see Supplementary Materials Section S4 for more discussions); but SSiE keeps this spike number distribution, while NRC makes this distribution broader (Figure 9D1, red upper panel). We found that *c*_{diffv} is larger under *P*_{TSiE+NRC} than under *P*_{TSiE+SSiE} (Figure 9D1, blue lower panel). This suggests that broader distribution of the spike numbers a neuron fires in different synchronous events increases DiffV.

We also found that P-D imbalance under these two patterns are almost the same (Figure 9D2), which suggests that AT_{SpikeNum} has no influence onto DriftV under heterogeneity of rates.

(5) To examine the effects of AT_{events}, we compared the efficacy variability under *P*_{SSiE} with that under *P*_{SSiE+ETS}. As ETS changes not only AT_{events} but also heterogeneity of cross-correlations (Figure 5), we first implemented SSiE to destroy heterogeneity of cross-correlations before further using ETS to study the effects of AT_{events}. ETS increases the burstiness of the occurrence of synchronous events (Figure 9E1), but we found that *c*_{diffv} is almost the same under these two patterns (Figure 9E2). However, this is possibly because that the irregularity of spike trains within synchronous events dominates DiffV, while the contribution 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. In Supplementary Material Section S5.1.4, we study the influence of *CV*_{events} onto DiffV in more details by comparing *P*_{TSiE+SSiE} with *P*_{TSiE+SSiE+ETS}, with TSiE being used to homogenize firing rate to make DriftV = 0. This makes the efficacy variability totally caused by DiffV, so that the estimated value of *c*_{diffv} becomes more precise. To increase the STDP interaction between spikes in different synchronous events, we also increase the STDP time scale τ_{STDP} (see Equation 5). We find that *c*_{diffv} under *P*_{TSiE+SSiE+ETS} is larger than that under *P*_{TSiE+SSiE} (Supplementary Figure 6), which suggests that burstier occurrence of synchronous events tends to increase DiffV.

How AT_{events} influences DriftV may be intriguing. Our previous paper (Bi and Zhou, 2016) suggests that if synchronous-event overlapping is absent (i.e., if a synchronous event ${S}$ happens during the interval [*t*_{1}, *t*_{2}], then no other synchronous event should happen within [*t*_{1} − τ_{delay}, *t*_{2}], with τ_{delay} being the axonal delay), then the burstiness of the occurrence of synchronous events tends to potentiate synapses when *A*_{p} > *A*_{d} and depress synapses when *A*_{p} < *A*_{d}. We carefully avoid synchronous-event overlapping during ETS (see Section 2.6.3.3). Consistently, we found that ETS tends to potentiate synapses when *A*_{p} > *A*_{d} and depress synapses when *A*_{p} < *A*_{d} (Figures 9E3,E4, red upper panels), which is consistent with our previous results. In our model, this weakens (or strengthens) P-D imbalance (i.e., the absolute value of mean synaptic changes under STDP) when *A*_{p} > *A*_{d} (or *A*_{p} < *A*_{d}), which decreases (or increases) DriftV under heterogeneity of rates (Figures 9E3,E4, blue lower panels).

**Figure 9. How spike shuffling methods influence the efficacy variability and spike pattern statistics under synchronous states**. **(A1,A2)** Comparing *P*_{SSiE} (solid line) with *P*_{SSiE+STR} (dashed line). **A1**: *c*_{diffv} is larger under *P*_{SSiE} than under *P*_{SSiE+STR}. **A2:** *c*_{driftv} (blue) and ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ (red) are both nonzero under *P*_{SS}, and are both almost zero under *P*_{SS+STR}. Note that there are two panels marked in different colors in this subplot and some subplots below. See the caption of Figure 8A2 for the meaning of ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$. **(B)** Comparing *P*_{SSiE} (solid line) with *P*_{SSiE+TSiE} (dashed line). *c*_{driftv} and ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ as functions of *A*_{d}, keeping *A*_{p} = 1 and ${\tau}_{d}^{I}=7\mathrm{\text{ms}}$. Note that under *P*_{SSiE+TSiE}, *c*_{driftv} ≈ 0; while under *P*_{SSiE}, *c*_{driftv} increases with $\left|{\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})\right|$, and *c*_{driftv} ≈ 0 only when $\left|{\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})\right|\approx 0$ (indicated by the vertical dashed line). **(C1,C2)** Comparing *P*_{TSiE+SSiE} (solid line) with *P*_{TSiE+NRC} (dashed line). **C1**: Both *c*_{diffv} and σ_{SpikeNum} are larger under *P*_{TSiE+NRC} than under *P*_{TSiE+SSiE}. Here, σ_{SpikeNum} represents the standard deviation of the distribution of spike number per neuron per synchronous event. **C2**: ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ under these two patterns strongly overlap. **(D1,D2)** Comparing *P*_{TSiE} (solid line) with *P*_{TSiE+SSiE} (dashed line). **D1**: Both *c*_{diffv} and *CV*_{WithinEvent} are larger under *P*_{TSiE} than under *P*_{TSiE+SSiE}. Here, *CV*_{WithinEvents} quantifies burstiness/regularity of the pieces of spike trains within synchronous events (See Section 2.5 for details). **D2**: ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ under these two patterns strongly overlap. **(E1–E4)** Comparing *P*_{TSiE} (solid line) and *P*_{TSiE+ETS} (dashed line). **E1**: The burstiness of the occurrence of synchronous events (quantified by *CV*_{events}, see Section 2.5) under these two patterns. **E2**: *c*_{diffv} under these two patterns both when *A*_{p} = 4.∕3 and *A*_{d} = 1 (blue), and when *A*_{p} = 1 and *A*_{d} = 4.∕3 (red). **E3**: ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ (red) and *c*_{driftv} (blue) under these two patterns when *A*_{p} = 4.∕3 and *A*_{d} = 1. **E4**: The same as **E3**, except that *A*_{p} = 1 and *A*_{d} = 4.∕3. In **(A1–E4)**, *A*_{p} = *A*_{d} = 1 by default. Simulations lasted for 41 s of biological time, and STDP and synaptic homeostasis started after 2 s of transient period. Error bars represent s.e.m. over 32 trials.

### 3.5. Understanding the Efficacy Variability under the Spike Patterns Generated by the LIF Network

The influence of synaptic time scales onto network dynamics has been discussed in theoretical studies (Brunel and Wang, 2003; Geisler et al., 2005). When the decay time scale of inhibitory conductance ${\tau}_{d}^{I}\le 6\mathrm{\text{ms}}$, our LIF network works in asynchronous state with the firing rate of each neuron fluctuating chaotically (Figure 6A). It is believed that the high-dimensional nature of the state trajectory of such networks provides substrate for complex computations such as classifying temporal signals (Ostojic, 2014) and learning complex spatio-temporal patterns (Sussillo and Abbott, 2009), also see the seminal papers by Maass et al. (2002) and Jaeger (2003) on reservoir computing. When ${\tau}_{d}^{I}=7$ or 8 ms, the network works in weak synchronous state (Figure 3A, middle panel). In this case, the oscillation frequency is larger than the firing frequency of most neurons, and the spike trains of most neurons look irregular despite the population coherent rhythm, which is consistent with the observations in hippocampus and cortex (Csicsvari et al., 1998; Fries et al., 2001). When ${\tau}_{d}^{I}$ continues to grow (especially when ${\tau}_{d}^{I}=13\mathrm{\text{ms}}$ or 14 ms), nearly every neuron fires burstly in a synchronous event (Figure 3A, right panel), just like in the spike patterns observed in epilepsy (Gulyás and Freund, 2015). In the following discussions, we will use spike shuffling methods to understand the contributions of different statistical features to the efficacy variability under these spike patterns generated by our LIF network, thereby gaining understanding on the efficacy variability under the spike patterns observed in these theoretical and experimental studies. For example, we would like to know which spike pattern statistics influences to the efficacy variability significantly, and which spike pattern statistics is the main reason for the change of efficacy variability with ${\tau}_{d}^{I}$. In our simulations, we kept the firing rate of the excitatory population at 20 Hz for different ${\tau}_{d}^{I}$s (see Section 2.3), so that the change of the efficacy variability with ${\tau}_{d}^{I}$ is totally caused by higher order pattern statistics.

#### 3.5.1. Asynchronous States

For asynchronous states, we sequentially shuffled the original spike pattern by WTS, STR, TM, SS, and NRC, and then observed the efficacy variability in the resulting spike patterns (Figure 10A). The change of the efficacy variability after each of these shuffling methods respectively manifests the contribution of pattern-network coupling (for WTS), population rate fluctuation with time (for STR), heterogeneity of cross-correlations (for TM), auto-correlation structure (for SS), and heterogeneity of rates (for NRC). Here, after WTS, we first use STR to assess the influence of population rate fluctuation onto efficacy variability. From Figure 5, STR not only flattens population rate, but also influence auto-correlation structure and heterogeneity of cross-correlations. However, as the fluctuation of population rate in asynchronous states is not so strong, the auto-correlation structure and heterogeneity of cross-correlation in the original patterns should largely preserve after STR. From Figure 10A, the efficacy variability is only slightly reduced by STR. After STR, both *CV* and index of heterogeneity of cross-correlation are slightly reduced (Supplementary Figure 8), both of which reduce the efficacy variability. This suggests that the reduction of the efficacy variability caused by population rate fluctuation alone is even smaller than the reduction by STR shown in Figure 10A. Afterwards, we implement TM onto patterns treated by STR for understanding the contribution of heterogeneity of cross-correlations, implement SS onto patterns treated by TM for auto-correlation structure, and compare the effects of SS and NRC on patterns treated by TM for heterogeneity of rates. All of these are consistent with the research design in the previous subsections (see the items about asynchronous states in Table 1).

**Figure 10. Understanding the contributions to the efficacy variability by different spike pattern structures under the spike patterns of the LIF network**. **(A)** The efficacy variance under the original asynchronous patterns and the patterns sequentially shuffled by WTS, STR, TM, SS, and NRC. Note that the variance under *P*_{STR+TM+SS} and that under *P*_{STR+TM+NRC} are strongly overlapped, because of the P-D *balance* caused by the absence of population rate fluctuation after STR and *A*_{p} = *A*_{d} in our model. **(B)** The efficacy variance under the original synchronous spike patterns and under *P*_{WTS}, *P*_{TSiE} and *P*_{SSiE}. Note that the results for the original pattern and *P*_{WTS} are strongly overlapped. **(C)** ${\mathrm{\text{E}}}_{(ab)}(\Delta {\stackrel{~}{w}}_{ab})$ (blue) and the standard deviation of firing rate σ_{rate} (red) as functions of ${\tau}_{d}^{I}$ under *P*_{SSiE}. **(D)** Index of HCC in *P*_{WTS} and *P*_{SSiE}. **(E)** Firing profiles within synchronous events under *P*_{WTS} for neurons with high (blue), middle (red) and low (black) firing rates as well as the whole population (dashed), when ${\tau}_{d}^{I}=7$ ms (left panel) and 14 ms (right panel). The zero point indicates the middle time of a synchronous event, and the function value at time *t* means the firing rate at time *t* relative to the middle time of the synchronous event, averaging over all synchronous events in the spike patterns. See Section 2.5 for details. **(F)** The unit cross-correlation between the neurons with low rate and middle rate (blue), and that between the neurons with high rate and middle rate (red), when ${\tau}_{d}^{I}=7$ ms (left panel) and 14 ms (right panel). The definitions of “high rate,” “middle rate,” and “low rate” neurons are the same as those in **(E)** (see Section 2.5 for details). In **(A–F)**, *A*_{p} = *A*_{d} = 1. Simulations lasted for 41 s of biological time, and STDP and synaptic homeostasis started after 2 s of transient period. Error bars represent s.e.m. over 32 trials.

We found that under our parameter values, heterogeneity of cross-correlations and auto-correlation structure are the two major sources of the efficacy variability; and auto-correlation structure contributes most to the increase of the efficacy variability with ${\tau}_{d}^{I}$ (Figure 10A). As heterogeneity of cross-correlations mainly contributes to DriftV (which is of ${O}$(*t*^{2}) order), and auto-correlation structure mainly contributes to DiffV (which is of only ${O}$(*t*) order), the significant contribution of auto-correlation structure here highlights its importance to understand the efficacy variability in asynchronous states.

#### 3.5.2. Synchronous States

For synchronous states, we first used WTS to destroy the pattern-network coupling, and found that this coupling hardly influences efficacy variability (Figure 10B). Then, we used TSiE to destroy heterogeneity of rates and heterogeneity of cross-correlations, which are the two sources of DriftV, and found that the efficacy variability is reduced by more than 10 times after TSiE (Figure 10B). This suggests that the efficacy variability under the synchronous states in our model is dominated by DriftV, and we only need to focus on DriftV to understand the change of the efficacy variability with ${\tau}_{d}^{I}$.

To compare the contribution of heterogeneity of rates and heterogeneity of cross-correlation to DriftV, we compared the efficacy variability under *P*_{WTS} with that under *P*_{SSiE} (Figure 10B). After SSiE, the heterogeneity of cross-correlations is destroyed. Therefore, the DriftV under *P*_{SSiE} is contributed by the interaction of heterogeneity of rates and P-D imbalance. We found that the efficacy variability under *P*_{SSiE} gets its maximum around ${\tau}_{d}^{I}=9\mathrm{\text{ms}}$, and decreases when ${\tau}_{d}^{I}$ continues to grow (Figure 10B). This can be understood as a result caused by the interaction of the increase of P-D imbalance and the decrease of rate heterogeneity with ${\tau}_{d}^{I}$ (Figure 10C). However, the efficacy variability under *P*_{WTS} is smaller than that under *P*_{SSiE} when ${\tau}_{d}^{I}$ is small, but tends to monotonically increase with ${\tau}_{d}^{I}$ (Figure 10B), which manifests the contribution of heterogeneity of cross-correlations. Indeed, we found that the index of HCC (which quantifies ${\mathrm{\text{Var}}}_{a}({\int}_{-\infty}^{\infty}\text{d}\tau H(\tau ){C}_{a}(-{\tau}_{delay}-\tau ))$, see Section 2.5) tends to increase with ${\tau}_{d}^{I}$, especially in the strong-bursting regime (i.e., ${\tau}_{d}^{I}=13\mathrm{\text{ms}},14\mathrm{\text{ms}}$, see Figure 10D). Together, these observations make us come to the understanding that: (1) in weak synchronous states (when ${\tau}_{d}^{I}$ is small), heterogeneity of rates contributes most to the DriftV, and heterogeneity of cross-correlations reduces the DriftV caused by heterogeneity of rates; (2) in synchronous bursting states (when ${\tau}_{d}^{I}$ is large), heterogeneity of cross-correlation overwhelms heterogeneity of rates to be the dominating factor to DriftV, which pushes DriftV to continuously increase with τ_{d, I}.

To understand the origin of the heterogeneity of cross-correlations in synchronous states, we plot the firing profiles of the neurons during synchronous events (Figure 10E). We found that under *P*_{WTS}, in weak synchronous states (when ${\tau}_{d}^{I}$ is small), the neurons with high firing rates (*high-rate neurons*) tends to start to fire earlier and stop to fire later than the neurons with low firing rates (*low-rate neurons*) in a synchronous event (Figure 10E, left panels); in synchronous bursting states (when ${\tau}_{d}^{I}$ is large), however, high-rate neurons still stop to fire later than low-rate neurons at the end of a synchronous event, but high-rate and low-rate neurons tend to start to fire at the same time with the same rate at the beginning of a synchronous event (Figure 10E, right panels). To understand this, note that the high-rate (low-rate) neurons tend to receive more (less) excitatory inputs and less (more) inhibitory inputs in the network, which makes the inputs of high-rate neurons increase above (decrease below) the firing threshold earlier (later) than those of the low-rate neurons at the beginning (end) of a synchronous event. This is the reason for the firing profiles when ${\tau}_{d}^{I}$ is smaller. When ${\tau}_{d}^{I}$ is large, however, the strengths of the inhibitory currents in our model (Equation 8) are accordingly scaled down. This enlarges the time window for the supra-threshold excitatory currents to quickly push the firing rates of all the neurons to near saturation at the beginning of a synchronous event, before the time-integrations of the inhibitory currents gradually turn off the neuronal activities, from low-rate to high-rate neurons. Now suppose that the 0th neuron in the network receives both from the *h*th neuron with higher firing rate and from the *l*th neuron with lower firing rate. When ${\tau}_{d}^{I}$ is small, both the unit cross-correlation between the *h*th and 0th neurons *C*_{h0}(τ) and that between the *l*th and 0th neurons *C*_{l0}(τ) are almost symmetric around τ = 0 (Figure 10F, left panel). But when ${\tau}_{d}^{I}$ is large, there is an apparent left-shift of *C*_{h0}(τ) comparing with *C*_{l0}(τ) (Figure 10F, right panel): this causes ${\int}_{-\infty}^{\infty}\text{d}\tau H(\tau ){C}_{h0}(-{\tau}_{delay}-\tau )$ very different from ${\int}_{-\infty}^{\infty}\text{d}\tau H(\tau ){C}_{l0}(-{\tau}_{delay}-\tau )$ when ${\tau}_{d}^{I}$ is large. This is the reason for the increase of heterogeneity of cross-correlations with ${\tau}_{d}^{I}$.

Another issue that has not been considered till now is the possible dependence (on, say, spike times or spike numbers) between the pieces of spike trains in adjacent synchronous events. To check the influence of this inter-event dependence onto the efficacy variability, we compared the efficacy variability under *P*_{WTS} with that under *P*_{EOS} (see Section 2.6.3.4 for the method of EOS). We found that EOS hardly influences the efficacy variability (Supplementary Figure 9), which suggests that this inter-event dependence hardly influences the efficacy variability in our model.

## 4. Discussion

In this paper, we developed a systematic spike shuffling approach to alter four types of statistics of spike patterns (i.e., synchronous firing, auto-correlation structure, heterogeneity of rates and heterogeneity of cross-correlations) under both asynchronous and synchronous states, and then applied this approach to systematically study the influence of the four aspects of pattern structures onto the efficacy variability under STDP and synaptic homeostasis in the spike patterns self-organized by a biologically plausible LIF neuronal network. The main results are shown in Table 1, which can be summarized as (1) synchronous firing and burstiness tend to increase DiffV, (2) synchronous firing influences P-D imbalance, which can induce DriftV together with heterogeneity of rates, and (3) heterogeneity of cross-correlations induces DriftV together with heterogeneity of rates. We compared our results with those of our previous paper (Bi and Zhou, 2016) in Supplementary Material Section S5, and found their consistency. We also examined the contributions of different pattern statistics to the efficacy variability under the spike patterns of the LIF network, and found that auto-correlation structure 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. We believe our method can contribute to the library of spike shuffling methods, and provide new angles for experimentalists to analyze their data; and our results can help to understand the efficacy variability under the spike patterns observed in theoretical and experimental studies.

When the time scales of inhibition and excitation are comparable, our LIF network works in asynchronous states; when the inhibition time scale starts to grow, the LIF network transites to weak synchronous state at some point, performing oscillation in the gamma frequency regime (Figure 3A). In our simulations, we fixed the excitatory decay time constant at 4 ms while changed the inhibitory decay time constant within 3–14 ms, both are around typical values of AMPA and GABA_{A} currents [2 ~ 5 ms for AMPA (Zhou and Hablitz, 1998; Angulo et al., 1999), 5 ~ 15 ms for GABA_{A} (Xiang et al., 1998; Gupta et al., 2000)]. In physiological situations, neuronal network dynamics is not only determined by synaptic time scales, but also by other factors such as the neuronal response properties and connections; so it is difficult to directly compare the spike patterns in our simulations with those observed in experiments. However, we can still gain insight onto the functional roles of the efficacy variability caused by the dynamics through the computational tasks that the network takes.

In the asynchronous states, the firing rate of each neuron in our LIF network fluctuates chaotically (Figure 6A, lower panels), and our study suggests that the irregular auto-correlation structure of spike trains in this chaotic asynchronous state may contribute significantly to the efficacy variability of the recurrent connections. Theoretically, it has been proposed that larger variance of recurrent synaptic weights facilitates chaoticity of a network (Sompolinsky et al., 1988; Toyoizumi and Abbott, 2011), which again promotes irregularity of spike trains. This mutual facilitation of chaos and synaptic variance may have important implications in the development of brain areas such as the primary olfactory system and the cerebellum, which may use the high-dimensional nature of the state trajectory in chaos to do complex computations such as odor discrimination (Mazor and Laurent, 2005) and motion control (Buonomano and Mauk, 1994; Yamazaki and Tanaka, 2007). One problem of such chaos-based computations is that the dynamics may be sensitive to the noises and the initial conditions of the network. This sensitivity may be suppressed using a feedback loop from the readout unit (Jaeger and Haas, 2004; Sussillo and Abbott, 2009) or perform suitable plastic changes on the recurrent weights (Laje and Buonomano, 2013), so that the innate trajectory becomes stabilized.

Gamma oscillations is believed important for memory formation under normal physiological conditions (Sederberg et al., 2007; Jutras et al., 2009; Yamamoto et al., 2014). Our intuition shown in Figure 2A suggests that successful memory embedding requires small efficacy variability, which may be a physiological function of gamma oscillation. However, in our simulations, we did not observe small efficacy variability when our LIF network works in weak synchronous state (Supplementary Figure 7). This may be because that we did not consider the shunting effect of perisomatic fast-spiking interneurons in our model. On the one hand, these interneurons may help to homogenize neuronal firing rates (Vida et al., 2006), which reduces DriftV. On the other hand, they are also able to entrain high temporal precision of spike time (Cobb et al., 1995; Salkoff et al., 2015), resulting in spike-to-spike synchrony of pyramidal neurons, which helps to reduce DiffV (see Section 3.3 of our previous paper Bi and Zhou, 2016).

When the time scale of feedback inhibition is much larger than that of excitation, our LIF network exhibits low-frequency oscillation burstiness (Figure 3A, left panel). This is similar to the low-frequency (about 3Hz) “spike-and-wave” EEG pattern commonly observed in absence seizure (Hughes, 2009), in which the “spike” component is associated with neuronal firing, while the “wave” is associated with hyperpolarization of neurons. In normal case, the thalamocortical (TC) cells are mainly inhibited by fast GABA_{A} currents from the thalamic reticular (RE) cells recurrently connected with them in the thalamus, and oscillates with frequency around 10 Hz. In pathological case, however, the cortical pyramidal neurons (PN) becomes hyperexcited by, say, lack of synaptic inhibition (Maheshwari and Noebels, 2014) or impaired hyperpolarization-activated current (*I*_{h}) (Strauss et al., 2004) in PN. In this case, the strong excitatory corticothalamic feedback will evoke the IPSPs in TC cells from RE cells dominated by the slow GABA_{B} component, which lowers down the oscillation frequency to about 3 Hz (Destexhe, 2007, 2008). Therefore, absence seizure is closely related to the prolonged inhibitory time scale in the thalamus, which shares a similar mechanism as the synchronous burstiness observed in our LIF network. Our study shows that heterogeneity of cross-correlation is the main reason for the large efficacy variability in this state, which provides possible understanding to the memory deficit in children with absence seizure (Nolan et al., 2004; Henkin et al., 2005).

However, although spike shuffling methods provide important insights into the influence of spike pattern structure to synaptic plasticity, they have their own limitations. Firstly, the spike shuffling methods we use always randomize spike patterns, so that the spike patterns look more like homogeneous Poisson processes with homogeneous firing rate after being treated by each method. For example, Spike Swap (SS) is able to change the CV of stationary spike trains. But this change is always toward the direction of “randomization”: it increases (decreases) CV if that of the original spike train is smaller (larger) than 1. This may be solved using parametric spike shuffling methods. For example, spike trains can be regularized if we constrain the shuffling process using spike history effect, such as refractory period (Berry and Meister, 1998), or correlations among adjacent inter-spike intervals (Brandman and Nelson, 2002). Secondly, from Figure 5, some spike shuffling methods simultaneously strongly alter more than one pattern statistics. As a result, before we study the influence of a pattern statistics using a shuffling method, we have to first treat the spike pattern using other shuffling methods to randomize other pattern statistics thereby nullifying their influences (see Section 3.2). This limits us from understanding how a pattern statistics may influence the efficacy variability when other pattern statistics remain unchanged. From this aspect, the advantage of the strategy of our previous paper (Bi and Zhou, 2016) becomes manifested, in which the statistical features of the spike patterns can be explicitly controlled by the statistical models, which facilitates a thorough and systematic study in a large parameter range. Therefore, both the strategies in this paper and in our previous paper have their pros and cons. They complement each other, helping us to gain a thorough and convincing understanding on the problem.

During plasticity, synaptic efficacies and network dynamics interact with each other. In both of our papers, we only studied the influence of network dynamics onto the efficacy variability under STDP and synaptic homeostasis, by supposing that the neurons fire spikes according to pre-specified spike patterns, irrelevant with the synaptic weights. Although this approach saved us from the complex co-evolution of the two (Figure 3B), it does not provide a thorough understanding on the dynamics of plastic networks. Studies on the synapse-dynamics coevolution usually use the assumption that the timescale of spiking covariance is far smaller than that of plasticity, thereby focusing on the evolution of the expectations of synaptic weights (Babadi and Abbott, 2013; Ocker et al., 2015), and neglecting the trial-to-trial variabilities; or use heavy simulations (Fiete et al., 2010; Litwin-Kumar and Doiron, 2014) to get phenomenological understandings. The key difficulty here lies in the efficacy variability during the plasticity process: initial variability of synaptic weights generates different network dynamics for different trials, which then further amplifies the trial difference of network structure during plasticity. This may make the synapse-dynamics coevolution sensitive to initial conditions and noises, thereby resulting in sharply different functional performances for different individuals. How do realistic neuronal networks develop robust brain functions through the jungle of noises? We believe our work is able to help to understand this intriguing problem in future research.

## Author Contributions

ZB and CZ conceived the idea and designed the research. ZB performed the research. ZB and CZ wrote the paper.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer CK and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

## Acknowledgments

This work was partially supported by Hong Kong Baptist University (HKBU) Strategic Development Fund, NSFC-RGC Joint Research Scheme (Grant No. HKUST/NSFC/12-13/01) and NSFC (Grant No. 11275027). ZB thanks his PhD supervisor Hai-Jun Zhou for persistent support to his study on computational neuroscience. ZB thanks Dongping Yang for kind helps on coding LIF networks at the beginning of the research. ZB also thanks David Hansel and Carl van Vreeswijk for helpful discussions during his visit in CNRS UMR 8119 in Université Paris Descartes in 2013.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fncom.2016.00083

## References

Allen, C., and Stevens, C. F. (1994). An evaluation of causes for unreliability of synaptic transmission. *Proc. Natl. Acad. Sci. U.S.A.* 91, 10380–10383. doi: 10.1073/pnas.91.22.10380

Amarasingham, A., Harrison, M. T., Hatsopoulos, N. G., and Geman, S. (2012). Conditional modeling and the jitter method of spike resampling. *J. Neurophysiol.* 107, 517–531. doi: 10.1152/jn.00633.2011

Angulo, M. C., Rossier, J., and Audinat, E. (1999). Postsynaptic glutamate receptors and integrative properties of fast-spiking interneurons in the rat neocortex. *J. Neurophysiol.* 82, 1295–1302.

Babadi, B., and Abbott, L. F. (2013). Pairwise analysis can account for network structures arising from spike-timing dependent plasticity. *PLoS Comput. Biol.* 9:e1002906. doi: 10.1371/journal.pcbi.1002906

Bailey, C. H., Giustetto, M., Huang, Y.-Y., Hawkins, R. D., and Kandel, E. R. (2000). Is heterosynaptic modulation essential for stabilizing Hebbian plasiticity and memory? *Nat. Rev. Neurosci.* 1, 11–20. doi: 10.1038/35036191

Bartos, M., Vida, I., and Jonas, P. (2007). Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks. *Nat. Rev. Neurosci.* 8, 45–56. doi: 10.1038/nrn2044

Berke, J. D., and Hyman, S. E. (2000). Addiction, dopamine, and the molecular mechanisms of memory. *Neuron* 25, 515–532. doi: 10.1016/S0896-6273(00)81056-9

Berry, M. J. II, and Meister, M. (1998). Refractoriness and neural precision. *J. Neurosci.* 18, 2200–2211.

Bi, Z., and Zhou, C. (2016). Spike pattern structure influences synaptic efficacy variability under STDP and synaptic homeostasis. I: spike generating models on dendritic motifs. *Front. Comput. Neurosci.* 10:14. doi: 10.3389/fncom.2016.00014

Brandman, R., and Nelson, M. E. (2002). A simple model of long-term spike train regularization. *Neural. Comput.* 14, 1575–1597. doi: 10.1162/08997660260028629

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. *J. Comput. Neurosci.* 8, 183–208. doi: 10.1023/A:1008925309027

Brunel, N., and Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. *Neural. Comput.* 11, 1621–1671. doi: 10.1162/089976699300016179

Brunel, N., and Wang, X.-J. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? I. synaptic dynamics and excitation-inhibition balance. *J. Neurophysiol.* 90, 415–430. doi: 10.1152/jn.01095.2002

Buonomano, D. V., and Maass, W. (2009). State-dependent computations: spatiotemporal processing in cortical networks. *Nat. Rev. Neurosci.* 10, 113–125. doi: 10.1038/nrn2558

Buonomano, D. V., and Mauk, M. D. (1994). Neural network model of the cerebellum: temporal discrimination and the timing of motor responses. *Neural Comput.* 6, 38–55. doi: 10.1162/neco.1994.6.1.38

Butler, C. R., and Zeman, A. Z. (2008). Recent insights into the impairment of memory in epilepsy: transient epileptic amnesia, accelerated long-term forgetting and remote memory impairment. *Brain* 131, 2243–2263. doi: 10.1093/brain/awn127

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. *Science* 304, 1926–1929. doi: 10.1126/science.1099745

Buzsáki, G., and Mizuseki, K. (2014). The log-dynamic brain: how skewed distributions affect network operations. *Nat. Rev. Neurosci.* 15, 264–278. doi: 10.1038/nrn3687

Buzsáki, G., and Wang, X.-J. (2012). Mechanisms of gamma oscillations. *Annu. Rev. Neurosci.* 35, 203–225. doi: 10.1146/annurev-neuro-062111-150444

Cancedda, L., and Poo, M. M. (2009). “Synapse formation and elimination: competition and the role of activity,” in *Encyclopedia of Neuroscience*, ed L. R. Squire (Oxford: Academic Press), 697–703.

Caporale, N., and Dan, Y. (2008). Spike timing-dependent plasticity: a Hebbian learning rule. *Annu. Rev. Neurosci.* 31, 25–46. doi: 10.1146/annurev.neuro.31.060407.125639

Clause, A., Kim, G., Sonntag, M., Weisz, C. J., Vetter, D. E., Rübsamen, R., et al. (2014). The precise temporal pattern of prehearing spontaneous activity is necessary for tonotopic map refinement. *Neuron* 82, 822–835. doi: 10.1016/j.neuron.2014.04.001

Cobb, S. R., Buhl, E. H., Halasy, K., Paulsen, O., and Somogyi, P. (1995). Synchronization of neuronal activity in hippocampus by individual gabaergic interneurons. *Nature* 378, 75–78. doi: 10.1038/378075a0

Cohn, R., Morantte, I., and Ruta, V. (2015). Coordinated and compartmentalized neuromodulation shapes sensory processing in Drosophila. *Cell* 163, 1742–1755. doi: 10.1016/j.cell.2015.11.019

Csicsvari, J., Hirase, H., Czurko, A., and Buzsáki, G. (1998). Reliability and state dependence of pyramidal cell-interneuron synapses in the hippocampus: an ensemble approach in the behaving rat. *Neuron* 21, 179–189. doi: 10.1016/S0896-6273(00)80525-5

Daie, K., Goldman, M. S., and Aksay, E. R. F. (2015). Spatial patterns of persistent neural activity vary with the behavioral context of short-term memory. *Neuron* 85, 847–860. doi: 10.1016/j.neuron.2015.01.006

Dan, Y., and Poo, M. M. (2006). Spike timing-dependent plasticity: from synapse to perception. *Physiol. Rev.* 86, 1033–1048. doi: 10.1152/physrev.00030.2005

Destexhe, A. (2007). Spike-and-wave oscillations. *Scholarpedia* 2:1402. doi: 10.4249/scholarpedia.1402

Destexhe, A. (2008). “Corticothalamic feedback: a key to explain absence seizures,” in *Computational Neuroscience in Epilepsy*, eds I. Soltesz and K. Staley (London: Academic Press), 184–214. doi: 10.1016/b978-012373649-9.50016-8

Dorval, A. D. Jr. and White, J. A. (2005). Channel noise is essential for perithreshold oscillations in entorhinal stellate neurons. *J. Neurosci.* 25, 10025–10028. doi: 10.1523/JNEUROSCI.3557-05.2005

Faisal, A. A., and Laughlin, S. B. (2007). Stochastic simulations on the reliability of action potential propagation in thin axons. *PLoS Comput. Biol.* 3:e79. doi: 10.1371/journal.pcbi.0030079

Faisal, A. A., Selen, L. P. J., and Wolpert, D. M. (2008). Noise in the nervous system. *Nat. Rev. Neurosci.* 9, 292–303. doi: 10.1038/nrn2258

Faisal, A. A., White, J. A., and Laughlin, S. B. (2005). Ion-channel noise places limits on the miniaturization of the brains wiring. *Curr. Biol.* 15, 1143–1149. doi: 10.1016/j.cub.2005.05.056

Fenno, L., Yizhar, O., and Deisseroth, K. (2011). The development and application of optogenetics. *Annu. Rev. Neurosci.* 34, 389–412. doi: 10.1146/annurev-neuro-061010-113817

Fiete, I. R., Senn, W., Wang, C., and Hahnloser, R. H. R. (2010). Spike time-dependent plasticity and heterosynaptic competition organize networks to produce long scale-free sequences of neural activity. *Neuron* 65, 563–576. doi: 10.1016/j.neuron.2010.02.003

Fries, P., Reynolds, J. H., Rorie, A. E., and Desimone, R. (2001). Modulation of oscillatory neuronal synchronization by selective visual attention. *Science* 291, 1560–1563. doi: 10.1126/science.1055465

Funahashi, S., and Inoue, M. (2000). Neuronal interactions related to working memory processes in the primate prefrontal cortex revealed by cross-correlation analysis. *Cereb. Cortex* 10, 535–551. doi: 10.1093/cercor/10.6.535

Ganmora, E., Segevb, R., and Schneidman, E. (2011). Sparse low-order interaction network underlies a highly correlated and learnable neural population code. *Proc. Natl. Acad. Sci. U.S.A.* 108, 9679–9684. doi: 10.1073/pnas.1019641108

Geisler, C., Brunel, N., and Wang, X.-J. (2005). Contributions of intrinsic membrane dynamics to fast network oscillations with irregular neuronal discharges. *J. Neurophysiol.* 94, 4344–4361. doi: 10.1152/jn.00510.2004

Gerstner, W., Kempter, R., van Hemmen, J. L., and Wagner, H. (1996). A neuronal learning rule for sub-millisecond temporal coding. *Nature* 384, 76–78. doi: 10.1038/383076a0

Graupner, M., and Brunel, N. (2010). Mechanisms of induction and maintenance of spike-timing dependent plasticity in biophysical synapse models. *Front. Comput. Neurosci.* 4:136. doi: 10.3389/fncom.2010.00136

Gulyás, A. I., and Freund, T. T. (2015). Generation of physiological and pathological high frequency oscillations: the role of perisomatic inhibition in sharp-wave ripple and interictal spike generation. *Curr. Opin. Neurobiol.* 31, 26–32. doi: 10.1016/j.conb.2014.07.020

Gupta, A., Wang, Y., and Markram, H. (2000). Organizing principles for a diversity of GABAergic interneurons and synapses in the neocortex. *Science* 287, 273–278. doi: 10.1126/science.287.5451.273

Hansel, D., Mato, G., Meunier, C., and Neltner, L. (1998). On numerical simulations of integrate-and-fire neural networks. *Neural Comput.* 10, 467–483. doi: 10.1162/089976698300017845

Henkin, Y., Sadeh, M., Kivity, S., Shabtai, E., Kishon-Rabin, L., and Gadoth, N. (2005). Cognitive function in idiopathic generalized epilepsy of childhood. *Dev. Med. Child Neurol.* 47, 126–132. doi: 10.1017/S0012162205000228

Hochbaum, D. R., Zhao, Y., Farhi, S. L., Klapoetke, N., Werley, C. A., Kapoor, V., et al. (2014). All-optical electrophysiology in mammalian neurons using engineered microbial rhodopsins. *Nat. Methods* 11, 825–833. doi: 10.1038/nmeth.3000

Hughes, J. R. (2009). Absence seizures: a review of recent reports with new concepts. *Epilepsy Behav.* 15, 404–412. doi: 10.1016/j.yebeh.2009.06.007

Jacob, V., Petreanu, L., Wright, N., Svoboda, K., and Fox, K. (2012). Regular spiking and intrinsic bursting pyramidal cells show orthogonal forms of experience-dependent plasticity in layer V of barrel cortex. *Neuron* 73, 391–404. doi: 10.1016/j.neuron.2011.11.034

Jaeger, H. (2003). “Adaptive nonlinear system identification with echo state networks,” in *Advances in Neural Information Processing System 15*, eds S. Becker, S. Thrun, and K. Obermayer (Cambridge: MIT Press), 593–600.

Jaeger, H., and Haas, H. (2004). Harnessing nonlinearity: predicting chaotic systems and saving energy in wireless communication. *Science* 304:78. doi: 10.1126/science.1091277

Ji, D., and Wilson, M. A. (2007). Coordinated memory replay in the visual cortex and hippocampus during sleep. *Nat. Neurosci.* 10, 100–107. doi: 10.1038/nn1825

Jutras, M. J., Fries, P., and Buffalo, E. A. (2009). Gamma-band synchronization in the macaque hippocampus and memory formation. *J. Neurosci.* 29, 12521–12531. doi: 10.1523/JNEUROSCI.0640-09.2009

Kamioka, H., Maeda, E., Jimbo, Y., Robinson, H. P. C., and Kawana, A. (1996). Spontaneous periodic synchronized bursting during formation of mature patterns of connections in cortical cultures. *Neurosci. Lett.* 206, 109–112. doi: 10.1016/S0304-3940(96)12448-4

Keck, T., Mrsic-Flogel, T. D., Afonso, M. V., Eysel, U. T., Bonhoeffer, T., and Hübener, M. (2008). Massive restructuring of neuronal circuits during functional reorganization of adult visual cortex. *Nat. Neurosci.* 11, 1162–1167. doi: 10.1038/nn.2181

Laje, R., and Buonomano, D. V. (2013). Robust timing and motor patterns by taming chaos in recurrent neural networks. *Nat. Neurosci.* 16, 925–933. doi: 10.1038/nn.3405

Lin, I.-C., Okun, M., Carandini, M., and Harris, K. D. (2015). The nature of shared cortical variability. *Neuron* 87, 644–656. doi: 10.1016/j.neuron.2015.06.035

Litwin-Kumar, A., and Doiron, B. (2014). Formation and maintenance of neuronal assemblies through synaptic plasticity. *Nat. Commun.* 5:5319. doi: 10.1038/ncomms6319

Maass, W., Natschläger, T., and Markram, H. (2002). Real-time computing without stable states: a new framework for neural computation based on perturbations. *Neural Comput.* 14, 2531–2560. doi: 10.1162/089976602760407955

Maheshwari, A., and Noebels, J. L. (2014). Monogenic models of absence epilepsy: windows into the complex balance between inhibition and excitation in thalamocortical microcircuits. *Prog. Brain Res.* 213, 223–252. doi: 10.1016/B978-0-444-63326-2.00012-0

Markram, H., Gerstner, W., and Sjöström, P. J. (2012). Spike-timing-dependent plasticity: a comprehensive overview. *Front. Synaptic Neurosci.* 4:2. doi: 10.3389/fnsyn.2012.00002

Masquelier, T. (2013). Neural variability, or lack thereof. *Front. Comput. Neurosci.* 7:7. doi: 10.3389/fncom.2013.00007

Mazor, O., and Laurent, G. (2005). Transient dynamics versus fixed points in odor representations by locust antennal lobe projection neurons. *Neuron* 48, 661–673. doi: 10.1016/j.neuron.2005.09.032

Mongillo, G., Barak, O., and Tsodyks, M. (2008). Synaptic theory of working memory. *Science* 319, 1543–1546. doi: 10.1126/science.1150769

Monteforte, M., and Wolf, F. (2010). Dynamical entropy production in spiking neuron networks in the balanced state. *Neural Comput.* 105:268104. doi: 10.1103/physrevlett.105.268104

Nádasdy, Z., Hirase, H., Czurkó, A., Csicsvari, J., and Buzsáki, G. (1999). Replay and time compression of recurring spike sequences in the hippocampus. *J. Neurosci.* 19, 9497–9507.

Narayanan, N. S., and Laubach, M. (2009). “Methods for studying functional interactions among neuronal populations,” in *Dynamic Brain Imaging*, ed F. Hyder (Totowa, NJ: Humana Press), 135–165. doi: 10.1007/978-1-59745-543-5_7

Nirenberg, S., and Latham, P. E. (2003). Decoding neuronal spike trains: how important are correlations? *Proc. Natl. Acad. Sci. U.S.A.* 100, 7348–7353. doi: 10.1073/pnas.1131895100

Nolan, M. A., Redoblado, M. A., Lah, S., Sabaz, M., Lawson, J. A., Cunningham, A. M., et al. (2004). Memory function in childhood epilepsy syndromes. *J. Paediatr. Child Health* 40, 22–27. doi: 10.1111/j.1440-1754.2004.00284.x

Ocker, G. K., Litwin-Kumar, A., and Doiron, B. (2015). Self-organization of microcircuits in networks of spiking neurons with plastic synapses. *PLoS Comput. Biol.* 11:e1004458. doi: 10.1371/journal.pcbi.1004458

O'Connor, D. H., Peron, S. P., Huber, D., and Svoboda, K. (2010). Neural activity in barrel cortex underlying vibrissa-based object localization in mice. *Neuron* 67, 1048–1061. doi: 10.1016/j.neuron.2010.08.026

Ostojic, S. (2014). Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons. *Nat. Neurosci.* 17, 594–600. doi: 10.1038/nn.3658

Ostojic, S., Brunel, N., and Hakim, V. (2009). How connectivity, background activity, and synaptic properties shape the cross-correlation between spike trains. *J. Neurosci.* 29, 10234–10253. doi: 10.1523/JNEUROSCI.1275-09.2009

Padmanabhan, K., and Urban, N. N. (2010). Intrinsic biophysical diversity decorrelates neuronal firing while increasing information content. *Nat. Neurosci.* 13, 1276–1282. doi: 10.1038/nn.2630

Panzeri, S., Golledge, H. D. R., Zheng, F., Pola, G., Blanche, T. J., Tovée, M. J., et al. (2002). The role of correlated firing and synchrony in coding information about single and separate objects in cat V1. *Neurocomputing* 44–46, 579–584. doi: 10.1016/S0925-2312(02)00443-5

Peron, S., and Svoboda, K. (2011). From cudgel to scalpel: toward precise neural control with optogenetics. *Nat. Methods* 8, 30–34. doi: 10.1038/nmeth.f.325

Roxin, A., Brunel, N., Hansel, D., Mongillo, G., and van Vreeswijk, C. (2011). On the distribution of firing rates in networks of cortical neurons. *J. Neurosci.* 31, 16217–16226. doi: 10.1523/JNEUROSCI.1677-11.2011

Salkoff, D. B., Zagha, E., Yüzgeç, Ö., and McCormick, D. A. (2015). Synaptic mechanisms of tight spike synchrony at gamma frequency in cerebral cortex. *J. Neurosci.* 35, 10236–10251. doi: 10.1523/JNEUROSCI.0828-15.2015

Schneidman, E., Berry, M. J. II, Sege, R., and Bialek, W. (2006). Weak pairwise correlations imply strongly correlated network states in a neural population. *Nature* 440, 1007–1012. doi: 10.1038/nature04701

Schulz, D. J., Goaillard, J.-M., and Marder, E. (2006). Variable channel expression in identified single and electrically coupled neurons in different animals. *Nat. Neurosci.* 9, 356–362. doi: 10.1038/nn1639

Schwindt, P., and Crill, W. (1999). Mechanisms underlying burst and regular spiking evoked by dendritic depolarization in layer 5 cortical pyramidal neurons. *J. Neurophysiol.* 81, 1341–1354.

Sederberg, P. B., Schulze-Bonhage, A., Madsen, J. R., Bromfield, E. B., McCarthy, D. C., Brandt, A., et al. (2007). Hippocampal and neocortical gamma oscillations predict memory formation in humans. *Cereb. Cortex* 17, 1190–1196. doi: 10.1093/cercor/bhl030

Shafi, M., Zhou, Y., Quintana, J., Chow, C., Fuster, J., and Bodner, M. (2007). Variability in neuronal activity in primate cortex during working memory tasks. *Neuroscience* 146, 1082–1108. doi: 10.1016/j.neuroscience.2006.12.072

Softky, W. R., and Koch, C. (1993). The highly irregular firing of cortical cells is inconsistent with temporal integration of random epsps. *J. Neurosci.* 13, 334–350.

Sompolinsky, H., Crisanti, A., and Sommers, H. J. (1988). Chaos in random neural networks. *Phys. Rev. Lett.* 61:259. doi: 10.1103/PhysRevLett.61.259

Song, S., Sjöström, P. J., Reigl, N., Nelson, S., and Chklovski, D. B. (2005). Highly nonrandom features of synaptic connectivity in local cortical circuits. *PLoS Biol.* 3:e68. doi: 10.1371/journal.pbio.0030068

Stark, E., Roux, L., Eichler, R., Senzai, Y., Royer, S., and Buzsáki, G. (2014). Pyramidal cell-interneuron interactions underlie hippocampal ripple oscillations. *Neuron* 83, 467–480. doi: 10.1016/j.neuron.2014.06.023

Strauss, U., Kole, M. H., Bräuer, A. U., Pahnke, J., Bajorat, R., Rolfs, A., et al. (2004). An impaired neocortical Ih is associated with enhanced excitability and absence epilepsy. *Eur. J. Neurosci.* 19, 3048–3058. doi: 10.1111/j.0953-816X.2004.03392.x

Sussillo, D., and Abbott, L. (2009). Generating coherent patterns of activity from chaotic neural networks. *Neuron* 63, 544–557. doi: 10.1111/j.0953-816X.2004.03392.x

Tiesinga, P., and Sejnowski, T. J. (2009). Cortical enlightenment: are attentional gamma oscillations driven by ING or PING? *Neuron* 63, 727–732. doi: 10.1016/j.neuron.2009.09.009

Toyoizumi, T., and Abbott, L. F. (2011). Beyond the edge of chaos: amplification and temporal integration by recurrent networks in the chaotic regime. *Phys. Rev. E* 84:051908. doi: 10.1103/PhysRevE.84.051908

Trousdale, J., Hu, Y., Shea-Brown, E., and Josić, K. (2012). Impact of network structure and cellular response on spike time correlations. *PLoS Comput. Biol.* 8:e1002408. doi: 10.1103/PhysRevE.84.051908

Turrigiano, G. (2011). Too many cooks? intrinsic and synaptic homeostatic mechanisms in cortical circuit refinement. *Annu. Rev. Neurosci.* 34, 89–103. doi: 10.1146/annurev-neuro-060909-153238

Turrigiano, G. G., and Nelson, S. (2004). Homeostatic plasticity in the developing nervous system. *Nat. Rev. Neurosci.* 5, 97–107. doi: 10.1038/nrn1327

van Vreeswijk, C., and Sompolinsky, H. (1998). Chaotic balanced state in a model of cortical circuits. *Neural Comput.* 10, 1321–1371. doi: 10.1162/089976698300017214

VanRullen, R., Busch, N. A., Drewes, J., and Dubois, J. (2011). Ongoing EEG phase as a trial-by-trial predictor of perceptual and attentional variability. *Front. Psychol.* 2:60. doi: 10.3389/fpsyg.2011.00060

Vida, I., Bartos, M., and Jonas, P. (2006). Shunting inhibition improves robustness of gamma oscillations in hippocampal interneuron networks by homogenizing firing rates. *Neuron* 49, 107–117. doi: 10.3389/fpsyg.2011.00060

Wang, X.-J., and Buzsáki, G. (1996). Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. *J. Neurosci.* 16, 6402–6413.

Xiang, Z., Huguenard, J. R., and Prince, D. A. (1998). GABAA receptor-mediated currents in interneurons and pyramidal cells of rat visual cortex. *J. Physiol.* 506, 715–730. doi: 10.1111/j.1469-7793.1998.715bv.x

Yamamoto, J., Suh, J., Takeuchi, D., and Tonegawa, S. (2014). Successful execution of working memory linked to synchronized high-frequency gamma oscillations. *Cell* 157, 845–857. doi: 10.1016/j.cell.2014.04.009

Yamazaki, T., and Tanaka, S. (2007). The cerebellum as a liquid state machine. *Neural Netw.* 20, 290–297. doi: 10.1016/j.neunet.2007.04.004

Yang, Y., and Calakos, N. (2013). Presynaptic long-term plasticity. *Front. Synaptic Neurosci.* 5:8. doi: 10.3389/fnsyn.2013.00008

Zador, A. (1998). Impact of synaptic unreliability on the information transmitted by spiking neurons. *J. Neurophysiol.* 79, 1219–1229.

Zhou, F.-M., and Hablitz, J. J. (1998). AMPA receptor-mediated EPSCs in rat neocortical layer II/III interneurons have rapid kinetics. *Brain Res.* 780, 166–169. doi: 10.1016/S0006-8993(97)01311-5

Keywords: spike pattern structure, synaptic plasticity, efficacy variability, STDP, synaptic homeostasis, spike shuffling methods

Citation: Bi Z and Zhou CS (2016) Spike Pattern Structure Influences Synaptic Efficacy Variability under STDP and Synaptic Homeostasis. II: Spike Shuffling Methods on LIF Networks. *Front. Comput. Neurosci.* 10:83. doi: 10.3389/fncom.2016.00083

Received: 19 March 2016; Accepted: 25 July 2016;

Published: 09 August 2016.

Edited by:

Pulin Gong, University of Sydney, AustraliaCopyright © 2016 Bi and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zedong Bi, zedong.bi@outlook.com

Changsong Zhou, cszhou@hkbu.edu.hk