In this section, we describe the two stages of our computational model. In the first stage, pulse trains of current evoke spikes in auditory nerve (AN) fibers through a stochastic model that accounts for stimulus and spike-history effects (Boulet, 2016; Tabibi et al., 2021). The pulses of current represent the electrical signals delivered to cochlear implant electrodes. AN spikes then produce excitatory synaptic conductances that drive voltage activity in the second stage of the model.

The second stage of the model describes an MSO neuron with two separate spatial domains (compartments). The first compartment represents the soma and dendrites, the regions where synaptic inputs are integrated. The second compartment represents the axon initial segment (AIS), the locus of spike-generation. The voltage dynamics of the MSO neuron model are governed by Hodgkin-Huxley type equations that have been modified for use in auditory brainstem neurons (Rothman & Manis, 2003; Khurana et al., 2011; Goldwyn et al., 2019). The two-compartment structure is adapted from previous models of MSO neurons (Khurana et al., 2011; Goldwyn et al., 2019).

We generated auditory nerve fiber spike trains using a phenomenological, stochastic model of responses to cochlear implant stimuli. The model, developed by Boulet (2016) and also described in Tabibi et al. (2021), randomly generates spikes in response to electrical current pulses that represent the CI stimulus. Our implementation of the model differs from theirs by considering spiking activity on a pulse-by-pulse basis, rather than as as a process that can produce spikes at any instant in time. We detail our approach below. All parameter values, function expressions, and other model features are identical to those used by Boulet (2016).

We describe AN activity using a state variable V(t) that represents (phenomenologically) the voltage response of a neuron to stimulation. Voltage is a dynamic variable that responds to input current I(t) according to the linear differential equation

$$\begin \frac = \frac + \frac \end$$

(1)

with parameters \(\tau _m = 0.1350~\text \) and \(C_m = 0.0714~\text \), as in Tabibi et al. (2021). We use the input function I(t) to represent a pulse train of current delivered to a cochlear implant electrode. A single pulse train (to one “ear” of the model) is specified by pulse duration \(\delta \) and interpulse interval \(\Delta \), see Fig. 1. Interpulse interval is the reciprocal of pulse rate, and we typically describe a pulse train in terms of the number of pulses per second (pps). For mathematical simplicity, we use monophasic pulses (positive-going current pulse only), whereas pulses in clinical devices are biphasic (positive and negative-going phases of current). This model is phenomenological in nature, and thus this equation does not represent a biophysically-resolved description of AN voltage dynamics. Among the simplifications this approach entails, the input current I(t) does not create a spatially-distributed electric field that impacts AN fibers through extracellular interactions. While others have modeled such interactions (Rattay et al., 2001; Cohen, 2009; Goldwyn et al., 2010; Kalkman et al., 2015; Malherbe et al., 2015, as examples), by following the approach of Boulet (2016) and Tabibi et al. (2021) we omit these details.

The parameter \(A_i\) is the amplitude of the \(i^\) pulse. In most simulations we take \(A_i\) to be constant across pulses. The exception to this is the final set of simulations in which we vary \(A_i\) values according to a sinusoidal function (to represent a signal with an amplitude-modulated envelope). In calculations described below that detail the model implementation, we denote the start time of the \(i^\) pulse as \(s_i\) and the end time of each pulse as \(e_i\). We refer to \(s_i\) as the times of pulse onsets and \(e_i\) as the times of pulse offsets. In simulations of the model, we use two pulse trains of input current that are identical except for a possible time difference between the two signals (Fig. 1). This time difference represents the interaural time difference (ITD) of a sound arriving at the two ears. Our method of quantifying sensitivity to time differences in the MSO neuron model is to simulate firing rate activity as a function of time differences between the input pulse trains.

We characterize the dynamics of V(t) in two distinct time intervals. The first interval is between the start of a pulse at time \(s_i\) and the end of that pulse at time \(e_i\). The pulse is on during this interval. The second interval is the time from the end of a pulse (\(e_i\)) until the time at the start of the next pulse (\(s_\)). Pulses are off during this time. The voltage variable increases during the first interval, when the pulse is on with \(I(t) = A_i\). Voltage decreases during the second interval (the time between pulses), since \(I(t) = 0\) during these times. We solve the differential equation in Eq. (1) over each of these two intervals using standard methods (Boyce et al., 2021) and obtain the expressions

$$\begin V(t ) \!=\! V(s_i) e^\!+\!\frac}\left( 1 - e^\right) & \text ~s_i \!\le \! t \le e_ \\ V(e_i) e^ & \text ~e_i \!<\! t\!<\! s_ \end\right. }. \end$$

