Axon initial segment plasticity caused by auditory deprivation degrades time difference sensitivity in a model of neural responses to cochlear implants

TL;DR

Auditory deprivation changes the axon initial segment in brainstem neurons, reducing their ability to detect time differences in sound, which is important for cochlear implant users.

Contribution

The study shows how AIS plasticity in MSO neurons reduces time difference sensitivity in a cochlear implant model.

Findings

AIS plasticity converts MSO neurons from phasic to tonic firing type.

AIS plasticity degrades time difference sensitivity in auditory deprived MSO neurons.

Changes in AIS biophysics may hinder sound source localization in cochlear implant users.

Abstract

Synaptic and neural properties can change during periods of auditory deprivation. These changes may disrupt the computations that neurons perform. In the brainstem of chickens, auditory deprivation can lead to changes in the size and biophysics of the axon initial segment (AIS) of neurons in the sound source localization circuit. This is the phenomenon of axon initial segment (AIS) plasticity. Individuals who use cochlear implants (CIs) experience periods of hearing loss, and so we ask whether AIS plasticity in neurons of the medial superior olive (MSO), a key stage of sound location processing, would impact time difference sensitivity in the scenario of hearing with cochlear implants. The biophysical changes that we implement in our model of AIS plasticity include enlargement of the AIS and replacement of low-threshold potassium conductance with the more slowly-activated M-type potassium conductance. AIS plasticity has been observed to have a homeostatic effect with respect to excitability. In our model, AIS plasticity has the additional effect of converting MSO neurons from phasic firing type to tonic firing type. Phasic firing is known to have greater temporal sensitivity to coincident inputs. Consistent with this, we find AIS plasticity degrades time difference sensitivity in the auditory deprived MSO neuron model across a range of stimulus parameters. Our study illustrates a possible mechanism of cellular plasticity in a non-peripheral stage of neural processing that could impose barriers to sound source localization by bilateral cochlear implant users.

Figures

Schematic of model for MSO responses to CI input. Cochlear implant stimuli are represented as pulses of current that drive spiking activity in a model for auditory neurons. Auditory neuron spikes produce excitatory conductance. The combined synaptic conductance from “right ear” and “left ear” inputs enter in the first compartment of the MSO neuron model. Parameters that define the current pulse input are time difference between “right ear” and “left ear” inputs, pulse duration ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb}

Spiking activity in the auditory nerve (AN) model. Responses to constant amplitude pulse trains. Raster plots show spike times in 50 repeated simulations of the AN model in response to constant amplitude pulse trains (A1) pulse rate is 250 pps and pulse amplitude is 55 nA, (B1) pulse rate is 500 pps and pulse amplitude is 60 nA, and (C1) pulse rate is 1000 pps and pulse amplitude is 75 nA. As specified in the model, spike times are aligned to the pulse offset with a small amount of spike time jitter, as evident in the post-stimulus histograms (PSTHs) shown in the second row for each stimulus s

Voltage traces in each compartment of the MSO neuron model in response to cochlear implant stimulation. Inputs are: (A1-4): 250 pps at 53 nA, (B1-4) 500 pps at 60 nA, and (C1-4) 1000 pps at 80 nA. These stimulus settings were selected so the firing rate is roughly 150 sp/s in the control model in all cases. Responses of the control model are shown in rows 1 and 2. Responses of the deprived model are shown in rows 3 and 4. The inputs have 0 ms time difference in rows 1 and 3, and have non-zero time differences in rows 2 and 4 (1 ms time difference for 250 pps and 500 pps inputs, and 0.5 ms time

Spike rate for different current strengths, using 300 ms inputs with 0 ms time difference and constant pulse amplitude (current values given on horizontal axis). Pulse rates in each panel are (A) 100 pps, (B) 250 pps, (C) 500 pps, (D) 750 pps, and (E) 1000 pps. Firing rates are mean of 100 repeated simulations calculated at intervals of 1 nA. Error bars are standard deviations around those means, plotted at intervals of 5 nA. Firing rates increase with current and saturate at the pulse rate for the lower pulse rates (panels A and B). The auditory deprived model is more excitable than the contr

Summary of firing rate sensitivity to time differences. Spike rates of MSO neuron model to constant amplitude pulse trains for various pulse rates: (A) 100 pps pulse rate, (B) 250 pps pulse rate, (C) 500 pps pulse rate, (D) 750 pps pulse rate, and (E) 1000 pps pulse rate. Maximal firing rates occur for inputs with 0 ms time difference (upper solid curves). Minimal firing rates occur for larger time differences (lower solid curves, using 1 ms time difference for lower pulse rates and half the interpulse intervals for 750 pps and 1000 pps). Values are averages over 100 repeated trials. The shade

Envelope time-difference sensitivity. Spike rates of MSO neuron model to 3000 pps pulse trains with sinusoidally amplitude-modulated pulse amplitudes. The sinusoidal oscillation of pulse amplitudes varies between a minimum level of 70 nA and a maximum level shown on the *x*-axis. The modulation frequencies are: (A) 0 Hz, unmodulated pulse train, (B) 100 Hz modulation frequency, (C) 200 Hz modulation frequency, and (D) 300 Hz modulation frequency. Maximal firing rates occur for inputs with 0 ms time difference (upper solid curves). Minimal firing rates occur for non-zero time difference (lower

Firing type is phasic for control model and tonic for auditory deprived model. AIS voltage response ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document}$ , voltage in the second compartment) in response to steps of current injections (50 ms duration), with current level ranging from 400 pA to 2200 pA in steps of 300 pA. In A (left column): control model responds with a single spike at stimulus ons

Full text

Introduction

Cochlear implants (CIs) are neural prostheses that can produce a sense of auditory perception by stimulating the auditory nerve (AN) with electrical signals. These devices can improve users’ speech understanding (Krueger et al., 2008; Gaylor et al., 2013, e.g.) and improve perception of other important aspects of sound (Wilson & Dorman, 2008, for review). For some individuals, CIs are implanted in both ears (bilateral CIs). The benefits of implanting a second CI include improved speech understanding in noise and improved sound source localization relative to users of a single CI (unilateral implant) (van Hoesel, 2004; Brown & Balkany, 2007). Compared to normal hearing listeners, however, users of bilateral CIs exhibit decreased sound source localization accuracy (Grantham et al., 2007).

Two primary cues for sound source location that listeners can use are the difference in sound level at the two ears (interaural level difference, ILD) and the difference in the arrival times of sounds at the two ears (interaural time difference, ITD). Bilateral CI users primarily rely on ILDs to determine sound source location (van Hoesel & Tyler, 2003; Seeber & Fastl, 2008). Factors that may limit the utility of ITD cues for bilateral CI users include mismatched electrode positions in the two ears (Kan et al., 2013), differences in AN temporal responses in response to CI and acoustic stimuli, and AN fiber degeneration that can vary depending on the etiology and duration of deafness. A number of modeling studies have explored how neuron death, demylenation, or other degenerative effects may impact AN activation by CIs (Briaire & Frijns, 2006; Resnick et al., 2018; Goldwyn et al., 2010). Recent modeling work also shows that demyelination of AN axons and synaptopathy of the inner hair cell to AN synapse can disrupt the sound localization circuit in the context of hidden hearing loss (Budak et al., 2022). Here we use modeling to study whether possible pathological changes to the axonal structure of neurons in the medial superior olive (MSO) may degrade the functioning of these neurons in the context of neural responses to bilateral CI inputs.

