At small areas the model fired spontaneously (**Figure 6A**), similar to the HH model above (**Figure 1**). Also, similar to before, our theory well predicts the variance of voltage fluctuations for large membrane areas, but diverges from the simulation

Spontaneous firing rate as function of membrane area for CA1 pyramidal cell model. Qualitative behavior is similar to Hodgkin–Huxley model, **Figure 1**. **(B)** Comparison of voltage noise variance from theory and simulation as a function of membrane area. As with Hodgkin–Huxley model (**Figure 4F**), theory is a good predictor for large membrane areas when there is no spontaneous spiking. **(C)** Contributions to voltage variance from the different ion channel types in the model, for both theory and simulation. Two types of potassium channels dominate the noise in this model. Membrane area of 1000µm2.

results for small areas < 500µm2 when the neuron spikes (**Figure 6B**).

In **Figure 6C** we plot the contributions to voltage variance from each of the three voltage-dependent ion channel types in the CA1 pyramidal neuron model. Similar to the HH model above, noise from K+ channels dominates over noise from Na+ channels. Of the two K+ channel types in the model, the delayed rectifier K+ channel contributes more voltage noise than the Atype K+ channel. Interestingly, Ka has a larger current noise variance than Kdr (not shown), but because Ka has faster kinetics than Kdr, its current noise is more heavily filtered by the membrane capacitance leading to a switch in the relative magnitudes of their contributions to voltage variance.

These results demonstrate two points: first, our method is readily applicable to any neuron model; and second, the dominant source of ion channel noise depends on the physiological details of the neuron. However, it should not be concluded from these results that K+ channels will always contribute more noise than Na+ channels. Substantial channel noise can arise from channels permeable to any ion: Na+, K+ or any other. The relative amplitudes and effects of channel noise simply depend on all of the earlier outlined factors and will need to be evaluated on a case-by-case basis.

## **DISCUSSION**

We have picked apart the various factors that determine a specific channel population's contribution to membrane noise. Although applied only to the Hodgkin Huxley and hippocampal pyramidal neuron models here, the method is applicable to any voltage gated channel model. In summary the factors are:

Another qualitative factor is the polarity of current flow. Open Na+ channels further depolarize the cell, hence increasing the probability for other Na+ channels to open and acting as a positive feedback loop. Hence, stochastic Na+ channels increase excitability of the cell through regenerative depolarizing excursions in membrane potential (Dudman and Nolan, 2009). Open K+ channels, in contrast, hyperpolarize the cell and act as negative feedback to changes in membrane potential. This negative feedback coupled with their relatively slow kinetics can, in some cases, enable stochastic K+channels to drive sub-threshold oscillations (Schneidman et al., 1998).

The combination of these factors yields an accurate prediction of the membrane voltage noise. While it is possible to obtain a coarse estimate of spontaneous firing rates, this is far from perfect and highlights two current hiatus in the theory, namely, the absence of an accurate phenomological model for channel noise when the noise can not be assumed to be small, and the lack of theory for colored noise driven spontaneous activity.

In the case of the stochastic HH model we have shown that the fluctuations from stochastic gating of potassium channels is the dominant source of noise by three different measures. First, a HH model where only K+ channels gate stochastically spontaneous fires at higher rates than a HH model where only Na+ channels gate stochastically (**Figure 1**) (Skaugen and Walløe, 1979; Schneidman et al., 1998; van Rossum et al., 2003). Second, examining the dynamics of Na+ and K+ currents in the milliseconds preceding a spontaneous action potential in the "all stochastic" HH model shows that, on average, spikes are generated by a drop in K+ current that precedes the increase in Na+ current (**Figure 2**). Third, direct calculation of the voltage noise spectra from each channel population at resting potential shows that K+ channel fluctuations contribute ∼ 75% of the total membrane noise (**Figure 4**). This finding, although consistent with results reported by Schneidman et al. (1998); van Rossum et al. (2003), is in contrast with other simulations (Chow and White, 1996; Faisal et al., 2005). We discuss these two studies separately.

Chow and White (1996) used approximate analytical methods to directly calculate the spontaneous firing rate in the stochastic HH model, and compared the predictions to numerical simulations to find apparently very good agreement. Our own simulations produce quantitatively similar results to their simulations (data not shown), so it is likely that their simulated data are correct. However, their analytical calculations were based on the assumption that spontaneous spiking is driven *solely* by stochastic activation of Na+ channels. No matter how elegant, their result can not be accurate as it ignores the K+ noise, which is the main cause for spontaneous firing (**Figure 1**). If anything, their analytical model should be a better approximation of our simulations when K+ channels are modeled deterministically. However, their calculations do not match this. The errors could have arisen in any of the multiple approximating steps necessary for their calculation. For example, they assume a static absolute voltage threshold when in reality the HH model has (1) no hard threshold for any type of stimulus (Izhikevich, 2007) and (2) different apparent spike thresholds for stimuli of different temporal structure (Koch, 1999). A more fruitful method for future studies could be to derive a higher-dimensional version of the spike threshold that incorporates both fast and slow channel variable states (Newby et al., 2013).

Faisal et al. (2005) find in cable axon HH models that Na+ channels contribute more to spontaneous spiking than K+ channels. We believe this to be a numerical simulation error as it is inconsistent with our simulations (when implemented with both PSICS and NEURON), and also those of Schneidman et al. (1998) and van Rossum et al. (2003)—implemented with the "NeuronC" simulator (Smith, 1992). Furthermore, we found the K+ channels to be dominant not only in single compartment models but in cable structures as well. Without access to their simulator, it is difficult to tell where the discrepancy lies. Nevertheless, it is possible that Na+ channel noise does drive spontaneous spiking in models other than the HH squid giant axon. Our general theory should help to quickly estimate such possibilities without resorting to full simulations.

#### **AUTHOR CONTRIBUTIONS**

Conception: Cian O'Donnell. Carried out analytical and computer calculations, wrote the paper: Cian O'Donnell and Mark C. W. van Rossum.

#### **FUNDING**

Cian O'Donnell was supported by EPSRC/BBSRC/MRC through the Doctoral Training Centre in Neuroinformatics.

#### **ACKNOWLEDGMENT**

Discussions with Matthew Nolan, Robert Cannon and Carson Chow are gratefully acknowledged.

#### **REFERENCES**

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

*Received: 14 March 2014; accepted: 17 August 2014; published online: 04 September 2014.*

*Citation: O'Donnell C and van Rossum MCW (2014) Systematic analysis of the contributions of stochastic voltage gated channels to neuronal noise. Front. Comput. Neurosci. 8:105. doi: 10.3389/fncom.2014.00105*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2014 O'Donnell and van Rossum. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Diffusion approximation-based simulation of stochastic ion channels: which method to use?

## *Danilo Pezo1, Daniel Soudry2 and Patricio Orio1,3\**

*1 Centro Interdisciplinario de Neurociencia de Valparaíso, Universidad de Valparaíso, Valparaíso, Chile*

*2 Department of Statistics and the Center for Theoretical Neuroscience at Columbia University, New York, NY, USA*

*3 Instituto de Neurociencia, Facultad de Ciencias, Universidad de Valparaíso, Valparaíso, Chile*

#### *Edited by:*

*Joshua H. Goldwyn, New York University, USA*

#### *Reviewed by:*

*Daniele Linaro, University of Antwerp, Belgium Roberto Fernández Galán, Case Western Reserve University, USA*

#### *\*Correspondence:*

*Patricio Orio, Centro Interdisciplinario de Neurociencia de Valparaíso, Pasaje Harrington 287, Valparaíso 2360103, Chile e-mail<EMAIL_ADDRESS>
To study the effects of stochastic ion channel fluctuations on neural dynamics, several numerical implementation methods have been proposed. Gillespie's method for Markov Chains (MC) simulation is highly accurate, yet it becomes computationally intensive in the regime of a high number of channels. Many recent works aim to speed simulation time using the Langevin-based Diffusion Approximation (DA). Under this common theoretical approach, each implementation differs in how it handles various numerical difficulties—such as bounding of state variables to [0,1]. Here we review and test a set of the most recently published DA implementations (Goldwyn et al., 2011; Linaro et al., 2011; Dangerfield et al., 2012; Orio and Soudry, 2012; Schmandt and Galán, 2012; Güler, 2013; Huang et al., 2013a), comparing all of them in a set of numerical simulations that assess numerical accuracy and computational efficiency on three different models: (1) the original Hodgkin and Huxley model, (2) a model with faster sodium channels, and (3) a multi-compartmental model inspired in granular cells. We conclude that for a low number of channels (usually below 1000 per simulated compartment) one should use MC—which is the fastest and most accurate method. For a high number of channels, we recommend using the method by Orio and Soudry (2012), possibly combined with the method by Schmandt and Galán (2012) for increased speed and slightly reduced accuracy. Consequently, MC modeling may be the best method for detailed multicompartment neuron models—in which a model neuron with many thousands of channels is segmented into many compartments with a few hundred channels.

**Keywords: channel noise, stochastic simulation, ion channel, Langevin, conductance based models**

## **INTRODUCTION**

Understanding the effect of stochastic phenomena on the behavior of the nervous system requires stochastic simulation algorithms that effectively and accurately capture the dynamics of the underlying modeled phenomena. Among the sources for variability, the stochastic opening and closing of ion channels has caught the attention of several works over the past years. The best description of stochastic gating of ion channels is attained with the use of continuous time, discrete states Markov Chain (MC) processes (Neher and Stevens, 1977; Colquhoun and Hawkes, 1981), however this approach can be very slow in simulations with a large number of channels.

As an alternative to the explicit MC simulation, the Diffusion Approximation (DA) calculates the trajectory of a population of independent MCs using a Stochastic Differential Equation (SDE), sometimes called the Chemical Langevin Equation (Gillespie, 2000, 2007). Its application to the simulation of stochastic ion channels was suggested almost 20 years ago (Fox and Lu, 1994; Fox, 1997), but in the beginning there were some errors in the application of the scheme. This led to the belief that the approximation was not good enough (Mino et al., 2002; Bruce, 2007, 2009). Later, revised implementations of the algorithms were published showing that indeed the DA can reproduce the statistical properties of a population of discrete ion channel fluctuating between open and closed states (Goldwyn et al., 2011; Goldwyn and Shea-Brown, 2011). Other works (Dangerfield et al., 2012; Orio and Soudry, 2012) also offered a simplified description of the algorithm, making it easy to apply to any given kinetic scheme.

