High Entrainment Constrains Synaptic Depression Levels of an In vivo Globular Bushy Cell Model

Globular bushy cells (GBCs) located in the ventral cochlear nucleus are an essential part of the sound localization pathway in the mammalian auditory system. They receive inputs directly from the auditory nerve and are particularly sensitive to temporal cues due to their synaptic and membrane specializations. GBCs act as coincidence detectors for incoming spikes through large synapses—endbulbs of Held—which connect to their soma. Since endbulbs of Held are an integral part of the auditory information conveying and processing pathway, they were extensively studied. Virtually all in vitro studies showed large synaptic depression, but on the other hand a few in vivo studies showed relatively small depression. It is also still not well understood how synaptic properties functionally influence firing properties of GBCs. Here we show how different levels of synaptic depression shape firing properties of GBCs in in vivo-like conditions using computer simulations. We analyzed how an interplay of synaptic depression (0–70%) and the number of auditory nerve fiber inputs (10–70) contributes to the variability of the experimental data from previous studies. We predict that the majority of synapses of GBCs with high characteristic frequencies (CF > 500 Hz) have a rate dependent depression of less than 20%. GBCs with lower CF (<500 Hz) work also with strong depressing synapses (up to 50% or more). We also showed that synapses explicitly fitted to in vitro experiments with paired-pulse stimuli did not operate properly in in vivo-like conditions and required further extension to capture the differences between in vitro and in vivo experimental conditions. Overall, this study helps to understand how synaptic properties shape temporal processing in the auditory system. It also integrates, compares, and reconciles results of various experimental studies.


