## Abstract

The strength of cortical synapses distributes lognormally, with a long tail of strong synapses. Various properties of neuronal activity, such as the average firing rates of neurons, the rate and magnitude of spike bursts, the magnitude of population synchrony, and the correlations between presynaptic and postsynaptic spikes, also obey lognormal-like distributions reported in the rodent hippocampal CA1 and CA3 areas. Theoretical models have demonstrated how such a firing rate distribution emerges from neural network dynamics. However, how the other properties also display lognormal patterns remain unknown. Because these features are likely to originate from neural dynamics in CA3, we model a recurrent neural network with the weights of recurrent excitatory connections distributed lognormally to explore the underlying mechanisms and their functional implications. Using multi-timescale adaptive threshold neurons, we construct a low-frequency spontaneous firing state of bursty neurons. This state well replicates the observed statistical properties of population synchrony in hippocampal pyramidal cells. Our results show that the lognormal distribution of synaptic weights consistently accounts for the observed long-tailed features of hippocampal activity. Furthermore, our model demonstrates that bursts spread over the lognormal network much more effectively than single spikes, implying an advantage of spike bursts in information transfer. This efficiency in burst propagation is not found in neural network models with Gaussian-weighted recurrent excitatory synapses. Our model proposes a potential network mechanism to generate sharp waves in CA3 and associated ripples in CA1 because bursts occur in CA3 pyramidal neurons most frequently during sharp waves.

**SIGNIFICANCE STATEMENT** The wiring structure of local cortical networks is known to be far from random. Here, we propose a recurrent neural network model with a long-tailed synaptic weight distribution to account for various properties of the synchronous firing of hippocampal neurons observed typically during sharp wave ripples. Furthermore, the model predicts that pathways of strong synapses route spike bursts much more efficiently than single spikes. Sharp wave ripples are crucial for memory encoding, but the underlying mechanism remains unknown. Our model suggests the crucial role of internal dynamics of nonrandom hippocampal circuits in generating and routing such activity patterns. Our results will have significant implications in understanding the mechanism of memory encoding by the hippocampus.

## Introduction

The advent of unbiased, large-scale recordings has started to reveal strikingly large differences in activities among neurons in local cortical areas. These differences are often seen in a skewed distribution of the average firing rates of individual neurons. In hippocampal CA1 and CA3 areas, firing rates as well as several other prominent features of population neural activity were shown to obey similar lognormal statistics (Buzsáki and Mizuseki, 2014). Because lognormal patterns are also seen in other regions of the brain, they may be part of the fundamental principles of computation by the brain. However, the origin of these lognormal patterns is poorly understood. In addition, little is known about the computational implications of the highly skewed statistical structure of cortical circuits and the resultant neural activity.

Here, we use computations to investigate whether the lognormal distribution of synaptic weights accounts for the various skewed distributions observed in hippocampal neuronal activities. The strengths of synapses, typically assessed as amplitudes of EPSPs among neurons, are distributed lognormally in the hippocampus (Ikegaya et al., 2013) and neocortex (Song et al., 2005; Sarid et al., 2007; Lefort et al., 2009). Spine sizes are also lognormally distributed (Yasumatsu et al., 2008; Loewenstein et al., 2011). In addition, the spontaneous firing rates of individual neurons are distributed lognormally in recurrent network models with lognormally weighted synapses (Koulakov et al., 2009; Teramae et al., 2012). Such a firing rate distribution can emerge from a nonlinear neuronal response curve, which also appears in a low-frequency regime of lognormally weighted recurrent networks (Roxin et al., 2011). Lognormal neural networks maintain sparse and highly irregular spontaneous firing by internally generating noise for stochastic resonance effects (Teramae et al., 2012). These results prompted us to examine whether the multiple lognormal patterns observed in the hippocampus originate from lognormally weighted synapses in a computational model of a CA3 recurrent network.

Hippocampal pyramidal neurons *in vivo* exhibit both single spikes and complex spike bursts (Ranck, 1973; Harris et al., 2001). To describe bursting activity, we used a multi-timescale adaptive threshold (MAT) neuron, known for its ability to mimic the irregular responses of cortical neurons (Jolivet et al., 2008; Kobayashi et al., 2009). We show that the low-frequency spontaneous activity in our network model well replicates the lognormal features of hippocampal activity. Moreover, the bursting activity exhibited by our model resembles the complex spike bursts of hippocampal pyramidal neurons (Harris et al., 2001) observed frequently, but not always, during sharp-wave ripples (SWRs) (Buzsáki et al., 2002; Lee and Wilson, 2002; Foster and Wilson, 2006), which are transient oscillatory activities in CA1 and play a crucial role in memory consolidation (Girardeau et al., 2009; Ego-Stengel and Wilson, 2010; Carr et al., 2011). The SWRs originate from the activation of CA3 neural ensembles and typically occur when sensory influences on the hippocampus decrease during slow-wave sleep and immobility. Although it is unknown how synchronized activity occurs spontaneously in CA3 (Csicsvari et al., 2000), our results suggest that a lognormal distribution of synaptic weights underlies the generation of the spontaneous synchronous firing.

## Materials and Methods

##### Network.

The structure of the network connections is essentially the same as that used previously (Teramae et al., 2012). The present network model consists of 12,000 (10,000 excitatory and 2000 inhibitory) MAT neurons connected randomly. Each neuron receives on average 1000 excitatory and 1000 inhibitory synaptic inputs. The MAT neuron model has two parallel dynamics for the membrane potential and spike threshold, with the characteristics of the spike trains tailored by the values of the parameters in the spike threshold dynamics. The membrane dynamic is the same as that in a leaky integrate-and-fire neuron, except that the membrane potential is not reset to the resting potential after spike generation as follows:
where the membrane time constant τ* _{m}* is 20 ms for excitatory neurons and 10 ms for inhibitory neurons, and

*RI*is the total input generated by excitatory and inhibitory reverberating synaptic inputs. The reversal potentials of leaky, excitatory, and IPSCs are

*V*= −70 mV,

_{L}*V*= 0 mV, and

_{E}*V*= −80 mV, respectively.

_{I}The adaptive threshold θ* _{i}*(

*t*) has fast (

*k*= 1) and slow (

*k*= 2) components and obeys the following: where

*t*is the

_{j}*j*-th spike time of the

*i*-th neuron, ω the resting value of the threshold, τ

*the*

_{k}*k*-th time constant, and α

*the weight of*