What the SDE system approach does is to approximate the trajectory in time of the fraction of channels at every state. For the simulation to remain physically meaningful, none of the fractions can be negative or greater than 1. We call this the "boundary constraint." This constraint would occasionally break in the numerical simulation of the SDE, if we use its naively discretized form (the Euler–Maruyama method). This is because stochastic fluctuations can make the variables leave the [0,1] interval. The problem amplifies when the number of channels is low and stochastic fluctuations increase.

If the boundary constraint is not maintained, this can generate additional technical problems in the simulation. Specifically, the calculation of the stochastic terms involves the square root of a term involving the fraction variables, which may yield complex values when the variables are negative. Such complex values must be avoided. Therefore, we get another constraint, which we call the "Real-valued Square Root" (RSR) constraint. This constraint is automatically fulfilled if the boundary constraint is maintained.

Finally, the sum of the fractions over all states must be equal to 1 at all times. We call this the "normalization constraint." Although this constraint is supposed to be guaranteed by the continuous-time equations (see below), machine rounding errors of the discretized equations can gradually break it. Additionally, any method that deals with the boundary constraint must also take normalization into account. For example, a naive truncation of any variable that leaves the [0,1] interval would break the normalization constraint.

To address these issues, a number of improvements have been proposed to the DA schemes.

Orio and Soudry (2012) proposed to allow the variables to freely change, ignoring the boundary constraint. In order to take care of the normalization constraint one of the variables was replaced with one minus the sum of the others. Additionally, an absolute value operation was added in the stochastic terms to maintain the RSR constraint. Previously, Goldwyn et al. (2011) also allowed the variables to freely change, but instead used a steady state approximation on the voltage to maintain the RSR constraint. However, such an approximation can be rather inaccurate (Dangerfield et al., 2012; Orio and Soudry, 2012; Huang et al., 2013a), even when the number of channels is relatively high.

Two methods aim to maintain the boundary constraint. Dangerfield et al. (2012) proposed that if variables break either the boundary or normalization constraint, they are "reflected" back into the valid region, in which all the constraints are kept. This is done using projection into a simplex (Chen and Ye, 2011). Huang et al. (2013a) proposed a different method. When the boundary constraint is broken, the variables should first be truncated. This breaks normalization, so the variables are then renormalized. Finally, in the next time step, the variables are incremented with the remainders from the truncation in the previous steps. As the boundary constraint is almost constantly broken, normalization is continuously corrected in both Dangerfield et al. (2012) and Huang et al. (2013a).

Other methods have been proposed with a different goal in mind. Schmandt and Galán (2012) aimed to reduce computational complexity and speeding up the simulation. They proposed to neglect stochastic noise terms in all state transitions, except those connecting the open state (or states), an approximation they call "Stochastic shielding." Güler (2013) introduced a stochastic HH model with colored noise in the conductance terms as well as in the current terms in order to capture the non-trivial crosscorrelation between the transmembrane voltage fluctuation and the component of open channel fluctuation attributed to multiple number of gates in individual ion channels. Another recently published method (Linaro et al., 2011), also used colored noise in the current terms (but not in the conductance terms). However, Linaro's method will not be examined here, because it involves a steady-state approximation in the stochastic terms (similarly to Goldwyn et al., 2011), which was shown to introduce inaccuracies (Orio and Soudry, 2012).

Despite the improvement in accuracy or computational efficiency that the new methods represent for the simulation of stochastic ion channel activity, we were concerned about the comparisons performed between them and the real benefit of implementing the numerical algorithms.

First, there is the computational cost issue. The initial motivation for developing DA methods was to make stochastic simulations faster than MC modeling. Therefore, if the extra computation needed to normalize and bound the variables makes it slower than MC then the purpose is defeated. Moreover, we already noticed that when the number of channels is low (when DA becomes more inaccurate) or with very small integration times, MC modeling can run faster than DA (see Figure 7 in Orio and Soudry, 2012). This, added to the fact that bounding and normalization of the DA requires more coding (and eventually, debugging), can render DA less attractive.

Second, we noted that the standard test employed to prove the accuracy of numerical methods for stochastic ion channels is the original Hodgkin and Huxley (HH) model. This was the only model used for testing in most previous papers, including a recent review (Rowat and Greenwood, 2014). This model, as standard and general as it is, reproduces the kinetics of ion channels of the squid axon at 6.3◦C, thus differing greatly from mammalian nerve excitable membranes. This difference can be very significant, as we noted (Orio and Soudry, 2012). There we found that the use of a steady-state approximation in the stochastic terms usually does not introduce severe inaccuracies in the context of the original HH model. However, deviations were detected in common current clamp-based simulations when the steady-state approximation is used in a model inspired in mammalian (therefore faster) ion channels. It is noteworthy that the difference between mammalian inspired and the squid axon model relies only in the parameters that describe the transition rate constants (and thus the time scale of the model), while the equations and the model framework are identical.

Thus, we see a necessity for testing the DA with and without the recently proposed corrections in a wider spectrum of simulation scenarios and taking into account other variables than simulation accuracy, namely:

In an attempt to test the real usability of the algorithms in the context of more complex neuronal models, we implemented them in one of the choice tools for biophysically-inspired modeling, the Neuron simulation environment (Hines and Carnevale, 2001; Carnevale and Hines, 2006). We conducted both singlecompartment and multi-compartment simulations using MC or DA algorithms and compared their performance as well as the ability of different DA implementations to reproduce the variability introduced by MC modeling.

In our results all DA algorithms deviate to some degree from the MC modeling when the number of channels falls below 1000, regardless of the attempts to deal with normalization and bounding of the variables. However, we see that in this condition MC modeling runs, in most of the scenarios, faster than DA implementations. Therefore, one of the most common motivations to use DA, which is to achieve faster computation times, is not accomplished when the number of channels is low. However, when the number of channels is high, DA algorithms can accurately reproduce MC, with improved speed. Specifically, in this regime, no inaccuracy was detected in both Orio and Soudry (2012) and Huang et al. (2013a); Schmandt and Galán (2012) was slightly inaccurate; Güler (2013) was somewhat inaccurate; and Dangerfield et al. (2012) was the least accurate. In terms of computational speed, the ranking is as follows (see **Figure 8**): (1) Stochastic Shielding (Schmandt and Galán, 2012) (2) Colored Noise (Güler, 2013) (3) Unbound DA (Orio and Soudry, 2012) (4) Reflected DA (Dangerfield et al., 2012) (5) Truncated and Restored DA (Huang et al., 2013a).

## **MATERIALS AND METHODS**

## **SIMULATIONS: MODELS EMPLOYED AND TESTS PERFORMED** *Original Hodgkin and Huxley model*

The original Hodgkin and Huxley (HH) model (Hodgkin and Huxley, 1952) was simulated with the equation:

$$C\_m \frac{dV(t)}{dt} = -\mathbf{g}\_{\rm Na}(t) \left( V(t) - E\_{\rm Na} \right) - \mathbf{g}\_{\rm K}(t) \left( V(t) - E\_{\rm K} \right)$$

$$- \lg \left( V(t) - E\_l \right) + I\_{\rm sim}(t) \tag{1}$$

With the exception of Güler's colored noise algorithm, sodium and potassium channels were treated as 8- and 5-state MCs, respectively. The corresponding kinetic schemes are:

$$\begin{array}{ccccccccc} m\_0 h\_0 & \stackrel{3\alpha\_m}{\underset{\beta\_m}{\rightleftharpoons}} & m\_1 h\_0 & \stackrel{2\alpha\_m}{\underset{2\beta\_m}{\rightleftharpoons}} & m\_2 h\_0 & \stackrel{\alpha\_m}{\underset{3\beta\_m}{\rightleftharpoons}} & m\_3 h\_0\\ \beta\_h \;\;\;\mid \;\alpha\_h & \beta\_h \;\;\mid \;\alpha\_h & \beta\_h \;\mid \mid \;\alpha\_h & \beta\_h \;\mid \mid \;\alpha\_h \; (\text{scheme1})\\ m\_0 h\_0 & \stackrel{3\alpha\_m}{\underset{\beta\_m}{\rightleftharpoons}} & m\_1 h\_0 & \stackrel{2\alpha\_m}{\underset{2\beta\_m}{\rightleftharpoons}} & m\_2 h\_0 & \stackrel{\alpha\_m}{\underset{3\beta\_m}{\rightleftharpoons}} & m\_3 h\_0\\ m\_0 & \stackrel{4\alpha\_n}{\underset{\beta\_n}{\rightleftharpoons}} & n\_1 & \stackrel{3\alpha\_n}{\underset{2\beta\_n}{\rightleftharpoons}} & n\_2 & \stackrel{2\alpha\_n}{\underset{4\beta\_n}{\rightleftharpoons}} & n\_3 & \stackrel{\alpha\_n}{\underset{4\beta\_n}{\rightleftharpoons}} & n\_4 \text{ (scheme2)} \end{array}$$

Sodium and potassium conductances at time *t* (*gNa*(*t*) and *gK*(*t*)) were calculated as the fraction of channels in the conducting states *m*3*h*1 and *n*4 multiplied by the maximum conductances *gNa* and *gK*, respectively. The kinetic rates α and β are given (in ms−1) by:

$$\alpha\_m(V) = \frac{0.1\left(V + 40\right)}{1 - \exp\left(-\frac{V + 40}{10}\right)}; \quad \beta\_m(V) = 4\exp\left(-\frac{V + 65}{18}\right)$$

$$\alpha\_h(V) = 0.07\exp\left(-\frac{V + 65}{20}\right); \beta\_h(V) = \frac{1}{1 + \exp\left(-\frac{V + 35}{10}\right)}$$

$$\alpha\_n(V) = \frac{0.01\left(V + 55\right)}{1 - \exp\left(-\frac{(V + 55)}{10}\right)}; \ \beta\_n(V) = 0.125\exp\left(-\frac{V + 65}{80}\right).$$