INTRODUCTION
Globular bushy cells (GBCs) are located in the ventral cochlear nucleus, which is the first processing station in the central auditory nervous system. Cochlear nucleus receives direct inputs from the inner ear through auditory nerve fibers (ANFs). In addition to GBCs, the population of neurons in ventral cochlear nucleus is made of several different types: spherical bushy cells, stellate cells, and octopus cells (Osen, 1969;Rhode et al., 1983). Each cell type processes different sound features, e.g., spectrum, temporal fine structure, or signal onset (Recio, 2000). The pre-processing of sound by the cochlear nucleus neurons is essential for higher auditory centers to perform their functions, e.g., sound localization or sound identification.
GBCs enhance temporal sound cues which are essential for sound localization. When stimulated with pure tones, their action potentials lock precisely to a certain phase of the signal. Some of the GBCs are experts in phase-locking and reach synchronization levels higher than any ANF. Neurons with such properties are sometimes called "high-sync" neurons (Joris et al., 1994). Additionally, they can fire an action potential every stimulus cycle up to 700 Hz, which is more than twice the maximum firing rate for ANFs. GBCs similarly show increased synchronization and entrainment in response to complex sounds (Rhode and Greenberg, 1994). Rothman et al. (1993) explained the mechanism of high synchronization (temporal precision of spikes) and high entrainment (one spike for each stimulus cycle) by the presence of many converging inputs. GBCs receive their main excitatory inputs from ANFs directly onto the soma through large synapses called modified endbulbs of Held (Spirou et al., 2005). To elicit an action potential, multiple subthreshold ANF input spikes must coincide. As a result, GBC spikes have greater temporal precision than the input action potentials from individual ANFs.
In most in vitro experiments, the strength of the excitatory inputs (excitatory postsynaptic currents, EPSC, synaptic strength) is strongly reduced after a presynaptic spike as shown in Figure 3. In the time between two input spikes, the synaptic strength recovers and if the pause between stimulations is long enough, it eventually reaches its original value (Wang and Manis, 2008;Xu-Friedman, 2009, 2015). This synaptic depression at the endbulb of Held is associated mostly with the depletion of neurotransmitter at the presynaptic site (Friauf et al., 2015) and receptor desensitization (Yang and Xu-Friedman, 2008). The level of depression depends strongly on the stimulation frequency (Yang and Xu-Friedman, 2009). A similar behavious was observed in a similar auditory synapse in medial nucleus of the trapezoid body, calyx of Held (Taschenberger et al., 2016). In contrast to in vitro experiments, in vivo extracellular recordings do not indicate strong variations of synaptic strength at the endbulbs of Held during stimulation (Young and Sachs, 2008;Borst, 2010;Kuenzel et al., 2011Kuenzel et al., , 2015. The few in vivo studies of endbulb of Held can be complemented with in vivo studies of calyx of Held, that shows slight decrease in synaptic strength only at high stimulation frequencies (Lorteije et al., 2009;Di Guilmi et al., 2014).
Currently, the difference between the in vitro and in vivo observations is not well understood. In particular, how the different levels of depression influence the firing properties in response to sounds is unclear. Both endbulb and calyx of Held are regarded as "model synapses." It is therefore important to obtain a precise understanding of their operation. In this study, we present a biophysically plausible computational model of GBCs, which is driven by realistic in vivo-like stimuli. The main goals of the study were: • To optimize the synaptic parameters of endbulbs of Held to accurately simulate GBCs including "high-sync" neurons (Sections 2.1 and 2.3).
• To better understand the behavior of synapses with different depression levels (from tonic, through weakly depressing, to strongly depressing as seen in in vitro experiments) in in vivo-like simulations (Figure 5, Sections 3.1 and 3.3). • To demonstrate how the higher number of ANF inputs increases the firing rate and the entrainment of bushy cells (Section 3.2). Experimental evidence shows similar variability (Joris et al., 1994;Spirou et al., 2005). • To provide a more comprehensive auditory model by combining an inner ear model with a GBC model, that can be used for further processing, e.g., for sound localization models. In this case, it is not sufficient that both models independently generate realistic outputs but also their combination, which provides further constraints.

Model Overview
The model consisted of three main components: (1) ANF inputs generated by an inner ear model (Zilany et al., 2014), (2) endbulbs of Held, and (3) a soma of a bushy cell. The general structure is shown in Figure 1 and was derived from the previous work of Rothman et al. (1993). First, the sound enters the cat inner ear model by Zilany et al. (2014) and is converted to ANF spike trains. The model of Zilany et al. (2014) is a state-of-the-art model capable of generating ANF responses with typical adaptation properties an with realistic phase locking (Zilany et al., 2009). The model includes several processing stages, such as a basilar membrane, inner hair cells, power-law dynamics at the inner hair cell synapses, and spike generators.
Next, we selected a number of ANF spike trains to drive synapses that were attached to the soma of the bushy cell. According to an electron microscopy study by Spirou et al. (2005), the number of ANF inputs shows large variations: between 9 and 69 ANF inputs across a population of 12 bushy cells. We simulated GBCs with the number of inputs varying from 10 to 70 per bushy cell. The inputs were coming solely from high-spontaneous rate fibers. We neglected low and medium spontaneous rates fibers after we verified (not shown) that the inclusion of a small number of them does not have significant influence on the firing properties in the simulated experiments Liberman (1991). The model synapses were instantly activated by ANF input spikes and had an exponential decay of the conductance with a short 0.2 ms time constant (Yang and Xu-Friedman, 2009;Cao and Oertel, 2010). The summation of multiple synaptic conductances was linear. With respect to synaptic depression, we examined three different synapse types: (1) a tonic synapse model, (2) single exponential recovery synapse with different levels of synaptic depression, and (3) double exponential recovery synapse tuned to in vitro experiments (Section 2.2).
The final component of the model was the soma of a GBC. We neglected dendritic trees, because the dendritic inputs are not well known and somatic ANF inputs are the major ones (Spirou et al., 2005). Each GBC was modeled as a single compartment with Hodgkin-Huxley-like channels (Rothman and Manis, 2003a,b,c) described by the following equation: where I is the membrane current, C is the membrane capacity, V denotes the membrane voltage, I Klt is the fast-activating slow-inactivating low-threshold K + current, I Kht is the highthreshold K + current, I h is the hyperpolarization-activated cation current, I leak is the leakage current, and I Na is the Na + current. The model was derived directly from Rothman and Manis (2003c) with the exception of sodium channels taken from Spirou et al. (2005). The original sodium channels produced unrealistically long refractory periods and needed to be replaced by channels with faster dynamics. The substitution resulted in the absolute refractory period of 0.66 ms. All channel parameters are summarized in Table 1 (Model Type II) of Rothman and Manis (2003c) and in the Appendix of Spirou et al. (2005) with the following adjustments. Time constants were corrected for the in vivo temperature of 37 • C using Q 10 = 2.5 for sodium and Q 10 = 3 for all other channels. To match measured firing thresholds (Yang and Xu-Friedman, 2009), the maximum Na + channel conductance was increased to 2,500 nS. Table 1 lists synapse models used in this study. The tonic synapse was the first and the simplest model that was examined. The synaptic weight was constant and resulted in constant EPSCs that did not depend on the ANF firing rates. Even though endbulbs of Held might not be perfectly tonic synapses, especially at higher stimulation rates (Young and Sachs, 2008;Borst, 2010), this model is equivalent to a 0%-depressing synapse and defines the bound of the parameter space.

Modeling of Synaptic Depression
The second model was a depressing synapse with a single exponential recovery (Tsodyks and Markram, 1997) which was referred as X%-depressing synapse throughout this study. The purpose of this model was to examine a wide range of synaptic depressions spanning from strong in vitro to weak depression observed in some in vivo experiments. The synapses were fitted only in the in vivo operational stimulation range. The peak synaptic conductance was proportional to the active internal synaptic resources. Every pre-synaptic spike caused an activation of some synaptic resources. Between the spikes, the synapse could recover exponentially. This process can be expressed iteratively with Equation (2): where w is the peak conductance of a completely recovered synapse, u is the fraction of available synaptic resources activated at each stimulation. As stated by Tsodyks and Markram (1997), depending on the physical mechanisms, u can be partially or completely determined by release probability of the neurotransmitter triggered by a presynaptic spike. t is the time from the last (n-th) synaptic event, and τ is the synaptic recovery time constant. τ was set to 90 ms based on in vitro experiments with in vivo-like stimuli of endbulbs (Yang and Xu-Friedman, 2015) and calyx of Held (Hermann et al., 2009). The initial value of g (g 0 ) was set to w/(1 − u) so that g 1 had the peak conductance of w. The depression level X% is defined as a relative decrease of steady-state synaptic strength s f between the spontaneous rate (SR) (f = 50 Hz) and the driven rate (f = 300 Hz). It can be calculated with Equation (3): The value of X% can be set by adjusting u in the model, because it directly changes the steady-state synaptic weights (s f ) as shown in Figure 2. The stimulation frequencies of 50 and 300 Hz were chosen as they correspond approximately to an input operational range bound by the spontaneous and driven rates of ANFs. In other studies, especially in vitro, X% can have a different meaning and refer to depression relative to an unconditioned synapse. As a result, the depression values are typically much higher in in vitro synapses. The third synapse model was a depressing synapse with a double exponential recovery. It is called yang2009mean mean in this study, because it is based on the mean in vitro data from Yang and Xu-Friedman (2009). The purpose of this model is to replicate behavior of synapses measured in vitro and validate them in simulated in vivo-like conditions. In their experiments, Yang and Xu-Friedman (2009) observed two exponential components in the recovery of EPSCs. Thus a model based on double exponential recovery, such as Cook et al. (2003), was a natural match to the experimental data. The model is described by Equation (4) and the parameters were set to the mean values of the population from Yang and Xu-Friedman (2009).
where w is the peak conductance of a completely recovered synapse, u is the fraction of the available synaptic resources being utilized at presynaptic event. We interpret this value as the vesicle release probability of 0.6, that was found for endbulbs of Held in vitro (Oleskevich et al., 2000). t is the time from the last (n-th) synaptic event, τ s (1,990 ms) and τ f (10.9 ms) are the slow and

Name Description
Tonic Constant synaptic weight X%-depressing Synapse with a single exponential recovery depressing by X% when stimulated at 300 Hz (Equation 3) yang2009mean Synapse with a double exponential recovery and mean parameters of an endbulb population from in vitro experiments by Yang and Xu-Friedman (2009) FIGURE 2 | Normalized peak synaptic weights of single exponential recovery synapses stimulated at two different frequencies.
The stimulation in the first interval (0-0.5 s) was 50 Hz (approximately SR of the input), and 300 Hz (approximately max rate of input) in the second one (0.5-1 s). The steady states in each interval (s f ) were used to calculate the effective depression level (X%) using Equation (3). All synapses share the same recovery time constant but have different values of u.
fast synaptic recovery time constants, respectively, and k (0.3) is the fraction of the fast recovery (relative to the slow recovery). The initial value of g (g 0 ) was set to w/(1 − u). The model exhibits different depression levels depending on the stimulation frequency. Figure 3 shows normalized EPSC of the model with different values of u for reference that was stimulated at 100, 200, and 333 Hz. They grayed area corresponds to the population data obtained in an analogous in vitro experiment from Yang and Xu-Friedman (2009). It is worth to note that Equation (2) describing X%-depressing models can be derived from Equation (4) of the yang2009mean model by setting k to 1. Also, the tonic synapse is a special case of X%-depressing model, where X = 0. However, because the models have different motivations and were fitted differently, they are presented as three distinct models in this study.

Fitting of Synaptic Weights
Modeling of GBCs, and "high-sync" neurons in particular, required adjustment of synaptic weights of the endbulbs of Held. The experiments of Joris et al. (1994) were used as a reference for this procedure. In our simulations, bushy cells were simulated with a in vivo-like stimulation protocol which consists of several steps: 1. A signal was generated. It consisted of 100 gated tones that were repeated every 100 ms and each tone lasted 25 ms. The frequency of the tones was equal to the CF of the stimulated neuron and their intensity was 50 dB SPL .
2. An inner ear model driven by that signal was used to generate the ANF spike trains. 3. The ANF spike trains were used to stimulate models of endbulbs of Held and GBCs. 4. The output spike trains of GBCs were recorded.
After each simulation, vector strength (VS), firing rate, and entrainment index (EI) were calculated from responses of GBCs.
VS is a measure that specifies the strength of phase locking in response to periodic stimuli (usually pure tones). The VS is calculated assuming that each spike represents a unit vector v with a phase equal to the phase of the stimulus at which it occurred. VS is equal to the magnitude of the sum of all such vectors normalized by the number of vectors n: The values of VS are between 0 and 1. VS of 0 indicates that spikes are uniformly distributed over the stimulus period and 1 indicates that all spikes occur at exactly the same phase of the stimulus. "High-sync" neurons are defined in terms of VS > 0.9 for CF < 700 Hz (Joris et al., 1994). EI specifies whether a neuron fired action potentials once for every period of the stimulus or less. The entrainment index p is calculated using the definition given by Joris et al. (1994). First an inter-spike interval histogram is constructed. Then the number of intervals k belonging to the first maximum (between 0.5 and 1.5 of the stimulus period) is divided by the total number of intervals n as given by: The values of EI are between 0 and 1. A value of 1 indicates that the neuron fired a minimum of one spike during each stimulus period. Entrainment falls to 0 when the neuron can no longer fire action potentials for adjacent stimulus periods. Generally, VS, EI, and SR were calculated for each combination of synapse type T (where T ∈ {tonic, X%-depressing, yang2009mean}), ANF number of inputs N (where N ∈ [10, 70]), and synaptic weight w. For simplicity, all synaptic weights had equal values for a given GBC. The parameter scan resulted in the following mapping: (T, N, w) → (VS, EI, SR). VS, EI, and SR can be seen as functions of w for all combinations of T and N: VS T,N (w), EI T,N (w), SR T,N (w). As an example, results with two synapse types and two different values of N are shown Figure 4.
On the one hand, an increase in w causes an increase in VS and EI, which is desirable for GBCs. On the other hand, increasing w can cause unrealistically high SR. As a result, the optimal weights were found by increasing the value of w up to the point where SR T,N (w) crossed the mean experimental SR from Spirou et al. (1990) which equals to 7.5 ± 13.8 spikes/s and is valid only for GBCs with CF < 6 kHz. The procedure of finding the best weight w o is equivalent to Equation (7): FIGURE 4 | SR, EI, and VS as a function of synaptic weight for a combination of two synapse models (10%-depressing and yang2009mean) and two different number of ANF inputs (20 and 35). SR and EI increase with increasing synaptic weight. VS also increases with synaptic weight for 10%-depressing synapses, but does not increase for synapses tuned to Yang and Xu-Friedman (2009). The thin dotted black line marks the desired SR of 7.5 spikes/s (Spirou et al., 1990).
where SR t is the target spontaneous rate (7.5 spikes/s). The procedure of finding optimal weights is further discussed in Section 4.1. Figure 5 shows postsynaptic currents of two GBCs with 10%depressing and yang2009mean synapses ( Table 1) having optimal weights. The GBCs were stimulated with a 50 ms pure tone of 650 Hz. Each plot has a corresponding curve in Figure 4 with 35 ANF inputs. The optimal synaptic weights can be found in Figure 4 at the intersection of the target SR (dotted line, 7.5 spikes/s) and the corresponding curves. The optimal synaptic weights for all analyzed synapse types and number of inputs are plotted in Figure 6.

Implementation Details
The model of GBCs was implemented in the NEURON simulation environment (Carnevale and Hines, 2006) and can be freely accessed, studied, and improved at https://github.com/ FIGURE 5 | Postsynaptic currents of two GBCs with 10%-depressing and yang2009mean synapses stimulated with a 50 ms pure tone of 650 Hz (green bar). The dashed line is plotted for reference and indicates a threshold for a single EPSC to initiate an action potential. Red stars indicate action potentials triggered by the stimulation. The less depressing synapses (10%-depressing) produced more action potentials than the more depressing synapses (yang2009mean, Table 1). The SR for both cells was 7.5 spikes/s. mrkrd/cochlear_nucleus (Rudnicki and Hemmert, 2017). Python programming language was used to manage the simulations and analyze the data. The model input was generated using the free cochlea Python library, which we made publicly available at Rudnicki and Hemmert (2014). Specifically, the library includes the inner ear model of Zilany et al. (2014).

RESULTS
This study mainly aimed to investigate the influence of different synapse types (tonic, with single, and with double exponential recovery), depression levels (0-70%), and the number of ANF inputs (10-70) on the firing properties of bushy cells. We generated in vivo-like stimulation patterns and analyzed VS and EI of the responses. The reference data came primarily from in vivo experiments of Joris et al. (1994).

Influence of Synaptic Depression on Synchronization and Entrainment
In the first set of simulations, we examined how different synapse types and depression levels influence synchronization and entrainment for cells with different CFs. We kept the number of ANF inputs constant, so every GBC was connected to 40 high-SR fibers. The stimulus was a train of 100 ramped pure tones (50 ms, 60 dB SPL ) at CF (Joris et al., 1994). Consistently with previous studies (Rothman et al., 1993;Xu-Friedman and Regehr, 2005;Joris and Smith, 2008), we found that GBCs improved phase locking relative to ANFs in the lowfrequency range up to approximately 1.5 kHz. Synchronization was similar for all synaptic models and depression levels, in good agreement with experimental data. However, there was a large difference between entrainment for synapses with different depression levels. The EI was unrealistically low for depression levels >50% between 400 and 800 Hz. The synapse with in vitro-like depression and double exponential recovery (Yang and Xu-Friedman, 2009) failed to drive GBCs above 400 spikes/s. Synapses with depression levels from 0 to 50% produced EI in the range of the reference data. Spirou et al. (2005) found a large variability (9-69) of ANF inputs converging onto the soma of individual GBCs. To examine the influence of the number of inputs on phase locking and entrainment, we analyzed different convergence patterns ranging from 20-50 inputs for a 10%-depressing synapse. Because "highsync" neurons are of the main interest for this study, we focused on bushy cells with average and high numbers of inputs. Figure 8 shows both the VS and EI for the simulated GBCs. We also included the original experimental data points (triangles and stars) from Joris et al. (1994) and simulated ANFs (dotted line) for reference. Figure 8 shows that the number of inputs has a strong effect on entrainment. High EI values for CFs up to 700 Hz could be reached only with a large number (50) of converging ANFs; smaller numbers of inputs caused a degradation of the EI. The change of entrainment due to a variation in the number of ANF inputs reflects the variability of the experimental data. The entire range of the measured data points was covered merely using a simple variation of the number of ANF inputs. It is further quantified in Section 3.3 where all data points were fitted to models by changing the number of ANF inputs. Compared to entrainment, the number of inputs had a smaller effect on phase locking, with a decreased number of ANF inputs leading to only a slight decrease of the VS.

Fitting Number of ANF Inputs and Synapse Types to Individual Data Points of EI
Previous simulations showed that VS does not change significantly for different synapse types (depression levels, Figure 7) or the number of ANF inputs (Figure 8). At the same time, EI strongly depends on both. To have better insight into the combined effect of the parameters, each synapse type was fitted to every experimental data point from the EI plot from Joris et al. (1994). The same data points are plotted in the lower panels of Figures 7, 8. The results of the fitting are shown in Figure 9 where colors encode the number of ANF inputs that was necessary to fit the data for each model. For clarity, only the frequencies in the transition region (CF > 500 Hz), in which EI degrades, was plotted. Data in the low frequency region (CF < 500 Hz) could be reproduced by any synapse model given enough ANF inputs, which is also visible in Figures 7, 8. The possible number of ANF inputs varied from 10-70 as observed by Spirou et al. (2005). They also discovered that the majority FIGURE 7 | Influence of synaptic depression on synchronization and entrainment. The points in the plots represent experimental data digitized from Joris et al. (1994). Each GBC was excited by 40 high SR ANFs. of GBCs receive 23 or less inputs. This constraint holds only for a few fitted GBCs in the transition region and is possible for synapses with small depression levels (0-20%). Other fitted GBCs required more inputs. Fields with white color mean that it was not possible to fit the model to the given EI at the given CF. Models with depression > 50% clould no longer reproduce data points with the highest CFs. The simulations also show that all (except one) EI data points could be reproduced by the less depressing models by varying the number of ANF inputs. SR of all modeled GBCs was set to 7.5 spikes/s. For completeness, Supplementary Material shows additional results for GBC models with higher SR.

DISCUSSION
GBCs play an important role in the sound localization pathway. They are highly specialized for that function because of fast channel kinetics and fast synapses. Together with multiple inputs from ANFs, such specialized features enable them to enhance the precision of temporal coding. Our investigation showed how to optimally adjust synaptic weights for endbulbs of Held with different levels of depression to produce realistic GBC responses including "high-sync" neurons. Additionally, simulations showed that synapses with depression greater than 50% can reproduce less data points than synapses with smaller depression. The depression is likely to be weaker (0-20%), FIGURE 8 | Influence of the number of ANF inputs on synchronization and entrainment. While the number of inputs only slightly influenced the VS, the effect on EI was strong. The variability of the experimental data can be explained in the model by the variation in the number of inputs to GBC. Reference data is from Joris et al. (1994).
FIGURE 9 | Fitting of number of ANF inputs for every synapse type to EI data point from Joris et al. (1994). Only data points with CF > 500 Hz are shown for clarity. Colors encode the number of ANF inputs that were necessary to fit the data. Each data point is represented as a tuple (CF, EI) on the horizontal axis. The EI data points are also shown in the lower panels of Figures 7, 8. A white field means that no fit was possible for a given synapse type.
because strongly depressing synapses require much more ANF inputs. Synapses with exponential recovery mechanisms fitted to in vitro experiments with paired-pulse stimuli did not operate properly with in vivo-like stimuli in the most of the transition region (CF > 500 Hz). It indicates that synaptic mechanisms which counteract vesicle depletion are activated in-vivo. It also shows that even though giant synapses have been studied in great detail before, it is still important to examine them as a part of a larger system, what provides further constrains. Finally, variation of the number of ANF inputs and synaptic depression can reproduce variability of entrainment observed in vivo. Last but not least, our modeled GBCs can generate realistic spike trains using a code that is freely available at https://github.com/mrkrd/ cochlear_nucleus (Rudnicki and Hemmert, 2017) for researchers who are interested in investigating neuronal processing at higher stages of the auditory pathway, especially in sound localization (Encke and Hemmert, 2015).

Adjustment of Synaptic Weights
Our model consisted of three main components: ANF inputs, endbulbs of Held and a GBC soma. The properties of ANF inputs and GBC soma had been well studied before. We used accurate biophysical models with fixed parameters from Zilany et al. (2014) and Rothman and Manis (2003c). The inner ear model of Zilany et al. (2014) allowed us to generate realistic in vivo-like ANF spike trains. GBC soma, despite being relatively realistic, still lacked features, such as multiple compartments and inhibitory inputs. In particular, inhibitory inputs might play an important role in GBCs as they do in spherical bushy cells (Nerlich et al., 2014;Keine and Rübsamen, 2015;Kuenzel et al., 2015). Since we wanted to examine a variety of endbulbs of Held, we decided to use phenomenological synapse models that allowed us to easily modify model parameters in a controllable way. Thus we could examine how synaptic properties influence input-output characteristic of the system, but did not examine the exact mechanisms of depression.
Synaptic weights were fitted separately for each combination of synapse type and the number of synapses (ANF inputs). The main objective was to achieve large VS and EI while keeping the SR realistic. While the SR of GBCs can vary a lot (Spirou et al., 1990;Smith et al., 1993;Rhode, 2008) observed correlation between SR and CF. They noticed that GBCs at low CFs tend to have low SR and GBCs at high CFs have higher spontaneous activity. Because our model was adjusted for GBCs with low CFs it might not be valid for GBCs with higher CFs and SRs. Additionally, Kopp-Scheinpflug et al. (2008) showed indirectly similar pattern, where principal neurons in MNTB with low CF tend to have lower SRs than neurons with high CFs for various species. Firing rates of MNTB neurons are related to GBCs, because they are directly driven by GBCs through calyxes of Held which work as relays (Borst and Soria van Hoeve, 2012). However, the SRs of MNTB neurons were higher than SRs of GBCs.
Because the VS and EI (as well as driven rate) are a nondecreasing function of synaptic weights (Figure 4), the SR was the limiting factor during the parameter scan. The SR was set to 7.5 spikes/s which was the average SR of GBCs with CF < 6 kHz from Spirou et al. (1990). The optimal synaptic weights for various numbers of ANF inputs and levels of synaptic depression are shown in Figure 6. The plot shows that the weights decrease with increasing number of ANF inputs to the asymptote of 0. Additionally, synapses with stronger depression have larger synaptic weights. The weights in the plot represent the initial synaptic conductance which is reduced during synaptic operation. Rothman et al. (1993) presented a study of a bushy cell model with tonic synapses that was strongly influential for the design of our model. The authors performed a detailed parameter space scan varying the number of inputs and synaptic strengths. Their model could reproduce all typical peristimulus time histograms of bushy cells (PL, PL N , On) and improve synchronization. The main differences with our study were lower number of ANF inputs (1-20) and consideration of tonic synapses only. The model of Rothman et al. (1993) most successfully captured the properties of spherical bushy cells, whereas our model focused on GBCs. Furthermore, Kuenzel et al. (2015) and Nerlich et al. (2014) showed advancements in modeling of spherical bushy cells which included inhibition and dynamic synapses to better explain in vitro and in vivo experiments. Spirou et al. (2005) presented a detailed biophysical model of GBCs. They combined structural data from electron microscopy with physiological data about vesicle release probability, synaptic depression, and receptor kinetics from Graham (2001). The model was tuned to closely reflect physiological data and successfully reproduced peristimulus time histograms of various GBCs. However, this model failed to improve the temporal precision for depressing synapses, with none of the simulated cells achieving a synchronization index larger than 0.9. In contrast, improved phase locking was one of the main goals of our model. We also examined various depression levels, but only for phenomenological synapse models.

Levels of Depression in In vivo-Like Simulations
We presented three different types of models of endbulb of Held: tonic, with single, and with double exponential recovery synapse. Tonic synapses were reported in some earlier in vivo studies of endbulbs (Young and Sachs, 2008;Borst, 2010;Kuenzel et al., 2011) and calyxes (Lorteije et al., 2009) of Held. They also represent a boundary case when analyzing depression, because they are equivalent to 0%-depressing synapses. Simple depressing synapses with single exponential recovery allowed us to study a whole range of depression levels (0-70%). Finally, depressing synapse with double exponential recovery represented an average synapse from in vitro experiments, where the data was fitted to a double exponential function, see Figure 3 (Yang and Xu-Friedman, 2009). All synapse models were purely phenomenological, i.e., they did not attempt to explain the mechanism of depression. It allowed us to easily manipulate synapse parameters and observe the influence of the depression on firing properties in in vivo-like stimulations.
A qualitative comparison of two synapses is shown in Figure 5. It shows synaptic currents of a weakly depressing (10%-depressing) and a strongly depressing (yang2009mean ,  Table 1) synapse during spontaneous and driven activity. First, the responses during the spontaneous activity look very similar. This is a direct result of the optimization where all synapses were fitted to have a specific SR. Second, the average amplitude of the synaptic currents is smaller compared to the 10%-depressing synapse. The peaks often do not reach the threshold necessary to generate a single action potential. As a result, there are many failures in generating action potentials. The strong difference in the driven region between both synapses is expected, because they represent the two extreme cases: one of the least and one of the most depressing synapses in our study. The other depression levels are expected to lie between them and are shown next. Figure 7 shows how different synapse types and levels of depression influence firing properties of GBCs in in vivo-like simulations. Interestingly, the depression level had little effect on phase locking. All GBC spike trains had VS greater than VS of ANF trains for frequencies <1 kHz, which is also observed in vivo. The results quantify how EI degrades with increasing synaptic depression. It is consistent with Figure 5 where the strongly depressing synapse was not able to drive GBCs to high firing rates. For the given number of ANF inputs (40), the maximum depression was approximately 50-60%. Stronger depression levels caused EI to be lower than the EI of the reference data. The interplay between the depression level and the number of ANF inputs is explained in the next paragraph. Figure 9 shows fitting results of every depression model to individual EI data points in the transition region (CF < 500 Hz) from Joris et al. (1994). The allowed number of ANF inputs was in the range of 10-70 (Spirou et al., 2005). If no fit was found, the corresponding field was left blank. Additional results including GBCs with higher SRs (20 spikes/s and 50 spikes/s) are shown in the Supplementary Material. First of all, in the low frequency region (not shown in Figure 9, but visible in Figure 7) all data could be reproduced by all examined synapses. It is consistent with in vitro observations of a large variability of the depressing synapses Yang and Xu-Friedman (2009). However, the double exponential recovery synapse, that was tuned explicitly to in vitro data, did not reach the desired entrainment for GBCs with CFs > 500 Hz, which could have two explanations: 1. The synaptic properties depend on the stimulus, e.g., a synapse might have different recovery time constants for paired-pulse (often used in slices) and ongoing (in vivo-like) stimuli as observed by Yang and Xu-Friedman (2015). 2. Physiological conditions in vitro and in vivo are responsible for different synaptic properties as suggested by Lorteije et al. (2009).
Therefore, for the in vitro synapse model to properly operate in in vivo conditions, requires either (a) adjustment of the synaptic parameters for the in vivo condition, or (b) a development of a more complex model (including the actual bio-physiological processes). This result also suggests that slice measurements of synaptic parameters cannot be easily adopted to in vivo conditions. The finding above is consistent with many in vivo observations of calyx of Held, where there was a small depression observed (Wang et al., 2013(Wang et al., , 2015Di Guilmi et al., 2014) or no evidence of synaptic depression (Lorteije et al., 2009). Also, Kuenzel et al. (2011) and Young and Sachs (2008) did not observe evidence of strong synaptic depression in vivo in endbulbs of Held. Hermann et al. (2009) observed reduction (but not abolition) of the synaptic depression in vitro using conditioning with spontaneous activity. They concluded that calyxes of Held operate in chronic synaptic depression in vivo. A similar explanation came from Friauf et al. (2015) where they proposed that in vivo synapses are in their normal stimulation range and the in vitro synapses are in a "manic" state due to lack of the spontaneous stimulation. Additionally, a possibly similar effect in vivo was shown by Di Guilmi et al. (2014) and Wang et al. (2015) where calyxes associated with low-SR (< 20 Hz) MNTB neurons had larger synaptic depression (23%) than high-SR neurons (8%). The difference between in vivo and in vitro synapses could be also partially explained by adjusting the Ca 2+ concentration. Yang and Xu-Friedman (2015) showed large changes in synaptic depression by varying Ca 2+ concentration from range 1 mmoll −1 to 2 mmoll −2 which resulted in depression varying from 0 to > 50%. Lorteije et al. (2009) observed reduction of the synaptic depression in vitro when varying Ca 2+ concentration. They also noted that deviations in Ca 2+ concentration was not enough to explain the difference between in vitro and in vivo synapses. Additionally, Yang and Xu-Friedman (2015) showed that not only the depression levels, but also time constant are different between in vitro and in vivo synapses, which suggest a different mode of operation in vitro and in vivo.

Number of ANF Inputs
In the simulations, we varied the number of ANF inputs (10-70) similarly to the observations of Spirou et al. (2005). Figure 8 show that on the one hand, the number of ANF inputs did not significantly influence VS in pure tone stimulations. On the other hand, EI was strongly influenced by the number of ANF inputs. Figure 9 shows most EI data points can be reproduced by varying the number of ANF inputs for synapses with depression smaller than 50%. In addition, Spirou et al. (2005) observed that the majority of GBC received less than 23 inputs, which could be reached only by some synapses in the transition region (CF > 500 Hz) with low synaptic depression (0-20%).

Physiological Consequences
The main finding of our study that the combination of strongly adapting inputs from ANFs with a biophysically plausible model of GBCs generates realistic output spike trains only when the endbulb of Held synapses show low depression levels has also important consequences on how this synapse works.
The synaptic terminals of the primary ANFs are remote from the soma. This distance limits and delays the delivery of new supplies to the synapse. For sustained operation terminals have to locally recycle their vesicles and neurotransmitter. GBCs are very specialized neurons, they can fire series of sustained and precisely phase-locked action potentials up to 700 Hz, despite decreasing (adapting) input from ANFs. The modeling results presented here show that this is only possible, if the endbulb of Held synapse has active mechanisms in place, which counteract depression. One mechanism is probably the high number of release sites (up to 150, Oleskevich et al., 2000;Nicol and Walmsley, 2002;Borst, 2010), in combination with a release probability, which has to be much smaller than the value of 0.6 determined in vitro. Only with low release probabilities vesicles can be spared such that the available vesicle pool size does not deplete too fast and too much during phases with high-frequency trains of presynaptic action potentials. Nevertheless, even small amount of vesicle depletion leads to depression and therefore compensating mechanisms have to be in place. Neher and T. Sakaba (2008) discuss that global intraterminal Ca 2+ concentration, which increases during phases of high activity, accelerates vesicle recruitment. A second mechanism was described by Hosoi et al. (2007) in the calyx of Held, where they concluded that during a 100 Hz stimulus train the fusion probability increased during the initial 5-10 APs by a factor of five. If both factors compensate pool depletion, synaptic depression can be more or less suppressed such that GBCs can fire fast sequences of action potentials with little fatigue. For GBCs, this behavior is essential to shape and relay sustained and precisely phase-locked high-frequency pulse trains, which are required to evaluate interaural time and level differences in the lateral and medial superior olives.

AUTHOR CONTRIBUTIONS
MR and WH are listed as authors according to the guidelines of Frontiers in Computational Neuroscience.

FUNDING
This work was supported by the German Research Foundation (DFG) within the Priority Program "Ultrafast and temporally precise information processing: normal and dysfunctional hearing" SPP 1608 (HE6713), the Technische Universität München within the funding programme Open Access Publishing, and the German Federal Ministry of Education and Research within the Munich Bernstein Center for Computational Neuroscience (reference number 01GQ1004B).