Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Minimal Conductance-Based Model of Auditory Coincidence Detector Neurons

Abstract

Sound localization is a fundamental sensory function of a wide variety of animals. The interaural time difference (ITD), an important cue for sound localization, is computed in the auditory brainstem. In our previous modeling study, we introduced a two-compartment Hodgkin-Huxley type model to investigate how cellular and synaptic specializations may contribute to precise ITD computation of the barn owl's auditory coincidence detector neuron. Although our model successfully reproduced fundamental physiological properties observed in vivo, it was unsuitable for mathematical analyses and large scale simulations because of a number of nonlinear variables. In the present study, we reduce our former model into three types of conductance-based integrate-and-fire (IF) models. We test their electrophysiological properties using data from published in vivo and in vitro studies. Their robustness to parameter changes and computational efficiencies are also examined. Our numerical results suggest that the single-compartment active IF model is superior to other reduced models in terms of physiological reproducibility and computational performance. This model will allow future theoretical studies that use more rigorous mathematical analysis and network simulations.

Introduction

Sound localization, or the ability to find the location of a sound source, is one of the fundamental functions of the auditory system. The interaural time difference (ITD), which is the difference in arrival timing of acoustic waves at the two ears, is an important cue for azimuthal sound localization (see [14] for reviews). In mammals, ITD is computed in the medial superior olive (MSO), where inputs from the bilateral cochlear nuclei converge [58]. In birds and reptiles, binaural neurons in the third order auditory station, nucleus laminaris (NL), show sensitivity to ITDs [912]. Recent advances in electrophysiological recording techniques in vivo allow us to examine how binaural synaptic inputs are computed in these neurons [13,14]. Since MSO and NL neurons sense precise timing of synaptic inputs to change their spike rates, they are often called auditory coincidence detectors.

Conductance-based Hodgkin-Huxley-type (HH) models have been widely used for investigating biophysical mechanisms of ITD coding in auditory coincidence detector neurons [1521]. These modeling studies have revealed that finely-tuned synaptic and cellular properties are essential for precise ITD computation. The HH-type model is useful to study detailed biophysical effects of channel kinetics, but it requires a large computational power [22]. Furthermore, mathematically rigorous analysis of HH-type models is generally difficult due to the large number of parameters and nonlinearities. To circumvent these difficulties, more abstract integrate-and-fire (IF) models with fewer parameters were also used in previous theoretical studies of sound localization [19,2325], although the biophysical basis of the IF model is much weaker than the HH model [22,26]. Combining the strengths of the both models, Svirskis and Rinzel [27] introduced an IF model with low-voltage-activated potassium (KLVA) conductance to study how KLVA channels contribute in the detection of weak signals in the model MSO neuron.

In our previous studies, we used HH-type two-compartment models that mimicked physiological characteristics of the barn owl's NL neuron [13,28,29]. Although our models were able to explain the roles of various biophysical factors in ITD computation, it was still unclear which model features were essential for binaural sound information processing. Thus, in the present study, we reduce our two-compartment NL model into three simplified models and compare simulation results with available physiological data. The main goal of this study is to find a minimal NL neuron model that reproduces known fundamental electrophysiological properties related to ITD coding. Numerical efficiency and reliability of the models are also examined.

Materials and Methods

Model construction

Anatomical [30,31] and physiological [13] evidence show that spikes in an avian NL neuron are initiated at a remote site from the cell body, presumably at the first node of Ranvier. Our previous NL neuron model [28,29] consisted of two compartments (Fig 1A). The somatic compartment receives synaptic inputs, whereas the nodal compartment serves as a spike generator. We refer to this model as the "Active Na" model, because action potentials are generated by the Na conductance in the nodal compartment. In this study, we reduced the active Na model in three steps. First, active conductances in the nodal compartment were replaced with a thresholding unit (shown by Θ in Fig 1E1; see next section for equations). Active and passive properties of the soma remain unchanged between these two models. This reduced model was named "Two-compartment active IF" model. The word "active" indicates that the model reserves low-voltage-activated (KLVA) conductance in the soma.

thumbnail
Fig 1. NL models.

(A) Circuit of the conductance-based two-compartment HH-type "active Na" model [13,29]. The somatic compartment contains leak and KLVA currents whereas the nodal compartment has high voltage activated potassium (KHVA) and Na currents (required for spike generation) in addition to leak and KLVA. (B) Circuit of the "non-spiking" conductance-based single-compartment model with leak and KLVA currents. This model was designed for simulating subthreshold membrane responses [13,32]. (C) Membrane impedance of the somatic compartment of the active Na model. (D) Membrane impedance of the nodal compartment of the active Na model. (C-D) In the "RC only" and "RC + KLVA" conditions, gaxon was fixed to zero (i.e., no axonal current) to isolate the compartment. In the "RC only" and "RC + axon" conditions, gKLVA of the compartment was fixed to zero. Membrane potential was fixed at -61 mV. (E1-3) Two-compartment active IF model: (E1) circuit diagram, (E2) average spike shape at the soma, and (E3) membrane responses to DC step current injection. Spike amplitude H = 9.7 mV. Half-amplitude spike width W = 0.3 ms. (F1-3) Single-compartment active IF model: (F1) circuit diagram, (F2) average spike shape, and (F3) membrane responses to DC step current injection. Spike amplitude H = 9.8 mV. Half-amplitude spike width W = 0.3 ms. (G1-3) Single-compartment active IF model: (G1) circuit diagram, (G2) average spike shape, and (G3) membrane responses to DC step current injection. Spike amplitude H = 9.8 mV. Half-amplitude spike width W = 0.3 ms. (E1, F1, G1) Θ denotes the IF spike generator. (E3, F3, G3) Holding potential = -60 mV.