(2)

The first equation applies during the pulses. The second expression applies between pulses.

To randomly generate spike times, we transform the value of the voltage variable to a spike probability value. We do this using the firing efficiency function, a concept introduced into cochlear implant modeling by Bruce et al. (1999a). Spike probability as a function of voltage is

$$\begin p(V) = \frac\left[ 1 + \textrm \left( \frac\sigma } \right) \right] . \end$$

(3)

The abbreviation \(\textrm\) is for the error function, the parameter \(\theta \) is the spike threshold (value of V at which spike probability is 1/2), and the parameter \(\sigma \) determines the variability of firing. Rescaling \(\sigma \) by the threshold parameter defines the quantity \(\rho = \sigma / \theta \). This quantity is known as relative spread (Bruce et al., 1999a) and is frequently used in modeling and neurophysiological studies of AN responses to CI stimulation.

Bruce et al. (1999b) introduced the technique of making threshold \(\theta \) and relative spread \(\rho \) depend on spike-history in order to describe a refractory period. Boulet, Tabibi, and colleagues expanded on this approach to include dynamics in \(\theta \) and \(\rho \) that account for observed phenomena of facilitation, accommodation, and spike rate adaptation in auditory nerve fiber responses to CI stimulation (Boulet, 2016; Tabibi et al., 2021). We detail the dynamics of these quantities below, in subsequent sections of Sect. 2. Our implementation of the model differs from Boulet, Tabibi, and colleagues by restricting spikes to occur only at the end of each pulse. To determine spike probability we must only evaluate V, \(\theta \), and \(\rho \) at the time of pulse offset within each interpulse interval. Our rationale for this choice is that pulse durations are brief (\(50~\mu \)s in our model) and V(t) rises rapidly during a pulse and decays rapidly after a pulse, so p(V) is largest at the end of each pulse.

We calculate the voltage at the end of each pulse using Eq. (2). In particular, evaluating the second expression in Eq. (2) at the end-time of a pulse gives

$$\begin V(e_i) = V(s_i) e^+\frac}\left( 1 - e^\right) . \end$$

(4)

Next, we observe that \(V(s_)\) depends on the voltage at the end-time of the previous pulse according to the exponential decay rule \(V(s_i) = V(})e^\). Thus, using the second expression in Eq. (2), we can replace \(V(s_i)\) with a term that depends on voltage at the ending time of a pulse. From these observations we create an iterative update rule for voltage at pulse offsets:

$$\begin V(e_) = V(e_) e^+\frac}\left( 1 - e^\right) . \end$$

(5)

Beginning from an initial value of \(V_0 =0\), we compute the voltages at the end of each pulse using this formula. We then place these voltage values into Eq. (3) to produce pulse-by-pulse spike probabilities, which we denote \(p_i\) (where \(i=1,2, \ldots \)) labels the pulse number). We reset the voltage variable to zero when spikes occur.

As mentioned above, the values of \(\theta \) and \(\sigma \) in the spike probability equation Eq. (3) are dynamic variables. The values of these parameters change from pulse-to-pulse due to spike history and stimulus history in order to incorporate the effects of refractoriness, facilitation, accommodation, and adaptation. We are calculating spike probabilities at the end of each pulse, thus we require expressions for the values of \(\theta \) and \(\sigma \) at the times \(e_i\) for \(i=1,2,\ldots \). We obtain these expressions by solving relevant linear differential equations, breaking the pulse-to-pulse intervals, in a manner similar to our iterative for \(V(e_)\) in Eq. (5). Following Boulet (2016), we use dynamic expressions for relative spread \(\rho \) and define \(\sigma \) through the expression \(\sigma = \rho \theta \).

We define \(\theta (t)\) and \(\rho (t)\) as the product of their baseline values (\(\theta _0 = 30\) and \(\rho _0=0.04\)) with dynamic multipliers that implement refractory period, facilitation, accommodation, and adaptation (Boulet, 2016; Tabibi et al., 2021):

$$\begin \theta (t) & = \theta _0 x_(t) x_(t) x_(t) x_(t)\\ \rho (t)& = \rho _0 y_(t) y_(t) y_(t) y_(t). \end} \end$$

(6)

We next describe the dynamics of each of these multiplier terms. The implementation we are presenting here follows the same approach of Boulet, Tabibi, and colleagues and uses their parameter values (Boulet, 2016; Tabibi et al., 2021). The adjustment we are making is to calculate spike probabilities only at pulse offset times.