where the terms were corrected to adjust the resting potential to −65 mV. Correspondingly, the rest of parameters are: *Cm* = 1μF/cm2, *gNa* = 120mS/cm2, *gK* = 36mS/cm2, *gl* = 0.1mS/cm2, *ENa* = 50mV, *EK* = −77 mV, *El* = −54.3 mV.

With this model, the following tests were conducted:

## *Schmidt-Hieber and Bischofberger model—single compartment*

The Schmidt-Hieber and Bischofberger (SB) model was proposed after the characterization of sodium channels both at the soma and at the axon initial segment of granule cells of the hippocampus (Schmidt-Hieber and Bischofberger, 2010). Sodium channels are described by the following kinetic scheme:

$$\begin{array}{ccccccccc} m\_0 h\_0 & \stackrel{\alpha\_1}{\underset{\beta\_1}{\rightleftharpoons}} & m\_1 h\_0 & \stackrel{\alpha\_2}{\underset{\beta\_2}{\rightleftharpoons}} & m\_2 h\_0 & \stackrel{\alpha\_3}{\underset{\beta\_3}{\rightleftharpoons}} & m\_3 h\_0\\ \beta\_h \mid \mid \alpha\_h & \beta\_h \mid \mid \alpha\_h & \beta\_h \mid \mid \alpha\_h & \beta\_h \mid \mid \alpha\_h \text{ (scheme3)}\\ m\_0 h\_0 & \stackrel{\alpha\_1}{\underset{\beta\_1}{\rightleftharpoons}} & m\_1 h\_0 & \stackrel{\alpha\_2}{\underset{\beta\_2}{\rightleftharpoons}} & m\_2 h\_0 & \stackrel{\alpha\_3}{\underset{\beta\_3}{\rightleftharpoons}} & m\_3 h\_0 \end{array}$$

where the kinetic rates are given by

$$\alpha\_i(V) = \alpha\_{i,0} \exp\left(\alpha\_{i,1} V\right)$$

$$\beta\_i(V) = \beta\_{i,0} \exp\left(-\beta\_{i,1} V\right)$$

$$\alpha\_h(V) = \frac{\alpha\_{h,0}}{1 + \alpha\_{h,1} \exp\left(\alpha\_{h,2} V\right)}$$

$$\beta\_h(V) = \frac{\beta\_{h,0}}{1 + \beta\_{h,1} \exp\left(\beta\_{h,2} V\right)}$$

Parameters α*i*,*j* and β*i*,*j* are given in **Table 1** for both axonal and somatic channels. In single-compartment simulations, somatic parameters were employed.

**Table 1 | Activation parameters for somatic and axonal sodium channels in the Schmidt-Hieber and Bischofberger model.**

Potassium channels in SB model are simulated by the same kinetic scheme as HH model (scheme 2) with the following voltage dependent transition rates:

$$\alpha\_n(V) = s\varepsilon \times 0.01 \frac{V + 55}{1 - \exp\left(-\frac{V + 55}{10}\right)}$$

$$\beta\_n(V) = s\varepsilon \times 0.125 \exp\left(-\frac{V + 65}{80}\right)$$

where *sc* is a scale parameter that adjusts the kinetic constants according to the compartments in which the channels are being simulated. For single compartment simulations, equation (1) was used with the following parameters: *Cm* = 1μF/cm2, *gNa* = 20 mS/cm2, *gK* = 4 mS/cm2 , *gl* = 0.1 mS/cm2, *ENa* = 75 mV, *EK* = −95mV, *El* = −70 mV. With this model, the number of channels was controlled by the membrane area, given a unitary conductance of 20 pS/cm2. We tested areas ranging from 15.7 to 628 μm2, resulting in 157 to 6283 sodium channels and 31 to 1257 potassium channels.

With this model, the following tests were conducted at different values of membrane area:

#### *Schmidt-Hieber and Bischofberger model—idealized multicompartment model*

An idealized model similar to the one described in Schmidt-Hieber and Bischofberger (2010, see **Figure 5A**) was simulated using stochastic algorithms for the ion channels. The parameters for ion channel densities and kinetic constants were used as described in the article with minor modifications, such as the absence of an axonal bleb and a longer axon for some simulations.

#### **SIMULATION ALGORITHMS**

The methods described and tested here are designed to simulate a number of independent and identical Markov Chains (MCs) with a discrete number of states, keeping track of the number of channels in any state at any given time. For the description of the algorithms, we denote S as the total number of states in a MC, *i* ∈ {1,..., *S*} are the individual states, and Ni is the number of MCs in state i.

#### *Markov chain simulations (MC)*

Markov chains were simulated using the Stochastic Simulation Algorithm (SSA) (Gillespie, 1976) with some modifications. Briefly, the method consists in:

At time *t*, calculate the effective transition rate λ(*t*) as

$$\lambda(t) = \sum\_{i}^{s} N\_i(t)\xi\_i(t)$$

where ζ*i*(*t*) is the sum of transition rates for transitions escaping from state *i*.

Calculate the time for the next transition *tn* as

$$t\_n = t\_{\mathcal{P}} - \frac{\log \xi\_1}{\lambda(t)}$$

where *tp* is the time of the previous transition (0 at the beginning of the simulation) and ξ1 is a random number uniformly distributed within [0,1], drawn after the previous transition.

If *tn* > *t*, continue integrating the time and the membrane voltage equation.

If *tn* < *t*, perform a transition:

Calculate the probability of all *j* transitions:

$$P\_{\dot{j}}(t) = \frac{N\_{\dot{i}}(t)\alpha\_{\dot{j}}(t)}{\sum\_{j} N\_{\dot{i}}(t)\alpha\_{\dot{j}}(t)}$$

where *i* is the state originating transition *j* and α*j* its rate.

Build a cumulative sum of all transition probabilities. Draw a random number ξ2 uniformly distributed in [0,1] and find the first term in the cumulative probability that is greater than ξ2.

Execute the transition indicated by the term found in the previous step, and draw a new random number ξ1 to be used for the time of the next transition.

## *Unbounded diffusion approximation (UA) (Orio and Soudry, 2012)*

The DA algorithm was implemented with SDEs described previously (Orio and Soudry, 2012; see also Mélykúti et al., 2010). In matrix form, the equations for the sodium and potassium channels are, respectively:

$$\frac{dX\_{\rm Na}}{dt} = A\_{\rm Na} X\_{\rm Na} + \frac{1}{\sqrt{N\_{\rm Na}}} S\_{\rm Na}(X\_{\rm Na}) \xi(t)\_{\rm Na}$$

$$\frac{dX\_K}{dt} = A\_K X\_K + \frac{1}{\sqrt{N\_K}} S\_K(X\_K) \xi(t)\_K \tag{2}$$

where *XNa* = *m*0*h*0 *m*1*h*0 *m*2*h*0 *m*3*h*0 *m*0*h*1 *m*1*h*1 *m*2*h*1 *m*3*h*1 *T* and *XK* = *n*0 *n*1 *n*2 *n*3 *n*4 *T* are column vectors with the fraction of channels at any given state, and ξ (*t*)*Na* and ξ (*t*)*K* are column vectors of independent normally distributed random variables (mean 0, variance 1) with length 10 and 4, respectively. *NNa* and *NK* are the number of sodium and potassium channels, respectively. The rate matrices *ANa* and *AK* and square root matrices *SNa*(*XNa*) and *SK*(*XK*) can be directly found from the state diagram of the corresponding ion channel type. This is explained in detail around equations 1 and 13 and the Supplemental Material in Orio and Soudry (2012). For example, in the case of Potassium channels we have

$$A\_K = \begin{bmatrix} -4\alpha\_n & \beta\_n & 0 & 0 & 0\\ 4\alpha\_n & -3\alpha\_n - \beta\_n & 2\beta\_n & 0 & 0\\ 0 & 3\alpha\_n & -2\alpha\_n - 2\beta\_n & 3\beta\_n & 0\\ 0 & 0 & 2\alpha\_n & -\alpha\_n - 3\beta\_n & 4\beta\_n\\ 0 & 0 & 0 & \alpha\_n & -4\beta\_n \end{bmatrix}$$

and

$$\begin{array}{c} \text{SK } (\text{X} \boldsymbol{\chi}) = \begin{bmatrix} \sqrt{4\alpha\_n n\_0 + \beta\_n n\_1} & 0\\ -\sqrt{4\alpha\_n n\_0 + \beta\_n n\_1} & \sqrt{3\alpha\_n n\_1 + 2\beta\_n n\_2} \\ 0 & -\sqrt{3\alpha\_n n\_1 + 2\beta\_n n\_2} \end{bmatrix} \\ \begin{bmatrix} 0 & 0\\ 0 & 0 \end{bmatrix} \\ \sqrt{2\alpha\_n n\_2 + 3\beta\_n n\_3} & 0 \\ -\sqrt{2\alpha\_n n\_2 + 3\beta\_n n\_3} & \sqrt{\alpha\_n n\_3 + 4\beta\_n n\_4} \\ 0 & -\sqrt{\alpha\_n n\_3 + 4\beta\_n n\_4} \end{bmatrix} . \end{array}$$

Note that in the case of the Schmidt-Hieber and Bischofberger model, the only difference from the HH model is that the rate constants for the sodium channel equations are different.

To take care of normalization, variables *m*1*h*0 ... *m*3*h*1 and *n*1 ... *n*4 were advanced by an Euler-Maruyama scheme and the remaining two were calculated as *m*0*h*0 = 1 − *m*1*h*0 − *m*2*h*0 − *m*3*h*0 − *m*0*h*1 − *m*1*h*1 − *m*2*h*1 − *m*3*h*1 and *n*0 = 1 − *n*1 − *n*2 − *n*3 − *n*4. As we do not control the bounding of the variables between 0 and 1, in order to ensure real-valued square roots we calculated the stochastic terms *S*(*X*) (and only those terms) taking the absolute value of the variables. Thus, we refer to this algorithm as UA – Unbounded with Absolute values in stochastic terms.

## *Reflected SDEs (Ref) (Dangerfield et al., 2012)*