Periods of hearing loss can induce degenerative and plastic changes in auditory brainstem neurons (Walmsley et al., 2006). In this work, we focus on changes that can occur at the site of spike initiation in neurons in the auditory pathway due to periods of auditory deprivation. Specifically, we draw on a number of studies that demonstrate changes to spike initiation regions of neurons during periods of auditory deprivation (Kuba et al., 2010; Kuba, 2012; Kuba et al., 2014, 2015; Kim et al., 2019). This phenomenon is known as axon initial segment (AIS) plasticity (Yamada & Kuba, 2016). AIS plasticity is conceptualized as a homeostatic response to changed synaptic inputs (Yamada & Kuba, 2016). During periods of decreased excitatory input rates to a neuron (by ablation of the cochlea, in experiments), the neuron responds by increasing its intrinsic excitability, thereby compensating for the diminished excitatory input. Changes to the neuron during AIS plasticity include enlarged AIS size, and alterations to both the density and composition of a variety of potassium channels.

We investigate how AIS plasticity may impact the dynamics of MSO neurons. MSO neurons are a critical component of the mammalian sound localization circuit that encodes ITDs. MSO neurons function as temporally-precise coincidence detectors that encode small time differences in sound arrival times at the two ears (Grothe et al., 2010; McAlpine & Grothe, 2003, for review). Although AIS plasticity has not been observed in human MSO neurons, we consider AIS plasticity in these neurons plausible in the context of human CI users. We believe this because CI users experience prolonged periods of hearing loss, during which time the inputs to their auditory systems are diminished. Additionally, binaural neurons in the ITD-detecting circuit of chickens (nucleus magnocellularis neurons, NM) exhibit AIS plasticity during auditory deprivation (Kuba et al., 2014) and binaural coincident detector neurons in chickens in the nucleus laminaris (NL) region (the homologue of MSO) have AIS regions with lengths tuned to the frequency of their synaptic inputs (Kuba et al., 2006; Kuba, 2012).

We hypothesized that if there is AIS plasticity in MSO neurons in a manner similar to what has been observed in NM neurons, then MSO neurons would become over-excitable and less able to perform their role as temporally-precise coincidence detectors of well-timed binaural excitatory inputs. We tested this hypothesis by simulating MSO neural activity driven by auditory nerve responses to bilateral CI input over a range ITD values. We compared responses of the model of MSO neurons with standard parameter values (Rothman & Manis, 2003; Khurana et al., 2011) to responses of a model with parameters adjusted to reflect AIS plasticity caused by auditory deprivation.

In our previous studies, we showed that temporally-precise coincidence detection by auditory brainstem neurons depends on optimally calibrated soma and axon regions (Goldwyn et al., 2019; Drucker & Goldwyn, 2023). Consistent with those studies, we show here that alteration to the soma-axon structure by AIS plasticity degrades ITD sensitivity. In particular we find that, although AIS plasticity is thought of as a homeostatic change with respect to firing rate, the nature of the neural dynamics (the “firing type”) changes in a consequential way. The parameter changes we make for AIS plasticity convert the model from phasic firing type in the control parameter setting to tonic firing type in the auditory deprived model. Principal cells of the MSO are phasic (Svirskis et al., 2002), as are the homologous ITD detectors in the nucleus laminaris of birds (Reyes et al., 1996). Phasic dynamics enhance ITD sensitivity in MSO and NL neurons (Meng et al., 2012; Huguet et al., 2017; Drucker & Goldwyn, 2023). Accordingly, we find a loss of time-difference sensitivity when AIS plasticity converts the model to tonic type from phasic type.

In Sect. 2 we describe our implementation of the model for AN spiking in response to CI stimulation. Our model is an implementation of an established model for AN responses to CI stimulation (Boulet, 2016; Tabibi et al., 2021). We also describe the MSO model and the changes we make to the AIS based on the observations of AIS plasticity by Kuba and colleagues (Kuba et al., 2010; Kuba, 2012; Kuba et al., 2014, 2015). In Sect. 3 we present how our model of AIS plasticity degrades time-difference processing in simulations of MSO responses to CI stimulation. We present our main results in terms of time-difference tuning curves (firing rates as a function of the difference in the timing of two streams of pulses that represent the CI stimulation at two different ears). We close in the Sect. 4 section with a review of our findings and their limitations, as well as comments on the implications of our results for ITD-based sound localization by bilateral CI users.Fig. 1. Schematic of model for MSO responses to CI input. Cochlear implant stimuli are represented as pulses of current that drive spiking activity in a model for auditory neurons. Auditory neuron spikes produce excitatory conductance. The combined synaptic conductance from “right ear” and “left ear” inputs enter in the first compartment of the MSO neuron model. Parameters that define the current pulse input are time difference between “right ear” and “left ear” inputs, pulse duration ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} ), and interpulse interval ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta $$\end{document} , the reciprocal of pulse rate)

Methods

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

Auditory nerve model

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dV(t)}{dt} = \frac{V(t)}{\tau _m} + \frac{I(t)}{C_m} \end{aligned}$$\end{document}

with parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _m = 0.1350~\text {ms}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_m = 0.0714~\text {pF}$$\end{document} , 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} and interpulse interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta $$\end{document} , 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_i$$\end{document} is the amplitude of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i^{th}$$\end{document} pulse. In most simulations we take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_i$$\end{document} to be constant across pulses. The exception to this is the final set of simulations in which we vary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_i$$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i^{th}$$\end{document} pulse as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} and the end time of each pulse as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i$$\end{document} . We refer to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} as the times of pulse onsets and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i$$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} and the end of that pulse at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i$$\end{document} . The pulse is on during this interval. The second interval is the time from the end of a pulse ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i$$\end{document} ) until the time at the start of the next pulse ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{i+1}$$\end{document} ). Pulses are off during this time. The voltage variable increases during the first interval, when the pulse is on with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I(t) = A_i$$\end{document} . Voltage decreases during the second interval (the time between pulses), since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I(t) = 0$$\end{document} 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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V(t ) \!=\! {\left\{ \begin{array}{ll} V(s_i) e^{-(t-s_i)/\tau _m}\!+\!\frac{A_i \tau _m }{C_{m}}\left( 1 - e^{-(t-s_i)/\tau _m}\right) & \text { for}~s_i \!\le \! t \le e_{i} \\ V(e_i) e^{-(t-e_i)/\tau _m} & \text { for}~e_i \!<\! t\!<\! s_{i+1} \end{array}\right. }. \end{aligned}$$\end{document}

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p(V) = \frac{1}{2}\left[ 1 + \textrm{erf} \left( \frac{V - \theta }{\sqrt{2}\sigma } \right) \right] . \end{aligned}$$\end{document}

The abbreviation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{erf}$$\end{document} is for the error function, the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} is the spike threshold (value of V at which spike probability is 1/2), and the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} determines the variability of firing. Rescaling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} by the threshold parameter defines the quantity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho = \sigma / \theta $$\end{document} . 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} and relative spread \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} depend on spike-history in order to describe a refractory period. Boulet, Tabibi, and colleagues expanded on this approach to include dynamics in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} 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, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} at the time of pulse offset within each interpulse interval. Our rationale for this choice is that pulse durations are brief ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$50~\mu $$\end{document} 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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V(e_i) = V(s_i) e^{-\delta /\tau _m}+\frac{A_i \tau _m }{C_{m}}\left( 1 - e^{-\delta /\tau _m}\right) . \end{aligned}$$\end{document}

Next, we observe that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(s_{i})$$\end{document} depends on the voltage at the end-time of the previous pulse according to the exponential decay rule \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(s_i) = V({e_{i-1}})e^{-(\Delta -\delta )/\tau _m}$$\end{document} . Thus, using the second expression in Eq. (2), we can replace \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(s_i)$$\end{document} 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:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V(e_{i}) = V(e_{i-1}) e^{-\Delta /\tau _m}+\frac{A_i \tau _m }{C_{m}}\left( 1 - e^{-\delta /\tau _m}\right) . \end{aligned}$$\end{document}

Beginning from an initial value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_0 =0$$\end{document} , 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_i$$\end{document} (where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2, \ldots $$\end{document} ) labels the pulse number). We reset the voltage variable to zero when spikes occur.