https://doi.org/10.1371/journal.pone.0122796.g001

As the second step, we replaced the entire nodal compartment and the axon with a single thresholding unit in the soma (Fig 1F1; see below for equations). This new model was labeled as "Single-compartment active IF model". This model can also be regarded as a simple expansion of the "Non-spiking" NL model (Fig 1B), which was introduced in our previous research to study subthreshold membrane responses [13,32].

As the last step, we dropped the KLVA conductance from the soma and obtained the "Single-compartment passive IF" model (Fig 1G1). Similar models to this passive IF model (with slightly different membrane parameters) were already used in previous theoretical studies of sound localization [23,24].

Two-compartment active IF model.

The two-compartment active IF model consists of the cell body (soma) and the node, which are connected with an axonal resistance (Fig 1E1). The membrane potentials of the soma (Vsoma) and node (Vnode) obey the standard membrane equations: and where Csoma and Cnode are the membrane capacitances, are the leak currents, is the KLVA current, is the axonal current, is the synaptic input to the soma, IIF is the current induced by the IF thresholding unit (defined below), and Iext is the external stimulus current (set to zero by default). The activation variable d(V,t) of the KLVA conductance in the soma obeys the first-order differential equation: The KLVA kinetics was taken from a previous slice study [33]: with Q10 = 2.5, T0 = 23 oC, T1 = 40 oC. The unit for αd and βd is 1/ms.

Parameter values are the same as our previous active Na model [29] and summarized in Table 1. The KLVA conductance in the soma was unchanged between the active Na model and the active IF model, because it plays an essential role in characterizing the membrane impedance (Fig 1C). The KLVA conductance in the node, however, was eliminated from the present active IF model, because it has only limited effects on the nodal membrane impedance, which is dominated by the axonal conductance (Fig 1D).

thumbnail
Table 1. Default parameter values for the membrane dynamics of the three reduced models.

https://doi.org/10.1371/journal.pone.0122796.t001

The current from the thresholding unit consists of two parts: IIF = Iconst + Ispike. The first term corresponds to the constant opening of ion channels in the nodal compartment. The second term Ispike denotes the spike-associated transient current: when the nodal membrane potential Vnode crosses the threshold Vθ at time t = Tθ, a spike current Ispike(t-Tθ) is initiated. The spike current is modeled as a sum of two exponential functions:

We chose exponential functions because they enable exact calculation at discrete time steps [34]. Once the membrane potential reaches the threshold Vθ, the threshold-crossing detector will be in an absolute refractory period of Tref. Namely, Vθ = ∞ for Tθ < t < Tθ+Tref. See Table 2 for the default parameter values of the thresholding unit, which were fixed throughout our series of simulations unless otherwise stated.

thumbnail
Table 2. Default parameter values for the thresholding unit of the three reduced models.

https://doi.org/10.1371/journal.pone.0122796.t002

Single-compartment active IF model.

In the single-compartment IF model, the spike-initiating node of the two-compartment model together with the axonal resistance was reduced into a thresholding unit in the soma (Fig 1F1), while other ionic conductances of the soma were unchanged. The membrane potential of the soma (Vsoma) satisfies the equation:

The current from the IF unit is the same as in the two compartment model: IIF = Iconst + Ispike. The equations and parameter values for Csoma, , IKLVA, Isyn, Iconst, and Iext are the same as those in the two-compartment active IF model introduced above (see Table 1 for parameters). Thus the subthreshold membrane property is essentially identical to that of the two-compartment model.

The spike initiation mechanism is similarly defined. When the membrane potential Vsoma reaches the threshold Vθ at time t = Tθ, a spike current Ispike(t-Tθ) is injected. Note that, in the single-compartment model, spike currents are applied directly to the soma. The spike current Ispike(t) is similarly defined as above with a set of slightly different parameters. The definition of the absolute refractory period Tref was the same as that of the two-compartment model (see Table 2 for parameters).

Single-compartment passive IF model.

By replacing the active KLVA conductance with passive leak conductance, we obtain the single-compartment passive IF mode (Fig 1G1). The membrane potential of the soma (Vsoma) obeys the passive equation:

The equations and parameter values for Csoma, , Isyn, Iconst, and Iext are the same as those in the above-defined models except for the leak conductance . The definitions for thresholding and spike initiation were also the same as in the single-compartment active IF model.

Membrane and threshold parameters.

We used the same parameter sets for the soma of the active IF models (Table 1) as used for our previous active Na model [29]. With the parameters shown in Table 1, somatic membrane resistance at -61 mV was 4.4 MΩ (membrane time constant: ~0.1 ms), which was comparable to the electrophysiological data from the owl NL [13]. In the passive IF model, the somatic leak conductance was increased so as to keep the membrane resistance. Parameters for the spike-associated current (Table 2) were determined so that simulated spike waveforms resembled what was observed in vivo [13]. To date, there is no report available that systematically studied the voltage threshold Vθ and the refractory period Tref of owl NL neurons. We determined the threshold of each model (Table 2) to make the "modulation depth" (defined later) of the simulated spike rates close to maximum. The refractory period was fixed to 0.9 ms. Possible effects of varying these parameters will be examined in Results.

