Abstract
In mammalian auditory systems, the spiking characteristics of each primary afferent (type I auditory-nerve fiber; ANF) are mainly determined by a single ribbon synapse in a single receptor cell (inner hair cell; IHC). ANF spike trains therefore provide a window into the operation of these synapses and cells. It was demonstrated previously (Heil et al., 2007) that the distribution of interspike intervals (ISIs) of cat ANFs during spontaneous activity can be modeled as resulting from refractoriness operating on a non-Poisson stochastic point process of excitation (transmitter release events from the IHC). Here, we investigate nonrenewal properties of these cat-ANF spontaneous spike trains, manifest as negative serial ISI correlations and reduced spike-count variability over short timescales. A previously discussed excitatory process, the constrained failure of events from a homogeneous Poisson point process, can account for these properties, but does not offer a parsimonious explanation for certain trends in the data. We then investigate a three-parameter model of vesicle-pool depletion and replenishment and find that it accounts for all experimental observations, including the ISI distributions, with only the release probability varying between spike trains. The maximum number of units (single vesicles or groups of simultaneously released vesicles) in the readily releasable pool and their replenishment time constant can be assumed to be constant (∼4 and 13.5 ms, respectively). We suggest that the organization of the IHC ribbon synapses not only enables sustained release of neurotransmitter but also imposes temporal regularity on the release process, particularly when operating at high rates.
Introduction
In several sensory systems, the conversion of stimulus representations from graded changes in membrane potential to stochastic trains of spikes is achieved by ribbon synapses (Matthews and Fuchs, 2010). Uniquely, in mammalian auditory systems, each primary afferent (a type I auditory-nerve fiber; ANF) contacts only a single receptor cell (an inner hair cell; IHC) and is excited by a single ribbon synapse. Furthermore, there is evidence that each transmitter release event from a given IHC ribbon synapse gives rise to a postsynaptic spike unless the ANF is refractory (Siegel, 1992; Rutherford et al., 2012). Therefore, the ANF spike patterns provide information about the release statistics of these synapses. These statistics can be studied in the absence of experimenter-controlled sound because ANFs are spontaneously active (Liberman, 1978).
In models of ANF spike patterns, the excitatory (release) events are commonly thought to be produced by a Poisson point process that is homogeneous during spontaneous activity (Kiang et al., 1965; Gaumond et al., 1983; Young and Barta, 1986; Carney, 1993; Li and Young, 1993; Miller and Wang, 1993; Johnson, 1996; Zhang et al., 2001). A fractal doubly stochastic point process has also been proposed to account for rate trends over longer timescales (>100 ms; Lowen and Teich, 1991, 1992; Delgutte, 1996; Jackson and Carney, 2005). Both processes lead to identical exponential distributions of interevent intervals (Teich et al., 1990). The deviations of empirical distributions of interspike intervals (ISIs) from the hypothesized exponential distribution of interevent intervals have therefore traditionally been attributed to refractory properties of the ANFs. However, the long refractory periods (tens of milliseconds) required by this assumption are at variance with the short refractory periods (<1–2 ms) obtained from direct measures (Brown, 1994; Cartee et al., 2000; Miller et al., 2001; Shepherd et al., 2004; Morsnowski et al., 2006). Heil et al. (2007) suggested a solution to this problem in the form of a non-Poisson excitatory point process and found that, when followed by refractoriness consistent with the experimental data, it describes the ISI distributions from spontaneous activity of cat ANFs nearly two orders of magnitude better than the Poisson model.
ISIs of ANFs during spontaneous activity can also exhibit nonrenewal properties that manifest, for example, as serial ISI correlations or a spike-count variability different from that of a memoryless renewal process (Lowen and Teich, 1992). These aspects have received little attention, and their origin has remained unclear. Here, we examine whether the previously discussed non-Poisson excitatory process can account for these nonrenewal properties. We find it does not provide a simple explanation for the dependence of nonrenewal measures on mean ISI, so we investigate an alternative model of excitation based on vesicle-pool depletion and replenishment. This alternative model is physiologically plausible and accounts for both the nonrenewal properties and the ISI distributions, with the variation of only a single parameter between spike trains. Our findings will help clarify the processes of vesicle release at ribbon synapses and information transfer in the auditory nerve.
Materials and Methods
Data
The data for this study include those data used in a previous study of the ISI distributions of ANFs during spontaneous activity (Heil et al., 2007). All details of data acquisition can be found there. Briefly, in five barbiturate-anesthetized adult cats (three females, two males), spikes of 171 individual ANFs were recorded extracellularly with microelectrodes from the left auditory nerve near its exit from the internal auditory meatus. Continuous samples of spontaneous activity (between 12.5 and 134.4 s long) were recorded, along with responses to various stimulus protocols for other purposes (Heil et al., 2008, 2011). Spike times were taken as the instants at which the amplified and filtered electrode signal crossed a Schmitt trigger level, and were stored on disk with a precision of 1 μs.
Data analysis
All data analysis was performed using MATLAB R2014a (The MathWorks) or C++ code compiled with Microsoft Windows SDK 7.1. The data presented here are always from complete samples of spontaneous activity, unlike the previous study (Heil et al., 2007) in which pruning was applied to remove initial periods of nonstationary spike rate.
Serial ISI correlation coefficient.
First-order ISIs (i.e., those between consecutive spikes) were computed from the spike times with 1 μs precision. In a renewal point process, the successive time intervals between adjacent events are drawn independently from a common probability distribution (Cox, 1962; Teich and Khanna, 1985). Consequently, there is no serial correlation between first-order or higher-order interevent intervals and the process is said to be memoryless. A homogeneous Poisson point process is a special case of a renewal process where the interevent intervals are exponentially distributed (Cox, 1962; Lowen and Teich, 1992; Ratnam and Nelson, 2000). To probe the spontaneous spike trains of ANFs for nonrenewal properties, the serial interspike interval correlation coefficient (SIICC) was computed for each sample of spontaneous activity using the equation proposed by Lowen and Teich (1992):
Here, ρ is the SIICC and N is the number of first-order ISIs. The mean ISI 〈τ〉 is computed from all N individual intervals τi. The equation is not defined for fewer than three intervals. For comparative purposes, we also computed SIICCs using similar equations provided by Ratnam and Nelson (2000), Chacron et al. (2001), Engel et al. (2008), Schwalger and Lindner (2010), and Urdapilleta (2011), but because the results were similar, we report only those obtained with Equation 1. Because the first-order SIICCs were of low magnitude (see Results), we refrained from analyzing higher-order interval correlations, shown by others to be even smaller (Ratnam and Nelson, 2000; Chacron et al., 2001; Avila-Akerberg and Chacron, 2011).
Conditional means of ISIs.
For comparative purposes, nonrenewal properties of ANF spike trains during spontaneous activity were also probed using conditional means of ISIs, similar to analyses performed by Tsuchitani and Johnson (1985) and Johnson et al. (1986). For each spike train, a criterion ISI was selected. The conditional mean of all ISIs directly following an ISI shorter than or equal to the criterion was then calculated. The conditional mean of the remaining ISIs, those directly following an ISI longer than the criterion, was also calculated. This was done for a range of criteria, from 1 to 25 ms in steps of 0.5 ms. For a renewal process, the two conditional means will not differ systematically and, on average, their ratio will be 1. In the presence of negative serial ISI correlations, however, the conditional mean of the ISIs following an ISI longer than the criterion will be shorter than the conditional mean of the ISIs following an ISI shorter than the criterion. On average, the ratio of the two conditional means will therefore be <1. This latter behavior was observed in the data, with 87% of the ratios <1 (from 180 spike trains with ≥500 ISIs). From these data, we derived a single ratio for each spike train, the geometric mean of the individual ratios across criteria. We found that the logarithm of this geometric mean ratio was closely and positively correlated with the SIICC (r2 = 0.83; n = 180). Therefore, we report only the results based on the SIICC analysis.
Spike-count distributions and Fano factors.
Additional information about spike-train variability was derived from spike-count distributions, also referred to as pulse-number distributions (PNDs; Teich and Khanna, 1985). Each spike train was first divided into nonoverlapping time windows of duration T. Then, the number of spikes in each window was counted. Division by the total number of spikes yields the probability P(n, T) of observing n spikes during the counting time T. A histogram of these probabilities yields the PND. The PNDs can be characterized by the Fano factor F(T) (Fano, 1947), which is defined for each counting time T as the ratio of the variance to the mean number of spikes per window:
The Fano factor has units of spikes. For a Poisson process, F(T) = 1 for all T (Cox and Lewis, 1966; Chacron et al., 2001; Amarasingham et al., 2006). Processes with F(T) <1 are less variable than a Poisson process (Gabbiani and Koch, 1998). For each of the 446 spontaneous spike trains with at least one spike, we computed Fano factors for a range of counting times in octave steps from ∼0.1 to 250 ms. Specifically, the values were T = 500/(2n) ms, where n represents the integers from 1 to 12. Nine spike trains were omitted because their total durations could not be divided into integer numbers of these counting times.
Simulations and models
Homogeneous and inhomogeneous Poisson point processes.
For various comparative purposes described below, events from a homogeneous Poisson point process were simulated by drawing random intervals from an exponential distribution with a specified mean parameter (using the MATLAB function exprnd). To simulate nonstationary (rate-varying) processes, event trains were generated for discrete time series with time steps of 0.01 ms or 0.05 ms. A nonhomogeneous Poisson point process was approximated by drawing a random number from the standard uniform distribution at each time step and comparing it to the instantaneous event probability per unit of time, which is the expected number of events per time step (the product of the rate in events per seconds and the duration of the time step in seconds). An event was said to occur at each time point for which the randomly drawn number was less than or equal to the instantaneous event probability.
Simulation of refractory effects.
Effects of refractoriness were also modeled. Each simulated event triggers a simulated spike unless the event falls into the refractory period that follows the previous spike. The refractory period, or dead time, consists of a portion with a fixed duration tD and a portion whose duration varies and was drawn after each spike from an exponential distribution with mean tR. The cumulative distribution function of dead times, p(t), is therefore given by the product of the Heaviside function H(t) and an exponential function, both shifted by tD:
where t is the time since the last spike and tD ≥ 0. This implementation of refractoriness is identical to that proposed by Young and Barta (1986) and to that used in the previous study of ISI distributions (Heil et al., 2007). Mathematically, it is described by the convolution of the distribution of interevent intervals with the refractory function of Equation 3 (Young and Barta, 1986; Heil et al., 2007; Neubauer et al., 2009).
In this implementation of refractoriness, the duration of the refractory period is a random variable. Notably, the outcome of this implementation differs from that of implementations of refractoriness in popular models of the auditory periphery (Carney, 1993; Zhang et al., 2001; Sumner et al., 2002, 2003; Zilany and Bruce, 2006; Zilany et al., 2009). These follow the (possibly more conventional) view that there is a nonrandom stereotypical refractory function composed of an absolute refractory period during which the probability that an event will trigger a spike is zero, followed by a relative refractory period during which the probability that an event will trigger a spike increases exponentially with time. It should be noted that unlike the model of Sumner et al. (2002), the model of Zilany et al. (2009) does not explicitly model synaptic release events and subject them to refractory effects to yield spikes. Instead, the refractory effects are incorporated directly into the calculation of the (nonhomogeneous) rate for the Poisson spike generator; in either case, the result is the same. Although the probability distribution of refractory periods given by Equation 3 has the same form as the conventional refractory function composed of an absolute refractory period with duration tD and a relative refractory period characterized by an exponential recovery with time constant tR ≥ 0, the two implementations are not equivalent. This difference can be seen by comparing the effects of each implementation on simulated Poisson event trains. The model by Sumner et al. (2003) assumes that in the relative refractory period the probability that an event will trigger a spike increases exponentially from a probability of c to 1 (where 0 ≤ c ≤ 1). The model of Zilany et al. (2009) is similar but uses a double-exponential relative refractory function. When the second time constant is set to zero and the remaining parameters are identical, the distributions of ISIs produced by the models of Sumner et al. (2003) and Zilany et al. (2009) are identical. However, with c = 0 to match Equation 3, the mean ISI produced by each of these models is systematically shorter than that produced by our approach and that of Young and Barta (1986), Heil et al. (2007), and Neubauer et al. (2009), which predicts a mean ISI of tD + tR + tE, where 1/tE is the rate of the Poisson process. The absolute and relative differences between implementations increase with increasing duration of the time constant tR of the relative refractory period.
The reason for this difference becomes evident from considerations of probability. Following a simulated spike, there may be a series of simulated release events that fail to trigger spikes due to refractoriness, followed by one event that finally succeeds in triggering a spike. In our implementation of refractoriness, the probability that the duration of the dead time lies between t = 0 and t = T is given by ℙ(t ≤ T) = p(T). Thus, the probability that the dead time exceeds t = T is given by ℙ(t > T) = 1 − p(T). All simulated events falling into this dead time fail to trigger spikes. In the implementation of Sumner et al. (2003) and Zilany et al. (2009), however, the probability that all events falling into the interval T fail to trigger spikes is given by the product of the individual probabilities that these events at intervals t1, t2,…tI fail to trigger spikes: ℙ(t > T) = (1 − p(t1))(1 − p(t2))…(1 − p(tI)). Because this probability is lower than ℙ(t > T) = 1 − p(T) when more than one event falls into this dead time, the number of events triggering a spike will be higher and the mean ISI will be shorter. In all simulations, both the absolute refractory period and the time constant of the relative refractory period were set to 0.6 ms, based on previous modeling of the ISI distributions of the same cat-ANF spontaneous spike trains (Heil et al., 2007).
The constrained-failure model.
Gamma-mixture distributions of the first-order intervals between simulated events (IEIs) were generated from simulated homogeneous Poisson event trains by removing specified fractions of events, but without ever removing more than one in a row. This “never-two-in-a-row rule” is required to generate gamma-mixture IEI distributions from Poisson event trains (Heil et al., 2007; Neubauer et al., 2009). Following this rule, the effects of different removal scenarios on serial correlation coefficients and Fano factors were explored. These removal scenarios are described in Results.
The depletion–replenishment model.
The model of vesicle-pool depletion and replenishment used here is equivalent to that proposed by Goldman et al. (2002) for central synapses, and contains three parameters. The model describes a pool of readily releasable units whose maximum size is specified by parameter nmax. Each unit has the same constant and independent probability per unit of time of being depleted, specified by parameter pdepl. Finally, each depleted unit has the same constant and independent probability per unit of time, prepl, of being replenished. On average, this yields an exponential replenishment with a time constant of τrepl = 1/prepl (Goldman et al., 2002). Figure 1A illustrates the model as it would apply to the auditory periphery, in the form of a synaptic ribbon with tethered vesicles, some of which are readily releasable. Although a single releasable unit is represented here as a single vesicle, it could also be a group of vesicles that are released simultaneously in a coordinated fashion (multivesicular release; Glowatzki and Fuchs, 2002; Grant et al., 2010). Following Goldman et al. (2002), depletion is constrained such that maximally one unit can be released during a given time step. The instantaneous probability that a synaptic release event occurs, pevent(ti), depends only on pdepl and the instantaneous number of available units, n(ti) (where 0 ≤ n(ti) ≤ nmax):
where ti is the current time step. When a unit is depleted, n is reduced by 1. If each depleted unit was replenished instantaneously, the pool of available units would always remain full (i.e., n would always equal nmax) and the release process would be a homogeneous Poisson point process with a mean interval of approximately 1/(pdepl · nmax) or a rate of approximately pdepl · nmax. This is so because Pevent(ti) can be approximated by pdepl · n(ti) for small values of pdepl and n(ti), as were used in our simulations. However, replenishment is not instantaneous in the model, so that the time-varying number of readily releasable units in the pool is given by the following:
where ti is the current time step, ti−1 is the prior time step, and ti − trel is the time since the last release event. Each release event in the depletion–replenishment model triggers a spike unless the release event falls into the refractory period following the previous spike, as described above. Negative SIICCs occur due to the fluctuation of the pool size, as captured by Equation 5, and its effect on release probability, as captured by Equation 4; for example, when two release events occur in close succession to form an interval shorter than the mean, the pool of readily releasable units will tend to need a long replenishment period and the next interval will therefore tend to be longer than the mean (Fig. 1B).
The depletion–replenishment model. A, Schematic ribbon synapse with tethered vesicles, some of which (green and blue) are ready for release. They constitute the pool of releasable units of maximal size nmax. In this scheme, a releasable unit is a single vesicle, but in the model, such a unit could also represent a group of vesicles released simultaneously. Each unit has a constant and independent probability per unit of time, pdepl, of being released. When a unit is released (blue vesicle), the pool size is depleted by 1 and replenished exponentially with time constant τrepl. If replenishment was instantaneous, the release process would be a homogeneous Poisson point process with a mean interval of approximately 1/(pdepl · nmax). B, A 50 ms simulation demonstrating the change in size of the readily releasable pool for a single combination of pdepl, τrepl, and nmax (0.05/ms, 15 ms, and 4, respectively). C, Mean SIICCs (color coded) resulting from the combinations of pdepl, τrepl, and nmax that produce a mean interval of 10 ms after refractoriness, as a function of τrepl and nmax. D, Mean SIICC (color coded) from a denser sampling of the relevant parameter space. A tenth-order polynomial was fitted to the pairs of τrepl and nmax capable of producing a mean ISI of 10 ms and a mean SIICC of −0.1 (cf. Fig. 2D; see text for rationale).
Procedure to determine pdepl in the depletion–replenishment model.
When simulating cat-ANF spike trains with the depletion–replenishment model, it is important that the model spike train has the same or a very similar mean ISI as the cat-ANF spike train to be simulated. To explore how to achieve this, event trains were simulated (with a time step of 0.05 ms and between 20 and 100 repetitions) for various combinations of the three model parameters (pdepl, τrepl, and nmax) and then subjected to refractoriness to yield simulated spike trains from which mean ISIs were calculated. For a given pair of τrepl and nmax, we found the mean ISI to increase approximately linearly with increasing 1/pdepl. The value of pdepl required to yield a particular mean ISI, given τrepl and nmax, was therefore calculated from a linear fit.
Procedure to obtain preliminary parameters for the depletion–replenishment model.
A multistep simulation procedure was used to obtain a set of preliminary parameter combinations that might account for the serial correlations observed in the cat data. This procedure was based on the simplifying assumption that only the depletion probability pdepl varies between synapses to account for differences in spontaneous rates between spike trains; therefore, we specifically sought to identify a set of τrepl and nmax pairs that might account for the data, with pdepl free to vary as needed to yield the appropriate mean ISI to match each sample of cat data. We first targeted parameter combinations that would produce a mean SIICC of −0.1 and a mean ISI of 10 ms (i.e., a mean rate of 100 spikes/s). This single point was chosen because it closely matches the mean SIICC in cat-ANF spike trains with mean ISIs of ∼10 ms and because of the high density of data points near it (cf. Fig. 2D); pairs of τrepl and nmax that could not reproduce this point were considered unable to account for the dataset, whereas those that could reproduce this point were included in the set of preliminary pairs to be explored further. The first step in the procedure was to determine the values of pdepl necessary to obtain a mean ISI of 10 ms after refractoriness. This was done for 10,000 different pairs of τrepl and nmax (100 values of τrepl, log-spaced between 1 and 50 ms, times 100 values of nmax, log-spaced between 1 and 100, forming the grid in Fig. 1C) using the linear-fit procedure described above. For certain pairs of long τrepl and small nmax, such a short mean ISI cannot be obtained (Fig. 1C, white area), and these pairs were not considered further. Next, the mean SIICCs (color coded in Fig. 1C) were calculated from 12.5 s long simulated spike trains for each of the remaining pairs (with a time step of 0.05 ms and 200 repetitions). The narrow yellow band marks the combinations of τrepl and nmax which, together with an appropriate pdepl, give rise to SIICCs near −0.1. For longer τrepl and smaller nmax, the resulting SIICCs are more negative and vice versa. Next, and as a first approximation, a power law was fitted to the pairs of τrepl and nmax that yielded SIICCs within the range of −0.1 ± 0.004; to improve the fit, short values of τrepl < 5 ms were excluded to omit the vertical portion below the bend in the yellow band in Figure 1C. The fitted power law was used to define a region of parameter space that was then resampled at a higher density (Fig. 1D). Specifically, the power law was used to calculate values of ncenter for 100 log-spaced values of τrepl between 5 and 50 ms. For each ncenter, 50 log-spaced values of nmax were specified between ncenter · 0.85 and ncenter/0.85. The steps described above were then repeated with this new set of τrepl and nmax pairs to determine the appropriate value of pdepl for each and the resulting mean SIICCs (color coded in Fig. 1D). This time, a more precise tenth-order polynomial was fitted to the pairs of τrepl and nmax that yielded SIICCs within the range of −0.1 ± 0.001. For 150 log-spaced values of τrepl between 5 and 50 ms, the corresponding values of nmax were calculated from this polynomial fit to yield, together with the appropriate pdepl, 150 preliminary parameter combinations, each of which produces a mean ISI of 10 ms and a mean SIICC of −0.1. Finally, for each of the 150 pairs of τrepl and nmax, 180 values of pdepl were determined (via the linear-fit procedure described above) to approximate the mean ISIs of all 180 ANF spike trains with ≥500 ISIs used to evaluate this model (see Results). These parameters were then used to simulate a complete set of 180 spike trains (with a time step of 0.01 ms), such that each cat-ANF spike train had a simulated counterpart of the same duration and with approximately the same number of ISIs. Twenty such sets of 180 simulations were made for each of the 150 pairs of τrepl and nmax and were compared with the cat data to identify which pairs of τrepl and nmax best reproduce the observed SIICCs and ISI distributions (described in Results).
Results
The results are organized into three parts. First, the nonrenewal properties of ANF spike trains are characterized using SIICCs and Fano factors. Second, the constrained-failure model of excitation is evaluated to determine whether it can adequately reproduce the nonrenewal properties. Finally, we show that the alternative depletion–replenishment model of synaptic vesicle pools can better account for the nonrenewal properties, while also accounting for the ISI distributions.
Characterization of the nonrenewal properties of ANF spike trains
Serial interspike intervals of ANF spike trains are predominantly negatively correlated
In this section, we demonstrate the existence of predominantly negative serial correlations in the spike trains of ANFs during spontaneous activity. For each recorded spike train with ≥4 spikes (≥3 ISIs), we computed the SIICC using the equation suggested by Lowen and Teich (1992) (Eq.1 in Materials and Methods). A negative value (anticorrelation) means that short ISIs tend to be followed by long ISIs and vice versa, whereas a positive value means that short ISIs tend to be followed by short ISIs and long ISIs by long ISIs. According to Lowen and Teich (1992), the SIICC ranges from −1 (perfect anticorrelation) to 1 (perfect correlation) and is zero for uncorrelated intervals.
Figure 2A shows the SIICCs obtained from 414 samples of spontaneous activity (each between 12.5 and 134.4 s long) of 156 cat ANFs as a function of the number of first-order ISIs in the sample. Most SIICCs (70%) are negative. For the smallest sample size of three ISIs, all SIICCs are ≤0 and some are even less than −1. With increasing sample size, the spread of the SIICCs decreases, as well as the magnitude of the mean SIICC (open circles). Still, for large sample sizes, most SIICCs are negative.
Spike trains of cat ANFs are nonrenewal and exhibit serial interspike interval correlations. A, SIICCs obtained from 414 samples of spontaneous activity as a function of the number of first-order ISIs in the sample (gray dots). The means, computed for 18 bins containing 23 entries each, are also shown (open circles). B, Plot of individual SIICCs (gray dots) and their means (open circles) obtained from simulated spike trains as a function of sample size. Spike trains were generated by a homogeneous Poisson process, so there is no systematic serial correlation of ISIs. The SIICC measure (Eq.1; adopted from Lowen and Teich, 1992) exhibits a bias toward negative values, particularly for small sample sizes. C, Cumulative probability of SIICCs obtained from 180 samples of ANF spontaneous activity with ≥500 ISIs (continuous line) and after duplicating and shuffling each train of ISIs 10 times (dashed line). Shuffling removes the systematic serial correlations. D, SIICCs obtained from the 180 cat-ANF spontaneous spike trains with ≥500 ISIs plotted as a function of mean ISI (two samples with mean ISIs longer than 50 ms are not shown).
We performed simulations to determine which SIICCs can be expected by chance from a renewal point process and to identify a possible bias in the SIICC measure. Spikes were generated from a homogeneous Poisson point process (see Materials and Methods) so that the intervals between them (ISIs) exhibit no systematic serial correlations. For each sample size (between 3 and 8000), ∼250,000 ISIs were partitioned into samples from which individual SIICCs were calculated. Figure 2B shows the individual SIICCs (gray dots) along with the means (open circles) as a function of sample size. Despite the absence of systematic serial correlations in these ISIs, the mean SIICC measure is strongly negative for small sample sizes and therefore has a systematic bias. With increasing sample size, the spread decreases and the mean asymptotically approaches zero. For large sample sizes, the small spread appears symmetrical around zero. Because of the bias toward negative values for small sample sizes, we restrict our analysis of SIICCs to the 180 samples with ≥500 ISIs (each between 12.5 and 52.5 s long). Here, the spread is small and the means from the simulated data are close to zero (Fig. 2B). This restriction, however, results in the exclusion of samples with low spontaneous rates (SRs), because they were not long enough to reach the criterion of at least 500 ISIs; the lowest rate present in our SIICC analysis is 13.7 spikes/s.
Figure 2C shows the cumulative probability of SIICCs obtained from the samples of ANF spontaneous activity with ≥500 ISIs (solid line; n = 180). The distribution has a median of −0.060 and an interquartile range from −0.083 to −0.028. The median differs significantly from zero (Wilcoxon signed rank test: z = −9.88, p = 4.97 · 10−23). Figure 2C also shows the cumulative probability of SIICCs obtained from the same spike trains after duplicating and shuffling the ISIs of each train 10 times (dashed line, n = 1800), which removes the serial correlations. The distribution of the SIICCs after shuffling has a median of −0.0015 and an interquartile range from −0.022 to 0.018. The median does not differ significantly from zero in any of the 10 shuffled sets (Wilcoxon signed rank test: −1.65 ≤ z ≤ 1.22, 0.10 ≤ p ≤ 0.98).
The spread of SIICCs in the original (unshuffled) spontaneous spike trains has a strong dependence on the mean ISI. Figure 2D shows the SIICCs from the 180 ANF spike trains with ≥500 ISIs as a function of the mean ISI. The SIICCs from two samples with mean ISIs >50 ms (SRs <20 Hz) are not shown to allow for better resolution. The most negative SIICCs tend to come from spike trains with short mean ISIs (high SRs). As the mean ISI increases (SR decreases), the SIICCs increase systematically and approach zero. The few positive SIICCs here are due mainly to nonstationary spike rates, as we briefly demonstrate in the next section.
Positive correlations in serial interspike intervals of ANFs are likely due to nonstationary spike rates
One factor contributing to the spread of SIICCs in the original (unshuffled) spike trains may be nonstationarities in the spike rate. Scrutiny of the samples with the most positive SIICCs revealed obvious (monotonic or non-monotonic) changes in spike rate over time, whereas this was not the case for samples with negative SIICCs. To examine the issue systematically, we investigated how SIICCs were affected when the ISIs of each spike train were shuffled within contiguous, nonoverlapping time bins. For a given spike train, the bins were chosen such that each contained the same number of ISIs and thus the same proportion of the total number of ISIs in the spike train. If positive serial correlations are caused by long-term nonstationarities in rate (as opposed to, for example, localized bursting behavior), then shuffling the ISIs locally within such bins would remove only the negative serial correlations, but leave the positive correlations intact. For each of the 180 ANF spike trains with ≥500 ISIs, we first removed the last few spikes (≤9 spikes), such that the remaining number of ISIs was a multiple of 10; this was done to increase, on average, the number of factors that divide the total number of ISIs into groups of equal integer size. ISIs within each group were then shuffled. The mean SIICCs calculated from 500 repetitions of this procedure are shown in Figure 3 as a function of the proportion of ISIs in each group relative to the total number of ISIs in the spike train. Each black line represents one spike train in which the original (unshuffled) intervals yielded an SIICC >0.025, as indicated by the leftmost dot. The rightmost dot indicates the case in which all ISIs were shuffled within a single group, and the dots in between represent the intermediate group sizes. The samples in which the original intervals yielded an SIICC <0.025 are shown as gray dots. Samples with SIICCs <0.025 tend to converge to zero even when the proportion of ISIs in each group is small (<∼0.02). In contrast, samples with SIICCs >0.025 tend not to converge to zero until the proportion of ISIs in each group is larger (>∼0.1). This indicates that nonstationary spike rates are responsible for the few clearly positive SIICCs that exist. However, due to the relatively small number of such conditions, we neither excluded them from further analyses nor attempted to incorporate nonstationarity into our modeling.
Positive SIICCs are due to nonstationary spike rates. The ISIs of each spike train were divided into nonoverlapping groups of equal size and then shuffled within each group. This was repeated for each spike train using all group sizes that fit into it an integer number of times. SIICCs were calculated for 500 shuffles and averaged to yield mean values, which are plotted as a function of the proportion of ISIs in each group. Each line connecting black dots represents one spike train in which the original intervals (leftmost black dots) yielded an SIICC >0.025, whereas gray dots represent the spike trains in which the original intervals (leftmost gray dots) yielded an SIICC <0.025. Samples with SIICCs <0.025 converge to zero after relatively localized shuffling, but samples with SIICCs >0.025 are only slightly affected until the proportion of ISIs in each group becomes much larger. The fact that localized shuffling removes negative but not positive serial correlations indicates that the latter are due to changes in spike rate over longer timescales.
Spike-count variability
In the following sections, we examine the spike-count variability of ANF spontaneous spike trains using the Fano factor. These analyses are not restricted to spike trains with ≥500 ISIs, but instead include trains with at least one spike (n = 446). Thus, unlike the SIICC analysis, the spike trains here have a broad range of spontaneous rates. For each of these spontaneous spike trains, the Fano factor F (the ratio of the variance to the mean number of spikes per window) was computed for a range of counting times T (see Materials and Methods).
Large Fano factors at long counting times are likely due to nonstationary spike rates
For short counting times, F(T) was nearly always <1, but for counting times in excess of a few hundred milliseconds, F(T) increased in most samples in a power-law fashion with increasing T. A representative example is shown in Figure 4A. Such behavior of F(T) has been described previously for ANFs (Lowen and Teich, 1992) and for electrosensory afferents (Ratnam and Nelson, 2000). Chacron et al. (2001) and Avila-Akerberg and Chacron (2011) have shown that the power-law-like increase of F(T) at long counting times is likely a consequence of nonstationarities in the data, because even subtle changes in the mean spike rate on long timescales can lead to weak positive serial correlations that extend out to higher-order intervals (long lags). This was also the case in our data. Figure 4B shows, for the same 50 s sample as in Figure 4A, the cumulative sum of the number of ISIs as a function of time since the first spike. A linear fit of these data through the origin is also shown (dashed line), as well as the differences between the fit and the data (residuals; gray curve, right ordinate). The residuals are pronounced and systematic in that they are positive for the initial 35 s and negative for the remaining time, reflecting a gradual increase in the spike rate over time. The RMS value of the residuals (∼49 ISIs for this example) served as a quantitative measure of nonstationarity.
Large Fano factors at long counting times are due to nonstationary spike rates. A, Fano factor F(T) plotted as a function of counting time T for a representative 50 s sample of spontaneous activity (A7-U5-R2, 49.74 spikes/s). As T increases, F(T) gradually decreases until T = 125 ms, after which F(T) increases in a power-law-like fashion. This rapid increase disappears after shuffling the ISIs (dashed line represents the geometric mean of 1000 shuffles). B, Cumulative sum of the number of ISIs versus time since the first spike (black continuous line) for the same sample used in A, and a linear fit of these data through the origin (dashed line). The differences between the linear fit and the data (residuals; gray line and right ordinate) emphasize the gradually increasing spike rate in this sample. C, Fano factor for T = 5 s, F(T = 5 s), versus the RMS of these residuals, for all samples with ≥100 spikes and a duration of at least 35 s (black circles). Note the positive relationship. The measure is also shown after 1000 shuffles of the ISIs of each sample (open symbols represent the geometric mean of F(T = 5 s) and the arithmetic mean of the RMS residuals).
Figure 4C shows the Fano factor obtained for a counting time of 5 s, F(T = 5 s), versus the RMS value of the residuals (filled symbols). The analysis was restricted to samples with at least 100 spikes and with durations of at least 35 s (n = 56), so that the computation of F(T = 5 s) was always based on at least seven counting windows. Figure 4C shows a correlation between the RMS values and F(T = 5 s). Shuffling the ISIs from these spike trains abolishes the nonstationarities (and the ISI correlations) and results in lower RMS and F(T = 5 s) values (Fig. 4C, open symbols, each averaged from 1000 shuffles). Because we are not primarily interested in the effects of slow changes in spike rate on Fano factors, we restrict further analyses to counting times ≤250 ms.
Fano factors at short counting times are determined by refractoriness and spontaneous rate
Figure 5A shows, for all 446 samples with at least one spike, F(T) as a function of counting times ≤250 ms. For counting times shorter than ∼1 ms, all F(T) are <1. For longer counting times, the majority (91%) of F(T) are <1. The spike trains therefore tend to be more regular than a homogeneous Poisson point process, for which F(T) = 1. The spread of F(T) increases with increasing T. To scrutinize the spread at a higher resolution, the 321 samples for which F(T) < 1 for all T ≤ 250 ms are shown in Figure 5B in the form 1 − F(T) on a logarithmic axis. This measure can be thought of as the excess regularity of the spike train compared with that of a homogeneous Poisson point process, for which 1 − F(T) = 0. Figure 5B shows that, in all samples, 1 − F(T) grows linearly with T (with a slope of 1 in the log-log plot) up to counting times of ∼1 ms. In some samples, this linear growth continues out to the longest value of T plotted, whereas in others, the growth becomes compressive at longer values of T before saturating. Samples with low values of 1 − F(T) at short counting times tend to have linear growth out to much longer counting times than samples with higher values of 1 − F(T) at short counting times. The spread of 1 − F(T) at a given value of T, i.e., the differences in the proportionality coefficient between 1 − F(T) and T, are related to differences in SR. This is illustrated in Figure 5C, which shows, in double logarithmic coordinates, 1 − F(T) as a function of SR, for selected values of T. Again, for short counting times (represented by T = 125/512 ms in Fig. 5C), 1 − F(T) grows linearly with SR, whereas for longer counting times, values of 1 − F(T) from high-SR fibers fall below the trend extrapolated from low-SR fibers.
Fano factors for short counting times are determined by refractoriness and spontaneous rate. A, Fano factors F(T) for all ANF spike trains with at least one spike (n = 446) computed for values of T = 500/(2n) ms, where n represents the integers from 1 to 12. Data from a given sample are connected. B, Plot of 1 − F(T) as a function of T for those 321 samples for which F(T) <1 at all T ≤ 250 ms. Note the linear increase for short T. C, Plot of 1 − F(T) as a function of SR for all samples and for four selected counting times. Note the linear increase of 1 − F(T) with SR for short T and for low SR for all T. D, Plot of 1 − F(T) as a function of the mean number of spikes per counting window. Note the equality up to ∼0.04 spikes/T.
The linear growth of 1 − F(T) with T and with SR can be understood from the effects of refractoriness (Rajdl and Lansky, 2014). The definition of the Fano factor (Eq. 2) can be reformulated to the following:
where ni is the number of spikes in counting window i and NT the total number of counting windows for a given T. When T is shorter than the shortest ISI in the spike train, each counting window can contain either no spike or one spike. Under this condition, it follows that ∑i=1Nni2 = ∑i=1Nni, so that F(T) = 1 − 〈ni〉, where 〈ni〉 = ∑i=1Nni/NT is the mean number of spikes per T. This number increases linearly with T for a given SR and linearly with SR for a given T. Hence, as long as each counting window contains either no spike or one spike, 1 − F(T) = 〈ni〉. In our data, we find this to be true for 〈ni〉 <∼0.04 spikes per counting window (Fig. 5D). Because the shortest ISI in a spike train is bounded by the absolute refractory period, that period and the spike rate uniquely determine the Fano factor for short values of T.
Measures of nonrenewal characteristics based on serial ISI correlations and spike-count variability are related
Next, we explore the relationship between two measures of the nonrenewal characteristics of ANF spike trains: the (predominantly negative) serial correlations of ISIs and a measure derived from the Fano factor analysis. The latter measure is the difference in the values of F(T) derived from the original spike trains and those from the same spike trains after the ISIs were shuffled. For a renewal process, both measures of nonrenewal should equal zero on average. A negative Fano factor difference means that the original spike train is more regular than the spike train with shuffled ISIs. It also represents an increase in the variance of the number of spikes per counting time after shuffling the ISIs, because the mean number of spikes per counting time is not affected by the shuffling. In other words, the PNDs from the original spike trains are narrower than those from the trains with shuffled ISIs. Figure 6 shows scatterplots of the Fano factor difference and the SIICC for all 180 samples with ≥500 ISIs and for selected counting times. As expected, the slope of the relationship between Fano factor difference and SIICC increases with the counting time T. The most negative Fano factor differences tend to be obtained from the spike trains also characterized by the most negative SIICCs, those with high SRs (Fig. 2D).
Serial ISI correlations and spike-count characteristics used to measure nonrenewal properties of ANF spike trains are related. The panels show, for each of 180 samples with ≥500 ISIs and for selected counting times, the difference in F(T) between the original spike train and that with shuffled ISIs as a function of the SIICC from the original spike train. Note that the two measures of nonrenewal characteristics are positively correlated and that the slope of the relationship increases with increasing counting time T.
Evaluation of the constrained-failure model
The ISI distributions of cat-ANF spontaneous spike trains have previously been described as resulting from refractoriness operating on trains of excitatory events with IEIs whose distribution is a mixture of a gamma distribution with shape factor 2 and an exponential distribution (a gamma distribution with shape factor 1), with a common rate factor (Heil et al., 2007). In the following sections, we explore whether processes that could generate such a “gamma-mixture” distribution provide a parsimonious explanation for the nonrenewal properties of ANF spike trains described above.
A gamma-mixture distribution of IEIs can be generated from a homogeneous Poisson point process generating “primary” events (Fig. 7A, filled circles in top row). The primary events can give rise to “secondary” events (Fig. 7A, squares, asterisks, and open circles) that, for the present purpose, can be conceived of as excitatory (release) events (each of which triggers a spike in the ANF unless the ANF is refractory). However, some of the primary events may fail to become secondary events. A single failure of a primary event results in an interval between secondary events, a secondary IEI, whose duration corresponds to the sum of two successive intervals from the exponential distribution of primary IEIs. The resulting secondary IEI will therefore be an interval from a gamma-2 distribution, with the same rate as the original exponential distribution. As the proportion of failures increases, the proportion of gamma-2 intervals between the secondary events in the gamma-mixture distribution increases and the proportion of gamma-1 intervals decreases. If every second primary event failed, the resulting distribution of secondary IEIs would be a pure gamma-2 distribution. To generate a gamma-mixture distribution of secondary IEIs, failures must not occur entirely at random, otherwise the resulting distribution of secondary IEIs would still be exponential, but with a lower rate. Specifically, the model has the constraint that two or more failures must not occur in direct succession (never-two-in-a-row rule; Heil et al., 2007).
Constrained failures of events from a homogeneous Poisson point process yield serial interevent interval correlations. A, Three scenarios in which the failure of primary events from a homogeneous Poisson point process (primary; filled circles) can lead to a gamma mixture of intervals between secondary events (squares, asterisks, and open circles). See description in text. B, SIICCs resulting from the different scenarios applied to simulated event times, plotted as a function of failure probability. The SIICCs after refractoriness are also shown for two simulated spike rates (50 and 100 spikes/s). The mean (filled triangle) and SD (error bars) of the SIICCs obtained from 180 cat-ANF spontaneous spike trains with ≥500 ISIs (cf. Fig. 2C) are also shown and plotted at the failure probability estimated from the ISI distributions (Heil et al., 2007). C, D, SIICCs of 180 simulated spike trains, with numbers of spikes and mean ISIs approximately matching those of the cat-ANF spike trains, generated according to the “regular” scenario (C) and the “irregular” scenario (D), both with a failure probability of 1/3. Note that neither scenario reproduces the dependence of SIICCs on the mean ISI in the cat data (cf. Fig. 2D).
Despite this rule, there are several ways in which failures could occur to yield a given gamma-mixture distribution of secondary IEIs. Three scenarios are illustrated in Figure 7A, using a failure probability of 1/3. This failure probability results in fractions of gamma-2 and gamma-1 intervals of 1/2 each, close to the fractions estimated from the ISI distributions of ANF spike trains (Heil et al., 2007). In the “regular” scenario, every third primary event fails to become a secondary event. In the “irregular” scenario, failures occur irregularly but without violating the never-two-in-a-row rule. In the “block” scenario, every second primary event fails during (the initial) two-thirds of the total time, whereas no events fail during the remaining one-third of the time.
The constrained-failure scenarios cause serial interevent interval correlations
Figure 7B shows the SIICCs for the secondary IEIs resulting from these scenarios for a range of failure probabilities. Each point represents the mean SIICC obtained from 400 simulated repetitions, each with >20,000 secondary IEIs. For a failure probability of 1/2, where every second primary event fails and the distribution of the secondary IEIs is a pure gamma-2 distribution, the mean SIICC is zero for all scenarios, as expected. For the other failure probabilities, however, serial correlations emerge from all scenarios. The regular and the irregular scenarios result in negative mean SIICCs and the block scenario in positive mean SIICCs. For all failure probabilities <1/2, the magnitudes of the mean SIICCs for the irregular scenario are smaller than those for the regular scenario. The most negative mean SIICC (−0.143) is obtained from the regular scenario when every third primary event fails.
To examine the effects that refractoriness might have on the serial correlations in such event trains, we simulated trains of secondary events according to each of the three failure scenarios and applied refractoriness, as described in Materials and Methods, to obtain simulated spike trains. Primary event rates were chosen such that the resulting spike rates had means of ∼50 and 100 spikes/s. The mean SIICCs obtained from these simulations are also plotted in Figure 7B (thin and thick gray lines). The simulations show that refractoriness, as implemented here, generally reduces the magnitude of the SIICCs. The magnitude of the reduction depends on the failure scenario, the failure probability, and the mean spike rate. The reduction becomes more pronounced as the spike rate increases and is largest for the regular scenario with a failure probability of 1/3. The mean SIICC does not change when the failure probability is zero, so refractoriness per se does not introduce serial correlations.
The constrained-failure scenarios do not provide a parsimonious explanation for the dependence of serial correlations in auditory-nerve fiber spike trains on mean interspike interval
The block scenario does not account for the cat data, because the mean SIICC and the associated SD obtained from the 180 cat-ANF spike trains with ≥500 ISIs (Fig. 7B, filled triangle and error bars) falls largely between the SIICCs obtained from simulations of the regular and irregular scenarios. To address whether the latter scenarios can account for the dependence of the SIICC on the mean ISI in the cat data (cf. Fig. 2D), a set of 180 event trains was generated with each scenario, with a failure probability of 1/3, and subjected to refractoriness. The simulated spike trains approximately matched the ANF spike trains with respect to numbers and means of ISIs. Figure 7, C and D, show the resulting SIICCs as a function of mean ISI. SIICCs from the regular scenario (Fig. 7C) are all negative, but, unlike those from the cat data (Fig. 2D), do not increase as mean ISI increases. Rather, the SIICCs tend to become more negative as mean ISI increases, a trend expected from the spike-rate-dependent effect of refractoriness on SIICC magnitude (Fig. 7B). SIICCs from the irregular scenario cluster near zero with a slight tendency to be negative (Fig. 7D) and show no obvious dependence on mean ISI, as expected (cf. Fig. 7B). Thus, to account for the increase in the SIICC with mean ISI in the cat data (Fig. 2D), one could assume that as the mean ISI increases, the failure scenario gradually changes from regular to irregular or the spike trains are increasingly affected by nonstationarities. A third possibility is that the failure probability decreases as mean ISI increases. Although ISI distributions from cat-ANF spike trains were found to be consistent with a constant failure probability (Heil et al., 2007), the probability could not be determined with high confidence in individual spike trains, so a dependence on mean ISI cannot be ruled out. In any case, none of these failure scenarios provides a parsimonious explanation for the dependence of the SIICC on the mean ISI as observed in the cat data.
Similarly, neither scenario provided a parsimonious account of the dependence of the Fano factor on mean ISI (i.e., spontaneous rate; data not shown). In the following, we therefore explore whether another model, termed the depletion–replenishment model, can do so. Such a model has been shown to produce negative SIICCs whose magnitudes decrease as mean ISIs increase (Goldman et al., 2002).
Reproducing cat-ANF spike-train statistics with the depletion–replenishment model
The depletion–replenishment model provides a parsimonious explanation for the dependence of serial correlations in auditory-nerve fiber spike trains on mean interspike interval
The set of 180 cat-ANF spike trains with ≥500 ISIs was simulated 20 times for each of the 150 preliminary pairs of replenishment time constant τrepl and maximum pool size nmax derived using the procedure described in Materials and Methods. By selecting the appropriate depletion probability pdepl separately for each spike train (see Materials and Methods), simulated spike trains had the same mean ISIs (and durations) as the corresponding ANF spike trains. The goal in this section is to determine which pair(s) of τrepl and nmax can best reproduce the relevant statistics for the complete dataset. Each set of simulated spike trains is therefore evaluated on the basis of how well it reproduces the SIICCs from the data.
The SIICCs from one set of simulations with τrepl = 13.7 ms and nmax = 4.0 are plotted in Figure 8A as a function of mean ISI (black dots). For comparison, the SIICCs from the cat data are replotted from Figure 2D (gray dots). The model clearly reproduces the trend in the data; only the positive SIICCs of some ANF spike trains, which are likely due to nonstationarities in the spike rate (cf. Fig. 3), are not reproduced by the simulations. The effect of nonstationarity is also reflected in the differences in the cumulative probabilities of the SIICCs derived from this set of simulations and the cat spike trains (Fig. 8B; cat data replotted from Fig. 2C). The two distributions are close together at more negative SIICCs (<−0.05), but diverge somewhat at less negative and at positive SIICCs. Nevertheless, a Kolmogorov–Smirnov test provided no reason to reject the null hypothesis that the two sets of values are drawn from the same continuous distribution (D = 0.088; p = 0.46). The other 19 sets of simulations for this pair of τrepl and nmax provided similar results.
The depletion–replenishment model provides a parsimonious account of the serial ISI correlations in cat-ANF spike trains and their dependence on mean ISI. A, Plot of the SIICCs from 180 cat-ANF spontaneous spike trains with ≥500 ISIs as a function of mean ISI (gray dots; replotted from Fig. 2D) and from 180 simulated spike trains (black dots). Simulated spike trains were generated by the depletion–replenishment model (using τrepl = 13.7 ms and nmax = 4.0) followed by refractoriness (tD = 0.6 ms and tR = 0.6 ms), with pdepl and train duration varied to match the cat-ANF spike trains with respect to numbers of spikes and mean ISIs. B, Cumulative probability of the SIICCs from these 180 cat-ANF spike trains (gray; replotted from Fig. 2C) and from the simulations (black). A Kolmogorov–Smirnov test provided no reason to reject the null hypothesis of equal distributions. C, Plot of the proportion of times the null hypothesis of equal distributions of SIICCs from the cat and the simulated spike trains had to be rejected as a function of the combinations of τrepl (lower abscissa) and nmax (upper abscissa) used in the simulations. The proportions shown are the means across 20 sets of simulations and after smoothing by a moving average over 21 points. Note the global minimum near τrepl = 13.3 ms and nmax = 3.9 (vertical dashed line).
Next, we examine how well the SIICCs are reproduced across all 150 pairs of τrepl and nmax tested. Figure 8C shows the proportion of times the null hypothesis of equal distributions of SIICCs from the simulated and cat spike trains had to be rejected (averaged across the 20 sets of simulations and smoothed by a moving average over 21 points). This function reveals a global minimum at a τrepl of ∼13.3 ms (Fig. 8C, lower abscissa) and a corresponding nmax of ∼3.9 (Fig. 8C, upper abscissa). With a criterion for a rejection fraction of 0.05, the ranges for τrepl and nmax capable of reproducing the SIICCs from the cat data are 10.8–17.2 ms and 3.2–4.9, respectively.
For comparative purposes, the same analysis was performed for 20 simulated datasets using the constrained-failure model (regular failure scenario with a failure probability of 1/3); the null hypothesis was rejected all 20 times (i.e., the proportion rejected equals 1), so the constrained-failure model does not account for the observed distribution of SIICCs.
The depletion–replenishment model also accounts for the distributions of ISIs from auditory-nerve fiber spike trains
Although the depletion–replenishment model can account for the serial ISI correlations and their dependence on mean ISI, it should also account for the ISI distributions from the cat-ANF spontaneous spike trains. These ISI distributions provide a separate view of the data, and therefore a separate opportunity to estimate the optimal pair(s) of τrepl and nmax for the depletion–replenishment model. The SIICC values assessed in the previous section neither determine nor follow trivially from the ISI distributions, so it is not necessarily the case that the same parameter pair would best describe both; the intervals from an ISI distribution could be arranged so as to yield an arbitrary SIICC value. In this section, each preliminary parameter pair is evaluated on the basis of how well it reproduces the mean ISI distribution from the data. First, we computed the probability density function (PDF) of each ANF spike train with ≥500 ISIs (using a bin width of 0.05 ms) and of its counterpart from the simulations, and computed the difference between the two PDFs in each bin. Figure 9A shows these differences as a function of ISI (gray dots) for the same set of simulations as used for Figure 8A and B (i.e., τrepl = 13.7 ms and nmax = 4.0). The mean difference in each bin is also shown (black). The mean difference scatters unsystematically around zero and its magnitude is small (≤0.001 in all bins). This correspondence is further illustrated in Figure 9B, which shows the mean PDF from the cat data (gray), computed by averaging the probabilities across the 180 spike trains with ≥500 ISIs, and the mean PDF from the corresponding simulations (black). The two mean PDFs are in close agreement. This demonstrates that the depletion–replenishment model can also reproduce the ISI distributions of ANF spontaneous spike trains. The other 19 sets of simulations for this parameter pair provided similar results.
The depletion–replenishment model provides a parsimonious account of the ISI distributions of cat-ANF spike trains. A, Plot of the differences of the PDFs (0.05 ms time steps) of ISIs computed from 180 cat-ANF spontaneous spike trains with ≥500 ISIs and from 180 simulated spike trains as a function of ISI (gray dots). The mean difference is also shown (black). Simulated spike trains were generated by the depletion–replenishment model (using τrepl = 13.7 ms and nmax = 4.0) followed by refractoriness (tD = 0.6 ms and tR = 0.6 ms), with pdepl and train duration varied to match the cat-ANF spike trains with respect to numbers of spikes and mean ISIs. B, Mean PDFs from these 180 cat-ANF spike trains (gray) and from the simulations (black). C, Plot of the root mean square differences (RMSD) between the PDFs from the 180 cat-ANF spike trains with ≥500 ISIs and those from the simulations as a function of the combinations of τrepl (lower abscissa) and nmax (upper abscissa) used in the simulations. The RMSDs are averaged across 20 sets of simulations, normalized to the minimum, and smoothed by a moving average over 21 points. Note the global minimum near τrepl = 13.6 ms and nmax = 4.0 (vertical dashed line).
As with the SIICCs above, we examined how well the PDFs from the cat data and the simulations agree across all 150 pairs of τrepl and nmax tested. The same 20 sets of simulated spike trains per parameter pair were used as for the analysis of the SIICC distributions. Figure 9C shows the root mean square differences between the PDFs from the cat spike trains and the corresponding simulated spike trains (averaged across the 20 sets of simulations, normalized to the minimum, and smoothed by a moving average over 21 points) as a function of τrepl (lower abscissa) and nmax (upper abscissa). This function reveals a global minimum at a τrepl of ∼13.6 ms and a corresponding nmax of ∼4.0, very close to the numbers obtained from the SIICC analyses (cf. Fig. 8C). With a criterion for a normalized root mean square difference of 1.025, the ranges for τrepl and nmax capable of reproducing the ISI distributions from the cat data are 10.7–19.4 ms and 3.2–5.3, respectively.
For comparative purposes, the same analysis was performed for 20 simulated datasets using the constrained-failure model (regular failure scenario with a failure probability of 1/3); the normalized root mean square difference was 1.03. Thus, although the constrained-failure model does not account for the dependence of SIICCs from the cat data on mean ISI (and is therefore not preferred), it does reproduce the ISI distributions well (Heil et al., 2007), almost as well as the depletion–replenishment model.
Fano factors from auditory-nerve fiber spike trains are consistent with the depletion–replenishment model
Finally, we examine whether the depletion–replenishment model also accounts for the Fano factors. For this purpose, we generated 446 spike trains with the depletion–replenishment model to approximately match the cat-ANF spike trains with respect to mean ISI and spike numbers. Parameters were fixed to τrepl = 13.4 ms and nmax = 4.0 to approximate the best values estimated above. Figure 10A summarizes the cat data (446 spike trains) by showing the geometric means of F(T) from each of three sample groups: the one-third of the samples with the lowest SRs, the one-third of the samples with the highest SRs, and the one-third of the samples with SRs in between. The geometric means of F(T) are also plotted for the same sample groups after shuffling the ISIs within each spike train to remove the serial correlations (dashed lines). For very short counting times, shuffling has no effect on F(T), due to the deterministic effect of refractoriness. As the counting time increases, the difference in F(T) between the original spike trains and those with shuffled ISIs initially increases. Generally, the values of F(T) from the original spike trains are lower than those from the trains with shuffled ISIs. This difference is most pronounced for the sample group with the highest SRs and least pronounced for the sample group with the lowest SRs. As the counting time increases further, the differences become smaller and then reverse sign, presumably due to nonstationary spike rates (cf. Fig. 3).
The depletion–replenishment model provides a parsimonious account of the Fano factors from ANF spike trains. A, Geometric mean of F(T) from each of three ANF sample groups (symbols): the one-third of the samples with the lowest SRs, the one-third of the samples with the highest SRs, and the one-third of the samples with SRs in between. The geometric mean of F(T) is also plotted for the same sample groups after shuffling the ISIs to remove the serial correlations (dashed lines). B, Geometric means of F(T) obtained from simulations. Simulated spike trains (n = 446) were generated by the depletion–replenishment model (using τrepl = 13.4 ms and nmax = 4.0) followed by refractoriness (tD = 0.6 ms and tR = 0.6 ms), with pdepl and train duration varied to match the cat-ANF spike trains with respect to numbers of spikes and mean ISIs. The geometric means of F(T) are also plotted for the same simulated sample groups after shuffling the ISIs to remove the serial correlations (dashed lines). Legend in A also applies to B. C, D, Plots of the geometric means of F(T) from the simulations in B against those from the cat data in A, from the original spike trains (C), and after shuffling the ISIs (D). The legend in C also applies to D. The continuous gray line represents the diagonal. Note the tight correlations.
The 446 simulated spike trains were likewise grouped and averaged. The results are shown in Figure 10B. Apart from the trends in the data at larger counting times, attributable to nonstationary spike rates (see above), the mean F(T) behaves similarly to that from the cat data, for all three sample groups, and both before and after shuffling the ISIs (cf. Figs. 10A,B). The close match of the Fano factors from the cat data and the model simulation is emphasized in Figure 10C, which shows the geometric mean of the Fano factors (for all counting times ≤250 ms) from the original spike trains in the cat data plotted against the corresponding values from the simulation. Figure 10D shows the result after shuffling the ISIs.
Discussion
We analyzed nonrenewal properties of cat ANFs during spontaneous activity and first examined whether a previously suggested excitatory process, constrained failures of events from a homogeneous Poisson point process (Heil et al., 2007), can account for them. While such failure scenarios can produce negative SIICCs and Fano factors of the magnitudes observed, they do not straightforwardly account for certain trends in the data. We therefore examined another parsimonious, physiologically more plausible, three-parameter model of vesicle-pool depletion and replenishment. We show that, combined with refractoriness, it accounts not only for the nonrenewal properties of cat-ANF spontaneous spike trains on short timescales but also for the dependence of SIICCs and Fano factors on mean rate and for the ISI distributions, with only release probability varying between spike trains.
Comparison with the literature
Several previous studies have addressed nonrenewal properties of ANF spike trains, particularly over long timescales (Teich, 1989; Teich et al., 1990; Lowen and Teich, 1991, 1992; Jackson and Carney, 2005). Nonrenewal properties over shorter timescales have attracted less attention (Kiang et al., 1965; Lowen and Teich, 1992). Kiang et al. (1965) regarded successive ISIs from cat-ANF spontaneous spike trains as independent, based on visual inspections of plots of each ISI against that which precedes it (also called ISI return maps; Chacron et al. (2000)). They did not compute slopes of best-fit lines, which would have yielded the SIICCs (Wang, 1998). Lowen and Teich (1992) computed SIICCs for 10 long (∼60–620 s) samples of spontaneous activity from different ANFs in the cat. Based on a statistical test and simulated spike trains with exponentially distributed ISIs, Lowen and Teich (1992) suggested significant negative SIICCs in four samples, but significant positive correlations in three, and no significant correlations in the remaining three. SIICCs ranged from approximately−0.045 to 0.094, with a median of 0.013 and interquartile range from −0.038 to 0.087. Their distribution is shifted toward more positive values than our data, which have predominantly negative serial ISI correlations (cf. Fig. 2C). Our analyses (Fig. 3) and simulations (Fig. 7, block scenario) show that positive SIICCs can emerge from nonstationary spike rates. We therefore hypothesize that the more positive SIICCs in the spike trains analyzed by Lowen and Teich (1992) resulted from nonstationary spike rates in the course of the long measurements, and that nonstationarity obscures the negative serial correlations present on shorter timescales.
The functions relating the Fano factor to counting time (Figs. 4A, 5A) generally agree with the only example from ANFs in the literature (Lowen and Teich, 1992; their Figs. 2 and 3) including the effects on F(T) of shuffling the ISIs. We have shown additionally that the power-law-like increase in F(T) with T at long counting times is likely due to nonstationary spike rates (Fig. 4), that the shape of the F(T) function depends on SR, and that ANF refractoriness and SR determine F(T) for short counting times (Fig. 5).
Origin of negative serial ISI correlations
Avila-Akerberg and Chacron (2011) summarize current hypotheses and models of how negative serial correlations may be produced. Generally, the mechanisms must prevent or delay the onset of a spike when the preceding ISI was short. Neuron-intrinsic and network mechanisms have been proposed. For example, Chacron et al. (2000, 2001) assume in their models that the voltage threshold for a spike is dynamic and carries memory. A spike is generated when the voltage equals the threshold. While the voltage is reset to its baseline value after each spike, the threshold is incremented by some constant amount, maintained at this value for the duration of the absolute refractory period, and then relaxed exponentially toward its equilibrium value until the next spike occurs. These properties result in cumulative refractory effects that produce negative serial ISI correlations. Geisler and Goldberg (1966) proposed a related model, in which the memory is carried in the membrane potential following a spike, rather than in the threshold. Specifically, to account for negative SIICCs, the model requires post-spike hyperpolarizations whose amplitude depends inversely on the previous ISI. Delayed inhibitory feedback may also produce nonrenewal properties of spike trains (Doiron et al., 2003).
We show here that the nonrenewal properties of the ANF spike trains could arise from a presynaptic process involving the depletion and replenishment of a pool of readily releasable vesicles. Our model is based on that proposed by Goldman et al. (2002) to account for serial ISI correlations of central neurons receiving input from spiking neurons. It assumes that each release event from the IHC triggers a spike unless the release event falls into the ANF refractory period following the previous spike. There is good evidence for this assumption (Siegel, 1992; Rutherford et al., 2012). Multivesicular release at the ribbon synapse (Glowatzki and Fuchs, 2002; Grant et al., 2010; Graydon et al., 2011; Kim et al., 2013) may ensure that each EPSP is large and fast enough to trigger a spike unless the ANF is refractory. When combined with physiologically plausible refractory periods (see references in Introduction), our model requires only a rate parameter pdepl to vary between spike trains. The maximum size of the readily releasable pool nmax and the replenishment time constant τrepl can be assumed to be constant. This leads to the prediction that, on average, the number of units in the readily releasable pool is smaller at synapses driving high-SR fibers than at those driving low-SR fibers, just as has been observed (Kantardzhieva et al., 2013). With these parsimonious assumptions, the model not only accounts for the nonrenewal properties of ANF spike trains on shorter timescales, but also for their dependence on the mean ISI, as well as for the ISI distributions. In the mammalian inner ear, each ANF is driven by a single ribbon synapse, but each IHC is likely contacted by ANFs of different SRs (Liberman et al., 1990; Moser et al., 2006; Wu et al., 2014). This unique arrangement suggests that the differences in pdepl likely originate from differences in release probability at individual ribbons. If the IHC is isopotential, such differences could arise, for example, from differences in the number or spatial arrangement of the voltage-gated Ca2+ channels required for exocytosis, as discussed elsewhere (Heil and Neubauer, 2001, 2010; Frank et al., 2009; Heil et al., 2011). The value of nmax was estimated here to be ∼4. If each release event comprises the coordinated simultaneous release of ∼8 vesicles (Grant et al., 2010), the readily releasable pool would consist of ∼32 vesicles, within the range estimated by other methods (Lenzi et al., 1999, 2002; Khimich et al., 2005; Nouvian et al., 2006; Kantardzhieva et al., 2013; Wong et al., 2014). With the further assumption that τrepl does not change during acoustic stimulation, the estimate of ∼13.5 ms would support a maximum release rate of 256 s−1 if units are released the moment they are replenished. This rate is high enough to account for the maximum sustained spike rates observed in these ANFs (Heil et al., 2011). Of course, we cannot exclude variation in nmax or τrepl between synapses or according to release history, but such assumptions are not needed to explain the data.
We also cannot exclude other mechanisms (e.g., those discussed above) contributing to the spiking properties analyzed here. Comparison of the statistics of EPSP intervals with those of ISIs may help identify possible contributions by spiking-related mechanisms.
Functional significance of the nonrenewal properties of ANF spike trains
As shown repeatedly (Ratnam and Nelson, 2000; Chacron et al., 2001; Avila-Akerberg and Chacron, 2011), the narrow spike-count distributions (PNDs) that result from negative serial ISI correlations have implications for information transmission and stimulus detection, at least in the ideal-observer framework of signal detection theory (Green and Swets, 1966). The theory assumes Gaussian-like distributions of an internal variable in the absence (“noise”) and presence (“signal plus noise”) of a stimulus, which differ in mean (and possibly variance). The d′ measure for discriminability of a signal from noise (i.e., for detection) is defined as the ratio of the absolute difference of the means of the two distributions divided by the square root of their summed variances. If a weak stimulus gives rise to some extra spikes, the d′ measure will therefore increase more when the variances of the two distributions are narrower. The negative serial ISI correlations and reduced Fano factors in the spontaneous spike trains of ANFs (relative to those after shuffling the ISIs) are thus beneficial for the detection of weak signals. It is noteworthy that the high-SR fibers, which have the lowest rate thresholds (Winter et al., 1990; Yates, 1991; Taberner and Liberman, 2005) and therefore have been seen as responsible for weak signal detection, are also those with the most negative SIICCs.
Conclusion
The ISI distributions, the negative serial ISI correlations and the spike-count variability, and their dependencies on mean ISI (or spike rate), as observed in spontaneous activity of cat primary auditory afferents, are all consistent with a physiologically plausible, parsimonious model assuming vesicle-pool depletion and replenishment. Only the release probability must be assumed to vary, whereas the maximum pool size and the mean time required for replenishment of a released unit can be constant. We suggest that the organization of the IHC ribbon synapse not only enables indefatigable neurotransmitter release, but also imposes regularity on the release process, particularly at high rates.
Footnotes
This work was supported by a grant of the Deutsche Forschungsgemeinschaft to P.H. (He1721/11-1) within the Priority Program 1608. Data collection was supported by other grants (He1721/5-1 and He1721/5-2). We are grateful to Mel Brown for help with data collection and to Artur Matysiak and Pawel Kordowski for their insights regarding the operation of refractoriness.
The authors declare no competing financial interests.
- Correspondence should be addressed to Peter Heil, Leibniz Institute for Neurobiology, Brenneckestrasse 6, 39118 Magdeburg, Germany. peter.heil{at}lin-magdeburg.de