_{k}*k*-th component. Each neuron generates a spike when its membrane potential reaches the instantaneous value of the spike threshold, and the threshold value is increased according to Equations 3 and 4 at every spike firing and decays exponentially to the resting value ω when the neuron does not fire. The total synaptic input to the neuron is interrupted during a refractory period of 1 ms. In excitatory neurons, the magnitudes of the fast component α

_{1}are drawn randomly from a normal distribution with the mean 1.5 mV, and the SD 0.25 mV, and the magnitude of the slow component is set as α

_{2}= 0.5 mV. In inhibitory neurons, the magnitudes of fast and slow components are set as α

_{1}= 3 mV and α

_{2}= 0. Excitatory and inhibitory neurons have identical time constants, τ

_{1}= 10 ms, for fast components, and excitatory neurons have τ

_{2}= 200 ms for slow components. For both neuronal types, ω = −55 mV.

The time course of the synaptic conductance g(*t*) of the *i*-th neuron is described as follows:
where the decay constant τ* _{S}* = 2 (ms),

*d*is the synaptic delay,

_{j}*G*the synaptic weight from the

_{ij}*j*-th to the

*i*-th neuron, and δ(

*t*) is Dirac's δ function. The delays of the excitatory-to-excitatory synaptic connections are uniformly distributed in the range of 1–3 ms, and the delays of the other connections have an identical value of 1 ms. Setting the different delays enhances the stability of spontaneous activity (Teramae et al., 2012). The weights of the excitatory-to-excitatory synaptic connections are drawn from the following lognormal distribution independently for individual neurons: where σ = 1.0 and μ = log(0.2) − σ

^{2}are the mean and SD of the variable's natural logarithm. Equation 6 mimics the typical EPSP amplitude distributions observed in the experiment (Song et al., 2005; Lefort et al., 2009). Here, any unrealistic value that is >20 mV is avoided by drawing a new value from the distribution. Thus, in our model, each excitatory neuron has terminals with similar heterogeneous synaptic connections weighted as lognormal. The weights of excitatory-to-inhibitory, inhibitory-to-excitatory, and inhibitory-to-inhibitory synaptic connections are uniform for mathematical simplicity as the corresponding experimental distributions are not known well, and their values are fixed at 0.018, 0.0035, and 0.0025, respectively, to ensure stable spontaneous activity. Synaptic transmission fails at excitatory-to-excitatory synapses at a rate depending on the amplitude of EPSP as follows: The network is randomly connected so that the innervation probabilities of excitatory and inhibitory neurons are 0.1 and 0.5, respectively. Brief external Poisson spike trains are applied to all neurons to initiate spontaneous activity in the network model.

##### Gaussian-weighted network of MAT model.

This approach inherits its structure from the original lognormal model, maintaining its parameters and characteristics but for the aspects we mention below. The weights of excitatory-to-excitatory synapses were drawn from a normal distribution as follows:
where μ* _{G}* and σ

*are the mean and SD. As this distribution allows negative values, we have truncated its results, allowing synaptic weights to take only non-negative and realistic numbers <20 mV. To make a fairer comparison between this approach and its lognormal counterpart, we ask that μ*

_{G}*should be so that the mean for the resulting truncated distribution approaches μ*

_{G}*(with an error of 0.001), and σ*

_{L}*= σ*

_{G}*, where μ*

_{L}*and σ*

_{L}*are the mean and SD of the lognormal counterpart, respectively. They are given as μ*

_{L}*=*

_{L}*e*

^{μ+σ2/2}and σ

*=*

_{L}*= −0.114082 and σ*

_{G}*= 0.896338.*

_{G}As changing the weight distribution function broke our previous spontaneous activity, we lowered the resting threshold value ω to −60 mV, a procedure that showed itself successful for achieving longer stability in the lognormal network. Such results, however, were not found in the Gaussian network, even though durations for dying activities were increased in general. Therefore, in the presence of the lowered ω, we further applied external background input represented by Poisson spike trains at 1 Hz to a randomly chosen subgroup of 5% of the excitatory neuron population. By analyzing the results for the duration of activities as shown in Figure 9*A*, we chose the weight values for inhibitory-to-excitatory and inhibitory-to-inhibitory synapses as 0.0023 and 0.0037, respectively, keeping the original value of the excitatory-to-inhibitory weight, 0.0018. This balance gives us 86% of chance of having a simulation with activity longer than 5 s, and only such activities were used to extract our results. Activities lasting 10 s were prioritized for any single analysis.

##### Synchrony and population synchrony.

The words “synchrony” and “population synchrony” are used herein to describe distinct phenomena. In Figure 1*F*, the cross-correlograms of each neuron pair were calculated from presynaptic and postsynaptic spikes separated by intervals within (−20, 20 ms). The cross-correlograms over all neuron pairs were then averaged, the mean value (over bins) from the peak value was subtracted, and the resultant value normalized using the mean to define synchrony. In Figures 4 and 6, the magnitude of the population synchrony was defined as the fraction of neurons that fired at least once in a sliding time window of 30 ms. We calculated the instantaneous firing rates of a neuron for population synchrony of given magnitude level by dividing by 0.03 s the average number of spikes generated by the neuron over repeated population synchrony of that level. The time step for the sliding time windows was 10 ms.

##### Correlation coefficient.

In each time window, a correlation coefficient between spike trains of two excitatory neurons was calculated as follows:
where *n*_{1} and *n*_{2} are the spike counts of neuron 1 and neuron 2 in the time window, Cov(*n*_{1}, *n*_{2}) is the covariance of spike counts between the two neurons, and Var is the variance in spike counts.

##### Selectivity indices.