As mentioned above, the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} at the times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,\ldots $$\end{document} . We obtain these expressions by solving relevant linear differential equations, breaking the pulse-to-pulse intervals, in a manner similar to our iterative for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(e_{i})$$\end{document} in Eq. (5). Following Boulet (2016), we use dynamic expressions for relative spread \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} and define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} through the expression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma = \rho \theta $$\end{document} .

We define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta (t)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (t)$$\end{document} as the product of their baseline values ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _0 = 30$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0=0.04$$\end{document} ) with dynamic multipliers that implement refractory period, facilitation, accommodation, and adaptation (Boulet, 2016; Tabibi et al., 2021):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} \theta (t) & = \theta _0 x_{ref}(t) x_{fac}(t) x_{acc}(t) x_{ad}(t)\\ \rho (t)& = \rho _0 y_{ref}(t) y_{fac}(t) y_{acc}(t) y_{ad}(t). \end{matrix}} \end{aligned}$$\end{document}

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.

Refractory period

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_{abs}=0.332$$\end{document}  ms, meaning spike probability is zero for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T\le 0.332$$\end{document}  ms and we only use Eq. (3) to calculate spike probabilities when pulse offset time is outside the absolute refractory period. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0.332$$\end{document}  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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} x_{ref} & = \left[ 1-\exp (-(T-t_{abs})/\tau _{x_{ref}}) \right] ^{-1} \\ y_{ref} & = 1 + \exp \left( -(T-t_{abs})/\tau _{y_{ref}} \right) \end{matrix}} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{x_{ref}} = 0.411~\textrm{ms}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{y_{ref}}= 0.2~\textrm{ms}$$\end{document} .

Spike rate adaptation

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 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ad}$$\end{document} ) relaxes back to its baseline value of one according to linear dynamics on a slow time scale given by time constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{x_{ad}}=50~\textrm{ms}$$\end{document} . This is roughly one hundred times slower than the time scales of refractoriness and facilitation. The dynamics of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ad}$$\end{document} are given by:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dx_{ad}}{dt} = (1 - x_{ad}) / \tau _{x_{ad}}. \end{aligned}$$\end{document}

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ad}$$\end{document} at any later time (but before the next spike) by initializing the value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ad}$$\end{document} at the time of the spike to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X+0.04$$\end{document} and we find

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} x_{ad} = 1 + (X-0.96)e^{-T/\tau _{x_{ad}}}. \end{aligned}$$\end{document}

This equation gives the value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ad}$$\end{document} at any time T after a previous spike and before the next spike. The same equations govern the multiplier for relative spread, denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{ad}$$\end{document} .

Facilitation

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} (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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{fac}$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{fac} = 1+z$$\end{document} where the dynamic variable z is governed by the linear differential equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dz}{dt} = -\frac{z}{\tau } + \frac{a}{\theta _0} V(t-\delta ). \end{aligned}$$\end{document}

The forcing term for z is the delayed value of voltage \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(t-\delta )$$\end{document} with rescaling of voltage by the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _0$$\end{document} (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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = -0.15$$\end{document} /ms with time constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau =0.5$$\end{document}  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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(t-\delta )$$\end{document} in Eq. (10)). As introduced previously, we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_i$$\end{document} to denote the onset time of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i^{th}$$\end{document} pulse and we let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i = s_i +\delta $$\end{document} be the time at pulse offset. During a short time immediately after pulse offset, that is for the time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i< t< e_i + \delta $$\end{document} 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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dz}{dt} = -\frac{z}{\tau } + \frac{a}{\theta _0}\left( V(s_i)e^{-(t-s_i)/\tau _m} +\frac{A \tau _m }{C_{m}} \left( 1 - e^{-(t-s_i)/\tau _m}\right) \right) , \quad \text { for } e_i< t< e_i + \delta . \end{aligned}$$\end{document}

For the remainder of the time until the end of the next pulse, that is for the time interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i + \delta< t < e_i + \Delta $$\end{document} , the dynamics of z are driven by the exponentially-decaying phase of V. Thus, z is governed by the equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dz}{dt} = -\frac{z}{\tau } + \frac{a}{\theta _0} \left( V(e_i) e^{-(t-e_i)/\tau _m}\right) , \quad \text { for } e_i + \delta< t < e_i + \Delta . \end{aligned}$$\end{document}

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 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i$$\end{document} ). 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z(e_i)$$\end{document} , the value of z at pulse offset, as the initial value for Eq. (11) (with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0$$\end{document} as the first initial value at the start of the stimulus). We solve this equation at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_i + \delta $$\end{document} and obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} z(e_i+\delta )&= z(e_i)e^{-\delta /\tau } + \frac{a}{\theta _0}\left[ \left( \frac{A\tau _m \tau }{C_m} \right) \left( 1-e^{-\delta /\tau }\right) \right. \nonumber \\&\left. \quad + \left( V(s_i) - \frac{A\tau _m}{C_m} \right) \mathcal {T}e^{-\delta /\tau _m} \left( 1 -e^{-\delta /\mathcal {T}}\right) \right] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal { T} = \frac{\tau _m \tau }{\tau _m - \tau }$$\end{document} .

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} z(e_{i+1})&= z(e_i + \delta )e^{-(\Delta -\delta )/\tau }\nonumber \\&\quad + \frac{a}{\theta _0} V(e_i) \mathcal {T}e^{-(\Delta -\delta )/\tau _m} \left( 1 - e^{ -(\Delta -\delta )/\mathcal {T}}\right) . \end{aligned}$$\end{document}

We add one to this value to obtain the value of the threshold multiplier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{fac}$$\end{document} at the pulse offset, which is the information we require to compute spike probability. The facilitation process resets ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0)$$\end{document} at every spike time and at the end of every pulse (Boulet, 2016). We use the same calculation scheme for the relative spread multiplier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{fac}$$\end{document} . Parameter values for the relative spread multiplier are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=0.75$$\end{document} /ms and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 0.3$$\end{document}  ms. As stated previously, we take all parameter values for the AN model from Tabibi et al. (2021).

Accommodation

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=0.5$$\end{document} /ms and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 1.5$$\end{document}  ms, for the slow accommodation threshold multiplier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 0.01$$\end{document} /ms and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 50$$\end{document}  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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 0.75$$\end{document} /ms and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau =0.3$$\end{document}  ms. There is no slow component for the relative spread multiplier ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a =0$$\end{document} /ms).

Spike time jitter

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

Synaptic input model

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g_{syn}(t) = 98.5 \left( e^{-t/0.18} - e^{-t/0.1} \right) . \end{aligned}$$\end{document}

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} I_{syn}(t) = \left( g_{contra}(t) + g_{ipsi}(t) \right) (V_1(t) - E_{syn}), \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{contra}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{ipsi}$$\end{document} are each composed from the sum of ten synaptic input streams, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} is the soma voltage of the MSO neurons (defined next), and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{syn} = 0$$\end{document}  mV is the reversal potential for synaptic excitation.