The refractory period refers to the time immediately after a spike during which a subsequent spike cannot be generated (absolute refractory period) or requires additional excitation (relative refractory period) (Miller et al., 2001). Let T represent the time since the last spike. A second spike cannot be produced for T smaller than the absolute refractory period of \(t_=0.332\) ms, meaning spike probability is zero for \(T\le 0.332\) ms and we only use Eq. (3) to calculate spike probabilities when pulse offset time is outside the absolute refractory period. For \(T>0.332\) ms, beyond the absolute refractory period, there is a relative refractory period during which spike probabilities are reduced (larger threshold) and more variable (larger relative spread). In particular, the multiplier terms in Eq. (6) for the refractory effects are

$$\begin x_ & = \left[ 1-\exp (-(T-t_)/\tau _}) \right] ^ \\ y_ & = 1 + \exp \left( -(T-t_)/\tau _} \right) \end} \end$$

(7)

where \(\tau _} = 0.411~\textrm\) and \(\tau _}= 0.2~\textrm\).

Spike probabilities also decrease over a longer time scale than the refractory period during periods of high neural activity. This is the phenomenon of spike rate adaptation (Zhang et al., 2007). Our model for adaptation, following Boulet (2016), is to increase the spike rate multiplier variable (and thereby increase the spike threshold) at the time of each spike by an increment of 0.04. In the periods between spikes, this multiplier variable (\(x_\)) relaxes back to its baseline value of one according to linear dynamics on a slow time scale given by time constant \(\tau _}=50~\textrm\). This is roughly one hundred times slower than the time scales of refractoriness and facilitation. The dynamics of \(x_\) are given by:

$$\begin \frac} = (1 - x_) / \tau _}. \end$$

(8)

Recall from above that we use T to denote the time since the last spike. We also introduce the symbol X to be the value of this adaptation multiplier at that previous spike time. We can then solve the above differential equation to find \(x_\) at any later time (but before the next spike) by initializing the value of \(x_\) at the time of the spike to be \(X+0.04\) and we find

$$\begin x_ = 1 + (X-0.96)e^}}. \end$$

(9)

This equation gives the value of \(x_\) at any time T after a previous spike and before the next spike. The same equations govern the multiplier for relative spread, denoted by \(y_\).

Facilitation refers to the increase in spike probability that occurs on short time scales if the neuron receives multiple subthreshold inputs that temporally-summate. Facilitation in auditory nerve fiber responses to CI stimulation has been characterized using pairs of pulses. A pulse is delivered at a subthreshold current level to decrease the threshold level of the neuron in response to a second pulse delivered a short time later (Cartee et al., 2000).

The facilitation multiplier is driven by changes in voltage, but with a delay of \(\delta \) (the pulse width) so that the multiplier begins to increase at pulse offset. Contrast this with voltage, which begins to increase at pulse onset. We define the facilitation multiplier \(x_\) as \(x_ = 1+z\) where the dynamic variable z is governed by the linear differential equation

$$\begin \frac = -\frac + \frac V(t-\delta ). \end$$

(10)

The forcing term for z is the delayed value of voltage \(V(t-\delta )\) with rescaling of voltage by the parameter \(\theta _0\) (the midpoint in the firing efficiency curve) and a, the parameter that controls the strength of the facilitation effect. The values for these parameters are \(a = -0.15\)/ms with time constant \(\tau =0.5\) ms. The value of a is negative in this case because facilitation increases excitability (lowers the threshold value). We solve this in two intervals using the expressions for V(t) reported in Eq. (2).

Recall from Eq. (2) that the dynamics of V are distinct in two intervals. During a pulse, voltage increases as it is driven by the input current. Between pulses, voltage decays exponentially since it receives no input. To facilitate analysis, we break up the dynamics of z into corresponding phases (but shifted, because of the delay term \(V(t-\delta )\) in Eq. (10)). As introduced previously, we use \(s_i\) to denote the onset time of the \(i^\) pulse and we let \(e_i = s_i +\delta \) be the time at pulse offset. During a short time immediately after pulse offset, that is for the time interval \(e_i< t< e_i + \delta \) the value of z is driven by the phase of V that is increasing in response to the pulse of current. Thus, in a short time interval after the end of a pulse, z is governed by the equation

$$\begin \frac = -\frac + \frac\left( V(s_i)e^ +\frac} \left( 1 - e^\right) \right) , \quad \text e_i< t< e_i + \delta . \end$$

