Modeling the influence of short term depression in vesicle release and stochastic calcium channel gating on auditory nerve spontaneous firing statistics

We propose several modifications to an existing computational model of stochastic vesicle release in inner hair cell ribbon synapses, with the aim of producing simulated auditory nerve fiber spiking data that more closely matches empirical data. Specifically, we studied the inter-spike-interval (ISI) distribution, and long and short term ISI correlations in spontaneous spiking in post-synaptic auditory nerve fibers. We introduced short term plasticity to the pre-synaptic release probability, in a manner analogous to standard stochastic models of cortical short term synaptic depression. This modification resulted in a similar distribution of vesicle release intervals to that estimated from empirical data. We also introduced a biophysical stochastic model of calcium channel opening and closing, but showed that this model is insufficient for generating a match with empirically observed spike correlations. However, by combining a phenomenological model of channel noise and our short term depression model, we generated short and long term correlations in auditory nerve spontaneous activity that qualitatively match empirical data.


INTRODUCTION
In the vertebrate auditory pathway, the inner hair cell and auditory nerve (IHC-AN) complex is the principal structure for the transduction of basilar membrane motion to stochastic trains of action potentials (Mulroy et al., 1974;Glowatzki and Fuchs, 2002;Johnson et al., 2009;Matthews and Fuchs, 2010). A computational model of the IHC-AN complex was proposed by Meddis (1986), and later modified by Sumner et al. (2002) to become a component in a larger computational model of the transformations of sounds by the middle ear. Unlike the Meddis (1986) model, in the Sumner et al. (2002) model, vesicle release from the IHC to the cleft was conceptualized as quantal and accruing with a probability that had a third power dependence on presynaptic calcium concentration. Later, the Sumner et al. (2002) model was modified by Meddis (2006) to take into account more physiological functions.
Here, we present a revised version of the Meddis (2006) model of the IHC-AN complex, with the aim of enhancing understanding of the biophysical sources of stochastic variability in the IHC-AN complex, by generating auditory nerve spontaneous spiking that provides an improved statistical match with empirical data.
The Meddis (2006) model includes a probabilistic "relative refractoriness" component, which is designed to replicate observed variation in the minimum time between spikes in AN fibers. Here we propose a pre-synaptic physiological explanation as the cause for what is attributed to post-synaptic relative refractoriness (note that we do not alter the original model's "absolute refractory" period, which models spike generation and membrane potential recovery). Specifically, we introduce a model of short term depression in pre-synaptic vesicle release, similar to short term plasticity models developed for cortical synapses (Tsodyks and Markram, 1997;Scott et al., 2012;Hennig, 2013;McDonnell et al., 2013). Unlike most such models, the conceptual model here is that there is a temporarily reduced probability of pre-synaptic vesicle release, following each actual release. Also unlike those models, the input to the synapse is not discrete spiking events, but instead the continuously valued membrane potential of the inner hair cell. The reason our model is suitable for capturing phenomena that have traditionally been attributed to post-synaptic relative refractoriness is that it introduces variability in the time between vesicle releases, which in turn leads to variability in the minimum time between post synaptic spikes. Our reasons for seeking this alternative conceptual model are given in the Discussion section.
We compare the resulting auditory nerve spontaneous firing statistics of our model with the firing statistics published by Heil et al. (2007). For spontaneous neural activity in auditory nerve fibers, inter-spike interval (ISI) distributions have been shown by Heil et al. (2007) to match empirical data better if the vesicle release inter-event interval (IEI) distribution was assumed to be a mixture of an exponential function and a gamma function with shape factor 2, both having the same scale parameters. We show that the probability density function (PDF) of ISI data obtained by Heil et al. (2007) fits PDF of ISI data obtained from our simulation if the time constant of short term depression is assumed to be around 2.5 ms.
Short and long term correlations have been observed in the spontaneous activity of auditory nerves (Teich, 1989;Lowen and Teich, 1992;Teich and Lowen, 1994). For individual auditory nerve fibers, it was shown that the Fano factor for spike counts increases for time scales from around 100 ms to tens of seconds indicating positive long term correlation and decreases slightly for time scales of around tens of milliseconds indicating short term negative correlation (Teich, 1989;Lowen and Teich, 1992;Teich and Lowen, 1994). Here we include a calcium channel noise model in the Meddis (2006) model. We show that for spontaneous activity, this biophysical noise model does not generate the short and long term correlations observed in the Teich and Lowen (1994) Fano factor curves.
However, we also modify the Meddis (2006) model to include a combination of a phenomenological model of IHC calcium channel noise and our model of short term depression in vesicle release. Using this model, for auditory nerve spontaneous activity, we generate Fano factor time curves that qualitatively match empirical Fano factor time curves of Teich and Lowen (1994); Teich (1989); Lowen and Teich (1992).

MATERIALS AND METHODS
Firstly, in Section 2.1, we review the previous models that our research is built upon: • The inner hair cell model of Meddis (2006).
Then in Section 2.2, we provide a review of previous statistical analysis of empirical auditory nerve spontaneous activity data including research published by Heil et al. (2007), Teich and Lowen (1994), Teich (1989) and Lowen and Teich (1992). The final models we describe in Section 2.3 are our modifications to the Meddis (2006) model. These are designed to enhance understanding of the biophysical origin of stochastic variability in AN spiking, and to generate auditory nerve spontaneous spiking that provides an improved statistical match with empirical results, as described in Section 2.2. Meddis (2006) describes a deterministic calcium-dependent model for converting the membrane potential of an inner hair cell, v(t), to a vesicle release rate, k(t). We use c(t) to describe the intra-cellular calcium concentration (relative to its rest concentration) as a function of time. In the model, the release-rate for available vesicles, k(t), is proportional to the cube of c(t). The calcium concentration depends on four constants, τ c , G c , E c , ν, on the membrane potential, v(t), and on an additional variable, m(t), where m 3 (t) represents the fraction of open channels at time t as well as the probability of a calcium channel to be open. This depends on three constants, γ , β and τ m , and on v(t). Note that m 3 (t) is bounded to the interval [0, 1], which is essential for it to physically represent a fraction of open channels. The maximum value of 1 occurs when v(t) is large and positive and the minimum value of 0 occurs when v(t) is large and negative.

Inner hair cell model
In summary, the model has the following parameters: • b is a parameter that can be varied to match data.
• E c is the calcium reversal potential.
• G c is the maximum calcium conductance.
• τ c is the time constant of calcium clearance.
• τ m , γ and β are constants that describe the voltage-dependent calcium current flow. • ν is the unit correction constant.
The values of these parameters are summarized in Table 1. The equations describing conversion from v(t) to k(t) are where k(t) has units of releases per second. We have modified the Meddis (2006) and Sumner et al. (2002) models by introducing a constant ν with units of MA −1 s −1 to ensure all terms in Equation (2) have units of Ms −1 , where M is the unit of molar concentration. By fitting to the saccular hair cells of the bull-frog data, it has been shown (Hudspeth and Lewis, 1988) that where F is Faraday constant, C v is the cell volume, ζ is the fraction of cell volume where calcium is accumulated to and L is the proportion of free calcium in the neuron. The values of these parameters are summarized in Table 2, with the result that ν = 2.3 × 10 9 MA −1 s −1 .
To confirm that our proposed model enhancements have no effect on previously established model features, in the Results section we compare the average vesicle release rates obtained from simulation of the proposed model to the average vesicle release  Values obtained from Hudspeth and Lewis (1988).
rate obtained from simulation of the Meddis (2006) model. We introduce the notation k as the simulated average vesicle release rate. We show that the changes that we make to Meddis (2006) model result in k that are close to k obtained from the original model of Meddis (2006). The parameter k for the various proposed models are summarized in the tables. A positive calcium current is required to increase the calcium concentration but in the Meddis (2006) and Sumner et al. (2002) models, calcium current is negative (i.e., inward) when v(t) < E c . Therefore, we have used (E c − v(t)) in Equation (2) instead of (v(t) − E c ) used in the Sumner et al. (2002) and Meddis (2006) models. The max( · ) function is included in Equation (1) since although it is possible for c(t) < 0 in the model (which represents calcium concentration less than its rest value), the rate k(t) cannot be negative. Note that the final term in Equation (2) has the form of the deterministic Hodgkin and Huxley (1952) voltagegated ion channel current model. Later, we replace this with a model of stochastically opening and closing ion channels.

Deterministic synapse model
The input to the deterministic synapse model of Meddis (1986) is the rate at which the neurotransmitter is released to the cleft, k(t). There are three continuous-time-dependent variables that describe transport between a vesicle "factory," an "immediate store," the synaptic cleft, and a vesicle "recycling pool": • the amount of releasable neurotransmitter, where M is the maximum amount of neurotransmitter in the immediate store. • the amount of neurotransmitter in the cleft, y(t).
• the amount of neurotransmitter being recycled, z(t).
There are four parameters that have units of rate: • r 1 is the rate of manufacture of neurotransmitter from the "factory." • r 2 is the rate of restoration of neurotransmitter from the recycling pool. • r 3 is the rate at which neurotransmitter is lost in the cleft.
• r 4 is the rate at which neurotransmitter is moved from the cleft to the recycling pool.
The values of these parameters are summarized in Table 3. The deterministic Meddis (1986) synapse model is of the following form where

Stochastic synapse model
Subsequently, Sumner et al. (2002) and Meddis (2006) modified Meddis (1986) to build a model where movement of neurotransmitter is stochastic rather than deterministic and neurotransmitter in the immediate store is quantal rather than continuous. The stochastic Meddis (2006) synapse model is of the following form, Stochastic movement of discrete vesicles of neurotransmitter is described by the binomial random variable, B(ρ, n): if there are n vesicles available during a small dt, each with equal probability of moving ρdt, then B(ρ, n) is the number of vesicles moving during dt. Vesicles in the immediate store are quantal so z(t) is mapped to the largest previous integer, z(t) .

Phenomenological synapse model
It has been shown that by using rate estimates from a fractional Gaussian noise driven Poisson process model, the shape of published histograms of spontaneous discharge rate (Liberman, 1978) can be reproduced (Jackson and Carney, 2005). This has been incorporated into a phenomenological model of the synapse in the IHC-AN complex by Zilany et al. (2009); Zilany and Carney (2010); Zilany et al. (2014). This synapse model has both exponential and a power-law adaptation functions. The exponential adaptation is implemented using the diffusion model of Westerman and Smith (1988). The exponential adaptation path is followed by two parallel fast and slow power-law adaptation function. The fractional Gaussian noise is incorporated in the slow power-law adaptation path. The input to the synapse model is the relative membrane potential of the inner hair cell.

Models for converting vesicle release to AN spikes
In the deterministic rate model of Meddis (1986), the amount of neurotransmitter in the cleft causes a post-synaptic spike at time t with probability, where h is a constant. An absolute refractory period of 1 ms during which no spike can occur is applied. A relative refractory period is not considered.
In the quantal stochastic model of Meddis (2006), each ejected vesicle to the cleft can generate a spike in the auditory nerve after an absolute refractory period (ARP) and relative refractory period (RRP) are considered. If a vesicle is released, a spike in the postsynaptic AN is generated if p conv (t) is greater than a uniformly distributed random number between 0 and 1.
where C r = 1, t R = 0.6 ms is the time constant of relative refractoriness, t A = 0.75 ms is the ARP, t is the current time, and t l is the time of the previous spike. The conversion model of Sumner et al. (2002) is very similar to the conversion model in the Meddis (2006) model. The differences are that in the Sumner et al. (2002), C r = 0.55 and t R = 0.8 ms.
In the Zilany et al. (2009), Zilany and Carney (2010) and Zilany et al. (2014) spike generator model, spike times in the auditory nerve are generated by a renewal process that simulate a non-homogeneous Poisson process driven by the output of the synapse model. Heil et al. (2007) has shown that the empirical ISI distribution for spontaneous neural activity in cat auditory nerve fibers is better described if the IEI distribution for vesicle release events is a mixture of an exponential distribution and a gamma distribution. The gamma distribution has a shape parameter equal to two, and both the gamma distribution and the exponential distribution have the same scale parameter.

Empirical vesicle release distribution
To calculate the ISI parameters, ARP and RRP in the form of Equation (11) are used. Two additional parameters are involved: • θ is the scale factor for both the exponential distribution and the gamma distribution; • ρ is the fraction of gamma distribution in the mixture. Heil et al. (2007) obtained the following equation describing the ISI probability density function (PDF),

Empirical firing correlations
The Fano-factor time curve is a measure of correlation over time.
Fano-factor is dispersion in a variable, as a function of an increasing time-window for obtaining data on which to estimate the dispersion. For a spike train, the Fano-factor is the variance of the number of spikes in a time window divided by the mean of number of spikes from a single spike train in that time window. We denote: • T as the size of a specific counting time window.
• F(T) as the Fano-factor for window size T. Teich and Lowen (1994), Kelly (1994), Teich (1989), Lowen and Teich (1992) plotted empirical Fano-factor time curves for neural activity in mammalian auditory nerve fibers as seen in Figures 1A,B. The Fano-factor is 1 for sufficiently small time windows. It slightly decreases to below 1 over time scales on the order of tens of ms after which it increases monotonically and reaches more than 10 for time windows of a few tens of seconds. It has been shown that negative short term correlation observed in the Fano factor curve of spontaneous activity of a simulated AN fiber model with second order refractory behavior matches the data of Lowen and Teich (1992) for time windows between 15 ms and 100 ms (Gaumond, 2002).

Short-term depression in vesicle release probability (STD v )
In AN spontaneous spike trains, the shortest ISIs occur much less frequently than the most likely ISIs (Heil et al., 2007). In Meddis (2006), this feature of ISI statistics is accounted for by ARP and RRP. Given this model includes variable relative refractory times in AN fibers, during which pre-synaptic vesicle release is unaffected, this would mean many vesicles are released that do not give rise to spikes. We therefore seek an alternative model in which what has been attributed to refractoriness is actually mainly due to pre-synaptic effects, due to vesicles not being released at all for durations longer than the absolute refractory period of the ANs. We return to this in Discussion.
Our hypothesis is that all vesicle releases, apart from any that occur during the absolute refractory period, cause action potentials, but that vesicle release is subject to short term depression. We introduce short term depression to pre-synaptic release probability in a manner analogous to standard stochastic models of cortical short term depression (Tsodyks and Markram, 1997;Wang, 1999;Hennig, 2013;McDonnell et al., 2013). In this model, immediately following release, the probability of release drops dramatically and then, increases back to a baseline level over a time frame that matches the spike data.
There are two additional parameters introduced in this model:  Teich and Lowen (1994), Kelly (1994), Teich (1989). (B) Time-window dependent Fano factor of (A) on linear axes for time windows shorter than 0.1 s.
• τ s is the time constant of short term depression.
• a is a fraction indicating an instantaneous decrease in release probability.
The model for the change of k(t) over time is (13) where t v i is the time of ith release.

Channel noise in inner hair cell calcium channels
Auditory nerve spike trains show positive long term correlation and usually negative short term correlation (Teich, 1989;Lowen and Teich, 1992;Kelly, 1994;Teich and Lowen, 1994). We hypothesize that a possible origin of the correlation is stochastic variability in the inner hair cell calcium channels. A biophysical model and a phenomenological model of calcium channel noise in the inner hair cell are built. (2006) (2006) model. Other possible origins of the observed long term correlation have been suggested, including fractal ion channel gating (Teich, 1989;Liebovitch and Toth, 1990), fractal behavior of the specialized proteins with direct role in exocytosis (Lowen et al., 1997), self-organized criticality in ion channel gating for example due to ion-conformational interaction (Kharkyanen et al., 1993;Brazhe and Maksimov, 2006), and fractal dynamics of transmitter diffusion in the synaptic junction (Teich, 1989).

Biophysical model. In the Meddis
An integrate and fire model with renewal point process input has been suggested to be capable of producing long term correlation that matches empirical data from spike trains of cortical neurons (Jackson, 2004). Unlike cortical neurons, inner hair cells encode graded input with a graded membrane potential (Van Steveninck and Laughlin, 1996). We aim to cast light on the possible biophysical mechanisms in the IHC-AN complex that can produce renewal point processes and hence long term correlation in the spike trains of auditory nerves. Meddis (2006) assumes the calcium concentration dependence of the release probability to be due to voltage dependent calcium channels. We hypothesize that a possible biophysical mechanism of the fractional Gaussian noise in Jackson and Carney (2005) We introduce to the Meddis (2006) model a four-state model of channel gating with standard transition rate formulae (Goldwyn and Shea-Brown, 2011;Schmerl and McDonnell, 2013), Equation (2) therefore changes to where n(t) is the number of open calcium channels out of total of N calcium channels and g c is the single calcium channel conductance.

Phenomenological model.
We consider a phenomenological model of calcium channel noise that we add to the Meddis (2006) model. Instead of modeling discrete channel noise, we add an Ornstein Uhlenbeck process to the mean fraction of open calcium channels, m 3 (t). Equation (2) changes to: where f ( · ) := max (0, min (1, ·)) ensures the fraction of open channels is restricted to the interval [0, 1] and X(t) is a noise driven from Ornstein Uhlenbeck process; i.e., where W t denotes the Wiener process and the mean (μ o ), time constant (τ o ) and variance (σ o ) of the noise are positive constants.

Noise in inner hair cell membrane potential
We also consider an alternative phenomenological model of noise where the IHC membrane potential is subject to an additive Ornstein Uhlenbeck process. Equation (2) changes to:

Combination of short term depression in vesicle release model and phenomenological calcium channel noise model
A possible origin of short term correlation in AN spike trains is a form of refractoriness (Teich and Lowen, 1994). We introduce a model that combines short term depression in vesicle release and phenomenological calcium channel noise as follows where k 1 (t) is the vesicle release rate when Ornstein Uhlenbeck noise is added as calcium channel noise in the Meddis (2006) model.

PARAMETERS
The parameters in Table 1 (except g c ), in Table 3, and for t A and t R (except in Table 5) were obtained from publicly available Matlab source code MAP BS at http://www.essexpsychology.macmate. me/HearingLab/modelling.html. The parameters in Table 2 and for g c were obtained from the literature (Hudspeth and Lewis, 1988;Zampini et al., 2013). The values of b and bc(t) 3 were chosen in order to produce the desired spontaneous rates in AN fibers. In Table 5, the parameters τ s , t R and a were obtained through parameter searches, in order to obtain a close quantitative fit to the data of Heil et al. (2007), while keeping estimated values of θ and ρ close to the results of Heil et al. (2007). The parameters τ o , σ o and μ o were chosen to produce spontaneous activity in the auditory nerve that qualitatively matches the empirical Fano factor data of Teich and Lowen (1994). The maximum number of readily releasable vesicles in the immediate store, M, in the Meddis (2006) model is considered to be 10. Moser and Beutner (2000) suggested the average number of vesicles in the immediate store to be about 14 vesicles per active zone. Khimich et al. (2005)   Empirical data from spontaneous activity in auditory nerve with SR∼ 65 spikes per second created from data in Teich and Lowen (1994). (B) Time-window dependent Fano factor of (A) on linear axes for time windows shorter than 0.1 s. Parameters are summarized in Table 4 (2006) model with a maximum number of readily available vesicles of M = 20 and a spontaneous rate of about 160 spikes per second. The Fano factor curves in red and green increase to higher values than the Fano factor curve in blue does for large time windows and are a better qualitative match to the empirical data of Figure 2A (Light blue). The Fano factor does not reach 10 for sufficiently large time windows. The magnitude of the decrease in Fano factor below one for shorter time windows is comparable to the empirical data of Figure 2A (Light blue.) As shown in Figure 2A, long term correlation in the Meddis (2006) model can be partially produced if either the maximum number of releasable vesicles is decreased or the firing rate is increased, both of which cause depletion of available vesicles in the immediate store. In this model, low spontaneous rate fibers are associated with smaller pools of vesicles, and high spontaneous rate fibers are associated with larger pools of vesicles. In the depleted model, the time scales of the correlation do not match the empirical data of Figure 2A (Light blue). Depletion of vesicles moves the onset of short and long term correlations to smaller time windows.
Depletion of available vesicles in the Meddis (2006) model (by decreasing the maximum number of available vesicles from 20 with SR of about 65 spikes per second to 6 with SR of about 65 spikes per second, or by increasing the spontaneous rate to around 160 spikes per second with maximum number of available vesicles of 20) produces an average vesicle release rate, k, of 107 and 55 (s −1 ) respectively that are both much larger than 5 (s −1 ), which is the k of the Meddis (2006) model with a maximum number of available vesicles of 20 and the spontaneous rate of about 65 spikes per second.

SHORT-TERM DEPRESSION MODEL
Here we consider the case where the relative refractoriness component of the Meddis (2006) model is removed and we use our alternative model of short term depression in vesicle release probability. That is, the release rate in the Meddis (2006) model , k(t), as given by Equation (1), was replaced by k(t) of Equation (13), and relative refractoriness in the auditory nerve was omitted. Using this model, ISI data for spontaneous activity in an AN fiber was simulated.
The effect of substituting relative refractoriness in the auditory nerve with short term depression in vesicle release in the Meddis (2006) model is more clearly observed in the simulated data when the absolute refractory period is (unrealistically) assumed to be zero. Under this assumption, Figure 3A shows that in the Meddis (2006) model, similar to including relative refractoriness in the auditory nerve, the alternative model of short term depression in vesicle release leads to the least probable ISIs being larger than they would otherwise be.
A distribution fitting application which returns maximum likelihood estimations of the model parameters was used to estimate the parameters that produce the best fit of the simulated ISIs to the empirical results. Figure 3B shows that the PDF of the simulated data for the Meddis (2006) model with AN relative refractoriness replaced by short term depression in vesicle release in blue and the best fit to Equation (12) in red. The refractory time constants, t A and t R , were kept at fixed values. The free parameters, θ and ρ, were estimated.
The models in Figures 3A,B were fitted to Equation (12), and the corresponding values of θ and ρ were estimated and summarized in Table 5. Parameters τ s and a were obtained through parameter search in order to obtain a good fit to data while keeping θ and ρ close to the result of Heil et al. (2007).
In two different neurons, Heil et al. (2007) obtained θ = 0.0988 (ms −1 ) and ρ = 0.39 for t A = 0.69 ms, and t R = 0.58 ms when SR = 65 spikes per second and θ = 0.0862 (ms −1 ) and ρ = 0.43 for t A = 0.73 ms, and t R = 0.41 ms when SR = 57.1 spikes per second. Using the short term depression in vesicle release model, we estimated θ and ρ to be 0.05 (ms −1 ) and 0.37, respectively. Thus, Heil et al. (2007) scaling factors, θ, and fraction of gamma distribution in the mix, ρ, are comparable to what we obtained with our model with comparable spontaneous rate.
However, while (Heil et al., 2007) assumed the post-synaptic refractory period to be less than 1 ms, we obtain our result with a post-synaptic refractory period of a few milliseconds. Despite Fitting PDF of ISI data for Meddis (2006) model with relative refractoriness in the auditory nerve substituted by short term depression in vesicle release where t A = 0.75 ms, τ s = 2.5 ms and bc 3 = 5 s −1 to Equation (12). In all traces, SR∼65 spikes per second. Parameters are summarized in Table 5 this difference, our model has introduced three features to the model's ISI distribution that are common with the data of Heil et al. (2007): an ISI PDF with a single maxima such that the PDF increases from zero to its peaks for small ISIs just above the absolute refractory period, a comparable scale factor and a comparable fraction of gamma distribution in the mix of exponential and gamma distributions.

Biophysical model
Here we consider the case where the biophysical model of calcium channel noise is added to the Meddis (2006) model. That is, in the Meddis (2006) model the state of each of N calcium channels is simulated, stochastic changes of states based on the state diagram of Figure 4 are permitted, and c(t) given in Equation (2) is replaced by c(t) given by Equation (16). Using this model, the time-window dependent Fano factor of spike counts in the auditory nerve model for different numbers of calcium channels were obtained and shown in Figures 5A,B. Unlike the empirical data of Figure 5A (Light blue), the Fano factor does not increase steadily to a value around 10 for long time windows. A slight decrease in Fano factor for shorter time windows is observed.
In the hair cells of a chick's cochlea, for each hair cell, around 100 calcium channels for short hair cells and 341 for tall hair cells are suggested (Martinez-Dunst et al., 1997), which in turn suggest quite small numbers of channels per synapse. In our model, no improvement was seen in long term correlation by decreasing the number of calcium channels from 200 (Black) to 50 (Green), 10 (Red) and 5 (Blue). We conclude that this calcium channel model fails to add a long term correlation to the spike trains of the auditory nerve in the Meddis (2006) model in a way that matches experimental observations shown in Figure 5A (Light blue).
Adding the biophysical calcium channel model with parameters summarized in Table 6 to the Meddis (2006) model produces k of 4, 6, 5 and 4 (s −1 ) respectively for 5, 10, 50, and 200 calcium channels which are all close to 5 (s −1 ), the k of the original Meddis (2006) model with a spontaneous rate around 65 spikes per second and a maximum number of available vesicles of M = 20.

Phenomenological model
Here we consider the case where the phenomenological model of calcium channel noise is added to the Meddis (2006) model. That is, in the Meddis (2006) model, Equation (2) is replaced by Equation (17).
Using this model, the time-window dependent Fano factor of spike counts in the auditory nerve model were obtained and shown in Figures 6A,B (Blue) . It can be seen in Figure 6A (Blue) that, like empirical Fano factor of Figure 6A (Light blue), the Fano factor increases to about 10 for large counting time windows. But, the Fano factor in Figure 6A (Blue) does not decrease below one for shorter time windows as much as the empirical Fano factor shown in Figure 6A (Light blue) does.
Adding the phenomenological channel noise with parameters summarized in Table 7 to the Meddis (2006) model produces k of 7 (s −1 ) which is close to 5 (s −1 ), the k of the original Meddis (2006) model with a spontaneous rate around 65 spikes per second and maximum number of available vesicles of M = 20.

COMBINING SHORT-TERM DEPRESSION AND CALCIUM CHANNEL NOISE
Here we consider a combination of short term depression in vesicle release with the phenomenological model of channel noise within the Meddis (2006) model. That is, in the Meddis (2006) model, k(t) from Equation (1) was replaced by k(t) from Equation (20) and relative refractoriness in the auditory nerve in the AN fiber was omitted. Figures 6A,B (Green) show the time-window dependent Fano factor for auditory nerve fiber spike counts for this model. The Fano factor for this model increases steadily to about 10 for large FIGURE 5 | (A) Time-window dependent Fano factor for spontaneous activity in an auditory nerve fiber model using the biophysical model of calcium channel noise in the IHC-AN complex applied to the Meddis (2006) model for different numbers of calcium channels, N. Light blue: Empirical data from spontaneous activity in auditory nerve with SR∼65 spikes per second created from data in Teich and Lowen (1994). (B) Time-window dependent Fano factor of (A) on linear axes for time windows shorter than 0.1 s. Parameters are summarized in Table 6 Meddis (2006) model. Light blue: Empirical data from spontaneous activity in an auditory nerve fiber with SR∼65 spikes per second created from data in Teich and Lowen (1994). counting time windows. It can be seen in Figure 6B (Green) that for counting time windows of a few milliseconds, Fano factor decrease is slightly more than that of Figure 6A (Blue) and hence a better match to the empirical data of Figure 6A (Light blue). Adding the combination model of phenomenological channel noise and short term depression in vesicle release with parameters summarized in Table 7 to the Meddis (2006) model produces k of 5 (s −1 ) which is the same as the k of the original Meddis (2006) model with a spontaneous rate around 65 spikes per second and maximum number of available vesicles of M = 20.
As the maximum number of available vesicles in the immediate store decreases, as shown in Figure 6C, the corresponding minima in the Fano factor curve for shorter time windows increases and the short and long term correlations compare quantitatively to the results from the Zilany et al. (2014) model. k = 10 remains close to the k = 5 from the original Meddis (2006) model with a spontaneous rate around 65 spikes per second and a maximum number of available vesicles of M = 20.
This combination model produces a release rate for which the baseline level is mainly controlled by Ornstein Uhlenbeck noise and the post release behavior is mainly controlled by short term depression in vesicle release as shown in Figure 6D (Green).

COMPARISON OF CALCIUM CHANNEL NOISE WITH MEMBRANE POTENTIAL NOISE
Here we consider the inclusion of the phenomenological model of noise in the inner hair cell membrane potential model in (Meddis, 2006) model. That is, in the Meddis (2006) model, Equation (2) is replaced by Equation (19). The time-window dependent Fano factor for AN spike counts in this model is shown in Figures 6A,B (Red). Like the situation of Figure 6A (Blue) where the Ornstein-Uhlenbeck noise is instead included as calcium channel noise, the Fano factor increases steadily to 10 for larger counting time windows, but it decreases below unity less than the empirical Fano factor of Figure 6A (Light blue) for smaller counting time windows.
Adding the phenomenological membrane potential noise with parameters summarized in Table 7 to the Meddis (2006) model produces k of 5 (s −1 ) which is the same as the k of the original Meddis (2006) model with a spontaneous rate around 65 spikes per second and a maximum number of available vesicles of M = 20.

DISCUSSION
We have shown that adding a combination of short term depression in vesicle release, and time-correlated channel noise, to the existing model of Meddis (2006) results in qualitatively similar results for spontaneous inter-spike interval correlations observed in empirical data. We make the case that it is the qualitative features of the Fano factor curve (namely the occurrence of positive and negative correlations, and the order of magnitude of the positive correlation) that are of most interest. Our model generates auditory nerve spontaneous spike trains for which the spike-count Fano-factor matches empirical data at short and long time scales qualitatively. The qualitative features of the Fano factor curve obtained from the proposed models are summarized in the last columns of Tables 4, 6-8. However, the time scales of maximum negative correlation, and the onset of positive correlation do not exactly match the data. Moreover, the long term correlation in the biophysical model of IHC calcium channel noise does not match empirical data. There are several reasons for these discrepancies. First, the simulation data is only as good as the overall model, which omits many details of the complex calcium channel dynamics of ribbon synapses. We have seen, for example, that a standard biophysical model of channel noise does not induce long-term correlations, while replacing that model with a phenomenological model based on Ornstein-Uhlenbeck noise does so. We suggest that a more biophysically detailed model of calcium channel noise can improve the long term correlation to match empirical data. For example, a model where a single calcium channel controls vesicle release at each docking site (Weber et al., 2010) could potentially lead to a more complicated release dynamics and might produce long term correlation in the auditory nerve spontaneous spiking activity. A second reason might be that the parameters we used (including parameters in the Meddis (2006) model) need to be better tuned to match the empirical data. We left this for future work.
There are several justifications for replacing auditory nerve relative refractoriness with short term depression in vesicle release probability in the model. First, extensive neurotransmitter release can be toxic to neural tissues and cleaning up the excessive transmitters by glia cells requires a large amount of energy (Glowatzki et al., 2006). Short term depression in vesicle release will reduce the number of vesicles released, which in turn will reduce potential for toxicity and energy usage. Moreover, since it is thought that single vesicle produces spikes in AN fibers, for energetics reasons it is wasteful to release vesicles during the refractory period when spikes cannot occur. A possible mechanism for short term depression in vesicle release could be the presence of auto-inhibitory metabotropic receptors called auto-receptors (Billups et al., 2005). To our knowledge, however, there is no evidence either for or against the presence of such auto-receptors in inner hair cells. Alternatively, it is possible that complex intra-cellular calcium dynamics and its relationship to vesicle exocytosis could cause such effects.
We hypothesize that observations of variable minimum time between spikes attributed to "relative refractoriness" above) in the IHC-AN complex is mainly due to pre-synaptic effects, namely that vesicle release sometimes doe not occur for a period longer than are the absolute refractory period. However, it is also possible that actual relative refractoriness in auditory nerve recovery following a spike (Cartee et al., 2000), and short term depression in vesicle release probability in the ribbon synapse could co-exist.
To obtain a fit close to the data of Heil et al. (2007), we have chosen the time constant of short term depression in the vesicle release to be 2.5 ms. Short term depression in vesicle release has been observed in synapses other than the ribbon synapse of inner hair cells (e.g., Stevens and Wang, 1995;Hjelmstad et al., 1997). Whole cell recordings from hippocampal pyramidal neurons showed that a 20 ms refractory period was required between    vesicle releases (Stevens and Wang, 1995). In a different experiment, Hjelmstad et al. (1997) observed a 6-7 ms period following release during which the synapse was incapable of transmission. Consequently, the time-scale of 2.5 ms is potentially biologically plausible.
In this paper we aimed to simulate auditory nerve spontaneous spiking patterns that provided an improved statistical match to empirical data. We modified a revised version of the Meddis (2006) model to develop a more biophysically detailed description of stochastic variability in the IHC-AN complex. It has been suggested (Morse and Evans, 1996;McDonnell et al., 2008) that significantly decreased stochastic variability in AN spiking generated by cochlear implants is a contributing factor to imperfect performance of these implants. A potential application of our model, therefore, is as a component in a larger model of the auditory system designed to predict differences in neural activity in higher brain regions, such as the cochlear nucleus, due to electrical stimulation by cochlear implants, in comparison with natural acoustic stimulation.
Based on our findings it will be interesting for future work to build on our study with a more detailed model of the calcium dynamics of the ribbon synapse in inner hair cells. Such a model might be capable of explaining both pre-synaptic shortterm depression in vesicle release, and long-term correlations due to calcium fluctuations.