MSO neuron model

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} for the control model and 15 M \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} for the deprived model, and the membrane time constant in the soma was 300 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} s for the control model and 500 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} s for the deprived model.Table 1. Parameters for MSO neuron modelparametercontroldeprivedcompartment 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_1 $$\end{document} 30.1 pF30.1 pF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{lk,1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-57.9$$\end{document}  mV \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-57.9$$\end{document}  mV \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{lk,1}$$\end{document} 16.5 nS16.5 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{h,1} $$\end{document} 76.9 nS38.4 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{KLT,1} $$\end{document} 208.6 nS104.3 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{KHT,1} $$\end{document} 0 nS0 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{Na,1} $$\end{document} 0 nS0 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{M,1} $$\end{document} 0 nS0 nScompartment 2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_2$$\end{document} 14.5 pF21.7 pF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{lk,2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-59.0$$\end{document}  mV \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-59.0$$\end{document}  mV \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{lk,2}$$\end{document} 3.95 nS5.93 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{h,2} $$\end{document} 19.2 nS19.2 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{KLT,2} $$\end{document} 52.2 nS0 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{KHT,2} $$\end{document} 150 nS225 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{Na,2} $$\end{document} 2000 nS3000 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{M,2} $$\end{document} 0 nS21.9 nS

The equations that govern the voltage dynamics in the soma/dendrite region ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} ) and the AIS region ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} ) are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} c_1 \frac{dV_1}{dt} & = -I_{lk,1} - I_{h,1} - I_{KLT,1} - I_{syn} - g_{ax}(V_1-V_2) \\ c_2 \frac{dV_2}{dt} & = -I_{lk,2} - I_{h,2} - I_{KLT,2} - I_{KHT} - I_{Na} - g_{ax}(V_2-V_1). \end{matrix}} \end{aligned}$$\end{document}

As in Khurana et al. (2011), the expressions for the voltage-dependent currents are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{lk,i} = g_{lk,i}(V_i - E_{lk,i})$$\end{document} for leak currents (we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=2$$\end{document} to label the compartment number), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{h,i} = g_{h,i} (0.65 r_{r,i} + 0.35 r_{s,i})(V_i - E_h)$$\end{document} for hyperpolarization-activated cation current, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{KLT,i} = g_{KLT,i}w_i^4z_i(V_i - E_K)$$\end{document} for low-threshold potassium, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{KHT} = g_{KHT}(0.85n^2 + 0.15 p)(V_2 - E_K)$$\end{document} for high-threshold potassium current, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{Na} = g_{Na}m^3 h (V_2 - E_{Na})$$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{syn}$$\end{document} previously in Eqs. (15) and (16).

All gating variables are governed by differential equations of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tau _x \frac{dx}{dt} = x_\infty - x. \end{aligned}$$\end{document}

Defining equations for time constants and equilibrium values are given in Khurana et al. (2011), we also include these in the Appendix.

Changes to MSO neuron model to describe AIS plasticity

Changes to MSO neurons due to auditory deprivation have not been reported, to our knowledge. We therefore draw on observations from the avian sound localization system to make changes to the MSO model that could represent possible effects of hearing loss on MSO structure and physiology. We refer to the model with parameters altered to reflect the effects of auditory deprivation as the deprived model. We contrast the deprived model with the control model. We construct the deprived model by changing the MSO model, as described below. All aspects of the AN model described above are identical for both control and deprived models.

The axon initial segment in neurons in the nucleus magnocellularis (NM) of chickens enlarges during periods of auditory deprivation. Specifically, the length of the axon region with dense sodium channels expands, although the distance from the soma to the AIS does not change (Kuba et al., 2010, 2014). Additional studies of AIS plasticity for NM neurons under conditions of auditory deprivation showed that, while there were no changes in the density of sodium channels, there were changes to the distribution and makeup of potassium channels in the soma and AIS (Kuba et al., 2010, 2015). In particular, low threshold potassium current (Kv1.1 subunit type) diminishes in the soma and AIS. There is a complementary increase in M-type potassium current (Kv7.2 subunit type) in the AIS during auditory deprivation, indicating a possible reconfiguration of potassium channels from low threshold type to M-type. Kuba and colleagues argue that these changes comprise a homeostatic adjustment that maintains neural activity when excitatory synaptic inputs are lost during auditory deprivation (Kuba et al., 2010; Kuba, 2012).

To model possible features of AIS plasticity induced by auditory deprivation, we changed several parameter values. To describe the elongation of the AIS during auditory deprivation, we increased passive parameters (capacitance, leak conductance) in the second compartment by 50%, consistent with an increase in the surface area of that compartment. Sodium channel density has not been observed to change due to auditory deprivation, so we also increased sodium conductance by 50% to retain the same sodium channel density for the larger-sized AIS compartment. To increase excitability – the homeostatic effect observed in experiments – we decreased low-threshold potassium conductance in the soma by 50%. We obtained similar time difference tuning curves, but with decreased excitability (non-homeostatic) using larger values of potassium conductance (simulations not shown). We matched the 50% decrease in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{KLT}$$\end{document} with a 50% decrease in the hyperpolarization-activated cation conductance ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_h$$\end{document} ) to preserve the balance of these conductances in the soma, as has been observed in vitro (Khurana et al., 2011). Lastly, and importantly, based on the experimental findings and modeling work in Kuba et al. (2015), we reduced potassium conductance in the AIS by 50% and converted the low-threshold potassium conductance to M-type potassium conductance. This results in a substantial slow-down of the potassium current. The activation time constant for low threshold potassium when voltage is at the resting potential of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-58~mV$$\end{document} is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _w = 1.24~ms$$\end{document} . This is ten times faster than the activation time constant for M-type potassium which is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _u = 12~ms$$\end{document} when voltage is at its resting value. Full expressions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _w$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _u$$\end{document} are given in the Appendix.

We list the resulting parameter set for the auditory deprived model in the right column of Table 1. Parameter values that are the same for both compartments and for both control and deprived models are reversal potentials ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_K = -106$$\end{document}  mV, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_h = -47$$\end{document}  mV, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{Na} = 55$$\end{document}  mV) and the axial conductance connecting the two compartments ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{ax} = 56.4$$\end{document}  nS).

Numerical methods

We used python code to perform the numerical calculations described above to simulate auditory nerve activity. We used a forward Euler method with time step of 0.002 ms written in C code to solve the differential equations that describe the dynamics of the MSO neuron model. Due to the stochastic nature of the model, most results are reported as mean and standard error of 100 repeated simulations of the model. To calculate firing rates, we counted spikes over the duration of a 300 ms long simulation, with spikes defined as upward crossings of -20 mV in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} variable (AIS voltage). Code is freely available for inspection and use at https://github.com/jhgoldwyn/auditoryDeprivedMSO.Fig. 2. Spiking activity in the auditory nerve (AN) model. Responses to constant amplitude pulse trains. Raster plots show spike times in 50 repeated simulations of the AN model in response to constant amplitude pulse trains (A1) pulse rate is 250 pps and pulse amplitude is 55 nA, (B1) pulse rate is 500 pps and pulse amplitude is 60 nA, and (C1) pulse rate is 1000 pps and pulse amplitude is 75 nA. As specified in the model, spike times are aligned to the pulse offset with a small amount of spike time jitter, as evident in the post-stimulus histograms (PSTHs) shown in the second row for each stimulus setting (PSTHs are computed from 1000 repeated trials using a bin width of 0.1 ms). In (D): average spike rate in response to a 300 ms constant amplitude pulse train, plotted as a function of pulse amplitude for various pulse rates. Values are averages from 1000 repeated simulations. AN firing rates saturate at one spike per pulse in responses to lower pulse rate stimuli (100 pps and 250 pps). In (E): the same firing rates as (D) but plotted as firing rate divided by the number of pulses to yield an estimate of average number of spikes per pulse. Spike probability saturates at one spike per pulse in the responses to the lower pulse rate stimuli and remains lower in the responses to the higher pulse rate stimuli

Results

AN responses to CI stimulation at varying pulse rates

We used the first stage of the model to create auditory nerve fibers (AN) spike trains in response to pulses of current that represent cochlear implant stimulation. In the model, AN spiking aligns closely with the ending times of pulses, as shown in the raster plots in the first row of Fig. 2. As we described in Sect. 2, to create temporal dispersion of spike times around the time of pulse offset, we added a small amount of random jitter to each spike time (standard deviation of 0.1 ms). This amount of jitter is not substantial relative to the interpulse intervals, particularly at lower pulse rates (apparent by the sharp peaks in the post-stimulus time histograms (PSTHs) in the second row of Fig. 2). At the lowest pulse rate shown (250 pps, panel A), the interpulse interval is 4 ms. This is sufficiently long to exceed the refractory period. The decrease in firing rate visible over the first 5 ms to 10 ms in this case, and also in response to higher pulse rates (visible in the decreasing amplitudes of peak heights in the PSTHs), is due to the phenomena of adaptation and accommodation. At the highest pulse rate in the figure (1000 pps, panel C), AN spike times can appear more disorganized across trials than the tightly phase locked responses to low pulse rates. This is particularly the case in the later portions of the responses where adaptation and accommodation take effect.