(11)

For the remainder of the time until the end of the next pulse, that is for the time interval \(e_i + \delta< t < e_i + \Delta \), the dynamics of z are driven by the exponentially-decaying phase of V. Thus, z is governed by the equation

$$\begin \frac = -\frac + \frac \left( V(e_i) e^\right) , \quad \text e_i + \delta< t < e_i + \Delta . \end$$

(12)

We always choose interpulse intervals to be larger than twice the pulse duration, and so we can first solve Eq. (11) and then use the appropriate value of z as the initial value for the next phase of z (governed by Eq. (12)). In this manner, we can iteratively solve these two differential equations to obtain values of z at the end-times of pulses (\(e_i\)). Recall, these are the only times at which we allow AN spike generation to occur in our implementation of the model. Our calculation proceeds as follows. First, we use \(z(e_i)\), the value of z at pulse offset, as the initial value for Eq. (11) (with \(z=0\) as the first initial value at the start of the stimulus). We solve this equation at time \(e_i + \delta \) and obtain

$$\begin z(e_i+\delta )&= z(e_i)e^ + \frac\left[ \left( \frac \right) \left( 1-e^\right) \right. \nonumber \\&\left. \quad + \left( V(s_i) - \frac \right) \mathcal e^ \left( 1 -e^}\right) \right] \end$$

(13)

where \(\mathcal = \frac\).

Next, we use this value of z as the initial value for Eq. (12). We solve this equation and evaluate z at the next pulse offset to find

$$\begin z(e_)&= z(e_i + \delta )e^\nonumber \\&\quad + \frac V(e_i) \mathcal e^ \left( 1 - e^}\right) . \end$$

(14)

We add one to this value to obtain the value of the threshold multiplier \(x_\) at the pulse offset, which is the information we require to compute spike probability. The facilitation process resets (\(z=0)\) at every spike time and at the end of every pulse (Boulet, 2016). We use the same calculation scheme for the relative spread multiplier \(y_\). Parameter values for the relative spread multiplier are \(a=0.75\)/ms and \(\tau = 0.3\) ms. As stated previously, we take all parameter values for the AN model from Tabibi et al. (2021).

On longer time scales than facilitation, multiple subthreshold inputs can have a cumulative effect to gradually increase spike threshold. This is the phenomenon of accommodation (Negm & Bruce, 2014). The multipliers for accommodation are split into two components, quick and slow accommodation. The dynamics of these variables are treated using the same scheme as facilitation. An auxiliary variable (denoted as z in the treatment of facilitation, see Eq. (10)) has dynamics driven by delayed voltage, and the accommodation multipliers are one plus the auxiliary variable. The accommodation effect makes the AN less excitable (increases spike threshold), thus the scaling parameter a is positive for both the quick and slow components of accommodation. For the quick accommodation threshold multiplier \(a=0.5\)/ms and \(\tau = 1.5\) ms, for the slow accommodation threshold multiplier \(a = 0.01\)/ms and \(\tau = 50\) ms. Notice that the dynamics of these components are distinguished by their constants, which differ by an order of magnitude. The quick component of the accommodation relative spread multiplier has \(a = 0.75\)/ms and \(\tau =0.3\) ms. There is no slow component for the relative spread multiplier (\(a =0\)/ms).

As constructed up to this point, our implementation of the AN model restricts spike times to occur only at pulse offsets. While AN spikes are well-timed to the pulses of CI stimulation, they do exhibit a small amount of temporal jitter (Miller et al., 2001). The neural computation we are simulating is coincidence detection between two streams of synaptic inputs, so incorporating this small amount of temporal jitter is an important element of our model. Our approach is to temporally displace each spike time by a random amount. In particular, we add to each spike time a normally-distributed random number with mean zero and standard deviation of 0.1 ms. This standard deviation value is based on physiological measurements of AN spike time jitter in response to CI stimulation (Miller et al., 2001).

We use auditory nerve fiber activity as inputs to a biophysically-based computational model of a principal cell in the medial superior olive (MSO). We omit additional, early stages of auditory pathway that are intermediate between the auditory nerve and MSO. These include cochlear nucleus neurons that pass excitatory inputs forward to MSO and neurons in the medial nucleus of the trapezoid body that inhibit MSO. Although these nuclei impact MSO tuning (cells in the anteroventral cochlear nucleus can enhance the temporal precision of inputs to MSO neurons Joris et al., 1994 and inhibitory inputs can adjust the positions of the peaks of ITD tuning curves Brand et al., 2002; Jercog et al., 2010; Myoga et al., 2014), our approach captures the essential architecture of the MSO circuit. Specifically, our model describes MSO neurons as neural coincidence detectors of well-timed excitatory inputs originating from both ears.