Synaptic inputs.

An NL neuron receives phase-locked inputs from converging axons of nucleus magnocellularis (NM) neurons on ipsi- and contralateral sides [9]. We used the same modeling procedure that was introduced in our previous study [11,32]. Briefly, phase-locked spikes of NM axons were modeled as an inhomogeneous Poisson process, whose intensity function λ(t) is periodic with the tonal stimulus frequency fstim. The degree of phase-locking was measured as vector strength [5], which was analytically related to the intensity function λ(t) [32]. Presynaptic spikes (modeled as delta functions) from each side were processed by a linear synaptic filter: to obtain synaptic conductances gipsi(t) and gcontra(t). We assumed that all inputs from each side are locked to the same phase at the same frequency, and thus the ITD only affects the relative phase δ of the model inputs from the two sides:

We used the same input parameters as in our previous modeling study [29,32]: fstim = 4000 (Hz), average intensity λ0 = 500 Hz [35], vector strength r = 0.6 [36], number of NM inputs from each side M = 150 fibers [37], half amplitude width of the unitary synaptic input Wα (= 2.446 τα) = 0.1 ms, amplitude of unitary synaptic input Hα = 1.3 nS. All synaptic inputs were assumed to converge onto the somatic compartment, because the owl's mid-to-high frequency NL neurons have only short and stubby dendrites [37].

Model evaluation

In order to evaluate the physiological reproducibility of the model neurons, we compare model outputs with published electrophysiological data. We examined spike waveforms, response types to step current injections, and spike rate modulations to simulated binaural inputs. To evaluate the numerical properties of the models, we examined their robustness to parameter changes and computational efficiency. For comparing membrane potential traces (Fig 2) we used the same data set from our previous in vivo intracellular recording study of the barn owl's NL [13].

thumbnail
Fig 2. Membrane potential traces.

(A) In vivo intracellular recording from a barn owl's NL neuron (data taken from [13]). Tonal stimulus at the best frequency (3600 Hz) was presented to the both ears by changing ITDs. (B) Simulated membrane potential of the two-compartment active IF model. (C) Simulated membrane potential of the single-compartment active IF model. (D) Simulated membrane potential of the single-compartment passive IF model. (B-D) Simulated binaural inputs are injected to the model neurons (see Materials and Methods for detail). "In-phase" means that the oscillation peaks of both inputs coincide resulted in the largest oscillation amplitudes. "Out-of-phase" means that the oscillating bilateral synaptic inputs are cancelled with each other to minimize the signal component at fstim = 4000 Hz [29].

https://doi.org/10.1371/journal.pone.0122796.g002

Spike shape.

For the calculation of average spike shapes (Fig 1E2, 1F2 and 1G2), we applied noise synaptic input to the model neurons. Random presynaptic inputs were modeled as non-phase-locked spike sequences (i.e., homogeneous Poisson process with an intensity λ0 of 500 Hz and a number M of 300 fibers), and were convolved with the synaptic filter α(t) to obtain a simulated synaptic input gsyn(t). Spikes are aligned to the threshold crossing point and then averaged. 1500 spikes were estimated to obtain an average spike waveform. Data from our previous in vivo intracellular recording [13] were used as a reference (typical values: spike amplitude = 4–13 mV, spike width = 0.3–0.5 ms).

Response to step currents.

In order to see the model membrane response to step current injection Istep (Fig 1E3, 1F3, and 1G3), we hold the membrane potential of each model at -60 mV by a constant current Ibase, and additionally applied step currents of varied amplitudes: Iext = Ibase + Istep. The model response was classified as "no spike" when the current injection did not induce an action potential. The response was categorized as "phasic spiking" when only a single spike at the onset of the step current was observed. The response was referred to as "tonic spiking" when repetitive spikes occurred during the step current injection. We examined whether the model neurons show a phasic-spiking property as is typical in auditory coincidence detectors (e.g., [38]).

Responses to binaural inputs and rate-phase difference curve.

We changed the phase difference δ of the bilateral model inputs (see above for definition) and calculated traces (Fig 2) and a rate-phase difference curve (Fig 3A1, 3B1, and 3C1) for each reduced model. The maximum spike rate is called the "in-phase rate", because it is attained when δ was zero or integer multiples of 2π [29]. The minimum rate, obtained with δ = (2N+1)π (N = 0, ±1,…), was referred to as the "out-of-phase" rate. Note that the simulated in-phase rate corresponds to the spiking rate with a "favorable ITD" in physiological experiments, whereas the out-of-phase rate is comparable to the rate with an "unfavorable ITD" [13]. Discharge rates of owls' NL neurons for favorable and unfavorable ITDs were systematically measured by Peña et al. [35]. Based on their results, we defined the "typical discharge range" as 79–522 spikes/sec, which was derived from the "mean-SD" of the discharge rate for unfavorable ITD and "mean+SD" for favorable ITD. We used this typical discharge range as a criterion to evaluate the binaural responses of the models, and targeted to fit in this range when determining the thresholding parameters.

thumbnail
Fig 3. Spiking properties of the reduced models.