An important distinction between responses to low-rate and high-rate inputs is the way in which AN firing rate increases with current level (Fig. 2D). In response to low pulse rates, there is a small range of currents over which the spike rate of the AN goes from zero (no spiking) to maximal firing rate (saturating at one spike per pulse). This range is determined by the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} in Eq. (3), the steepness of the firing efficiency curve. Thus, at high current levels and lower pulse rates, the AN can fire reliably in response to every pulse. In response to high pulse rate stimuli, firing rates increase more gradually over a larger range of current values. One way to view this distinction is to notice that, in response to high pulse rate stimuli, the AN can achieve high firing rates at relatively lower average number of spikes per pulse. (Fig. 2E).

MSO responses to CI stimulation at varying pulse rates

In the second stage of the model, we use simulated AN spike trains as inputs to a two-compartment model of an MSO neuron. Excitatory synaptic conductances drive voltage responses in the input compartment ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} in Eq. (17)). This variable represents voltage in the soma (recall we omit dendritic processing from our model). Traces of soma voltage responses to constant amplitude pulse trains of various pulse rates are shown in Fig. 3 (thick curves, black in control model and gray in deprived model). The second compartment of the model represents the site of spike initiation in the axon initial segment (AIS). Voltages in this compartment are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} in Eq. (17) and portrayed in Fig. 3 as thinner, orange curves. MSO spikes appear as large amplitude events in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} traces. We constructed the model to have strong forward coupling and weaker backward coupling (Goldwyn et al., 2019). Consequently, spikes back-propagate weakly into the first compartment. Spike amplitudes are smaller and less visible in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} than in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} .Fig. 3. Voltage traces in each compartment of the MSO neuron model in response to cochlear implant stimulation. Inputs are: (A1-4): 250 pps at 53 nA, (B1-4) 500 pps at 60 nA, and (C1-4) 1000 pps at 80 nA. These stimulus settings were selected so the firing rate is roughly 150 sp/s in the control model in all cases. Responses of the control model are shown in rows 1 and 2. Responses of the deprived model are shown in rows 3 and 4. The inputs have 0 ms time difference in rows 1 and 3, and have non-zero time differences in rows 2 and 4 (1 ms time difference for 250 pps and 500 pps inputs, and 0.5 ms time difference for 1000 pps inputs). Thick black or gray curves are voltages in the soma/dendrite compartment ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_1$$\end{document} in Eq. (17)). Thin orange curves are voltages in the axon initial segment compartment ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} in Eq. (17)). Spikes are defined as upward crossings by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} of the threshold \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-20$$\end{document}  mV

For the pulse amplitude used in the example simulations shown in Fig. 3, the MSO neuron model entrains reliably to low pulse rate stimuli with 0 ms time differences in both control (Fig. 3A1) and auditory deprived (Fig. 3A3) parameter settings. For higher pulse rates (500 pps in column B, and 1000 pps in column C), the model cannot generate spikes rapidly enough to entrain to the input on a pulse-by-pulse basis. The firing rates of both MSO models (control and deprived versions) decrease in response to very high-rate inputs, so we used substantially larger-amplitude input currents for the 1000 pps pulse train when creating Fig. 3C (pulse amplitudes were 53 nA and 60 nA for the 250 pps and 500 pps stimuli, respectively; but 80 nA for the 1000 pps stimulus).

To assess the effect of AIS plasticity on the sound localization function of MSO neurons, we calculated firing rates for various time differences (the lag time) between the two pulse trains that represent binaural inputs. In rows 2 and 4 of Fig. 3, we show responses to inputs with non-zero time difference (1 ms in the case of 250 pps and 500 pps stimuli, and 0.5 ms in the case of 1000 pps). Spiking activity in the control model decreases for larger time differences for lower pulse rate stimuli (compare Fig. 3A1 and A2, or Fig. 3B1 and B2), as is expected for a temporally-precise neural coincidence detector. Intriguingly, whereas there is a relatively large difference between spiking activity in response to inputs with zero time difference and non-zero time difference in the control model for the 500 pps input (Fig. 3B1 and B2), there is no difference in the spiking activities of the deprived model to these same inputs (Fig. 3B3 and B4). This indicates that time difference sensitivity in the deprived model is degraded relative to the control model, an observation we explore in more detail in the remainder of this work.

In Fig. 4, we show firing rate responses of the control model (black curves) and the auditory deprived model (gray curves) to constant amplitude, 0 ms time difference pulse trains over a range of current levels and pulse rates. Firing rate responses to the lower pulse rate inputs saturate at the pulse rate (one spike per pulse) at high current levels (Fig. 4A and B). There is no such saturation effect at higher pulse rates for the range of current levels we tested. The MSO neuron control model is most excitable at intermediate pulse rates (500 pps input in Fig. 4) and higher pulse rate stimuli evoke lower firing rates over the same range of pulse amplitudes (Fig. 4D and E).Fig. 4. Spike rate for different current strengths, using 300 ms inputs with 0 ms time difference and constant pulse amplitude (current values given on horizontal axis). Pulse rates in each panel are (A) 100 pps, (B) 250 pps, (C) 500 pps, (D) 750 pps, and (E) 1000 pps. Firing rates are mean of 100 repeated simulations calculated at intervals of 1 nA. Error bars are standard deviations around those means, plotted at intervals of 5 nA. Firing rates increase with current and saturate at the pulse rate for the lower pulse rates (panels A and B). The auditory deprived model is more excitable than the control model (firing begins at lower current levels). Firing rates in the deprived model exceed that of the control model at the lowest and highest pulse rates and at low current levels, but the firing rate of the control model exceeds the firing rate of the auditory deprived model at intermediate pulse rates for large input current levelsFig. 5Time difference tuning curves for control model and auditory deprived model for varying pulse rates: (A) 100 pps at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 51~nA$$\end{document} , (B) 250 pps at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 53~nA$$\end{document} , (C) 500 pps at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 60~nA$$\end{document} , (D) 750 pps at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 70~nA$$\end{document} , (E) 1000 pps at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 80~nA$$\end{document} . Spike rate values on the vertical axis are mean of 100 repeated simulations and error bars are standard deviation of those values. Time difference on the horizontal axis varies from 0 (in phase) to 1 ms, except in the case of higher pulse rates (500 pps, 1000 pps). In these cases, the largest time difference possible is half the interpulse interval. Spike rates decrease more as time difference increases for the control model than the auditory deprived model, indicating greater time difference sensitivity for the control model