Selectivity indices were defined to assess the change in dominance among excitatory neurons with an increased magnitude of population synchrony. Two values from each Lorenz curve *y* = *L*(*x*) of the distributions of firing rates are exploited; *y*_{1/2} = *L*(1/2) gives the fraction of spikes emitted by the lower 50% of excitatory neurons to the overall spikes emitted by all excitatory neurons. Similarly, *x*_{1/2} = 1 − *L*^{−1}(1/2) represents the most active fraction of excitatory neurons generating half of all spikes to all excitatory neurons. The two quantities for the long-term averaged firing rates of individual neurons (*x*_{1/2}^{(0)}, *y*_{1/2}^{(0)}) and the firing rates for each class of the ranked time windows (*x*_{1/2}^{(r)}, *y*_{1/2}^{(r)}) (see Fig. 4*D*) are calculated, where the index *r* refers to the six ranks from A to F. The value of *y*_{1/2}^{(r)} is smaller for a stronger population synchrony, as is the value of *x*_{1/2}^{(r)}. Now, the selectivity index of each class *r* is defined by *x*_{1/2}^{(r)}/*x*_{1/2}^{(0)}, or similarly by *y*_{1/2}^{(r)}/*y*_{1/2}^{(0)}, which is unity if spikes from each neuron contribute equally to the count of all spikes. In Figure 6, the former index is shown in orange, whereas the latter index is shown in yellow. The selectivity indices are small in the time windows for highly synchronous events.

##### Statistical test for sequence propagation.

Sequences propagating through the pathways formed by extremely strong synapses are expected to be statistically significant. We examined this by the following three methods. In the first method, we calculated the cumulative distributions of sequences following bursts of a reference excitatory neuron and having length *n* ≥ 11. We tracked only such sequences of single spikes and bursts that propagate along a cascade of the strongest 25 (instead of 5) synapses from each participating neuron and defined the temporal order of participating neurons using the first spikes within bursts if they generate bursts. We used the large number of strong synapses to pool many different sequences. The observed sequences were arranged in the decreasing order of their frequencies, and a cumulative distribution was generated. Then, as a control, we calculated a similar cumulative distribution for sequences propagating in arbitrary time periods temporally far from any activity from the reference neuron. We may term these sequences “imaginary sequences.” We examined by Kolmogorov–Smirnov (KS) test whether the distributions of actual and imaginary sequences are statistically different, which clarifies whether some sequences following bursts of a reference neuron occur statistically more often when the neuron is active than in other instances.

In the second method, we used different control sequences, in which the temporal order of the (*n*−1) neurons following a reference neuron was randomly permuted. We constructed cumulative distributions for actual and permuted sequences, and performed the KS test to examine whether they are statistically different. If this is indeed the case, the sequences following the firing of a reference neuron occur in specific temporal orders more often than expected by chance.

In the third method, we changed our control sequences to those following single spikes of a reference neuron. We analyzed cumulative distributions for burst-following and single-spike-following sequences by KS test to examine whether the firing pattern can influence the occurrence of statistically significant sequences.

## Results

### Dynamics of bursty neurons in the recurrent network

We examined whether the lognormal statistics observed in single-cell and neural ensemble activities in the hippocampus originated from interactions between the lognormal distributions of synaptic weights and bursting properties of hippocampal neurons. To this end, we constructed a network model of the CA3 region with recurrent excitatory connections in which the EPSP amplitudes were lognormally distributed in individual neurons (Fig. 1*A*). We used the MAT model to describe both single spikes and spike bursts of CA3 neurons, where a spike burst was defined as a series of two or more spikes with <6 ms intervals (Ranck, 1973; Harris et al., 2001), unless otherwise stated. Thus, the characteristic feature of the present network model was that each neuron demonstrated a bursting property. The MAT model is mathematically simple and enabled us to perform numerical simulations in large-scale networks. Because the proportion of very strong synapses in a lognormal network is small, individual neurons may acquire very strong connections only when they receive sufficient synaptic inputs, that is, when the network is sufficiently large.

Except for the bursty property of neurons, the present network model shares various dynamic properties with the lognormal network model of integrate-and-fire neurons (Teramae et al., 2012). As shown in that study, in this latter type of network, recurrent synaptic inputs to many weak synapses generate subthreshold noise useful for maximizing the signal-to-noise ratio of spike transmissions at extremely strong recurrent synapses. Through this mechanism, a brief application of external stimuli elicited a stable spontaneous activity in the neural network examined in the present study in which time-varying excitatory and inhibitory activities were balanced (Fig. 1*B*). During a spike burst, the spike threshold of a MAT neuron was incremented by a small amount whenever the membrane potential hit the threshold from below; otherwise, the spike threshold declined exponentially to a baseline level. A burst eventually stopped because the threshold rapidly increased during the burst (Fig. 1*C*). Only occasional strong excitatory inputs evoked single spikes or spike bursts from postsynaptic neurons, keeping neuronal firing very sparse and irregular (Fig. 1*D*,*E*). Hippocampal pyramidal neurons exhibit a refractory period of ∼100 ms between successive spike bursts due to the slow inactivation of sodium channels (Mickus et al., 1999; Su et al., 2001), and such burst refractoriness emerged in the MAT model from the slower decay constant of the spike threshold τ_{2} (see Materials and Methods). The network activity was self-sustained without external background input, so the neural dynamics were deterministic. Our model robustly generated stable spontaneous activity across a wide range of parameter space (Fig. 1*F*).

The quantitative features of the spike bursts in our model were in good agreement with the experimental observations of the hippocampus (Harris et al., 2001). Individual pyramidal cells exhibit both single spikes and spike bursts, and typical interspike intervals (ISIs) during a burst are ∼2–4 ms. The ISI histogram in our model was bimodal, reflecting the generation of single spikes and spike bursts, although the ISIs in a burst were typically 1–2 ms, slightly shorter than the experimentally derived hippocampal values (Fig. 2*A*). In the hippocampus, the number of spikes within a burst varies from event to event in a single neuron as well as across neurons, and bursts with more spikes occur less frequently. Typically, the probability that a pyramidal cell generates a burst of length *n* (i.e., a burst contains *n* spikes) decays exponentially as a function of *n* for the majority of pyramidal cells. However, some pyramidal cells are exceptionally bursty, and their distributions of burst length are supraexponential. Similar to experimental results, the majority of our model neurons showed an exponential distribution of the burst length (Fig. 2*B*), whereas exceptionally bursty neurons showed a clear trend away from exponential behavior (Fig. 2*C*). However, the correlations between the generation probability and burst length were continuously distributed over the entire excitatory neuron population, indicating no clear border between the two categories of bursty model neurons (Figs. 2*D*, 3*B*,*C*). In hippocampal pyramidal cells, the ISIs between successive spikes are progressively prolonged in a spike burst, which was also the case in our model (Fig. 2*E*). In the spontaneous activity of the model, each excitatory neuron fired when the balance between the excitatory and IPSCs was biased slightly toward excitation (Fig. 2*F*). This result is consistent with a similar finding in spontaneous sharp wave-like events that occur in slice preparations of the hippocampus (Behrens et al., 2005). Pyramidal cells have been shown to exhibit sharp wave-locked spiking when excitation briefly dominates over inhibition in CA3 (Hájos et al., 2013) and CA1 (Mizunuma et al., 2014).

### Firing rate distribution and burstiness of neurons

In the hippocampus, the distribution of the mean firing rate and burstiness (i.e., the burst-event rate and burst index, where the latter refers to the ratio of spikes within a burst to all generated spikes) in a neural population displays lognormal patterns regardless of brain states, such as slow-wave or rapid eye movement sleep, running, and immobility (Mizuseki and Buzsáki, 2013). As shown in Figure 3*A–C*, excitatory neurons in the present model displayed lognormal patterns in the distributions of those quantities. The distributions for inhibitory neurons were narrow and not regarded as lognormal. For excitatory neurons, the median, upper, and lower quartiles for the average firing rates were 0.87, 1.48, and 0.48 Hz, respectively; for burst-event rates, they were 0.027, 0.079, and 0.008 Hz, respectively; and for the burst index, they were 0.10, 0.18, and 0.05 Hz, respectively. In Figure 3*D*, we constructed Lorenz plots for the average firing rates of excitatory and inhibitory neurons to calculate the Gini coefficient (Gini, 1921), which measures inequality among the values of a frequency distribution, becoming zero if all neurons equally contribute to the distribution. The coefficient was calculated as 0.41 for excitatory neurons (Fig. 3*E*). All of these values obtained for our model were in good agreement with those measured for hippocampal pyramidal cells. In addition, the burst-event rate and burst index were positively correlated with the firing rate, as experimentally observed for the burst index (Fig. 3*F*). These results showed that population bursts are dominated by a relatively small number of highly active neurons in the lognormal network. However, unlike results from experimental models, the present model lacked low-frequency firing neurons displaying strong burstiness, suggesting that a mechanism additional to the lognormal pattern of excitatory recurrent synapses boosts burst generation in hippocampal neurons. For instance, although most hippocampal CA3 pyramidal cells produce bursts, the patterns of spike firing exhibit a considerable cell-to-cell variability (Brown and Randall, 2009). This diverse heterogeneity in excitatory neuron population was not incorporated into our model. Furthermore, our model also did not take into account the highly nonrandom wiring topology of hippocampal circuits (Takahashi et al., 2010). The nonrandom connectivity, especially that of inhibitory circuits, may further enhance the heterogeneity of firing patterns across the excitatory cell population.

As shown above, our model predicts that some activity indices also exhibit lognormal patterns in inhibitory neurons. However, the quantitative results obtained for inhibitory neurons significantly depended on how we modeled the unknown details of hippocampal inhibitory circuits. Therefore, the predictions for the inhibitory neurons should be considered qualitative rather than quantitative.

Although a rigorous treatment is difficult, we may heuristically argue the underlying mechanism of lognormal distributions in the network of MAT model. Leaky-integrate-fire (LIF) neurons show a lognormal firing rate distribution if the average firing rate of spontaneous activity is sufficiently low (e.g., Roxin et al., 2011; Teramae and Fukai, 2014). This follows from the facts that the sum of many independent identical spike inputs obeys a Gaussian distribution and that the firing rate versus input fluctuation (*f*-Δ*I* curve) of LIF neuron is an exponential function in the low-frequency regime (Brunel, 2000). Although the *f*-Δ*I* curve is not analytically known for MAT model, we can expect that the event rate for single spikes and bursts should obey lognormal distributions. Unlike in LIF neuron, spiking does not reset the membrane potential in MAT model. Rather, threshold is reset and then evolves with fast and slow time scales (τ_{1} ≪ τ_{2}). Because the typical interburst interval *T* is much larger than the fast time scale responsible for burst generation (Fig. 2*A*), we may treat multiple spikes within a burst effectively as a single firing event. Furthermore, we may regard threshold value at every burst onset as approximately constant if variations from burst to burst are negligible (Fig. 1*C*). Now, except for boundary conditions, MAT model and LIF neuron obey identical Fokker-Planck equations for the membrane potential (Brunel, 2000). Therefore, the first passage time to threshold, from which the exponential response curve is derived for LIF neuron (Roxin et al., 2011; Teramae and Fukai, 2014), may take approximately identical values in MAT model and LIF neuron. Thus, we expect that the rate *R*_{e} of events (including single spikes and bursts) obeys an approximate lognormal distribution in MAT model. Now, numerical simulations show that ∼90% of the events are single spikes and the remaining 10% are bursts (Fig. 2*B*). If all neurons generate single spikes in the fixed fraction of events, or if this fraction also varies from neuron to neuron lognormally, firing rate *f* and burst event rate *R _{b}* obey approximate lognormal distributions with different scale parameters. Then, because spike threshold increases by α at every postsynaptic spike during a burst, the number of spikes within a burst is approximately calculated as the amplitude of EPSP divided by α, and burst index is estimated as follows:
Equation 10 suggests that burst index also obeys a lognormal distribution, as the EPSP amplitude,

*f*and

*R*all obey lognormal distributions. If variables

_{b}*x*and

*y*are lognormal, variables

*xy*and 1/

*x*also obey lognormal.

The above arguments are far from rigorous. However, they suggest that different mechanisms underlie the various lognormal distributions in the present model. In particular, Equation 10 suggests that a lognormal fit for burst index requires the lognormal property of synaptic weights unlike for firing rate and burst event rate.

### Temporal and spatial distributions of population dynamics

The SWR is a transient activity observed in CA1 and is primarily triggered by synchronized firing of neural ensembles in CA3 (Csicsvari et al., 2000; Behrens et al., 2005; Nakashiba et al., 2009). Population synchrony refers to the synchronous activation of neurons in a specific time window, which is typically in the range of 100 ms (see Materials and Methods). Various indices measuring the magnitude of the population synchrony, for instance, the proportion of neurons activated during the time window, obey lognormal patterns in areas CA1 and CA3, both inside and outside of SWR events (Mizuseki and Buzsáki, 2013). We next examined whether our model generated such distributions for the population synchrony.

To characterize the magnitude of the population synchrony, we counted the numbers of excitatory and inhibitory neurons that fired within a time bin of 30 ms (Fig. 4*A*). The population synchrony exhibited significant temporal fluctuations in its magnitude, and the instantaneous firing rates and the magnitude of the synchrony were strongly correlated between excitatory and inhibitory neurons, implying that an excitatory–inhibitory balance was maintained during the population synchrony (Fig. 4*B*). The population firing rate and magnitude of the population synchrony were also positively correlated, and their distributions were skewed over excitatory and inhibitory neurons (Fig. 4*C*). Similar results have been reported for pyramidal cells (Mizuseki and Buzsáki, 2013).

Population bursts frequently occur during SWR events, although they do not always co-occur (Buzsáki et al., 2002; Lee and Wilson, 2002; Foster and Wilson, 2006). Although the precise relationship between SWR and population burst remains uncertain, it will be intriguing to compare the distribution of population bursts in this model with experimental observations for SWR and population burst. As shown above, the magnitude of the population synchrony varied significantly from event to event in our model. We categorized this magnitude using the standard normal deviate (Fig. 4*D*). Because we did not model neural circuits in CA1, we needed to define a reasonable method to compare the population synchrony in the model with the SWR events in the hippocampus. Figure 4*E* shows a typical temporal sequence of highly synchronous events in the model, which have a standard normal deviate >2. The distribution of interevent intervals was exponential, and the intervals were uncorrelated between the successive events, suggesting that the sporadic outbreaks of high excitation emerging among neural populations follow a Poisson process. The average interevent interval was ∼1.25 Hz, which approximately coincides with the intervals of spontaneous sharp waves measured in area CA3 *in vitro* (∼1.1–1.6 Hz) (Ellender et al., 2010; Hájos et al., 2013) and *in vivo* (0.01–2.0 Hz) (Buzsáki, 1986). The SWR events in CA1, which reflect the population synchrony in CA3, occur at average frequencies of 0.3–1.0 Hz during sleep (Eschenko et al., 2008; Girardeau et al., 2014). These values were also consistent with the prediction of our model. Furthermore, excitatory neurons dominantly generated single spikes in our model (average spike number = 1.43; Fig. 4*F*) similar to those observed during *in vitro* SWRs (spike number = 1.5 ± 0.2) (Behrens et al., 2005). All of these results support our selection of SWR-associated population synchrony in the model.

We next calculated the proportion of neurons participating in each highly synchronous event among the ∼250 randomly chosen neurons and found that the distribution of the proportion obeyed a lognormal pattern (Fig. 4*G*). This result is consistent with the experimental finding that the proportion of pyramidal cells participating in the population synchrony during SWRs shows a lognormal pattern (Mizuseki and Buzsáki, 2013). We also calculated spike correlations between excitatory neuron pairs (see Materials and Methods). As shown experimentally, the correlation coefficients were also distributed lognormally (Fig. 4*H*). As intuitively anticipated, neuron pairs having large correlation coefficients were found more often among those pairs connected by extremely strong synapses than among pairs connected by weaker synapses (Fig. 4*I*).

### Population activity patterns in highly synchronous events

Firing patterns of individual neurons during population synchrony exhibit skewed patterns in experiments (Mizuseki and Buzsáki, 2013). We investigated these patterns in highly synchronous events by calculating the same three indices that were used in the experimental analyses for every excitatory neuron: (1) the proportion of spikes within highly synchronous events to all spikes; (2) the proportion of highly synchronous events in which a neuron fires relative to all events; and (3) the mean number of spikes in a highly synchronous event (Fig. 5*A*). The distributions of these indices displayed skewed patterns over a broad range of values (Fig. 5*B*), with profiles and value ranges of the distributions similar to those observed experimentally. The three indices exhibited positive correlations with the firing rates of individual neurons (Fig. 5*C–E*), consistent with previous experimental findings, except for a positive correlation obtained for the proportion of spikes in highly synchronous events (*r* = 0.17; Fig. 5*C*). In the experiments, this index measured during sleep and awake SWRs is negatively correlated with firing rates (*r* = −0.12 to ∼−0.56), showing a rare qualitative discrepancy between our model and the experiments. The degrees of the correlations for the other two indices were consistent with those in experimental observations.

### Neural selection during highly synchronous events

Whether a highly synchronous event represents the overall activation of the entire neural population or the selective activation of a subpopulation has important functional implications. During a highly synchronous event, both excitatory and inhibitory neurons are strongly modulated by reverberating synaptic input. Accordingly, some neurons are more depolarized by increased excitatory activity, whereas other neurons are more hyperpolarized by increased inhibitory activity (Fig. 6*A*). To achieve further insight into these competitive neural dynamics, we categorized population synchrony into six levels of magnitude ranked from A (the lowest) to F (the highest) according to the standard normal deviate of the distribution of the magnitude, as illustrated in Figure 4*D*. Highly synchronous events are categorized as rank F. We calculated the distributions of instantaneous firing rates over the entire neural population for the various magnitude levels (Fig. 6*B*). It is clear from the figure that the distributions acquired longer tails as the magnitude of the population synchrony increased, implying that the neurons were not uniformly activated during a highly synchronous event. The highly activated neurons were not necessarily those neurons already firing at high average frequencies during spontaneous activity because some of these highly activated neurons were actually inhibited during highly synchronous events (Fig. 6*C*). Indeed, whereas the average firing rate of the neural population increased monotonically with the magnitude of the population synchrony, the proportion of neurons with increased firing rates to all neurons decreased during the epochs of highest synchrony (Fig. 6*D*), indicating that the activity of only a small number of neurons was selectively enhanced in these epochs.

We assessed the contributions of individual neurons to the overall spike count in the network. For the six ranked states of population synchrony, we calculated the proportion of the spikes *r _{s}* generated by the 50% of excitatory neurons that were less active to all spikes and the proportion of the excitatory neurons

*r*that produced 50% of all spikes. For instance, during the strongest population synchrony (rank F), the least active half contributed a mere 17% to all spikes (

_{n}*r*= 0.17), and only 15% of all neurons contributed to half of all spikes (

_{s}*r*= 0.15). We also calculated these proportions for the long-term averaged network activity (

_{n}*R*and

_{s}*R*) and defined the “selectivity index” for each ranked state as

_{n}*r*/