(A1-3) Two-compartment active IF model. (B1-3) Single-compartment active IF model. (C1-3) Single-compartment passive IF model. First row (A1,B1,C1): simulated spike rates plotted against the phase difference δ of the bilateral inputs. Peaks appear at δ = 0 or ±2π (in-phase rate) whereas troughs appear at ±π (out-of-phase rate). Second row (A2,B2,C2): simulated spike rates plotted against the threshold of the IF unit. Refractory period was fixed to 0.9 (ms). Modulation depth (gray lines) is defined as the difference between the in-phase and out-of-phase rates (black lines). Dotted lines show the criterion of 180 spikes/sec for the modulation depth. Third row (A3,B3,C3): simulated spike rates plotted against the length of the refractory period. Threshold was fixed to (A3) -56.7, (B3) -58.3, and (C3) -58.6 mV, respectively. Shaded area of each panel in the first and third rows shows the range of typical discharge rates of the owl's NL neuron (see Materials and Methods for the definition).

https://doi.org/10.1371/journal.pone.0122796.g003

Robustness to parameter changes.

To examine the robustness to changes in thresholding parameters, we changed the threshold Vθ and the refractory period Tref of each IF model around the default values. We tested whether the in-phase and out-of-phase rates were within the typical discharge range. We additionally tested whether the "modulation depth" (defined as the difference between the in-phase and out-of-phase rates) exceeded our criterion of 180 spikes/sec, which is derived from the spike rate changes between favorable and unfavorable ITDs [35].

Numerical efficiency and reliability.

As in our previous studies [28,29,32], we used the explicit (forward) Euler method for the numerical integration of the model equations. We varied the time step Δt from 0.1 μs to 30.0 μs. To examine the model efficiency, we calculated an average integration time of ten 1-second traces (i.e., 10 s in total) for each model. To test the numerical reliability, we compared simulated membrane potentials with different values of Δt. A numerical result was regarded as "unreliable" when the maximum discrepancy in the simulated potential differed by more than 0.5 mV from the potential calculated with the smallest time step of Δt = 0.1 μs. Note that this unreliability usually occurs at each spike initiation (see Fig 4). Numerical algorithms were implemented in D [39] and simulations were performed on a desktop computer (Dell Precision T1700) with 64 bit Windows 7 Professional Operating System, Intel Xeon CPU E3-1270 v3 (4 core, 3.5 GHz), and a 16 GB memory.

thumbnail
Fig 4. Model efficiency and reliability.

(A) Relative integration time of the active Na (upward triangles), two-compartment active IF (downward triangles), single-compartment active IF (circles), non-spiking (diamonds), and single-compartment passive IF (squares) models. Single-compartment active IF model with Δt = 10.0 μs (small arrow) was used as the reference point (i.e., relative integration time = 1.0). Open symbols denote numerical unreliability with corresponding panels indicated by the small letters. (B) Unreliability of the two-compartment active IF model. A large time step (Δt = 3 μs, thin black line) leads to an erroneous oscillation in the nodal potential induced by spike currents. Inset shows a magnified view of the spike-induced unreliability, which does not occur with a sufficiently small time step (e.g., Δt = 1 μs, gray line). (C) Unreliability of the single-compartment active IF model. Too large a time step (Δt = 20 μs in this example) leads to an imprecise computation of the somatic shape. (D) Unreliability of the single-compartment passive IF model. Too large a time step (Δt = 20 μs in this example) leads to an imprecise computation of the spike shape. Small arrows in C and D show the timings when the simulated traces with Δt = 20 μs started to diverge from the traces with Δt = 1 μs by more than 0.5 mV.

https://doi.org/10.1371/journal.pone.0122796.g004

Results

Basic properties of reduced models

As introduced in Materials and Methods, we stepwise reduced the active Na model (Fig 1A) into three simpler models: the two-compartment active IF model (Fig 1E1), the single-compartment active IF model (Fig 1F1), and the single-compartment passive IF model (Fig 1G1). Using a similar approach as Svirskis and Rinzel [27], spike generating conductances in the active models were replaced by a threshold-crossing detector, while the KLVA conductance of the cell body was kept intact. The values for KLVA and leak conductances were chosen so that the membrane resistance matched previous in vivo measurements [13].

Spike shapes of the reduced models are shown in Fig 1E2, 1F2 and 1G2. We selected the parameters for spike-associated currents Ispike so that the simulated spike amplitude (9.7–9.8 mV) and half-peak width (0.3 ms) of all three models corresponded to in vivo recording data [13]. For the passive model, a small effect of the refractory period was seen in the spike trace (arrow in Fig 1G2), which was not evident for the active models.

Due to the large KLVA conductance [40], the active models showed phasic spiking to step current injections (Fig 1E3 and 1F3), as was found in auditory brainstem neurons in vitro (chicken NL [41]; chicken NM [33,42]; mammalian MSO and octopus cells [38]). As shown previously [40,43], removing the KLVA current resulted in the lack of phasic spiking mode with a direct transition from no spike to tonic spiking (Fig 1G3).

Binaural phase-coding

Our previous in vivo intracellular recordings from owls NL [13] showed that phase-locked synaptic input induces an oscillating membrane potential. The oscillation amplitude depends on ITD (Fig 2A; see [13] for detailed discussion). Similar oscillations were found in the non-spiking model (Fig 1B; see [29,32] for theoretical examinations). All three reduced models, which have similar subthreshold response properties to the non-spiking model, also show these oscillations (Fig 2B–2D). The amplitude of the oscillation was maximal for δ = 0 (i.e., binaural inputs arrive in-phase) and minimal for δ = π (i.e., binaural inputs arrive out-of-phase) [29,32].