The parameter changes we made to the MSO neuron model (Table 1) reflect changes to axon initial segment observed in auditory brainstem neurons during periods of auditory deprivation (Kuba et al., 2010; Kuba, 2012; Kuba et al., 2014, 2015). These changes have been interpreted as a homeostatic response that compensates for the loss of synaptic drive with an increase in intrinsic excitability (Yamada & Kuba, 2016). Consistent with this notion, we find our auditory deprived model to be more excitable than the control model over a range of pulse rates (gray curves in Fig. 4). This is true in two senses. First, the current level at which spiking activity initiates is lower in the auditory deprived model than the control model across all pulse rates. Second, at all current levels for the lowest pulse rate (Fig. 4A) and the highest pulse rate stimuli (Fig. 4E), and at lower current levels for the intermediate pulse rates (Fig. 4B-D), the auditory deprived model fires at higher rates than the control model. In response to the highest current values and intermediate pulse rates there appears to be some firing rate saturation in the deprived model that allows the firing rate of the control model to exceed the firing rate of the deprived model (Fig. 4B-D, at high current levels).Fig. 6. Summary of firing rate sensitivity to time differences. Spike rates of MSO neuron model to constant amplitude pulse trains for various pulse rates: **(A) ** 100 pps pulse rate, **(B) ** 250 pps pulse rate, **(C) ** 500 pps pulse rate, **(D) ** 750 pps pulse rate, and **(E) ** 1000 pps pulse rate. Maximal firing rates occur for inputs with 0 ms time difference (upper solid curves). Minimal firing rates occur for larger time differences (lower solid curves, using 1 ms time difference for lower pulse rates and half the interpulse intervals for 750 pps and 1000 pps). Values are averages over 100 repeated trials. The shaded regions show the range of possible firing rates between these upper and lower firing rates. Time difference sensitivity, characterized here as the range of firing rates (vertical height of shaded region), is greatest at intermediate pulse rates and greater for the control model than the auditory deprived model

Time-difference tuning is degraded in the auditory deprived model

To assess time-difference sensitivity in our MSO neuron model, we simulate responses to inputs with time delays between the pulse train inputs to the AN model. This time-shift between the two pulse trains represents an interaural time difference between two sounds arriving at the two ears of a listener using bilateral implants. MSO neurons are coincidence detectors, known to fire maximally to inputs that are aligned in time (Svirskis et al., 2004; van der Heijden et al., 2013). The extent to which firing rates decrease as the time differences between the two inputs increases is a straightforward measure of the time difference sensitivity of the model neuron.

We find time difference sensitivity, indicated by firing rates that depend on time difference, in the responses of the control model to pulse rates below 1000 pps (black curves in Fig. 5). The decrease in firing rate with increasing time difference is most pronounced at the intermediate pulse rate (500 pps, panel C in the figure).

Overall, firing rates in the auditory deprived model are flatter with respect to changes in time difference, indicating that the parameter changes we made to create the auditory deprived model degraded time difference sensitivity. Indeed, firing rate responses of the deprived model to the lowest pulse rate (100 pps, Fig. 5A) and to the two higher pulse rates (750 pps in Fig. 5D and 1000 pps in Fig. 5E) are constant with respect to changes in the relative timing of the two pulse train inputs. These tuning curves were computed at specific input levels (see caption of Fig. 5). Next, we show that the same result holds over a large range of current levels.

In Fig. 6 we display time difference sensitivity across a large range of stimulus conditions by plotting firing rates at 0 ms time difference with respect to current level (same curves as in Fig. 4), along with firing rate responses to stimuli with larger time difference. We used 1 ms for this larger time difference in cases for which this did not exceed half the interpulse interval. For the pulse rates that exceeded this limit, we set the maximum possible time difference to be half the interpulse interval (Fig. 6D, E). We plotted firing rates in response to pulse trains with these larger time differences on the same axes, and shaded the region in between to display the range of firing rates produced by the model neuron between these two extreme time difference values.Fig. 7. Envelope time-difference sensitivity. Spike rates of MSO neuron model to 3000 pps pulse trains with sinusoidally amplitude-modulated pulse amplitudes. The sinusoidal oscillation of pulse amplitudes varies between a minimum level of 70 nA and a maximum level shown on the x-axis. The modulation frequencies are: (A) 0 Hz, unmodulated pulse train, (B) 100 Hz modulation frequency, (C) 200 Hz modulation frequency, and (D) 300 Hz modulation frequency. Maximal firing rates occur for inputs with 0 ms time difference (upper solid curves). Minimal firing rates occur for non-zero time difference (lower solid curves, 0.16 ms for the unmodulate pulse train in (A), and 1 ms time difference in all other cases). Values are averages over 100 repeated trials. The shaded regions show the range of possible firing rates between these upper and lower firing rates. Time difference sensitivity, characterized here as the range of firing rates (vertical height of shaded region) is greater for the control model than the auditory deprived model for the higher modulation frequencies (C, D)

The control model (dark gray regions) exhibits a small amount of time difference sensitivity in its responses to low pulse rate stimuli (100 pps in Fig. 6A) and high pulse rate stimuli (1000 pps in Fig. 6E). Responses to the intermediate pulse rate, however, show large differences between inputs with no time difference and inputs with a large time difference (1 ms, in this case). The response regions of the auditory deprived model (light gray regions) are substantially narrower than the response regions of the control model (dark gray regions) across all pulse rates. Thus, by this measure of firing rate differences, the auditory deprived model is less sensitive to time differences than the control model.

Envelope time-difference sensitivity in responses to high-rate pulse trains with time-varying pulse amplitude

In typical listening conditions for CI users, the amplitude of each pulse changes over time due to the dynamically changing intensity of sounds. It is also of interest, then, to study time difference sensitivity in the case of pulse amplitudes varying with time. To this end, we simulated responses to high-rate pulse trains with sinusoidally-varying pulse amplitudes. We exhibit time difference tuning as above, shading the firing rate regions between firing rates to inputs with and without a time shift. We used a high-rate pulse train that, if unmodulated, evokes responses that are constant with time difference (Fig. 7A, note the maximum time difference in this case is 0.16 ms, half the interpulse interval).

If this high-rate pulse train carries a signal with time-varying pulse amplitudes, however, the model exhibits time difference sensitivity (Fig. 7C and D). We used 1 ms as the maximum time difference in the case of the amplitude-modulated pulse trains. We find that the control model exhibits time difference tuning, with increasing time difference sensitivity (larger shaded response regions) at the higher modulation frequencies that we tested. The deprived model, in contrast, does not exhibit similar levels of time difference sensitivity to amplitude-modulated pulse trains. Even at the highest modulation frequency we tested (300 Hz) the auditory deprived model does not exhibit substantial firing rate differences between inputs with 0 ms time difference and inputs with larger time differences (Fig. 7D).Fig. 8. Firing type is phasic for control model and tonic for auditory deprived model. AIS voltage response ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_2$$\end{document} , voltage in the second compartment) in response to steps of current injections (50 ms duration), with current level ranging from 400 pA to 2200 pA in steps of 300 pA. In A (left column): control model responds with a single spike at stimulus onset (phasic firing type) to inputs above the threshold level of approximately 1360 pA. In B (right column): deprived model responds with repetitive firing (tonic firing type) to inputs above approximately 1250 pA. The threshold for a single onset spike in the auditory deprived model is 500 pA

Firing type differences in responses to constant current inputs

We demonstrated that time difference sensitivity in the deprived model is degraded relative to the control model. By this we mean that time difference sensitivity, as measured by the difference between maximum and minimum firing rates (the vertical height of the shaded regions in Fig. 6), is greater for the control model across a large range of input current values. A possible question regarding this result is whether it is simply due to over-excitability in the deprived model. In other words, perhaps the deprived model fires at high rates, generally, including in response to inputs with non-zero time difference. While the deprived model is more excitable, in line with the physiological reports on which we based our model, this alone does not explain why time difference sensitivity is diminished in the deprived model. In particular, if we were to compare tuning curves at fixed firing rates, rather than comparisons at fixed current levels (as in Fig. 5), we would still see greater time difference sensitivity in the control model. This can be deduced by comparing the widths of the response regions in Fig. 6 at locations at which the firing rates are matched, even if the current level is not.

What, then, underlies the degraded time difference sensitivity in the auditory deprived model? We answer this question by pointing out that the type of spiking dynamics changes due to AIS plasticity, in addition to the degree of excitability or average firing rate. By this we mean that control model is phasic firing – it responds to a constant-level current injection with (at most) a single spike at the stimulus onset (Fig. 8A). In contrast, the deprived model fires repetitively to constant-level current inputs (Fig. 8B). This type of response is known as tonic firing.