_{s}*R*or

_{s}*r*/

_{n}*R*(see Materials and Methods). We found that

_{n}*R*= 0.22 and

_{s}*R*= 0.21, and thus the selectivity indices were 0.77 and 0.71 for the F-ranked synchrony. As the population synchrony grew stronger, the selectivity indices rapidly decreased (Fig. 6

_{n}*E*), indicating that the excitatory neurons underwent stronger competition during the epochs of stronger population synchrony.

### Efficient propagation of spike bursts through the lognormal network

We showed that our lognormal neural network generated population neural activity that replicated the statistical properties of single spikes and complex spike bursts of hippocampal CA3 neurons. A question arises regarding the functional role of the spike bursts generated by the lognormal network. Spike transmission in a lognormal network is highly probabilistic even at strong synapses (Teramae et al., 2012). Similarly, single spikes and spike bursts showed only a finite “lifetime” in the present lognormal network. Spike bursts, however, exert a stronger impact than single spikes on a postsynaptic neuron (Lisman, 1997), and hence may travel longer distances and more robustly than single spikes in a recurrent network. This possibility is of particular interest from a functional viewpoint because a lognormal network involves a tremendous number of synaptic pathways linked by strong connections.

We therefore investigated how spike bursts and single spikes were distinctly propagated through such synaptic pathways. We selected an arbitrary excitatory neuron from the network and observed how signals propagated through the cascade of neurons starting from this “parent” neuron. To simplify the analysis, we tracked only spikes and bursts that propagated through the five strongest links from any neuron in the cascade. Figure 7*A* depicts an example of such cascade networks up to the fourth generation of neurons. Although in reality spikes and bursts may propagate through weaker connections, we did not take those pathways into account because the probability of such a transmission in the model was small. Figure 7*B* displays typical examples of spike propagation and burst propagation in the cascade network. A single spike from the parent neuron typically traveled only a few links along a few pathways. By contrast, a spike burst traveled a much longer distance along many more pathways to a significantly greater number of “children” neurons. Burst propagation originating from a parent neuron is stereotyped and repeatable, and may recruit thousands of children neurons within several steps of propagation (Fig. 7*C*). In Figure 7*D*, we summarize the statistics of burst propagation for initial bursts from the parent neuron comprising different numbers of spikes. For example, the transmission length (the number of links a burst propagates) could be as large as 30 if the initial burst contained more than three spikes. These results showed that spike bursts are much more effective than single spikes for routing signals in a lognormal network, and this may have significant implications for the processing of memory traces in CA3.