This method aims to normalize the variables *m*1*h*0 ... *m*3*h*1 and *n*1 ... *n*4 and to keep them bounded in the interval [0, 1] using the reflected stochastic equation approach, described in section V of Dangerfield et al. (2012). We have used this method together with the DA equation system (Equations 2).

## *Truncated and restored DA (TR) (Huang et al., 2013a)*

We used the DA equations (Equations 2) with an additional residual term. This residual term was introduced to ensure the boundary and normalization constraints, as explained in section II.F in Huang et al. (2013a).

## *Stochastic shielding approximation (SSmc) (Schmandt and Galan, 2012)*

In this method, transitions not connecting to the conducting states are approximated to be deterministic and solved as ODEs. Transitions connecting conducting with non-conducting states are solved as Markov Chains with the already mentioned algorithm. As there is a mixture of continuous (ODEs) and discrete (MC) treatment of variables, violations of the constraints occur. In our implementation, and inspired by the code by Schmandt and Galan, variables going off the [0,N] boundary are manually corrected and normalization was performed as in the UA algorithm.

## *Stochastic shielding approximation with DA (SSda)*

We modified the Schmandt and Galán (2012) approach by calculating the stochastic transitions with a DA approach (Equations 2) rather than using MCs. Therefore, SSda similarly uses the same DA equations, but the stochastic terms related to transitions not connecting to the conducting states were neglected.

For example, for potassium channels, we now used

$$S\_K(X\_K) = \begin{bmatrix} 0 \ 0 \ 0 & 0 \\ 0 \ 0 \ 0 & 0 \\ 0 \ 0 \ 0 & 0 \\ 0 \ 0 \ 0 & \sqrt{\alpha\_n n\_3 + 4\beta\_n n\_4} \\ 0 \ 0 \ 0 & -\sqrt{\alpha\_n n\_3 + 4\beta\_n n\_4} \end{bmatrix},$$

Thus, we needed only 2 Brownian terms for sodium channels and 1 for potassium channels. Boundary and normalization constraints were again treated as in the UA algorithm, that is: there was no bounding of the variables; the absolute value of the variables was used in the square roots of the stochastic terms; and normalization was applied by calculating *m*0*h*0 = 1 − *m*1*h*0 − *m*2*h*0 − *m*3*h*0 − *m*0*h*1 − *m*1*h*1 − *m*2*h*1 − *m*3*h*1 and *n*0 = 1 − *n*1 − *n*2 − *n*3 − *n*4.

#### *HH with colored noise terms (CN) (Güler, 2013)*

We simulated the Güler's Brownian harmonic oscillator using the system of stochastic differential equations 6.1, 6.2 in Güler (2013), using the constant parameters given in Table 2 in Güler (2013). Note that these equations are very different from the standard DA equations.

To take care of the normalization constraint at each time step, if any of the variables *m*, *h* or *n*, left the [0,1] interval, then the stochastic term η was redrawn until the variable fulfilled the boundary constraint.

## **SOFTWARE**

All the models and simulations algorithms presented here were implemented and run in the Neuron simulation environment (Hines and Carnevale, 2001; Carnevale and Hines, 2006). The different algorithms were written inside the MOD files for each channel. With the exception of TR, all algorithms run fine regardless of the numeric integrator specified (cnexp or euler). They also produce the same results in Python using the Euler-Maruyama integration method (tested in some selected cases). The TR algorithm required the specification of the Euler integrator within the MOD file to produce the results presented here, otherwise a much lower firing rate was obtained. Simulation control and the recording of variables were specified with Python scripts (Hines et al., 2009). Sample codes and .mod files can be found in ModelDB http://senselab.med.yale.edu/ModelDB/ Accession 167772 .

Data analysis and plotting was performed using the Python libraries numpy, scipy, and matplotlib.

## **RESULTS**

To test the accuracy of the methods we performed a series of simulations, comparing the variability of the results to that obtained with explicit MCs solved by the exact Gillespie algorithm (Gillespie, 1976, 2007). We employed the original Hodgkin and Huxley model (Hodgkin and Huxley, 1952), in order to reproduce previously published comparisons of the algorithms. In addition we performed some tests with a faster model, based on sodium channels from granular cells in the hippocampus (Schmidt-Hieber and Bischofberger, 2010). Finally, we simulated a model neuron with multiple compartments (Schmidt-Hieber and Bischofberger, 2010) and measured the variability in the generation *and conduction* of action potentials.

## **SIMULATIONS WITH THE ORIGINAL HODGKIN AND HUXLEY MODEL** *Firing variability—15 ms simulation with stimulus*

A widely used test to compare stochastic simulation algorithms (Mino et al., 2002; Bruce, 2007; Orio and Soudry, 2012 and others) consists of a short simulation (15 ms) in which a 2-ms current stimulus is given after a 1-ms delay (**Figure 1A**). Depending on the amplitude of the stimulus the probability of eliciting an action potential increases, and this relationship depends on the number of channels. **Figure 1B** shows the probability of firing an action potential in 10,000 trials at different stimulus amplitudes, for the algorithms tested with *NNa* = 5000. **Figure 1C** plots the variance of action potential timing, a measure of jitter. The Reflection method produces a higher firing probability at all stimulus amplitudes. This entails a lower variability in action potential timing. Additionally, Güler's CN method produces a higher variability than MC and other methods. This result is repeated with higher number of channels, however at a lower number of channels the spontaneous firing of action potentials makes the comparison unreliable. To compare the behavior of the models with *NNa* ≤ 1600, we modified the protocol so the stimulus is sustained during the simulation (**Figure 1D**, inset) and explored negative values of current amplitude. As can be seen in **Figures 1D,E**, the behavior of all the DA algorithms (as well as of SSmc) deviates from MC considerably for *NNa* = 50 and to a minor degree for *NNa* = 500. Again, the Reflection method produces a higher firing probability than the other methods.

## *Spontaneous firing rate*

The original HH model with stochastic ion channel produces spontaneous firing activity that increases as the number of channels is decreased. With each simulation algorithm and with sodium channel number (N*Na*) ranging from 50 to 50,000, we simulated 500 s and recorded the occurrence of action potentials. **Figure 1F** shows the mean frequency of spikes that were detected in the simulations. With *NNa* ≥ 1600, almost all the methods reproduce the behavior of MC modeling. The sole exception is the Reflection method, that showed a higher firing probability at *NNa* = 1600 and below. To discard some incompatibility of the Neuron simulation environment with the reflection procedure, we repeated this simulation using an Euler-Maruyama integration procedure written in Python and obtained the same result. Huang's truncated and restored DA method seems to be the one that more closely follows MC modeling at extremely low number of channels, only slightly underestimating the firing rate. At *NNa* = 160, the Unbound DA and Stochastic Shielding methods overestimate the firing rate, dropping abruptly when *NNa* = 50. This latter behavior is actually due to numeric overflows that made the simulations run without producing action potentials. This was corrected using a smaller time step (**Figure 1F**, segmented lines). Thus, the Unbound DA (note that SSda is also unbound) becomes numerically unstable when the number of channels is too low.

#### *Voltage clamp—noisy voltage trace with action potential*

A third test, to check how the different DA methods can reproduce the variability of channel openings obtained with MC modeling, consists on recording the response of the model channels to a fixed voltage trajectory obtained from a stochastic simulation. The voltage trace is shown in **Figure 2A** and it contains an action potential as well as a noisy background (zoomed in **Figure 2B**). With each model and condition, 2000 independent simulations were run and the time evolution of open channels was recorded. At each point in time, the mean and variance of the open channels was calculated. In addition to the comparison with the behavior of MC simulation, we compared to the expected mean of open channels which is calculated by applying the same voltage clamp simulation to deterministic HH channels. Moreover, we can compute the expected variance as explained in Methods.

The results are shown in **Figure 2**. During the subthreshold regime, Reflection method overestimates the mean of open channels, both for sodium and potassium (**Figures 2C,D**, left). Huang's TR algorithm also overestimates the mean of open sodium channels (**Figure 2C**, left) to a minor extent. However, during the action potential any difference between the DA methods and the MC modeling or the exact solution appears to be negligible (**Figures 2C,D**, right). Regarding the variance of

pulse, for the same simulations shown in **(B)**. **(D,E)** Firing efficiencies obtained with a constant pulse (inset) of the indicated amplitudes with *NNa* = 50, *NK* = 15 **(D)** and *NNa* = 500, *NK* = 150 **(E)**. **(F)** Mean number

SSmc, stochastic shielding approximation (Schmandt and Galán, 2012); SSda, stochastic shielding with DA approximation; CN, colored noise (Güler, 2013).

the open channels (**Figures 2E,F**), the main deviation seems to occur with Güler's CN algorithm, which overestimates the variance of both open sodium and open potassium channels, during the subthreshold regime (left) and the action potential (right). Schmandt's stochastic shielding approximation (both SSmc and SSda) underestimates the variance of open channels during the action potential, when the voltage changes more rapidly. The results shown here are for *NNa* = 500, *NK* = 160; with higher number of channels (*NNa* = 5000, *NK* = 1600) we naturally observed less fluctuations but the results maintained: Reflected DA overestimates the mean of open channels in subthreshold regime; Colored Noise overestimates the variance of open channels in subthreshold regime; and Stochastic Shielding underestimates the variance of open channels during the action potential (not shown).

#### **SCHMIDT-HIEBER AND BISCHOFBERGER MODEL—SINGLE COMPARTMENT**

We decided to use a model with faster sodium channels, resembling mammalian ion channels, to test the accuracy of the DA methods when the transitions between states occur at faster rates. We chose a recently published model that focuses on the fast opening of sodium channel in the axon initial segment of granule cells from the hippocampus (Schmidt-Hieber and Bischofberger, 2010). We will refer to this model as the "SB" model. We noted that the SB model does not show spontaneous firing when simulated stochastically. Therefore, the 500-ms simulation test was not performed.

#### *Firing variability—20 ms with 1 ms stimulus*