MSO neurons are phasic firing cells (Svirskis et al., 2002). Phasic firing is commonly associated with temporally-precise coincidence detection (Meng et al., 2012; Huguet et al., 2017) and therefore would be beneficial for sound localization. By changing the firing type of the model from phasic to tonic, the model loses this dynamical specialization for temporal precision. As a result, time difference sensitivity in the control model exceeds that of the auditory deprived model, as we have shown throughout. The parameter manipulations that caused this change in firing type include reductions to the low-threshold potassium and h-type conductances in the soma and transformation of low-threshold potassium conductance in the AIS compartment to M-type potassium (see Table 1). We made these changes based on the experimental and modeling study of AIS plasticity in NM neurons under conditions of auditory deprivation in Kuba et al. (2015).

Discussion

The use of cochlear implants by individuals with severe to profound sensorineural hearing loss can increase their speech understanding, particularly in noise-free environments (Krueger et al., 2008; Gaylor et al., 2013). Individuals with hearing loss in both ears who use a second implant (bilateral CI) can achieve higher rates of speech understanding than users of a single implant (Litovsky et al., 2006). Bilateral implants, by providing auditory inputs to both ears, have potential to improve perception of the aspects of sound that require binaural processing. Sound source location is one such percept. Two of the primary cues for sound source location are interaural level differences (ILDs) and interaural time differences (ITDs). While bilateral implant users can experience improved sound source localization relative to users with a single CI (Brown & Balkany, 2007), sound source localization and spatial hearing remain a deficit for CI users relative to normal hearing listeners (Grantham et al., 2007).

When evaluating the limitations of CIs, it is essential to consider the relative extents to which CI users’ hearing outcomes are limited by the auditory periphery and central stages of the auditory system. The peripheral factors that hearing outcomes with CIs are numerous. These include the signal processing algorithms that convert acoustic inputs to electrical stimuli, the engineering design of the devices, damage to the cochlea during surgery, and pathologies in the cochlea and auditory nerve that can depend on the etiology of hearing loss (what has been termed the electrode-neuron interface Bierer, 2010). While these factors certainly affect hearing outcomes for CI users, it is also possible that hearing outcomes are impaired by more central stages of neural processing in the auditory pathway.

Hearing development passes through critical periods of development that may be disrupted by hearing loss early in life (Kral, 2013; Persic et al., 2020). Hearing loss later in life and even into adulthood can also lead to neural plasticity that reshapes auditory cortex (Eggermont, 2008). In addition to these wide-ranging, circuit-level changes to the auditory system during hearing loss, periods of auditory deprivation in birds (as a model for hearing loss) can induce cellular-level physiological changes in individual neurons. Kuba et al. (2010) found axon initial segments in the nucleus magnocellularis of chickens grew during periods of auditory deprivation, and Kim et al. (2019) found the AIS in the medial nucleus of the trapezoid body in deaf mice to be larger and more distant from the soma than in normal-hearing mice.

AIS plasticity has not been studied in the MSO of humans, to our knowledge. We used Kuba’s studies of NM and NL (nuclei in the sound source localization circuit of birds), as motivation to model similar changes in MSO neurons. Specifically, we compared a control model of MSO neural dynamics to a modified model, which we called the (auditory) deprived model. Based on the experimental findings and computational modeling by Kuba and colleagues, we constructed the auditory deprived model by reducing the amount of low-threshold potassium conductance in the soma (Kuba et al., 2015) (first compartment of our two-compartment model), transforming the low-threshold Kv1.1 potassium channel into to M-type Kv7.2 potassium channels in the AIS region (Kuba et al., 2015) (second compartment in our two-compartment model), and increasing the size of the AIS (Kuba et al., 2010).

We consistently found the deprived model to be more excitable than the control model (Figs. 4 and 8) with poorer time difference sensitivity (as measured by the range of maximum to minimum firing rates in time difference tuning curves, see Figs. 5 and 6). Although neither model was sensitive to small differences in the timing of high-rate, constant-amplitude pulse trains (Fig. 6E), the control model exhibited greater time difference sensitivity to amplitude-modulation envelopes of high pulse rate carriers than the deprived model (Fig. 7).

Relation to previous models

Our approach to simulating MSO response to CI stimuli is similar to a previous modeling study by Chung et al. (2015). That paper builds on earlier work (Colburn et al., 2008) by incorporating a more detailed auditory nerve model. Similar to the model we used (Boulet, 2016; Tabibi et al., 2021), Chung et al. model auditory nerve responses to CI stimuli as a stochastic process governed by a voltage-like variable that evolves with linear dynamics and a probabilistic spike generation rule based on a dynamic threshold with additive Gaussian noise. Spike threshold in their model is dynamic to describe a refractory period (see also (Bruce et al., 1999a, b)) and also includes spike rate adaption (Nourski et al., 2006). The model we present follows this same general approach and includes refractory and adaptation effects. In addition, we include two additional effects: accommodation (a long time-scale decrease in spike probability) and facilitation (temporal summation of subthreshold inputs to increase spike probability on a short time-scale). These modifications, developed by Boulet, Tabibi and colleagues (Boulet, 2016; Tabibi et al., 2021), enhance the model by capturing additional nuances of neural responses to CI stimuli.

Chung and colleagues found that their MSO neuron model produced little to no spiking activity in response to unmodulated pulse trains at rates of 500 pps and higher. They modified their model to produce higher firing rates by strengthening synaptic inputs and increasing the speed of membrane dynamics. No such manipulation was needed in our model to drive MSO activity in response to high-pulse rate stimuli. That being said, time difference sensitivity in our model degrades at higher pulse rates, with no significant time difference sensitivity in responses to 1000 pps stimuli. This limitation at high pulse rates is consistent with a recent experimental study showing ITD sensitivity in adult rats implanted with bilateral CIs for pulse rates as high as 900 pps, but not for 1800 pps pulse trains (Buck et al., 2023). It is possible that by using fewer, stronger synaptic inputs and faster membrane dynamics, we could obtain time difference sensitivity in response to higher-pulse rate stimuli, but we have not tested that in our simulations. In our conversion of the control model to the auditory deprived model, we did alter the speed of membrane dynamics, but made no changes to synaptic inputs (Kuba et al., 2010). Following observations by Kuba and colleagues and their modeling study (Kuba et al., 2015), we removed a portion of the low-threshold potassium current (Kv1.1-type channels) from the soma to create the auditory deprived model. This has the effect of increasing the input resistance at the soma and slowing the soma dynamics. Our results, therefore, are consistent with typical findings that faster membrane dynamics improve time difference sensitivity and neural coincidence detection (Kempter et al., 1998).

The changes we made to the second compartment, to reflect the processes of AIS plasticity during sound deprivation, enlarged the surface area of this compartment and replaced low-threshold potassium with M-type potassium conductance. By enlarging the spike-initiation region in the deprived model, we increased backpropagation of voltage from the AIS into the soma. In a previous work, we described this as an increase in backward coupling (Goldwyn et al., 2019). In that prior study, we found that the greatest time difference sensitivity occurs with strong forward coupling and weak backward coupling between the soma and axon. The effect of auditory deprivation in our model is to move the soma-axon coupling out of this optimal configuration. Removing the low-threshold potassium conductance from the spike-generating region, where it functions as a source of dynamic negative feedback (Huguet et al., 2017), also diminishes time difference sensitivity in models of neural coincidence detection (Goldwyn et al., 2019).

Limitations of the model

In the first stage of our model, we used an up-to-date model of temporal response properties of auditory nerve fibers to CI stimulation (Boulet, 2016; Tabibi et al., 2021). The accommodation component of this model distinguishes it from previous models constructed along similar principles (Bruce et al., 1999b; Goldwyn et al., 2012). Accommodation is particularly relevant for modeling high-rate pulse trains as it strongly suppresses spiking under these conditions (Fig. 2). The auditory nerve stage of our model was identical in the control and deprived conditions. Degeneration of auditory nerve fibers certainly impacts neural responses to CI stimulation and listening outcomes for CI users. In order to isolate the possible effects of AIS plasticity in the MSO due to sound deprivation, we did not model any peripheral changes.