### Statistical test of sequences

We then tested whether sequences propagating through the present network model are statistically significant by two methods that use imaginary or permuted sequences as control groups (see Materials and Methods). Figure 8, *A* and *B*, exemplifies the distributions of actual and imaginary sequences with length 3 or 7, respectively, following bursts of a reference neuron, which was highly bursty. The distributions of actual and imaginary sequences were statistically different (KS test: *p* < 0.01), implying that some sequences following the bursts of the reference neuron occur statistically more often when the neuron is active than in other instances. Figure 8, *C* and *D*, shows similar cumulative distributions of real and permuted sequences for the same reference neuron as before. The distributions for actual and permuted sequences were statistically different (KS test: *p* < 0.01), meaning that the sequences following the bursting of a reference neuron prefer specific temporal orders to randomly generated orders. Because of practical limitations of computation time and memory, we repeated similar analyses in the five most bursty neurons (including the previous neuron, listed as the fourth most bursty neuron). We found that longer sequences are generally less likely to be statistically significant. Although the critical length varied from neuron to neuron, actual sequences were statistically significant compared with permuted sequences if the length was typically shorter than 5 or 6. By contrast, sequences were statistically significant compared with imaginary sequences only for lengths up to 3 or 4.

We then analyzed whether sequences following bursts of a reference neuron and those following single spikes of the same neuron are statistically different. As shown in Figure 8, *E* and *F*, for the same reference neuron as used in the previous tests, the differences are large, especially for short lengths. The cumulative distribution for single spikes increases linearly with the rank order, implying that all sequences following single spikes occur with almost equal probabilities. In three other neurons, the differences remained to be statistically significant for sequence lengths equal to or shorter than 3. In one reference neuron, single spikes did not produce any sequence. Thus, the occurrence of bursts favors particular sequences in detriment of others, altering an otherwise uniform frequency distribution of sequences generated by single spikes. Together, our results imply that burst sequences propagating through the pathways of strong synapses are statistically significant.