Spikes in auditory coincidence detector neurons are in general small at the cell body (gerbil MSO [44]; chicken NL [31]; owl NL [13]) probably because the actual spike initiation site is located away from the soma [28,30]. Reflecting this property, both measured (Fig 2A) and simulated (Fig 2B–2D) waveforms showed sinusoidal oscillations at the falling phase of spikes. It should be noted that, in the single compartment models, the spike current mimics the backpropagating spikes from the spike initiation site.

We changed the relative timing of ipsi- and contralateral inputs and calculated the spike rates of the model neurons. The maximum (in-phase) and minimum (out-of-phase) spike rates of the two active models were within the typical discharge rates of 79–522 spikes/sec (gray areas of Fig 3A1 and 3B1; see Materials and Methods for the definition) derived from previous owl NL recordings in vivo [35]. The single compartment passive IF model, however, showed a steeper rate-phase difference curve with the in-phase rate exceeding the typical range (Fig 3C1).

Effect of parameter changes

In order to examine the parameter dependence of the model, we varied the parameter values of the thresholding unit (Fig 3A2, 3B2 and 3C2). First we varied the threshold Vθ. Both the in-phase and out-of-phase rates decreased with increasing threshold (Fig 3A2 and 3B2, black lines). The spike modulation depth, which was defined as the difference between the in-phase and out-of-phase rates, showed a mild peak around the default threshold of the active models (Fig 3A2 and 3B2, gray lines). This qualitative tendency was observed in all the reduced models, but the single-compartment passive IF model was most strongly affected by threshold changes, showing steepest changes in the modulation depth (Fig 3C2). The range of thresholds for which the modulation depth exceeds our criterion of 180 spikes/sec was very similar for all the models (about 2.2 mV).

Next we changed the absolutely refractory period Tref (Fig 3A3, 3B3 and 3C3). The length of the refractory period generally affects the maximum spike rate of an IF model (e.g., [45]). Since owls' NL neurons sometimes show discharge rates exceeding 600 spikes/sec [9] or even higher [13], we restricted the value of Tref below 1.6 ms. When the refractory period exceeded 0.8 ms, it affected only moderately the in-phase-rate of the two-compartment IF model and almost negligibly the out-of-rate rate (Fig 3A3). If the period was shortened below 0.6 ms, both rates rapidly blew up to an unreasonable range. The refractory effect was even milder for the single-compartment active model (Fig 3B3). The out-of-phase spike rate of the passive model, however, was more sensitive to refractory periods than the other two active models (Fig 3C3). For example, varying Tref from 0.8 to 1.6 ms lead to changes in the in-phase rate of 172, 124, and 310 spikes/sec in two-compartment active, single-compartment active, and single-compartment passive models, respectively.

Efficiency and reliability

In order to examine the numerical efficiency of the models, we calculated their relative integration times (Fig 4A). The two-compartment active IF model can be computed roughly three times faster than the active Na model, because of its much smaller number of active currents. The allowable time step Δt for the two-compartment active IF model, below which numerical calculation was reliable, was ten times larger than that for the active Na model. For Δt ≥ 3 μs, numerical results show oscillations at the spike initiation leading to a computational unreliability (Fig 4B; see Materials and Methods for the definition of the reliability).

Numerical integration of the single-compartment active IF model was about 10% faster than the two-compartment model (Fig 4A). A substantially larger time step Δt than the two-compartment model was accepted, because the single-compartment model lacks the node which is small and more vulnerable to spike-induced currents. The maximum allowable time step was about 15 μs. Exceeding this criterion resulted in unreliable calculation of spike waveforms (Fig 4C). The computational efficiency of the passive IF model was roughly three times better than the other two active models and the non-spiking model (Fig 4A). However, the maximum time step, which was determined by the reliability of the numerical calculation (Fig 4D), was similar to that of the single-compartment active IF model.

Comparison of three reduced models

Starting from the active Na model (Fig 1A), we constructed three reduced models step-by-step. The two-compartment active IF model (Fig 1E1) and the single-compartment active IF model (Fig 1F1) showed similar physiological properties to previous experimental results. The single-compartment active IF model was slightly more robust to parameter changes, and allowed five times larger time step with a subtly better computational speed. Our results are summarized in Table 3.

thumbnail
Table 3. Summary of tested physiological reproducibility and computational performance of the three reduced models.

https://doi.org/10.1371/journal.pone.0122796.t003

The active KLVA conductance was found to reduce low frequency components of the input [28,32] and to promote coincidence detection [16,40]. Since the single-compartment passive IF model (Fig 1G1) lacks this conductance, its response properties differed substantially from the two active models. Changes in inputs (Fig 3C1) and parameters (Fig 3C2 and 3C3) more severely affect the coding property of the passive IF model.

In conclusion, we suggest the single-compartment active IF model as the minimal conductance-based model that satisfies fundamental electrophysiological properties of auditory coincidence detectors and that showed the best robustness to threshold parameter changes.

Discussion

Technical limitations