An absence in our CI and AN modeling is that we do not describe how the spatial spread of electric fields in the cochlea activates a distributed population of AN fibers, as has been done in other models (Rattay et al., 2001; Cohen, 2009; Goldwyn et al., 2010; Kalkman et al., 2015; Malherbe et al., 2015, for examples). Understanding how population-level auditory nerve activity flows into the MSO circuit remains an important question. On this point, some have suggested that MSO neurons may receive inputs from tonotopically misaligned input (Shamma et al., 1989; Bonham & Lewis, 1999; Day & Semple, 2011). For bilateral CI users, there is almost certainly some tonotopic mismatch in the inputs to binaural circuits because the surgical placements of electrode arrays in two different ears are not perfectly aligned. Mismatched place-of-stimulation does impair lateralization of sounds by binaural CI users (Kan et al., 2013). Future modeling studies of MSO responses to CI stimulation could consider these spatial and population-level effects.

In the second stage of our model, we used a two-compartment model of MSO neural dynamics, similar to previous MSO studies (Khurana et al., 2011; Goldwyn et al., 2019). With the two-compartment structure we can separate the integrative region of the cell where inputs arrive (soma and dendrite regions, represented by the first compartment) from the spike-generating region (in the AIS, represented by the second compartment). Our goal was to understand how AIS plasticity affects time difference sensitivity in the MSO. With this model structure, we represented the AIS and modified its parameters in accordance with AIS plasticity in a straightforward way. Morphological features that we did not include in our model also affect ITD tuning in MSO neurons, including their bipolar dendritic structure (Carr et al., 1998; Dasika et al., 2007; Mathews et al., 2010), and the distribution of multiple spike initiation sites along its axon (Lehnert et al., 2014). If further evidence regarding AIS plasticity in MSO neurons becomes available, our model could be elaborated to incorporate additional morphological detail.

The connection from the auditory nerve stage to the MSO stage is direct in our model. We omitted nuclei intermediate to these two stages (anteroventral cochlear nucleus, AVCN; and medial nucleus of the trapezoid body, MNTB). One possible consequence of this omission is that afferent inputs to MSO may be even more precisely-timed than what we used in our model. In response to acoustic stimulation, it is known that bushy cells in AVCN can enhance phase-locking so that spike timing of inputs to MSO is more precise relative to spike timing in the auditory nerve (Joris et al., 1994). In the AVCN and MNTB of gerbils, the jitter spike times in responses to pulse trains of CI stimuli is very small, with median values reported to be roughly one-fifth the size the spike time jitter of 0.1 ms that we use in our model (a value based on recordings in cat auditory nerve in Miller et al., 2001). These studies were done under conditions of acute deafening. It is possible that mechanisms of phase-locking enhancement in auditory brainstem may be degraded during periods of prolonged auditory deprivation. For instance, if transmission along axons or at synapses were severely disrupted during periods of hearing loss, then this could substantially impact MSO activity, as illustrated in the study by Budak et al. (2022). We also note that we did not include inhibitory inputs to MSO neurons, instead focusing on an “excitatory-excitatory” input motif for MSO neurons. Given the great interest in how inhibitory inputs impact sound localization processing in the MSO (Brand et al., 2002; Grothe, 2003; Jercog et al., 2010; Roberts et al., 2013; Myoga et al., 2014), future work could include this pathway.

Elements of the excitatory and inhibitory pathways exhibit activity-dependent refinement (Kapfer et al., 2002; Walmsley et al., 2006). Some of these changes have been described as homeostatic, such as increased excitatory post-synaptic currents entering the AVCN through the endbulb of Held in congenitally deafened mice (Oleskevich & Walmsley, 2002; Walmsley et al., 2006). MNTB neurons in congenitally deafened mice exhibit greater excitability due to changes in their complement of potassium channels (Walmsley et al., 2006; Leão et al., 2010). Changes in sodium channel development have been characterized in the MNTB of congenitally deaf mice (Leão et al., 2006) and AIS plasticity occurs in MNTB neurons during auditory deprivation (Kim et al., 2019). AIS plasticity occurs in the avian nucleus magnocellularis during sound deprivation (Kuba, 2012), so it is possible that similar changes occur in bushy cells in AVCN (the mammalian homologue to nucleus magnocellularis). A recent modeling study showed that demyelination of AN axons and synaptopathy of the inner hair cell to AN synapse can disrupt the sound localization circuit (Budak et al., 2022). In particular, they observed that these pathological changes could cause substantial reductions in MSO firing activity. Their work was done in the context of hidden hearing loss and acoustic stimulation. Similar work in the context of CI stimulation could add to understanding of deficits in sound localization by bilateral CI users.

Implication for sound localization by bilateral cochlear implant users

Cochlear implants can improve speech understanding for individuals with profound to severe hearing loss, and yet challenges remain. At the periphery of the auditory pathway, there are factors that can limit CI effectiveness. These include the signal processing performed by the devices, the coarse nature of neural activation by CI stimulation (spatially and temporally), and degradation of the electrode-neuron interface (death or degeneration of AN fibers, for instance). In this study, we used parameter variations in an MSO neuron model to explore how pathological changes to downstream stages of the auditory pathway may also limit how sound location information can be provided by CIs. In particular, we found that the changes we made to create our “auditory-deprived model” disrupted the temporally-precise neural computations that underlie sound source localization, in the sense that firing rate in the control model typically decreased with larger time difference more than in the deprived model (more sharply-peaked time difference tuning curves).

A point that we emphasize is that processes that may appear homeostatic in nature may alter the dynamics and function of neural systems in consequential ways. Changes to a neuron’s properties that increase excitability (to compensate for diminished excitatory drive under conditions of auditory deprivation), can transform the firing type of a neuron. Consequently, changes that appear homeostatic with respect to firing threshold or average firing rate can alter the functioning of neurons and neural circuits.

Many contemporary CI speech processing strategies use high-rate pulse trains. The Nucleus cochlear implant (Cochlear Ltd.), for instance, has a default stimulation rate of 900 pps and manufacturer-proposed settings that range from 500 pps to 2,400 pps (Wathour et al., 2021). Two features of our simulation results have direct bearing on the high-rate pulse train stimulation strategies used in contemporary CIs. First, we observed in the model a decrease in time-difference sensitivity at high pulse rates (1000 pps and higher). The range of available ITDs becomes increasingly limited at high pulse rates. For a 1000 pps pulse train, for instance, the maximum possible ITD is 0.5 ms. More generally, the loss of ITD sensitivity at high pulse rates is consistent with the typical understanding of the MSO as a nucleus that extracts ITD information from low-frequency sounds. Second, MSO neurons can be sensitive to time-differences conveyed by high pulse rate stimuli when those time differences are contained in slowly-varying pulse amplitudes (envelope ITDs, as opposed to fine-structure ITDs). This result is consistent with a previous modeling study of envelope ITD sensitivity in MSO responses to CI stimuli (Colburn et al., 2008). It also supports physiological evidence that the MSO neuron, although typically thought of as a low-frequency nucleus, is sensitive to ITDs in the envelopes of high-frequency inputs (Gai et al., 2014).

Lastly, it is important to point out that we have considered the possibility that there is AIS plasticity in MSO neurons during periods of auditory deprivation (envisioned here as the period of hearing loss that precede implantation), but we have not considered the possibility that AIS plasticity may occur post-implantation. After CIs are turned on and evoke neural activity, can they initiate homeostatic processes that reduce neural excitability and reverse the AIS plasticity processes that take place during auditory deprivation? Understanding the plasticity that occurs post-implantation, and the extent to which such processes can improve hearing outcomes for CI users, is an intriguing line of future work.