### Comparison with other network models

We have shown that lognormal connectivity is sufficient for producing various properties of spike bursts observed in the hippocampus during highly synchronous events. We examined whether other connectivity structures, typically Gaussian random connections, are also able to replicate these properties. We constructed a recurrent network model of MAT neurons connected randomly with Gaussian-weighted excitatory synapses. The weights of each excitatory synapse were drawn from a Gaussian distribution in which the negative portion was truncated. The mean and SD of the truncated Gaussian were the same as those of the lognormal model. The weights of other connection types were uniform and parameter values for MAT model were unchanged, except ω (see Materials and Methods).

It was previously shown that the generation of sparse (typically <10 Hz) spontaneous activity is difficult in truncated Gaussian recurrent networks of leaky integrate-and-fire neurons (Teramae et al., 2012). Similarly, we could not find out parameter values for which the present truncated Gaussian network generates stable sparse spontaneous activity. Therefore, we added weak background input to randomly chosen 5% of excitatory neurons. Then, this network could generate long-lasting sparse activity (at ∼1–2 Hz) in a narrow area of the parameter space (Fig. 9*A*). For the following results, we set *G _{EI}* = 0.0023 and

*G*= 0.0037, at which the sparse activity continued most robustly. Excitatory neurons displayed a skewed firing rate distribution (Fig. 9

_{II}*B*), as expected from the low firing rates of these neurons (Roxin et al., 2011). By contrast, the distributions of burst event rate and burst index were not lognormal-like (Fig. 9

*C*,

*D*). Indeed, the great majority of excitatory neurons frequently showed single spikes and only occasionally showed consecutive double spikes (Fig. 9

*E*). Here, we define “bursty” neuron as a neuron that achieves all activity sizes up to a burst of size 3, a group that represents ∼2% of the total excitatory population. As evidenced by the biggest burst size neuron (Fig. 9

*F*), analogous examples to exceptionally bursty neurons in the lognormal network (and in experiment) were not found (compare Fig. 2

*C*). In sum, we found that the lognormal network is much more susceptible to bursts of all sizes, whereas the truncated Gaussian network shows a higher average probability for single spikes (Fig. 9

*G*). Finally, we examined burst propagation and found that bursts could propagate only short distances in the truncated Gaussian network regardless of burst size (Fig. 9

*H*; compare Fig. 7

*D*). This result suggests that efficient routing of information by spike bursts is a virtue of the lognormal network.