We performed the test in which 20-ms were simulated with a 1-ms current stimulus (similar to **Figure 1A**). For each stimulus amplitude, 2000 simulations were run and the Firing Efficiency, mean firing time, and firing time variance were calculated. **Figure 3** shows that the results are similar to that obtained with HH. At high number of channels, most methods perform reasonably similar to MC with a higher excitability of the Reflection method (**Figures 3A,B**). However, when the number of channels is low all DA methods fail to approximate the results of MC, showing a much higher probability of firing at all amplitudes of the stimulus

scales. **(C)** Mean of open sodium channels during the subthreshold (left) and action potential (right) regimes, for the simulation algorithms tested. **(D)** Mean of open potassium channels. **(E)** Variance of open sodium channels. **(F)** Variance of open potassium channels. *NNa* = 500, *NK* = 160. Very similar results were obtained with *NNa* = 5000, *NK* = 1600 (see text).

(**Figure 3C**). Also, the Reflection method shows a lower firing time variability (**Figure 3D**).

#### *Voltage clamp*

We applied the voltage clamp test with the same voltage trace as the HH model to the SB model stochastic channels. Similar to what we observed for HH model, the Reflection method overestimates the mean of open channels (both sodium and potassium) during the subthreshold regime. Huang's truncated and restored method also overestimates it to a minor degree. During the action potential, the greatest deviation occurs with the stochastic shielding approximation, which underestimates the variance for both channels. With a higher number of channels, we observed similar results, with the exception of Huang's TR method performing better in the mean of open channels (not shown).

## **SCHMIDT-HIEBER AND BISCHOFBERGER MODEL—MULTI-COMPARTMENT SIMULATIONS**

To test the applicability of DA methods in more complex simulations of physiological relevance, we set up a multi-compartmental model of a neuron. We chose the idealized neuron described in Schmidt-Hieber and Bischofberger (2010) and shown in **Figure 4A**. Moreover, we kept the particular inhomogeneous sodium channel density for the axon that causes the action potentials to be initiated in the axon initial segment (AIS), about 10 μm from the soma (**Figure 4B**, bottom). The neuronal sections were spatially discretized according to their spatial constant λ, with a further increase in the number of segments in the AIS area. In total, the number of segments simulated were 895 with a 1500 μm axon and 2239 when the axon was extended to 7500μm. The distributions of segment areas is shown in **Figure 4C**. Together with the different ion channel densities and a unit conductance of 20 pS, the resulting distributions of number of channels per segment are shown in **Figures 4D,E**. It is noteworthy that the number of channels to be simulated in any given segment is rarely higher than 500 for sodium channels and never higher than 120 for potassium channels. In the model with the long axon, this adds up to 171189 sodium channels and 50670 potassium channels.

The model neuron with different stochastic channels was subject to a current clamp stimulus applied to the soma. The stimulus consisted in a 2 s noisy stimulus (**Figure 5A**) which in a deterministic simulation elicited 8 action potentials (represented by stars in the Figure). 400 independent simulations were performed with the same stimulus and **Figures 5B,C** show the raster plots (100 simulations) of spikes detected at the soma and at the tip of the axon, respectively. **Figure 5D** depicts a normalized firing probability calculated for the spikes at the tip of the axon. Both raster and firing probability plots show that the simulations with the Reflection method displayed a greater excitability, as several action potentials were only elicited with this algorithm and were not seen with the other DA methods or were seen with a much lower probability (i.e., around *t* = 1300 ms, *t* = 1700 ms and near the end of the trace). On the other hand, Stochastic Shielding with MC produced a lower excitability, firing near half of the action potentials per sweep than the other methods. The mean of spikes per sweep (**Figure 6A**) was significantly different to MC for all the algorithms, not only the Reflection and SSmc methods. With the Reflection method, however, the deviations from the other methods go beyond a higher excitability. Some spikes fired with high probability with all DA methods except for Reflection (see for instance around 800 ms), and some spikes had a slightly different timing with Reflection (700, 1500 ms). Therefore, the Reflection method in this test actually introduced a bias, producing spikes with different timings than the other DA method. The variability of the number of spikes elicited per trial varied with some DA algorithms compared to MC (**Figure 6B**) but only in the case of UA and SSmc a significant difference was observed.

We took the raster plots and searched for spikes that were repeated in at least 50% of the sweeps with the same timing ± 5 ms. These were called "common spikes" and for most algorithms 10–11 common spikes were found (**Figure 6C**), with the exception of Reflection and SSmc methods. Then, we measured how variable the timing of these spikes was at two axon locations, one near the initiation site and the other at the tip (**Figure 6D**). Although the TR method reproduced more closely the variability obtained with MCs, none of the observed differences was significant.

As a measure of variability with functional consequences, we examined how the duration of the action potential (duration measured at the detection threshold level of 0 mV) evolved as it propagates along the axon. The relevance of this measure is that the duration of action potentials at the release zone of a synapse will impact the amount of neurotransmitter released. **Figure 7A** plots the duration of all the action potentials recorded in a MC simulation at several sites of the axon, plotted against the duration at a site near the initiation. As a first observation, in the reference site there is a wide distribution of action potential durations, which gets narrower as the measurement site moves along the axon. Also, action potentials are shorter in the distal axon than in the initiation site and the Reflection method produces the shorter action potentials of all the methods (**Figure 7B**).

Besides getting shorter, action potential duration at the distal axon is completely uncorrelated to the duration at the initiation site (**Figure 7A**, right). We looked at how the correlation of action potential duration decays along the axon with the different stochastic simulation algorithms. Results are shown in **Figure 7C**. It is apparent that Reflection method produces a faster decay in the correlation, while the Unbound DA produces a longer propagation of correlation. To test for similarity, we fitted an exponential decay to the data points, obtaining a space constant λ. An extra sum of squares *F* test was performed to test the null hypothesis that the data points of each set could be fitted with the same λ as the MC data, showing that Unbound DA, Reflection and SSda methods produced a behavior significantly different to that obtained with MC (**Figure 7C**, inset). When the test was repeated with a fit to a double exponential decay, the same result was obtained.

## **COMPUTATION TIME**

To account for the usefulness of the simulation algorithms, we found important to compare the computational cost of each of them. **Figures 8A,B** show the time required to simulate 500 ms of the HH model with an integration time step (*dt*) of 5 and 0.5μs, respectively, for each algorithm used. As we reported previously (Orio and Soudry, 2012), DA methods are highly sensitive to *dt* but mostly insensitive to the number of channels. On the other hand, MC simulations are sensitive to both, but its sensitivity to the number of channels approaches a linear relationship as the number of channel increases. Importantly, as both *dt* and/or number of channel decrease, MC outperforms all of the DA methods, giving the best simulation times precisely in the condition where the DA methods showed to be more problematic.

In the case of the multi-compartment simulation, MC is faster than all DA methods except for Schmandt's stochastic shielding approximation (**Figure 8C**). This result is not surprising, given that the number of channels per simulated compartment was almost always below 1000 (**Figures 6C,D**).

## **DISCUSSION**