The synaptic inputs to the MSO neuron model are in the form of excitatory post-synaptic conductances (EPSGs). For each auditory nerve fiber spike, we create a unitary EPSG of the form

$$\begin g_(t) = 98.5 \left( e^ - e^ \right) . \end$$

(15)

The rise and decay time scales of this double exponential function are in units of milliseconds and are taken from reports of in vivo MSO recordings (Franken et al., 2015). The amplitude of this unitary EPSG is in units of nanosiemens and is selected so that a excitatory input depolarizes the soma voltage of the MSO neuron (described below) by roughly 6 mV. We selected this value based on MSO recordings reported in Roberts et al. (2013).

MSO neurons have been described as having few inputs (Couchman et al., 2010), so in our model we use ten independent synaptic inputs for each “side” of the MSO neuron. We vary the time delay of the two input streams (contralateral and ipsilateral, so to speak) to represent sounds that arrive at the two ears with different ITDs. In total, the synaptic input to the MSO neuron is in the form a current

$$\begin I_(t) = \left( g_(t) + g_(t) \right) (V_1(t) - E_), \end$$

(16)

where \(g_\) and \(g_\) are each composed from the sum of ten synaptic input streams, \(V_1\) is the soma voltage of the MSO neurons (defined next), and \(E_ = 0\) mV is the reversal potential for synaptic excitation.

Principal cells in the medial superior olive are the first nucleus in which bilateral inputs are combined to make fine-timing comparisons to encode sound source location (Grothe et al., 2010, for review). We model these neurons using a two-compartment model that describes separately the dynamics of the soma/dendrite region (first compartment) and the axon initial segment (second compartment). We adapt the model of Khurana et al. (2011), a Hodgkin-Huxley type model with low and high-threshold potassium currents, spike-generating sodium current (in the second compartment only), and hyperpolarization-activated cation current; all based on the Rothman-Manis model (Rothman & Manis, 2003). In earlier work, we discussed how two-compartment models can be parameterized by the relative strengths of forward and backward coupling (Goldwyn et al., 2019). Based on that work, we adjusted parameter values to increase the strength of forward coupling (less voltage attenuation from soma to axon), since stronger forward connection from soma to axon are beneficial for ITD sensitivity (Goldwyn et al., 2019).

We list the parameter values for our model in Table 1 in the column labelled “control.” These parameters result in a forward coupling strength of 0.8 (meaning the passive, steady state voltage in the AIS is 80% of the steady state voltage in the soma for constant current input to the soma) and a backward coupling strength of 0.5 (meaning the passive, steady state voltage in the soma is 50% of the steady state voltage in the AIS for constant current input to the AIS). See Goldwyn et al. (2019) for details and definition of coupling coefficients in two-compartment models. We also adjusted leak reversal potentials in both compartments to achieve a resting potential of -58 mV. For these parameter values the input resistance at the soma was 10 M\(\Omega \) for the control model and 15 M\(\Omega \) for the deprived model, and the membrane time constant in the soma was 300\(\mu \)s for the control model and 500\(\mu \)s for the deprived model.

The equations that govern the voltage dynamics in the soma/dendrite region (\(V_1\)) and the AIS region (\(V_2\)) are

$$\begin c_1 \frac & = -I_ - I_ - I_ - I_ - g_(V_1-V_2) \\ c_2 \frac & = -I_ - I_ - I_ - I_ - I_ - g_(V_2-V_1). \end} \end$$

(17)

As in Khurana et al. (2011), the expressions for the voltage-dependent currents are \(I_ = g_(V_i - E_)\) for leak currents (we use \(i=1\) or \(i=2\) to label the compartment number), \(I_ = g_ (0.65 r_ + 0.35 r_)(V_i - E_h)\) for hyperpolarization-activated cation current, \(I_ = g_w_i^4z_i(V_i - E_K)\) for low-threshold potassium, \(I_ = g_(0.85n^2 + 0.15 p)(V_2 - E_K)\) for high-threshold potassium current, and \(I_ = g_m^3 h (V_2 - E_)\) for sodium current. Notice that sodium and high-threshold potassium currents are only present in the spike-initiation zone (compartment two). We described the synaptic current term \(I_\) previously in Eqs. (15) and (16).

All gating variables are governed by differential equations