We also examined another model of Gaussian-weighted network in which only the positive portion of a Gaussian weight distribution having a zero mean was used to generate the weights of recurrent excitatory connections. We performed similar numerical simulations while changing the size of SD of this distribution. We found a narrowly restricted region of the parameter space in which sparse activity (at ∼1 Hz) can last for at least a few tens of seconds. All properties of this network, including difficulties in burst propagation, were similar to those of the previous Gaussian-weighted network (data not shown). Thus, properties and implications of spike bursts are different between the lognormal network and Gaussian-weighted network.

## Discussion

### Network mechanism of lognormal activity patterns

The lognormal distribution of synaptic weights implies that synaptic strength is distributed over a broad range of magnitudes. Theoretical models suggest that such distributions underlie the lognormal distributions of firing rates in cortical and hippocampal neurons *in vivo* (Koulakov et al., 2009; Teramae et al., 2012; Ikegaya et al., 2013) and significantly improve the pattern retrieval performance of associative memory networks in spiking neurons (Hiratani et al., 2013). For other weight distributions without long tails, such as Gaussian distributions, recurrent neural networks contain very few of the low-frequency regimes (typically <10 Hz) of spontaneous firing (Teramae et al., 2012) that are necessary for a lognormal pattern of firing rates (Roxin et al., 2011). As shown recently, lognormal distributions also appear in other statistical features of hippocampal activity, such as the burstiness of individual neurons (Mizuseki and Buzsáki, 2013). Here, we demonstrated that lognormal distributions naturally emerged in firing rates, burstiness, and the magnitude of the population synchrony in a recurrent network of bursty neurons that were connected lognormally. The model predicted exponential distributions of burst length for the majority of excitatory neurons and supraexponential distributions for a minority of neurons. These results were consistent with previous observations in the hippocampus (Harris et al., 2001). Our model also predicts that the temporal fluctuation of population firing rates and population synchrony display lognormal patterns in CA1 and CA3.

A class of neural networks possessing lognormal EPSP amplitude distributions has been shown to generate a lognormal distribution in the mean firing rate (Koulakov et al., 2009; Roxin et al., 2011; Teramae et al., 2012; Ikegaya et al., 2013). Yet, it is not obvious why similar lognormal statistics appear in other statistical indices of neuronal activity because several potential origins exist for these lognormal patterns. For example, a specific intrinsic property of neurons or a multiplicative stochastic process may give rise to a lognormal distribution of burstiness. These possibilities cannot be excluded because the intrinsic properties of a cell crucially influence the statistical properties of burst generation. Our results, however, suggest that a lognormal EPSP amplitude distribution is sufficient for explaining the various lognormal patterns observed in neuronal activity.

### Lognormal neural networks spontaneously activate neural ensembles

To describe the bursting activity in CA3, we used the MAT neuron model, known for its ability to mimic the stochastic spiking responses of *in vivo* cortical neurons (Jolivet et al., 2008; Kobayashi et al., 2009). The resultant network model with lognormally weighted synapses accounted for several known features of single spikes and spike bursts in hippocampal neurons (Harris et al., 2001). Consistent with the results of recent experiments (Hájos et al., 2013; Mizunuma et al., 2014), spike bursts occurred in our model because of a transient deviation from the balanced state of excitation and inhibition. It is widely thought that the transient synchronous activation of ensembles of CA3 pyramidal cells propagates to CA1 to cause SWRs (Buzsáki et al., 1983; Buzsáki, 1989; Csicsvari et al., 2000). Our model showed that similar temporal fluctuations in population firing rates and population synchrony organized spontaneously in recurrent neural networks connected by lognormally weighted synapses. The proportion of neurons activated in highly synchronous events obeyed a lognormal distribution (Fig. 4*G*), and the distribution of the correlation coefficients between connected excitatory neuron pairs was also lognormal (Fig. 4*H*), consistent with experimental observations (Buzsáki and Mizuseki, 2014). Thus, we propose that the long tails of the synaptic weight distribution in the neural network of CA3 consistently account for the sporadic outbreaks of the population activity observed in CA1 during SWR events.

The recurrent architecture among excitatory and inhibitory connections induces neurons to compete with one another during highly synchronous events. The instantaneous firing rates of the individual neurons fluctuate with an increased magnitude of synchronization, with winner (depolarized) and loser (hyperpolarized) neurons emerging in the network, while maintaining lognormal patterns for the overall distribution of firing rates. A similar selection mechanism potentially underlies the selective reactivation of neuron groups during SWR events that engage in the hippocampal representations of episodes (Wilson and McNaughton, 1994; Skaggs and McNaughton, 1996; Nakashiba et al., 2009; Carr et al., 2011). Neurons in CA3 generate more bursts during slow-wave sleep than in other brain states (Mizuseki and Buzsáki, 2013). Our results suggest that spike bursts propagating through the lognormal CA3 network generate an excess amount of sequence information during slow-wave sleep (Fig. 7). Because a single spike burst in Schaffer collaterals can induce long-term potentiation at synapses on CA1 pyramidal cells (Remy and Spruston, 2007), the selective reactivation of bursting CA3 pyramidal cells likely enhances the process of memory formation during SWR events in immobile or sleep states of an animal. Voltage-dependent spike-timing-dependent plasticity may serve to strengthen reciprocal connections between the co-activated neurons (Clopath et al., 2010), or, alternatively, logarithmic spike-timing-dependent plasticity (Gilson and Fukai, 2011) may organize these neurons into mutually connected cell assemblies (Hiratani and Fukai, 2014). Future studies will be necessary to reveal how particular neural ensembles are selected for this process based on past behavioral experiences.

Although the present model described several fundamental features of burst generation in hippocampal neurons during population synchrony, some quantitative aspects of the modeled activity deviated from those obtained during experiments. For instance, the ISI during a burst of 1–2 ms described in the modeled activity was shorter than that for the experimentally observed value of 2–4 ms. This may be due in part to our use of a simple neuron model, which presumably lacked some important details for the intracellular mechanism of burst generation. In the hippocampus, pyramidal cells with lower spontaneous firing rates tend to be activated more frequently during SWRs (Mizuseki and Buzsáki, 2013), implying a negative correlation between the spontaneous firing rate and the proportion of spikes during bursts to all spikes (correlation coefficients of ∼−0.12 to −0.56). In the model, however, the correlation coefficient was ∼0.17. We speculate that inhibitory neurons play a role in producing the counterintuitive negative correlation, although the cause of the above discrepancy remains to be further clarified.