In this work we numerically tested five Diffusion Approximation algorithms proposed to reproduce the behavior of a number of simultaneous Markov Chains, in the context of stochastic ion channel modeling. Most of these algorithms (Dangerfield et al., 2012; Orio and Soudry, 2012; Schmandt and Galán, 2012; Huang et al., 2013a) are based on a Langevin Equation proposed for stochastic modeling of the Hodgkin and Huxley model (Fox and Lu, 1994; Fox, 1997). However, they deal in different ways with numerical issues that appear in the simulation of stochastic trajectories: the requirement that the variables be bounded in [0,1] ("boundary constraint") and the requirement that at any given time the sum of variables must be equal to 1 ("normalization constraint"). The boundary constraint breaks more frequently when the number of channels being simulated is low or when the integration time step increases, because the stochastic terms scale with *dt*/*N*. Therefore, special attention should be put to the performance of these simulation algorithms—and their capacity of reproducing MC behavior—with a low number of channels.

are elicited in a deterministic simulation.**(B)** Raster plots of action potentials detected at the soma during 100 of the simulations performed

Our tests were primarily aimed with a practical question in mind: what method should be used to study the contribution of channel stochasticity to neural excitability in any given specific context? To answer this question, the simple test of the numerical accuracy is not enough. Modeling algorithms should also be examined for applicability and simplicity of implementation in different contexts and for the computational efficiency for the intended model to be solved.

potential detected in a bin. The resulting vectors were then smoothed by

convolving with a Blackmann filter function of length 20.

### **ACCURACY**

Regarding accuracy, in brief we found that all DA algorithms fail in the reproduction of MC behavior when the number of channels is low (generally speaking, below 1000), with no clear "winner." For example, the TR procedure improved the results in some current clamp simulations (**Figures 1D–F**), but introduced higher deviations than the UA in others (**Figure 3**). The Reflection method in our simulations performed the worse, introducing a higher firing probability in current clamp simulations. Most likely, this is related to the higher mean of open sodium channels observed in voltage clamp (**Figure 2C**).

When the number of channels is 5000 or higher, the TR, UA and SS methods perform well in reproducing MC behavior in current clamp tests. This means that in the high number of channels regime, bounding the variables to [0,1] is not essential, as the UA and SSda implementations (both unbounded) give the same results than TR. Inaccuracies were observed (**Figure 1**) in the Reflection method, and, to a lesser extent, the CN method.

**FIGURE 6 | Firing statistics and variability in the multi-compartment model. (A)** Mean of spikes per sweep detected at the axon tip in each of the stochastic simulation algorithms. Data was compared to MC using a *t*-test. **(B)** Variance of the spikes per trial in each stochastic simulation. Data was compared to MC using a Bartlett's test for equal variance. **(C)** Number of

common spikes found in each simulation. A common spike is a spike present in more than 50% of the trials with the same timing, considering a window of ±5 ms. **(D)** For all the common spikes detected, the variance in their timing was calculated near the initiation site (blue bars) and at the axon tip (green bars). In **(A–B)**, ∗∗∗*p* < 0.001, ∗∗*p* < 0.01.

**FIGURE 7 | Variability and correlations in the duration of action potentials. (A)** Duration of all action potentials (measured at 0 mV) at different positions of the axon, plotted against the duration at a site near the initiation site. This reference location is designated as "0 μm." Data presented correspond to simulation with MC. *r* is the Pearson correlation coefficient and p is the associated *p*-value when testing for a correlation different to 0. **(B)** Evolution of action potential duration along the axon, for the different simulation algorithms. Data is mean ±

*SD*. **(C)** Correlation of action potential duration at different points in the axon, with the duration at the reference "0μm" location. Lines represent the fit to a single exponential decay. Inset: length constants λ obtained in the fit of the different data sets. Error bars represent the *SD* of the parameter estimation. The fits were compared to MC using an extra sum of squares *F* test, to test the hypothesis that each data set and the MC set could be fit with the same parameter λ. ∗∗∗*p* < 0.001, ∗*p* < 0.05.

Additionally, both methods showed inaccuracy in the voltage clamp tests. Also on these tests, the SS methods showed some inaccuracy in the variance of open channels during the action potential, but the results in current clamp simulations suggest that this may not be relevant for the neural excitability.

The results with the multicompartment model deserve special attention, because some results were different to what was observed in a single compartment. Although the SSmc method introduced only minor inaccuracies in the single compartment test, this method severely altered the excitability of the multicompartment model. We could not identify the reason for this effect, and tested several alternatives of bounding and normalization which did not improve the results. Nevertheless, the overall effect is a reduced excitability, which is consistent with the deviation observed in single compartment (**Figure 1D**). We confirmed that this effect does not arise from our implementation of the algorithm, by repeating the single compartment test with the Matlab code published by Schmandt and Galán (ModelDB acc. 144468) and obtaining similar results (not shown). So, although it is an attractive method for increasing simulation speed in multicompartmental models, its use is not recommended until further testing is performed. The Reflection method altered the excitability of the neuron but the results were in agreement with the single compartment results. Curiously, the use of Stochastic Shielding in

a DA framework (SSda) improves the behavior of the simulations, bringing it closer to the behavior of MCs, like the TR does.

#### **APPLICABILITY AND SIMPLICITY**

In terms of applicability, a first observation is that the colored noise approach (Güler, 2013) can be used only if the ion channel is composed of independent subunits. This is sometimes true (e.g., the original HH model), but not always (e.g., the SB model). For example, channels with non-identical voltage sensors (Vandenberg and Bezanilla, 1991; Horn et al., 2000), cooperativity in the movement of voltage sensors (Bezanilla et al., 1994; Schoppa and Sigworth, 1998) or complex allosteric gating mechanisms (Horrigan and Aldrich, 2002) do not have independent subunits. Therefore, they cannot be modeled with the colored noise approach. Besides, Güler's equations have a number of constant parameters (γ*K*, γ*Na*, ω2 *K*, ω2 *Na*, T*K*, T*Na*) that were estimated empirically to obtain an adequate level of channel noise (Güler, 2013). It is not clear how these parameters can be derived for other ion channels, even if they are composed of identical and independent gating subunits.

The other DA-based algorithms can be applied to any given kinetic scheme but first require to obtain the corresponding system of SDEs. Until recently, it seemed quite complicated to implement, since the original descriptions involved the calculation of a matrix root square (Fox and Lu, 1994; Goldwyn et al., 2011). However, alternative derivations of the Langevin equation (Mélykúti et al., 2010; Orio and Soudry, 2012) yield an explicit form that does not use complex matrix operations. This method can be derived for any given kinetic scheme using simple and intuitive rules without using matrix notation (see Supplemental Material in Orio and Soudry, 2012). Importantly, any kinetic scheme can be translated to an SDE system and the equations can be written explicitly. This makes it simpler to employ low-level or limited languages such as C or Neuron's NMODL, which are compiled prior to execution code and therefore run faster.

The DA algorithms differ in the treatment of boundary and normalization constraints. In this regard, the Unbound DA (UA) is the simplest, not taking care of the boundaries issue and doing a simple normalization by making one variable to depend on the others (alternative normalization procedures can be implemented, for instance dividing all the variable values by their sum). Finally, to avoid non-real square roots in the stochastic terms, a simple absolute value operation is performed. The Stochastic Shielding approximation, when used with DA equations (SSda), can simplify the code even further because it uses less stochastic terms. Huang's truncation and restoration is also rather simple to follow and implement. However, it requires several lines of code and a series of nested *if* and *for* blocks when written in simple languages. As a side note, we found that the restoration procedure requires the specification of the Euler integration method within Neuron's NMODL files. Failing to do so, using instead the default Crank-Nicholson integrator, results in severely distorted results such as a much lower firing rate. Dangerfield's reflection method takes similarly amount of lines of code as TR but it was more complicated to follow and implement.

#### **COMPUTATIONAL EFFICIENCY**

Our comparison of simulation time shows that DA-based methods are not the best choice for all the situations. It has already been noted that when the number of channels is low, MC modeling runs faster than DA (Orio and Soudry, 2012). The limit (number of channels) at which this happens is variable, depending not only on the time step used for numerical integration (as shown in **Figure 8A** vs. **Figure 8B**) but also on the kinetics of the channels being simulated, as this determines the number of transitions occurring in MCs (see below). Interestingly, the Stochastic Shielding in the context of MC (SSmc) behaves like the DA methods: it is faster than MC only with large numbers of channels and its performance depends mainly on the integration time step. This results from the increase in the number of (non-stochastic) ODEs that have to be solved at each integration step, regardless of the actual transitions that occur.

At first glance, it seems that a similar comparison carried out in another work (see Figure 8 in Huang et al., 2013b) produced very different results. In Huang et al. the MC method is faster than DA only if there are less than 20 ion channels. In contrast, we found that MC is more efficient if there are less than ∼500 channels (*dt* = 5μs) or ∼5000 channels (*dt* = 0.5μs). The main reason for this large difference is that in the MC simulations of Huang et al. the state dynamics of each gating particle were updated individually. This method is highly inefficient and should be avoided in stochastic simulations of neurons. In contrast, we used the standard efficient MC method by Gillespie (1976), which tracks only the total number of channels in each state.

Interestingly, in another work (Rowat and Greenwood, 2014) it was found that Güler's CN method is much faster than Unbound DA. However, we have found that it is only a little faster (**Figures 8A,B**). The difference might be attributed to the language (Neuron vs. Python), or some other difference in the implementation.

Within all DA methods, Stochastic Shielding is the fastest in all circumstances. Next comes the colored noise approach (as noted, not applicable to all kinetic schemes) closely followed by Unbounded DA and Reflection. Recall that Stochastic Shielding (SSda) produces a minor loss in accuracy when used with a large number of channels, being in almost all cases indistinguishable from MC modeling. Finally, in our simulations the less efficient (slowest) algorithm was the Truncated and Restored DA (Huang's).

The numerical stability issue also deserves to be considered. Although we did not perform a systematical assessment of numerical stability in our simulations, we noted that simulations with the UA and SSda algorithms produced unreliable results with *NNa* = 50, even at the lowest *dt* of 0.5μs. As the other algorithms did not show this problem, it is most likely due to the lack of variable bounding. We did not pursue in finding a fix for these fails as the low channel number is a condition where MC modeling becomes the fastest and most accurate method.

#### **THEORETICAL ESTIMATES OF ACCURACY AND SPEED**

For a given model, when is it better to use DA instead of MC? Specifically, we would like to know in advance when a DA approach will be accurate, and also faster than MC. A definite answer usually requires some preliminary simulations. However, as we explain next, a rough estimate could be obtained based on the following numbers: the simulation timestep, the number of channels to be simulated, and their typical constant rates.

Suppose we have N ion channels (of some particular type), with X of these ion channels in some state A, and α being a kinetic rate from state A to another B. In each simulation timestep dt, let be the number of channels switching from state A to state B. As different channels are independent, is distributed according to a binomial distriubtion with *n* = *X*, and *p* = α*dt*. Therefore, the average number of channels switching from A to B in that timestep is *np* = *X*α*dt*.

These quantities can be used to estimate the expected accuracy of DA. As explained in Orio and Soudry (2012), the key idea in DA is to use the central limit theorem and approximate the distribution of to be Gaussian. This approximation becomes accurate when *np* = *X*α*dt* 1. This also means that *N*α*dt* 1 since *X* < *N*. For example, in the HH model the slowest kinetic rates (in the relevant voltage range), are about α ∼ 0.1*ms*−1. Then if *dt* = 5μs DA should be expected to be accurate only when *N* > 2000, which is comparable to what we found in our simulations.

Using the same quantities we can also estimate the relative speed of the MC and DA algorithms. Simulation time is roughly proportional to the number of times the simulation variables (the fraction of channels in each) are being updated at each *dt* timestep (in which the voltage is updated). On the one hand, the MC algorithm performs updates each time a single channel switches between states. The number of these updates in each timestep, for each type of switch, is proportional to . Recall that the mean of is equal to *np* = *X*α*dt*. Therefore, in total, about *N*α*dt* updates are performed on each timestep (where α is the appropriate average over all the kinetic rates). On the other hand, the DA algorithms perform a single update at each timestep. Therefore, DA should become more efficient than MC only when *N*α*dt* > 1. For example, in the HH model most state transitions occur near rest voltage in the (fast) *m* kinetics, and so we get approximately α ∼ 2ms−1. This condition yields results comparable to what we found in our simulations (**Figure 8**): if *dt* = 5μs then N > 100, and if *dt* = 0.5μs then N > 1000.

### **CONCLUSIONS**

Suppose *dt* is the simulation timestep, *N* is the number of ion channels, and α is the "typical" transition rate of the channel. Our results suggest that, as a rule of a thumb,

#### **ACKNOWLEDGMENTS**

This work was supported by grants Fondecyt 1130862, ACT-1104, and ACT-1113 to Patricio Orio (Consejo Nacional de Investigación Científica y Tecnológica CONICyT, Chile). The Centro Interdisciplinario de Neurociencia de Valparaíso is a Millenium Science Institute supported by funds P09-022-F, Ministerio de Economía, Chile.

#### **REFERENCES**

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

*Received: 26 May 2014; accepted: 13 October 2014; published online: 03 November 2014.*

*Citation: Pezo D, Soudry D and Orio P (2014) Diffusion approximation-based simulation of stochastic ion channels: which method to use? Front. Comput. Neurosci. 8:139. doi: 10.3389/fncom.2014.00139*

*This article was submitted to the journal Frontiers in Computational Neuroscience. Copyright © 2014 Pezo, Soudry and Orio. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

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

## *Bahar Moezzi\*, Nicolangelo Iannella and Mark D. McDonnell*

*Computational and Theoretical Neuroscience Laboratory, Institute for Telecommunications Research, University of South Australia, Mawson Lakes, SA, Australia*

#### *Edited by:*

*Joshua H. Goldwyn, New York University, USA*

#### *Reviewed by:*

*Ian Bruce, McMaster University, Canada John Wittig Jr., National Institutes of Health, USA*

#### *\*Correspondence:*

*Bahar Moezzi, Computational and Theoretical Neuroscience Laboratory, Institute for Telecommunications Research, University of South Australia, Building W, Mawson Lakes Boulevard, Mawson Lakes, SA 5095, Australia e-mail: bahar.moezzi@ mymail.unisa.edu.au*

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.

**Keywords: calcium dynamics, stochastic synapse, inner hair cell, auditory nerve, short term depression, neural variability, channel noise**

## **1. 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).

## **2. MATERIALS AND METHODS**

Firstly, in Section 2.1, we review the previous models that our research is built upon:

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.

## **2.1. PREVIOUS MODELS**

#### *2.1.1. Inner hair cell model*

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*, *Gc*, *Ec*, ν, 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.

In summary, the model has the following parameters:

The values of these parameters are summarized in **Table 1**. The equations describing conversion from *v*(*t*) to *k*(*t*) are

$$k(t) = \max(0, bc^3(t)),\tag{1}$$

$$\frac{dc(t)}{dt} = -\frac{c(t)}{\tau\_c} + \nu G\_c m^3(t)(E\_c - \nu(t)),\tag{2}$$

$$\frac{dm(t)}{dt} = -\frac{m(t)}{\tau\_m} + \frac{1}{\tau\_m \left(1 + \frac{e^{-\gamma\nu(t)}}{\beta}\right)},\tag{3}$$

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−1s −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

$$\nu = \frac{L}{2\mathcal{F}C\_{\nu}\xi},\tag{4}$$

where *F* is Faraday constant, *Cv* 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 × 109 MA−1s −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

**Table 1 | Parameters for inner hair cell calcium levels.**

*The value of gc was obtained from Zampini et al. (2013). The value of v was obtained by running MAP BS with no stimulus present. All other values are identical to those used in publicly available Matlab source code MAP BS at http:// www*.*essexpsychology*.*macmate*.*me/ HearingLab/ modelling*.*html.*

*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*) < *Ec*. Therefore, we have used (*Ec* − *v*(*t*)) in Equation (2) instead of (*v*(*t*) − *Ec*) 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.

### *2.1.2. 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":

There are four parameters that have units of rate:

The values of these parameters are summarized in **Table 3**. The deterministic Meddis (1986) synapse model is of the following form

$$\frac{d\mathbf{x}(t)}{dt} = \mathbf{A}(t)\mathbf{x}(t) + \mathbf{B},\tag{5}$$

where

**Table 3 | Parameters for neurotransmitter release with values identical to those used in publicly available Matlab source code MAP BS at http://www.essexpsychology.macmate.me/HearingLab/ modelling.html.**

$$\mathbf{A}(t) = \begin{bmatrix} -r\_1 - k(t) & 0 & r\_2 \\ k(t) & -r\_3 - r\_4 & 0 \\ 0 & r\_4 & -r\_2 \end{bmatrix},$$

$$\mathbf{B} = \begin{bmatrix} r\_1 M \\ 0 \\ 0 \end{bmatrix} \text{and } \mathbf{x}(t) = \begin{bmatrix} x(t) \\ y(t) \\ z(t) \end{bmatrix}. \tag{6}$$

#### *2.1.3. 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,

$$\frac{d\mathbf{x}}{dt} = \mathcal{B}(r\_1, (M - \mathbf{x}(t))) + \mathcal{B}(r\_2, \lfloor z(t) \rfloor) - \mathcal{B}(k(t), \mathbf{x}(t)), \tag{7}$$

$$\frac{dy}{dt} = \mathcal{B}(k(t), \mathbf{x}(t)) - r\_3 \mathbf{y}(t) - r\_4 \mathbf{y}(t),\tag{8}$$

$$\frac{dz}{dt} = r\_4 \mathbf{y}(t) - \mathcal{B}(r\_2, \lfloor z(t) \rfloor). \tag{9}$$

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

### *2.1.4. 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.

### *2.1.5. 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,

$$p\_{\rm conv}(t) = hy(t)dt,\tag{10}$$

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.

$$p\_{\rm conv}(t) = \begin{cases} 0 & \text{if } \quad t - t\_l < t\_A, \\ 1 - C\_r e^{-(\frac{t - t\_l - t\_A}{t\_R})} & \text{if } \quad t - t\_l \ge t\_A, \end{cases}, \tag{11}$$

where *Cr* = 1, *tR* = 0.6 ms is the time constant of relative refractoriness, *tA* = 0.75 ms is the ARP, *t* is the current time, and *tl* 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), *Cr* = 0.55 and *tR* = 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.

#### **2.2. PREVIOUS STATISTICAL ANALYSIS**

#### *2.2.1. Empirical vesicle release distribution*

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.

To calculate the ISI parameters, ARP and RRP in the form of Equation (11) are used. Two additional parameters are involved:

Heil et al. (2007) obtained the following equation describing the ISI probability density function (PDF),

$$D(t) = \begin{cases} \frac{\theta}{t\_R(\theta - \frac{1}{t\_R})} ( (e^{-(\frac{t - t\_A}{t\_R})} - e^{-\theta(t - t\_A)} )(1 - \rho + \frac{\rho \theta}{\theta - \frac{1}{t\_R}}) \\ -\rho \theta (t - t\_A) e^{-\theta(t - t\_A)} ) & \text{for } t \ge t\_A, \\ 0 & \text{for } t < t\_A, \end{cases} \tag{12}$$

#### *2.2.2. 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:

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

#### **2.3. NEW MODELS**

## *2.3.1. Short-term depression in vesicle release probability* **(STDv)**

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:

The model for the change of *k*(*t*) over time is

$$\frac{dk(t)}{dt} = \frac{\max(0, bc^3(t)) - k(t)}{\tau\_t} + ak(t) \sum\_{i} \delta(t - t\_{\mathbb{V}\_i}), \quad (13)$$

where *tvi* is the time of *i*th release.

#### *2.3.2. 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.

*<IP_ADDRESS>. Biophysical model.* In the Meddis (2006) model, long term correlation observed in AN fibers can be partially explained by depletion of readily available vesicles, as explained in the Results section. In the Results section, we show that the Meddis (2006) depleted model with readily available vesicles depleted by decreasing the maximum number of vesicles in the immediate store, or by increasing the spontaneous rate, both require much higher vesicle release rate than the non-depleted Meddis (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).

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); Zilany et al. (2009); Zilany and Carney (2010); Zilany et al. (2014) is the random fluctuations in the number of open and closed calcium ion channels as they are expected to cause variability in vesicle release probabilities.

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

$$\alpha(\nu) = \frac{1}{\mathfrak{r}\_m \left(1 + \frac{\varepsilon^{-\gamma \cdot \nu(t)}}{\beta}\right)},\tag{14}$$

$$
\beta(\nu) = \frac{1}{\tau\_m} - \alpha(\nu). \tag{15}
$$

Equation (2) therefore changes to

$$\frac{dc(t)}{dt} = -\frac{c(t)}{\mathbf{r}\_c} + \nu \mathbf{g}\_c n(t)(E\_c - \nu(t)),\tag{16}$$

where *n*(*t*) is the number of open calcium channels out of total of *N* calcium channels and *gc* is the single calcium channel conductance.

*<IP_ADDRESS>. 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:

$$\frac{dc(t)}{dt} = -\frac{c(t)}{\tau\_c} + \nu G\_c(f(m^3(t) + X(t)))(E\_c - \nu(t)), \quad (17)$$

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

$$dX(t) = \frac{dt}{\mathfrak{r}\_o}(\mu\_o - X(t)) + \sigma\_o dW\_t,\tag{18}$$

where *Wt* denotes the Wiener process and the mean (μ*o*), time constant (τ*o*) and variance (σ*o*) of the noise are positive constants.

#### *2.3.3. 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:

$$\frac{dc(t)}{dt} = -\frac{c(t)}{\mathfrak{r}\_{\varepsilon}} + \nu G\_{\varepsilon} m^3(t) (E\_{\varepsilon} - \nu(t) - X(t)),\qquad(19)$$

#### *2.3.4. 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

$$\frac{dk(t)}{dt} = \frac{k\_1(t) - k(t)}{\mathbf{r}\_s} + ak(t) \sum\_i \delta(t - \mathbf{t}\_{\mathbf{v}\_i}),\tag{20}$$

where *k*1(*t*) is the vesicle release rate when Ornstein Uhlenbeck noise is added as calcium channel noise in the Meddis (2006) model.

#### **2.4. PARAMETERS**

The parameters in **Table 1** (except *gc*), in **Table 3**, and for*tA* and *tR* (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 *gc* 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*, *tR* 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) suggested a readily release pool of about 22 docked vesicles in the IHC of mouse. Pangršic et al. ˇ (2010) estimated a readily releasable pool of 12 vesicles per active zone in the pachanga mouse. We assumed the maximum readily available pool size, *M*, to be 20 vesicles per active zone.

### **3. RESULTS**

#### **3.1. PREVIOUS MODELS**

**Figures 2A** (Gray) and **2B** (Gray) show the Fano factor time curve of a spike train generated by the Zilany et al. (2014) synapse model. These figures are obtained by running the model code available at http://www.urmc.rochester.edu/labs/ Carney-Lab/publications/auditory-models.cfm, with a relative membrane voltage input of 0 V. Like the empirical Fano factor of **Figure 2A** (Light blue), the Fano factor increases to about 10 for large counting time windows. The Fano factor in **Figure 2A** (Gray) does not decrease below one for shorter time windows as much as the empirical Fano factor shown in **Figure 2A** (Light blue) does. The time scales of the correlation do not match empirical data of **Figure 2A** (Light blue).

In the Meddis (2006) model, long term correlation observed in the auditory nerves can be partially explained by depletion of readily available vesicles in the immediate store. In **Figures 2A,B**, the blue trace is the Fano factor time curve for the original Meddis (2006) model with a maximum number of readily available

**FIGURE 2 | Time window dependent Fano factor for spontaneous activity in the auditory nerve obtained from previous models. (A)** Gray: Fano factor time curve for a spike train generated by the Zilany et al. (2014) model with SR∼100 spikes per second. Blue, green, red: in the original Meddis (2006) model, by decreasing the maximum number of available vesicles from M = 20 with SR∼65 spikes per second (Blue) to M = 6 with SR∼65 spikes per second (Red) or increasing the spontaneous firing

rate from SR∼65 spikes per second with M = 20 to SR∼160 spikes per second with M = 20 (Green), short term negative correlation starts at smaller time windows and long term positive correlation is increased. 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 4**.

vesicles of *M* = 20 and a spontaneous rate of around 65 spikes per second. The red trace is the Fano factor time curve for the original Meddis (2006) model with a maximum number of readily available vesicles of *M* = 6 and a spontaneous rate of around 65 spikes per second. The green trace is the Fano factor time curve for the original Meddis (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.

#### **3.2. 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, *tA* and *tR*, 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 *tA* = 0.69 ms, and *tR* = 0.58 ms when SR = 65 spikes per second and θ = 0.0862 (ms−1) and ρ = 0.43 for *tA* = 0.73 ms, and *tR* = 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 *tA* = 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.

#### **3.3. CALCIUM CHANNEL NOISE**

#### *3.3.1. 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.

#### *3.3.2. 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.

#### **3.4. 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**.

model. Green: adding a combination model of the phenomenological channel noise and short term depression in vesicle release to 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). **(B)** Time-window dependent Fano-factor of **(A)** on linear axes for time windows shorter than 0.1 s. **(C)** Fano factor for the Meddis (2006) model with

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

depression in vesicle release for different maximum numbers of vesicles in the immediate store on linear axes for time windows shorter than 0.1 s. **(D)** Vesicle release rate for **(A)**. The color representations are the same as in **(A)**. subfigure is a zoom into the top subfigure. The rapid decreases in the green trace in the **(D)** bottom plot for the combination model of phenomenological channel noise and short term depression in vesicle release are due to actual vesicle release while remaining fluctuations are due to channel noise. Parameter are summarized in **Tables 7**, **8**.

#### **3.5. 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.

## **4. 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

**Table 6 | Parameters of Meddis (2006) model with biophysical calcium channel noise.**

*M* = *20 and SR*∼*65 spikes per second.*

**Table 5 | Comparison of the original Meddis (2006) model and Meddis (2006) model with relative refractoriness in the auditoy nerve substituted by short term depression in vesicle release.**

*M* = *20 and SR*∼*65 spikes per second. For fitting to Equation (12), in the Equation (12) tA* = *0.75 ms and tR* = *3.5 ms were used. Log likelihood was used as a measure of goodness of fit.*

**Table 7 | Parameters for the phenomenological models of stochastic variability in the IHC-AN complex.**

*M* = *20, and SR*∼*65 spikes per second.*

**Table 8 | Parameters of the combination model of phenomenological channel noise and short term depression in vesicle release probability with various maximum numbers of vesicles in the available store.**

τ*o* = *1.2 ms,* σ*o* = *0.3,* μ*o* = *0.38,* τ*s* = *2.5 ms, a* = *0.001 and SR*∼*65 spikes per second.*

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.

#### **FUNDING**

Mark D. McDonnell's contribution was supported by the Australian Research Council under ARC grant DP1093425 (including an Australian Research Fellowship).

#### **ACKNOWLEDGMENTS**

We would like to thank Ray Meddis of the University of Essex, Bruce Graham of the University of Stirling, Nigel G. Stocks of The University of Warwick, Anthony N. Burkitt and David B. Grayden of The University of Melbourne, and Christian Stricker of Australian National University for helpful discussions.

#### **REFERENCES**

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

*Received: 05 June 2014; accepted: 26 November 2014; published online: 23 December 2014.*

*Citation: Moezzi B, Iannella N and McDonnell MD (2014) Modeling the influence of short term depression in vesicle release and stochastic calcium channel gating on auditory nerve spontaneous firing statistics. Front. Comput. Neurosci. 8:163. doi: 10.3389/ fncom.2014.00163*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

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

## On how correlations between excitatory and inhibitory synaptic inputs maximize the information rate of neuronal firing

## *Pavel A. Puzerey and Roberto F. Galán\**

*Department of Neurosciences, School of Medicine, Case Western Reserve University, Cleveland, OH, USA*

#### *Edited by:*

*Joshua H. Goldwyn, New York University, USA*

#### *Reviewed by:*

*Robert Rosenbaum, University of Pittsburgh, USA Woodrow Shew, University of Arkansas, USA*

#### *\*Correspondence:*

*Roberto F. Galán, Department of Neurosciences, School of Medicine, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106-4975, USA e-mail<EMAIL_ADDRESS>
Cortical neurons receive barrages of excitatory and inhibitory inputs which are not independent, as network structure and synaptic kinetics impose statistical correlations. Experiments *in vitro* and *in vivo* have demonstrated correlations between inhibitory and excitatory synaptic inputs in which inhibition lags behind excitation in cortical neurons. This delay arises in feed-forward inhibition (FFI) circuits and ensures that coincident excitation and inhibition do not preclude neuronal firing. Conversely, inhibition that is too delayed broadens neuronal integration times, thereby diminishing spike-time precision and increasing the firing frequency. This led us to hypothesize that the correlation between excitatory and inhibitory synaptic inputs modulates the encoding of information of neural spike trains. We tested this hypothesis by investigating the effect of such correlations on the information rate (IR) of spike trains using the Hodgkin-Huxley model in which both synaptic and membrane conductances are stochastic. We investigated two different synaptic input regimes: balanced synaptic conductances and balanced currents. Our results show that correlations arising from the synaptic kinetics, τ , and millisecond lags, δ, of inhibition relative to excitation strongly affect the IR of spike trains. In the regime of balanced synaptic currents, for short time lags (δ ∼ 1 ms) there is an optimal τ that maximizes the IR of the postsynaptic spike train. Given the short time scales for monosynaptic inhibitory lags and synaptic decay kinetics reported in cortical neurons under physiological contexts, we propose that FFI in cortical circuits is poised to maximize the rate of information transfer between cortical neurons. Our results also provide a possible explanation for how certain drugs and genetic mutations affecting the synaptic kinetics can deteriorate information processing in the brain.

**Keywords: stochastic Hodgkin-Huxley model, synaptic kinetics, input correlation, information, feed-forward inhibition**

## **INTRODUCTION**

The rate and timing of firing in cortical neurons is strongly affected by the interaction between synaptic excitation and inhibition (Salinas and Sejnowski, 2000). The architecture of cortical circuits ensures that the magnitude of excitatory and inhibitory synaptic inputs is approximately balanced on average and temporally correlated (Shu et al., 2003b; Haider et al., 2006), albeit with a small time delay for inhibition of ∼1–10 ms (Wehr and Zador, 2003; Okun and Lampl, 2008; Wu et al., 2008). This correlation in amplitude and timing presumably arises in feedforward inhibition (FFI) circuits, an anatomical motif present ubiquitously throughout the cortex which drives monosynaptic excitation and disynaptic inhibition onto target neurons (Porter et al., 2001; Sun et al., 2006; Cruikshank et al., 2007). The functional consequences of the correlations imposed by such a layout are far-reaching, encompassing a range of functions such as gain modulation for rapidly fluctuating synaptic inputs (Salinas and Sejnowski, 2000; Chance et al., 2002; Shu et al., 2003a; Pouille et al., 2009), shaping of neuronal tuning properties and stimulus selectivity (Wehr and Zador, 2003; Marino et al., 2005; Wu et al., 2006), directing the propagation of activity by selectively gating firing in neuronal ensembles (Kremkow et al., 2010a,b), and creating "windows of integration" during which excitatory inputs can temporally summate to promote spike generation before being rapidly suppressed by inhibition (Pouille and Scanziani, 2001; Pouille et al., 2009). Furthermore, Marsalek and colleagues demonstrated that small differences in the timing between presynaptic excitatory and inhibitory inputs (i.e., input correlation) is directly correlated with temporal jitter in postsynaptic spikes, i.e., output precision (Marsalek et al., 1997). Since neurons may represent information through the precise timing of spikes (decharms and Merzenich, 1996; Dan et al., 1998; Strong et al., 1998; Liu et al., 2001; Nemenman et al., 2008), it stands to reason that the control of spike timing by correlated excitation and inhibition is likely to govern the transfer of information between cortical neurons. Previous investigations using realistic simulations of cortical neurons have shown that, indeed, balanced excitatory and inhibitory synaptic currents maximize both coding and metabolic efficiency of neuronal spikes (Sengupta et al., 2013). In that study, however, the excitatory and inhibitory synaptic conductances were *uncorrelated* and as a result did not exhibit the correlation characteristic of cortical dynamics under experimental contexts (Wehr and Zador, 2003; Okun and Lampl, 2008; Wu et al., 2008). Furthermore, Kawaguchi and colleagues have shown that the relative balance between excitation and inhibition of a random synaptic input to simulated pyramidal neurons controls the maximal information content of spike trains in the presence of background synaptic noise (Kawaguchi et al., 2011), yet the timing of excitation and inhibition with respect to each other were not considered. To our best knowledge, the relevance of statistical correlations between balanced excitatory and inhibitory synaptic inputs for the information rates of neural spike trains have not been investigated. This may in part be due to the high computational cost of realistic simulations of stochastic neuronal dynamics. Here, we have overcome this limitation by using the stochastic-shielding approximation, which was recently introduced by our lab, accelerating stochastic simulations by up to two orders of magnitude while preserving accuracy (Schmandt and Galán, 2012).