In any neuronal modeling, there is always a compromise between specification and abstraction [26]. In this study we focused on finding the minimal configuration required for simulating known physiological functions of the owl's coincidence detector neuron. We introduced three different IF models that were reduced from the HH-type conductance based model. Despite its name "integrate", an IF model with a very short time constant behaves not like an integrator but rather like a coincidence detector [46]. As far as we tested, the single-compartment active IF model showed the best results in terms of physiological reproducibility and computational performance. The active KLVA conductance in our models contributed in stabilizing the model responses to parameter changes.

Because we replaced Na conductance with a thresholding unit, the reduced IF models do not show Na-related phenomena such as an improvement of small signal detection [16,27]. Reducing the number of parameters, however, may allow for more rigorous mathematical treatments (e.g., analysis on IF models receiving periodic inputs [24,4750]; two-dimensional phase-plane analysis [19,27]). Applying these mathematical techniques to our simplified models will be a future research subject.

Implementation of a neuronal model always requires a careful selection of the time step Δt, because it affects the numerical accuracy. Reduced models generally allow larger time steps (Fig 4A). The single-compartment active IF model was compatible with a maximum time step of 15 μs, when the explicit Euler method was used. Our previous study suggested that, to quantify phase-locking, the sampling frequency should be at least 10 (ideally 20) times larger than the locking frequency [51]. Therefore, a time step of 10–15 μs would be the lowest precision allowed for systems like the owls' auditory brainstem, which shows phase-locking up to about 8 kHz [36]. Nevertheless, if a larger time step is required, an implicit or exponential method may need to be considered [52].

Possible applications

Our simulations showed that it is not always necessary to include all active ionic currents (such as fast Na and delayed rectifier K) to properly simulate auditory coincidence detector neurons. This suggests that functional roles of various ion channels (e.g., [53]) may be separately investigated by restricting the number of ion channels combined with the IF-based model. Further physiological investigations would be necessary for evaluating the applicability of the reduced models for these purposes.

Because of the computational portability, IF models were often used for simulating a network of neurons (e.g., [5456]). In the auditory brainstem of mammals and birds, a huge extracellular field potential called the neurophonic is observed [57,58]. The amplitude of the neurophonic often exceeds 1 mV, but its functional significance (if any) is mostly unknown. Our reduced models may be useful for examining these collective effects of a large number of brainstem neurons. For such cases, the two-compartment active IF model, which showed less computational performance than the single-compartment model, might also be of some use, because the shape of a neuron affects the formation of the extracellular potential [58,59].

Both in birds and mammals, the auditory brainstem response (ABR) has been widely used as non-invasive diagnostic and research tools. Various auditory stations differentially contribute to the peaks of ABR waveforms [60,61], but how the underlying neural activity is related to the formation of ABR remains to be investigated. Future studies using a network of simplified brainstem neuron models may contribute in identifying the mechanisms of these sound-induced electric responses.

In summary, our results suggest that a single-compartment active IF model with relatively small number of parameters can simulate auditory coincidence detectors that sense submillisecond differences of binaural inputs. Our reduced models would serve as a useful tool for a theoretical investigation on the collective functions of these neurons.

Acknowledgments

The authors thank Catherine E. Carr for stimulating discussion.

Author Contributions

Conceived and designed the experiments: GA KF JK. Performed the experiments: GA. Analyzed the data: GA. Contributed reagents/materials/analysis tools: GA KF. Wrote the paper: GA KF JK.