Our model well describes spontaneous activity in CA3. However, the random recurrent connectivity makes it difficult to localize the activity in a restricted part of the network. Therefore, the model does not describe spatially localized task-related activity, such as place fields. The encoding of localized activity may require cell assemblies characterized by highly nonrandom recurrent connectivity (Takahashi et al., 2010), such as modeled previously (Klinshov et al., 2014), or a specific configuration of inhibitory connections.

### Relationship to other models of sharp waves and ripples

Computational models of SWR oscillations in CA1 have been proposed in which transient synchronous input from CA3 drives pyramidal cells and GABAergic interneurons in CA1 to cause high-frequency oscillations with (Traub and Bibig, 2000) and without axo-axonic gap junctions between pyramidal cells (Taxidis et al., 2012). Some models include the influences of different types of interneurons on pyramidal cell responses (Cutsuridis and Hasselmo, 2010). Analysis of the distinct roles of inhibition is important because recent experiments revealed that a transient inhibition of interneurons in CA3 selectively activates the axon initial segments of pyramidal cells during SWR events (Viney et al., 2013). Stimulus-evoked synchronous bursting in CA3 was previously studied in a recurrent network of multicompartmental hippocampal neuron models (Traub and Wong, 1982). Our present results demonstrated that a lognormal recurrent network spontaneously generated and coordinated the synchronous activation of ensembles of bursting CA3 neurons without external input. Sharp waves are also likely to be spontaneous events because they preferentially occur when sensory influences on the hippocampus are weak during sleep or immobile states (Carr et al., 2011). Population synchrony was also modeled previously in terms of frequency-dependent recurrent synapses (Tsodyks et al., 2000). In this model population, synchrony tends to occur periodically due to the inherent rhythm of dynamic synapses. By contrast, the population synchrony generated in a lognormal network obeys the Poisson process.

The most important finding in our model is that a lognormal network efficiently broadcasts spike bursts, but not single spikes, to distant neurons through pathways consisting of strong synaptic connections (Fig. 7). For instance, if a burst from one neuron in a pathway activates on average 3 recipient neurons, the initial burst can propagate to 6561 downstream neurons after 8 steps (3^{8}) and 531,441 downstream neurons after 12 steps (3^{12}) of synaptic transmissions in a recurrent network. By contrast, single spikes cannot travel a long distance. Although this estimation is oversimplified, we suggest that only spike bursts, and not single spikes, are efficient information carriers in the lognormal networks of the hippocampus and neocortex. Our results may account in part for the recent experimental findings indicating that the transmission of isolated spikes in the hippocampus is dispensable for the acquisition of contextual fear memory and that information transfer by hippocampal neurons relies solely on spike bursts during memory encoding (Xu et al., 2012). Thus, overall, our results suggest that lognormal networks are useful for sequence-based information processing in the brain.

## Footnotes

This work was supported in part by the Japanese Society for the Promotion of Science KAKENHI, Ministry of Education, Culture, Sports, Science and Technology Grant-in-Aid for Scientific Research Grants 22115013 and 15H04265 to T.F., and Japanese Society for the Promotion of Science KAKENHI Grant 23220009 to K.I., Ministry of Education, Culture, Sports, Science and Technology Grant-in-Aid for Scientific Research on Innovative Areas “Memory dynamism” 25115002 to K.I., the Mitsubishi Foundation to K.I., the Uehara Memorial Foundation to K.I., and the Takeda Science Foundation to K.I.; M.M.C. was supported by the International Mobility Program for Undergraduate Students of the São Carlos Institute of Physics. We thank Hajime Hirase and Tom McHugh for fruitful discussions about SWR events; Hajime Hirase for valuable comments on the manuscript; and Hideaki Shimazaki for suggestions on the statistical analysis of sequence propagation.

The authors declare no competing financial interests.

- Correspondence should be addressed to either of the following: Dr. Tomoki Fukai, Laboratory for Neural Circuit Theory, RIKEN Brain Science Institute, Hirosawa 2-1, Wako, Saitama 351-0198, Japan, tfukai{at}riken.jp; or Dr. Kaoru Inokuchi, Department of Biochemistry, Faculty of Medicine, Graduate School of Medicine and Pharmaceutical Sciences, University of Toyama, 2630 Sugitani, Toyama 930-0194, Japan, inokuchi{at}med.u-toyama.ac.jp

## References

- Behrens et al., 2005.↵
- Brown and Randall, 2009.↵
- Brunel, 2000.↵
- Buzsáki, 1986.↵
- Buzsáki, 1989.↵
- Buzsáki and Mizuseki, 2014.↵
- Buzsáki et al., 1983.↵
- Buzsáki et al., 2002.↵
- Carr et al., 2011.↵
- Clopath et al., 2010.↵
- Csicsvari et al., 2000.↵
- Cutsuridis and Hasselmo, 2010.↵
- Ego-Stengel and Wilson, 2010.↵
- Ellender et al., 2010.↵
- Eschenko et al., 2008.↵
- Foster and Wilson, 2006.↵
- Gilson and Fukai, 2011.↵
- Gini, 1921.↵
- Girardeau et al., 2009.↵
- Girardeau et al., 2014.↵
- Hájos et al., 2013.↵
- Harris et al., 2001.↵
- Hiratani and Fukai, 2014.↵
- Hiratani et al., 2013.↵
- Ikegaya et al., 2013.↵
- Jolivet et al., 2008.↵
- Klinshov et al., 2014.↵
- Kobayashi et al., 2009.↵
- Koulakov et al., 2009.↵
- Lee and Wilson, 2002.↵
- Lefort et al., 2009.↵
- Lisman, 1997.↵
- Loewenstein et al., 2011.↵
- Mickus et al., 1999.↵
- Mizunuma et al., 2014.↵
- Mizuseki and Buzsáki, 2013.↵
- Nakashiba et al., 2009.↵
- Ranck, 1973.↵
- Remy and Spruston, 2007.↵
- Roxin et al., 2011.↵
- Sarid et al., 2007.↵
- Skaggs and McNaughton, 1996.↵
- Song et al., 2005.↵
- Su et al., 2001.↵
- Takahashi et al., 2010.↵
- Taxidis et al., 2012.↵
- Teramae and Fukai, 2014.↵
- Teramae et al., 2012.↵
- Traub and Bibbig, 2000.↵
- Traub and Wong, 1982.↵
- Tsodyks et al., 2000.↵
- Viney et al., 2013.↵
- Wilson and McNaughton, 1994.↵
- Xu et al., 2012.↵
- Yasumatsu et al., 2008.↵