References

  1. 1. Konishi M. Listening with two ears. Sci Am. 1993; 268: 66–73. pmid:8446882
  2. 2. Joris P, Yin TCT. A matter of time: internal delays in binaural processing. Trends Neurosci. 2007; 30: 70–78. pmid:17188761
  3. 3. Grothe B, Pecka M, McAlpine D. Mechanisms of sound localization in mammals. Physiol Rev. 2010; 90: 983–1012. pmid:20664077
  4. 4. Ashida G, Carr CE. Sound localization: Jeffress and beyond. Curr Opin Neurobiol. 2011; 21: 745–751. pmid:21646012
  5. 5. Goldberg JM, Brown PB. Response of binaural neurons of dog superior olivary complex to dichotic tonal stimuli: Some physiological mechanisms of sound localization. J Neurophysiol. 1969; 32: 613–636. pmid:5810617
  6. 6. Yin TCT, Chan, JCK. Interaural time sensitivity in medial superior olive of cat. J Neurophysiol. 1990; 64: 465–488. pmid:2213127
  7. 7. Brand A, Behrend O, Marquardt T, McAlpine D, Grothe B. Precise inhibition is essential for microsecond interaural time difference coding. Nature. 2002; 417: 543–547. pmid:12037566
  8. 8. Pecka M, Brand A, Behrend O, Grothe B. Interaural time difference processing in the mammalian medial superior olive: the role of glycinergic inhibition. J Neurosci. 2008; 28: 6914–6925. pmid:18596166
  9. 9. Carr CE, Konishi M. A circuit for detection of interaural time differences in the brain stem of the barn owl. J Neurosci. 1990; 10: 3227–3246. pmid:2213141
  10. 10. Klump GM. Sound Localization in Birds. In: Dooling et al. editors. Comparative Hearing: Birds and Reptiles. New York: Springer; 2000. pp 249–307.
  11. 11. Köppl C, Carr CE. Maps of interaural time difference in the chicken’s brainstem nucleus laminaris. Biol Cybern. 2008; 98: 541–559. pmid:18491165
  12. 12. Carr CE, Soares D, Smolders J, Simon JZ. Detection of interaural time differences in the alligator. J Neurosci. 2009; 29: 7978–7982. pmid:19553438
  13. 13. Funabiki K, Ashida G, Konishi M. Computation of interaural time difference in the owl’s coincidence detector neurons. J Neurosci. 2011; 31: 15245–15256. pmid:22031870
  14. 14. van der Heijden M, Lorteije JA, Plauška A, Roberts MT, Golding NL, Borst JGG. Directional hearing by linear summation of binaural inputs at the medial superior olive. Neuron. 2013; 78: 936–48. pmid:23764292
  15. 15. Brughera AR, Stutmen ER, Carney LH, Colburn HS. A model with excitation and inhibition for cells in the medial superior olive. Aud Neurosci. 1996; 2: 219–233. pmid:9160626
  16. 16. Svirskis G, Kotak V, Sanes DH, Rinzel J. Enhancement of signal-to-noise ratio and phase locking for small inputs by a low-threshold outward current in auditory neurons. J Neurosci. 2002; 22: 11019–11025. pmid:12486197
  17. 17. Grau-Serrat V, Carr CE, Simon JZ. Modeling coincidence detection in nucleus laminaris. Biol Cybern. 2003; 89: 388–396. pmid:14669019
  18. 18. Zhou Y, Carney LH, Colburn HS. A model for interaural time difference sensitivity in the medial superior olive: interaction of excitatory and inhibitory synaptic inputs, channel dynamics, and cellular morphology. J Neurosci. 2005; 25: 3046–3058. pmid:15788761
  19. 19. Dodla R, Svirskis G, Rinzel J. Well-timed, brief inhibition can promote spiking: postinhibitory facilitation. J Neurophysiol. 2006; 95: 2664–2677. pmid:16551843
  20. 20. Mathews PJ, Jercog PE, Rinzel J, Scott LL, Golding NL. Control of submillisecond synaptic timing in binaural coincidence detectors by KV1 channels. Nature Neurosci. 2010; 13: 601–611. pmid:20364143
  21. 21. Khurana S, Remme MWH, Rinzel J, Golding NL. Dynamic interaction of Ih and IK-LVA during trains of synaptic potential in principal neurons of the medial superior olive. J Neurosci. 2011; 31: 8936–8947. pmid:21677177
  22. 22. Izhikevich EM. Which model to use for cortical spiking neurons? IEEE Trans Neural Netw. 2004; 15: 1063–1070. pmid:15484883
  23. 23. Gerstner W, Kempter R, van Hemmen JL, Wagner H. A neuronal learning rule for sub-millisecond temporal coding. Nature. 1996; 383: 76–78. pmid:8779718
  24. 24. Kempter R, Gerstner W, van Hemmen JL, Wagner H. Extracting oscillations: neuronal coincidence detection with noisy periodic spike input. Neural Comput. 1998; 10: 1987–2017. pmid:9804669
  25. 25. Fischer BJ, Steinberg LJ, Fontaine B, Brette R, Peña JL. Effect of instantaneous frequency glides on interaural time difference processing by auditory coincidence detectors. Proc Natl Acad Sci. 2011; 108: 18138–18143. pmid:22006305
  26. 26. Herz AVM, Gollisch T, Machens CK, Jaeger D. Modeling single-neuron dynamics and computations: a balance of detail and abstraction. Science. 2006; 314: 80–85. pmid:17023649
  27. 27. Svirskis G, Rinzel J. Influence of subthreshold nonlinearities on signal-to-noise ratio and timing precision for small signals in neurons: minimal model analysis. Network: Comput Neural Syst. 2003; 14: 137–150.
  28. 28. Ashida G, Abe K, Funabiki K., Konishi M. Passive soma facilitates submillisecond coincidence detection in the owl's auditory system. J Neurophysiol. 2007; 97: 2267–82. pmid:17135480
  29. 29. Ashida G, Funabiki K, Carr CE. Biophysical basis of the sound analogue membrane potential that underlies coincidence detection in the barn owl. Front Comput Neurosci. 2013; 7: 102. pmid:24265615
  30. 30. Carr CE, Boudreau RE. An axon with a myelinated initial segment in the bird auditory system. Brain Res. 1993; 628: 330–334. pmid:8313166
  31. 31. Kuba H, Ishii TM, Ohmori H. Axonal site of spike initiation enhances auditory coincidence detection. Nature. 2006; 444: 1069–1072. pmid:17136099
  32. 32. Ashida G, Funabiki K, Carr CE. Theoretical foundations of the sound analogue membrane potential that underlies coincidence detection in the barn owl. Front Comput Neurosci. 2013; 7: 151. pmid:24265616
  33. 33. Rathouz M, Trussell LO. Characterization of outward currents in neurons of the avian nucleus magnocellularis. J Neurophysiol. 1998; 80: 2824–2835. pmid:9862887
  34. 34. Rotter S, Diesmann M. Exact digital simulation of time-invariant linear systems with applications to neuronal modeling. Biol Cybern. 1999; 81: 381–402. pmid:10592015
  35. 35. Peña JL, Viete S, Albeck Y, Konishi M. Tolerance to sound intensity of binaural coincidence detection in the nucleus laminaris of the owl. J Neurosci. 1996; 16: 7046–7054. pmid:8824340
  36. 36. Köppl C. Phase locking to high frequencies in the auditory nerve and cochlear nucleus magnocellularis of the barn owl, Tyto alba. J Neurosci. 1997; 17: 3312–3321. pmid:9096164
  37. 37. Carr CE, Boudreau RE. Organization of the nucleus magnocellularis and the nucleus laminaris in the barn owl: encoding and measuring interaural time differences. J Comp Neurol. 1993; 334: 337–355. pmid:8376623
  38. 38. Golding NL, Oertel D. Synaptic integration in dendrites: exceptional need for speed. J Physiol. 2012; 590: 5563–5569. pmid:22930273
  39. 39. Alexandrescu A. The D programming language. Boston, MA: Addison-Wesley; 2010.
  40. 40. Rothman JS, Manis PB. The roles potassium currents play in regulating the electrical activity of ventral cochlear nucleus neurons. J Neurophysiol. 2003; 89: 3097–3113. pmid:12783953
  41. 41. Reyes AD, Rubel EW, Spain WJ. In vitro analysis of optimal stimuli for phase-locking and time-delayed modulation of firing in avian nucleus laminaris neurons. J Neurosci. 1996; 16: 993–1007. pmid:8558268
  42. 42. Reyes AD, Rubel EW, Spain WJ. Membrane properties underlying the firing of neurons in the avian cochlear nucleus. J Neurosci. 1994; 14: 5352–5364. pmid:8083740
  43. 43. Gai Y, Doiron B, Kotak V, Rinzel J. Noise-gated encoding of slow inputs by auditory brainstem neurons with a low-threshold K+ current. J Neurophysiol. 2009; 102: 3447–3460. pmid:19812289
  44. 44. Scott LL, Mathews PJ, Golding NL. Posthearing developmental refinement of temporal processing in principal neurons of the medial superior olive. J Neurosci. 2005; 25: 7887–7895. pmid:16135745
  45. 45. Gerstner W, Kistler W. Spiking Neuron Models. Cambridge UK: Cambridge University Press; 2002.
  46. 46. Bugmann G. Summation and multiplication: two distinct operation domains of leaky integrate-and-fire neurons. Network. 1991; 2: 489–509.
  47. 47. Keener JP, Hoppensteadt FC, Rinzel J. Integrate-and-fire models of nerve membrane response to oscillatory input. SIAM J Appl Math, 1981; 41: 503–517.
  48. 48. Fourcaud-Trocmé F, Hansel D, van Vreeswijk C, Brunel N. How spike generation mechanisms determine the neuronal response to fluctuating inputs. J Neurosci, 2003; 23: 11628–11640. pmid:14684865
  49. 49. Gedeon T, Holzer M. Phase locking in integrate-and-fire models with refractory periods and modulation. J Math Biol. 2004; 49: 577–603. pmid:15565447
  50. 50. Burkitt AN. A review of the integrate-and-fire neuron model: II. Inhomogeneous synaptic input and network properties. Biol Cybern. 2006; 95: 97–112. pmid:16821035
  51. 51. Ashida G, Carr CE. Effect of sampling frequency on the measurement of phase-locked action potentials. Front Cell Neurosci. 2010; 4: 100.
  52. 52. Börgers C, Nectow AR. Exponential time differencing for Hodgkin-Huxley-like ODEs. SIAM J Sci Comput. 2013; 35: B623–B643. pmid:24058276
  53. 53. Johnston J, Forsythe ID, Kopp-Scheinpflug C. Going native: voltage-gated potassium channels controlling neuronal excitability J Physiol. 2010; 588: 3187–3200. pmid:20519310
  54. 54. Golomb D, Ermentrout GB. Continuous and lurching traveling pulses in neuronal networks with delay and spatially decaying connectivity. Proc Natl Acad Sci. 1999; 96: 13480–13485. pmid:10557346
  55. 55. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci. 2000; 8: 183–208. pmid:10809012
  56. 56. Aviel Y, Mehring C, Abeles M, Horn D. On embedding synfire chains in a balanced network. Neural Comput. 2003; 15: 1321–1340. pmid:12816575
  57. 57. Kuokkanen PT, Wagner H, Ashida G, Carr CE, Kempter R. On the origin of the extracellular field potential in nucleus laminaris of the barn owl (Tyto alba). J Neurophysiol. 2010; 104: 2274–2290. pmid:20685926
  58. 58. Mc Laughlin M, Verschooten E, Joris PX. Oscillatory dipoles as a source of phase shifts in field potentials in the mammalian auditory brainstem. J Neurosci. 2010; 30: 13472–13487. pmid:20926673
  59. 59. Gold C, Henze DA, Koch C, Buzsáki G. On the origin of the extracellular action potential waveform: a modeling study. J Neurophysiol. 2006; 95: 3113–3128. pmid:16467426
  60. 60. Boettcher FA. Presbyacusis and the auditory brainstem response. J Speech Lang Hear Res. 2002; 45: 1249–1261. pmid:12546491
  61. 61. Lucas JR, Freeberg TM, Krishnan A, Long GR. A comparative study of avian auditory brainstem responses: correlations with phylogeny and vocal complexity, and seasonal effects. J Comp Physiol A. 2002; 188: 981–992. pmid:12471495