Phase Locking of Auditory-Nerve Fibers Reveals Stereotyped Distortions and an Exponential Transfer Function with a Level-Dependent Slope

Refractoriness affects phase locking in a frequency-dependent manner

Refractoriness limits the spike rate of ANFs and also slightly alters the shapes of the period histograms relative to those that would be observed in the absence of refractoriness (Avissar et al., 2013). Refractoriness has also been proposed to be responsible for the temporal asymmetry of period histograms (Horst et al., 2018). We removed the effects of refractoriness from each spike train to yield period histograms showing the recovered instantaneous rate of presumed synaptic release events rather than the rate of spikes (Fig. 2C). Figure 3A shows that the recovered mean event rate q_mean exceeds the mean spike rate r_mean by a factor that increases with increasing mean spike rate from just above 1 to ∼2.3. Figure 3B shows, as a function of stimulus frequency, the ratio of the vector strength computed from the recovered instantaneous event rate to the vector strength computed from the instantaneous spike rate. Ratios >1 represent a decrease of vector strength by refractoriness, and ratios <1 represent an increase of vector strength by refractoriness. Refractoriness typically decreased vector strength for frequencies <∼400 Hz and increased it for frequencies between ∼400 and ∼800 Hz. In either case, however, the effect is small (typically <5%). For frequencies >800 Hz, refractoriness has virtually no effect on vector strength. These findings are consistent with observations from chicken ANFs (Avissar et al., 2013, their Fig. 7b) when expressed with respect to the duration of the refractory period: vector strength is decreased by refractoriness when frequencies are low enough that the stimulus period is considerably longer than the mean refractory period; vector strength is increased by refractoriness when the stimulus period is only somewhat longer than the mean refractory period; and vector strength is unaffected by refractoriness when the stimulus period is shorter than the mean refractory period. Figure 3C shows, as a function of stimulus frequency, the difference between the phase angles of the mean vectors computed from the recovered instantaneous event rate and the instantaneous spike rate, expressed in units of time. This difference is always positive, meaning that refractoriness advances the phase angle of the mean vector (Fig. 2C). The advancement is typically small (<0.1 ms) and decreases with increasing stimulus frequency. On this double logarithmic plot, a slope of −1 would correspond to a constant difference in the phase angles when expressed in radians. Because the slope of the decrease is steeper than −1, the difference in the phase angles decreases with increasing frequency. The difference also decreases as the rate of release events decreases (gray scale).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Recovery of the rate of release events from spike trains. A, Raster plot (top) and poststimulus time histogram (bottom) for an example spike train in response to 50 repetitions of 100 ms tones of 211 Hz and 70 dB SPL, presented at 4 Hz (A3-U8-R3-C10; CF, 353 Hz; SR, 89.8 spikes/s). Only the first 50 ms of the 150 ms interstimulus intervals are shown. Gray vertical lines in the poststimulus time histogram mark complete stimulus cycles between 10 and 100 ms. Note that, over this time range, the spacing of the peaks matches the stimulus period, but that the interval following the peak before 10 ms is shorter. B, The excitability function (i.e., the probability that the ANF is not refractory), shown for the first 30 ms of the example spike train (gray line; top), and collapsed into a single period for all complete cycles between 10 and 100 ms (gray lines; bottom). Spike times and corresponding excitabilities are marked by white dots. To aid visibility, the vertical portions of the gray lines following each spike are omitted from the bottom panel. The mean excitability was computed for each time point in the cycle (black line). C, Smoothed period histograms showing the instantaneous spike rate (black) and recovered instantaneous event rate (gray) for the example spike train. D–G, Simulations were performed to examine the validity of the recovery method. D, The simulated (true) instantaneous event rate (dotted line), the instantaneous spike rate resulting from application of refractoriness (black), the recovered instantaneous event rate (gray), and the absolute error of the recovery method (light gray line), for one example simulation. E, Relative error as functions of the maximum and mean event rates. The inset shows that the recovered event rates match the true event rates up to high values. F, Relative error as a function of the stimulus frequency. G, Relative error as functions of the absolute and mean relative refractory periods assumed in the recovery method (true values were 0.6 ms each).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Refractoriness affects phase locking in a frequency-dependent manner. A, Ratio of the recovered mean event rate to the mean spike rate plotted as a function of mean spike rate. B, Ratio of the vector strength computed from the recovered instantaneous event rate to the vector strength computed from the instantaneous spike rate, plotted as a function of stimulus frequency. Ratios >1 represent a decrease of vector strength by refractoriness. Ratios <1 represent an increase of vector strength by refractoriness. Five ratios >1.1 are not shown. C, Difference between the phase angles of the mean vectors computed from the recovered instantaneous event rate and the instantaneous spike rate, expressed in units of time by dividing by 2π and the stimulus frequency f₁. The difference is plotted as a function of stimulus frequency (horizontal axis) and the mean rate of release events (gray scale). This difference is always positive, meaning that refractoriness advances the phase angle of the mean vector. In B and C, the frequency was randomly jittered to aid visualization. Dashed vertical lines mark the inverse of the mean refractory period (1.2 ms) assumed in the recovery method. T denotes the stimulus period. The white dot marks the example shown in Figure 2A–C.

Together, these analyses show that the dominant effect of ANF refractoriness is to scale down the rate (Fig. 3A). The small changes in vector strength (Fig. 3B) imply only small changes in the shapes of the histograms, regardless of any phase shifts (Fig. 3C). Refractoriness therefore is not the mechanism responsible for skewing the period histograms, in contrast to what has been proposed (Horst et al., 2018), or for preventing peak clipping at high SPLs.

Phenomenological model of phase locking

We seek to relate the period histogram to the instantaneous pressure during a cycle of the tone stimulus having peak amplitude P₁ and frequency f₁ (Fig. 4A). It is not our intention to present a biophysical model of ANF phase locking that incorporates detailed cochlear mechanics, IHC physiology, synaptic transmission, and spike-generating mechanisms. Instead, we attempt to use the simplest and fewest assumptions possible. The first necessary assumption is that the response driven by the stimulus is delayed relative to the AC voltage driving the speaker (Fig. 4B, top). Such a delay is accounted for by a phase-shift parameter (ϕ₁). The second necessary assumption is that the stimulus driving the response acquires harmonic distortion components, each with its own frequency (f₂ = 2·f₁, f₃ = 3·f₁, etc.), amplitude (P₂, P₃, etc.), and phase (ϕ₂, ϕ₃, etc.; Fig. 4B, middle and bottom). These distortion components are not present in the acoustic stimulus (Johnson, 1980), but are instead generated by nonlinearities in the auditory periphery (Kim and Molnar, 1979; Sachs and Young, 1980; Ruggero and Rich, 1983; Zwislocki, 1986; Cody and Mountain, 1989; Dallos and Cheatham, 1989; Kiang, 1990; Rhode and Cooper, 1996; Cheatham and Dallos, 1998; Cooper, 1998; Khanna and Hao, 1999; Mountain and Cody, 1999; Guinan, 2012). Consistent with the reported peripheral nonlinearities, we assume that distortion components are present already in the forces deflecting the IHC stereocilia. For convenience, we describe these forces using the units of the acoustic stimulus [pascals (Pa)] and refer to them as the “mechanical drive.” The waveform of the mechanical drive, P(t) (Fig. 4C), is therefore as follows (Eq. 14):

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Phenomenological model of phase locking, illustrated using results obtained from a single example spike train (A6-U14-R1-C20; CF, 1300 Hz; SR, 48.3 spikes/s; f₁, 1300 Hz; 76 dB SPL). A, Acoustic stimulus. B, Phase-shifted stimulus of frequency f₁ (top) and harmonic distortion components of frequencies f₂ = 2·f₁ and f₃ = 3·f₁ (middle and bottom). C, Mechanical drive obtained by summing the individual frequency components in B. D, IHC transfer function relating the mechanical drive to the rate of release events. Although the model includes a first-order Boltzmann function (gray, shown here with α = 1), the system only uses a portion in the lower tail, such that it can be well described by an exponential function (black). The vertical dotted line marks zero. E, Model event rate (black) resulting from the exponential transfer function in D. The gray bars show the period histogram of events recovered by removing the effects of refractoriness (F) from the empirical period histogram of spikes (G).

Kiang (1990) proposed the existence of a second fundamental component to account for a notch in the rate-level function and an ∼180° phase shift of the response to the fundamental at very high stimulus levels (Kiang, 1984; Liberman and Kiang, 1984). We do not include a second fundamental component because the sum of two sine waves having the same frequency yields, regardless of their amplitudes and phases, a single sine wave with that frequency. Therefore, for the purpose of fitting the model to period histograms, a single fundamental component suffices. In addition, we find no clear notches in the rate-level functions in our data, perhaps because stimulus levels were not high enough. The mechanical drive deflects the IHC stereocilia to give rise to a mechanoelectrical transducer (MET) current. We assume that the relationship between the mechanical drive and the MET current is described by a first-order Boltzmann function, as expected for a two-state channel that transitions between a closed and an open state (Fettiplace and Kim, 2014). Several alternative functions were also examined (see below Initial model fits reveal that saturation is not reached and that the model's implementation can be simplified). The asymmetry of the MET function around its operating point creates a DC component. The MET current activates a variety of voltage-gated K⁺ channels and changes the membrane potential, which is also characterized by AC and DC components (Russell and Sellick, 1978, 1983; Dallos, 1985; Dallos and Cheatham, 1989; Johnson, 2015; Altoè et al., 2018). For simplicity, we assume that changes in the membrane potential are proportional to changes in the MET current. Low-pass filtering by the IHC membrane (not modeled here) will attenuate the AC component in a frequency-dependent way and thereby affect the AC/DC ratio (Palmer and Russell, 1986). The change in the membrane potential gives rise to a change in the Ca²⁺ signal that initiates transmitter release at the ribbon synapse. For simplicity, we assume that changes in the Ca²⁺ signal are proportional to changes in the membrane potential (e.g., consistent with the finding of Li et al., 2014; see their Fig. 5). Given the assumed relationships, the mechanical drive and the synaptic Ca²⁺ signal are also related by a first-order Boltzmann function. The rate of synaptic release events is assumed to be proportional to some power of the Ca²⁺ signal (Beutner et al., 2001; Brandt et al., 2005; Johnson et al., 2005, 2008, 2009, 2010; Goutman and Glowatzki, 2007; Beurg et al., 2008; Heil and Neubauer, 2010), such that the instantaneous event rate of the model, R_model(t), is given by a first-order Boltzmann function raised to the power α (Fig. 4D, gray function, shown for α = 1), expressed as follows (Eq. 15):

Here, S_max is the maximum possible Ca²⁺ signal (in events ^−α/s^−α), P_0.5 is the mechanical drive (in Pa) at which the Ca²⁺ signal is half the maximum value (i.e., the inflection point of the Boltzmann function), and b is the slope factor (in Pa⁻¹). The maximum possible instantaneous event rate is R_max = S_max^α.

Figure 4E shows the resulting instantaneous event rate (black line) as a function of time in the stimulus cycle. It also shows the period histogram of the events (gray bars) recovered by removing the effects of refractoriness (Fig. 4F) from the period histogram of the spikes (Fig. 4G). Because we have an analytical function for the instantaneous event rate (Eq. 15), but not for the instantaneous spike rate, the model was fitted to the period histogram of the recovered instantaneous event rate. Nonetheless, we also fitted the model equation to the instantaneous spike rates, ignoring the effects of refractoriness. All results were qualitatively very similar to those obtained from fits to the recovered instantaneous event rates (reported below), which is consistent with the finding that the major effect of refractoriness is to scale down the rate (Fig. 3).

Initial model fits reveal that saturation is not reached and that the model's implementation can be simplified

Initial fitting of the model to the recovered instantaneous event rates revealed that it was impossible to reliably estimate S_max, α, or R_max (which is derived from S_max and α; see above Phenomenological model of phase locking). To investigate why this was the case, we systematically fitted the model with α fixed to 1 and R_max (equal to S_max when α = 1) fixed to several different values. Two harmonic distortion components were included in the model, resulting in seven free parameters (P₂, P₃, ϕ₁, ϕ₂, ϕ₃, b, and P_0.5). Figure 5A shows the resulting transfer functions for a single example histogram (gray lines; assuming two upper harmonic components). The vertical dotted lines mark the portion of each transfer function covered by the acoustic stimulus. Notably, aside from the four or five functions that saturate at the lowest values of R_max, the shapes of the transfer functions are virtually identical over this portion. Furthermore, this portion is well described by an exponential function (black line). A similar exercise was repeated using seven log-spaced values of R_max from 250 to 16,000 events/s for all 2675 period histograms that satisfied our inclusion criteria. For each fit, the negative log likelihood was computed and then normalized by dividing it by the minimum obtained across all R_max. Figure 5B shows these normalized values as functions of R_max (gray lines). If there were an optimal value for R_max, these functions would be U-shaped. But this is not the case. Instead, they tend to decrease rapidly and monotonically toward 1, and the negative log likelihood hardly changes once R_max is large enough, which explains why R_max could not be estimated reliably. This shows that the data are well described by an exponential transfer function for all SPLs examined (≤100 dB SPL). Replacing the Boltzmann function in Equation 15 with the exponential function yields the following (Eq. 16): where C is a scale factor and D is a slope factor. The normalized negative log likelihoods for this exponential transfer function (Fig. 5B, exp) are virtually identical to those for the Boltzmann transfer function and large R_max. The log likelihoods obtained with Equation 16 are identical for all values of α, because C^α represents a single scaling factor and α·D represents a single slope factor, such that C and D can perfectly compensate α. This explains why α (and S_max) could not be estimated reliably. The model therefore reduces to a single exponential transfer step (Fig. 4D, black function) between the instantaneous mechanical drive P(t) and the instantaneous event rate R_model(t), expressed as follows (Eq. 17): which has two fewer parameters than the full model (Eq. 15). Here, A is the scale factor that specifies the rate (in events/s) at P(t) = 0 Pa (i.e., the operating point) and accounts for any DC component of the response, and B is the slope factor (in Pa⁻¹) such that R_model changes by a factor of e^B when the pressure changes by 1 Pa. Note that a DC shift in the input to the exponential transfer function would also be included in the scale factor A.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Justification for the use of an exponential transfer function and two distortion components. A, First-order Boltzmann transfer functions (gray lines) estimated for a single example histogram (A6–U14-R1–C20; CF, 1300 Hz; SR, 48.3 spikes/s; f₁, 1300 Hz; 76 dB SPL), with α fixed to 1 and R_max fixed to several different values. Aside from the four or five functions having the lowest values of R_max, the shapes of the functions are virtually identical over the portion used in the fits (between vertical dotted lines) and are well described by an exponential function (black line). B, Normalized negative log likelihoods for all period histograms (gray lines) and their geometric mean (circles and black line), plotted as a function of R_max. These functions tend to decrease rapidly and converge to the result for an exponential transfer function (exp). C–E, Residuals (data minus model) for spike trains recorded from the same ANF as in A, in response to CF stimuli presented at 4–76 dB SPL in 4 dB steps. Residuals are shown with no (C), one (D), or two (E) distortion components included. The mean number of events per cycle ranged from 0.042 at low stimulus levels to 0.12 at high stimulus levels. F, Exponential transfer functions obtained from the model fits (shown for 76 dB SPL) depend little on the number of distortion components included. G, Normalized negative log likelihoods for all period histograms (gray lines) and their geometric mean (circles and black line), plotted as a function of the number of distortion components.

We also examined alternative functions used in other models, specifically the second-order Boltzmann function expected for a three-state channel that transitions between two closed states and an open state (Sumner et al., 2002; for review, see Fettiplace and Kim, 2014; Altoè et al., 2018) and a logarithm-based function (Zhang et al., 2001; Zilany et al., 2009; Bruce et al., 2018). Despite introducing additional free parameters, these functions generally did not perform better than the exponential function (data not shown). The tangent hyperbolic function used by Carney (1993) can be scaled and shifted to be mathematically equivalent to a first-order Boltzmann function and would therefore yield identical fits and parameter values.

Two harmonic distortion components are sufficient to describe most responses

We next systematically examined how many harmonic distortion components would be needed to satisfactorily account for the data using the exponential transfer function. We fitted Equation 17 to each period histogram of recovered events, assuming zero, one, two, or three distortion components in the mechanical drive (Eq. 14). Figure 5C–E shows results for a representative set of data recorded from one ANF in response to CF stimuli of 1300 Hz with levels from 4 to 76 dB SPL. Each function shows the residuals for one stimulus level. Residuals were computed by subtracting the number of events obtained from the model from the number of events obtained from the data. Without distortion components included (Fig. 5C), the residuals are pronounced, particularly for high SPLs, and oscillate around zero with a frequency of two times the stimulus frequency. When one distortion component is included (Fig. 5D), the residuals are much smaller and oscillate around zero with a frequency of three times the stimulus frequency. When two (Fig. 5E) or three (data not shown) distortion components are included, the residuals are smaller still and without a clear oscillatory pattern. These observations emphasize the need for including distortion components in the model. The transfer functions obtained from the model fits for this example depend little on the number of distortion components included (Fig. 5F, shown for 76 dB SPL). Figure 5G shows, for all period histograms, the normalized negative log likelihoods as functions of the number of distortion components (gray lines). These functions tend to decrease markedly with the inclusion of the first distortion component, and then less so with each subsequent component. Although a third component meaningfully improved the fits for a few histograms, we included only two distortion components to avoid overfitting. The relative improvement of the fits when adding distortion components tends to be greatest for responses to low stimulus frequencies (data not shown). Note that although the first and second distortion components would contribute to the Fourier spectrum of the period histogram, the spectrum will also reflect additional components introduced by the exponential transfer function, including additional first and second distortion components differing in phase from those in Equation 14. This means that the distortion components estimated from the model fits are not identical to those seen in Fourier spectra computed from period histograms (Sachs and Young, 1980). Detailed examination of the distortion components estimated from the model might therefore provide new insights into the mechanical drive to the IHC stereocilia.

Harmonic distortions in the mechanical drive are stereotypical, both in the absence and presence of peak splitting

All period histograms were fitted by the exponential model (Eq. 17) with the fundamental (f₁) and two distortion components (f₂ and f₃) in the mechanical drive (Eq. 14). This model has seven free parameters (P₂, P₃, ϕ₁, ϕ₂, ϕ₃, A, and B). Figure 6 shows the resulting parameter estimates of the f₂ and f₃ components (P₂, P₃, ϕ₂, ϕ₃) for all examples having a significant vector strength for the respective component. Figure 6A,B shows the estimated amplitudes of these components, P₂ and P₃, relative to the stimulus amplitude P₁ [in dB, computed as 20·log(P_i/P₁)] and plotted as functions of P₁ (in dB SPL). The color indicates the vector strength for f₂ and f₃, not for f₁ (Fig. 6A, scale). Period histograms characterized by clear peak splitting (Fig. 1F) yield parameter values circled in black. Figure 6A shows that P₂ tends to be ∼15 dB below P₁ at low stimulus levels and increases to ∼10 dB below P₁ at high stimulus levels. Figure 6B shows that P₃ tends to be ∼20 dB below P₁ at all stimulus levels. The vector strengths for f₂ and f₃ also increase with increasing stimulus level and with increasing relative amplitude of these components. Data points from period histograms characterized by peak splitting do not form separate clusters, but instead tend to lie near the upper edges of the bands of data points. In these cases, P₂ approaches or even exceeds P₁ and the vector strengths for f₂ and f₃ are rather high. Figure 6C shows P₂/P₁ and P₃/P₁ plotted as functions of the stimulus frequency. These relative amplitudes decrease only minimally with increasing frequency. Although significant phase locking is sometimes obtained for distortion components having relatively high frequencies (Fig. 6C, rightmost pink dots, for which f₁ is ∼2–4 kHz and f₃ is therefore ∼6–12 kHz), a visual inspection of the corresponding period histograms leads us to suspect these results are spurious. Figure 6D shows P₂/P₁ and P₃/P₁ plotted as functions of the distance in octaves of the stimulus frequency from CF. Relative amplitudes decrease only minimally over the range of distances examined. Figure 6E,F shows the estimated phases of the f₂ and f₃ components, ϕ₂ and ϕ₃ (horizontal axis), relative to the phase of the fundamental, ϕ₁, and plotted as functions of tone frequency (vertical axis). The relative phases are given respectively by ϕ₂/2 − ϕ₁ and ϕ₃/3 − ϕ₁. They are fairly similar across ANFs and conditions, such that the data points form a narrow band in each plot. The horizontal axis spans one period of f₁ and therefore spans two periods of f₂ and three periods of f₃, such that the same bands are shown two and three times, respectively. Color indicates the vector strength for f₁, and not for f₂ and f₃ as in Figure 6A,B. For both the f₂ and f₃ components, there is a subtle but clear tendency for the relative phase to increase with increasing frequency. These effects might be due to frequency-dependent phase shifts introduced by low-pass filtering (e.g., due to the IHC membrane capacitance). Most of the low vector strengths obtained for frequencies <1 kHz are from period histograms with peak splitting, marked by black circles in one band of data points in Figure 6E,F. For the f₂ component, the relative phases from these histograms align with the left edge of the band. This is because such a phase relationship produces two prominent peaks in the mechanical drive at values of P₂/P₁ lower than those required with other phase relationships. Figure 6G shows a polar histogram of the relative phases from Figure 6E,F. The relative phase of the f₂ component (blue) leads the fundamental (black) by ∼π/12 radians, and the relative phase of the f₃ component (pink) leads the fundamental by ∼π/7 radians. Figure 6H shows the typical waveform (black) that results from summing the fundamental (gray) and the typical f₂ (blue) and f₃ (pink) components (the typical values obtained from the model fits were as follows: P₂/P₁ = −13.1 dB and P₃/P₁ = −20.9 dB, ϕ₂/2 − ϕ₁ = 2.88 radians and ϕ₃/3 − ϕ₁ = 1.64 radians). The instantaneous amplitude of the typical waveform has a steeper leading edge and a shallower trailing edge than a sinusoid of the same frequency, which is similar to most period histograms (Fig. 1A,C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Parameter estimates of the harmonic distortion components. Results are shown for period histograms having significant vector strengths for both the fundamental (f₁) and the respective (f₂ or f₃) distortion component. A, B, Estimated amplitudes, P₂ and P₃, of the distortion components relative to the amplitude P₁ of the fundamental (the stimulus amplitude). The ratios are expressed in dB and plotted as functions of the stimulus amplitude in dB SPL. Color indicates the vector strength for the respective distortion component (see scale in A). Parameter values circled in black indicate period histograms characterized by peak splitting. C, D, Relative amplitudes of the f₂ and f₃ components plotted as functions of the stimulus frequency f₁ (C) and the distance in octaves of f₁ from CF (D). In A–D, six values of P₂/P₁ and nine values of P₃/P₁ are <−40 dB and are not shown. E, F, Estimated phases of the f₂ and f₃ components (horizontal axis), relative to the phase of the fundamental and plotted as functions of tone frequency (vertical axis). The horizontal axis spans one period of f₁ and therefore spans two periods of f₂ and three periods of f₃. Color indicates the vector strength for the fundamental. Parameter values circled in black in one of the bands indicate period histograms characterized by peak splitting. G, Polar histogram of the relative phases from E (blue) and F (pink). H, The typical waveform (black) that results from summing the fundamental component (gray) and the typical f₂ (blue) and f₃ (pink) components.

Harmonic distortions in the mechanical drive eliminate the need for assuming hysteresis in the transfer function

Figure 7 shows the instantaneous event rate plotted as a function of the instantaneous mechanical drive estimated from the model fits without (rows A, C, E) and with (rows B, D, F) harmonic distortion components included. All examples are from the same three ANFs as in Figure 1. Data are shown for the majority of SPLs for which there were ≥125 spikes and for which the vector strength was significant. Most panels include results for two adjacent SPLs. Model fits are shown as solid lines, whereas data are shown as symbols connected by dashed lines to form Lissajous figures. For display purposes, data are plotted using 30 bins per cycle. At each SPL, the data are well described by the exponential transfer function. Note that the steepness of the function (and the scaling of the horizontal axis) varies with SPL, an issue we return to in the next section. With only the fundamental included (three free parameters: ϕ₁, A, and B), the data can form loops that systematically deviate from the model function. This effect is most obvious for the ANF showing peak splitting (Figs. 1F, 7E), but is also visible in the other examples at high SPLs. Such looping behavior in the data could be interpreted as hysteresis in the transfer function (Patuzzi and Moleirinho, 1998; Brown et al., 2009; Meenderink and van der Heijden, 2011; Horst et al., 2018). However, the loops are reduced or eliminated with f₂ and f₃ components included in the model (Fig. 7, compare rows B, D, F, rows A, C, E).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

A–F, Instantaneous event rate as a function of instantaneous mechanical drive estimated from model fits without (rows A, C, E) and with (rows B, D, F) two harmonic distortion components. Examples are the same three as in Figure 1. Model fits are shown as solid lines, whereas data are shown as symbols connected by dashed lines in the form of Lissajous figures. Data are shown using 30 bins per cycle. Without distortion components included in the model, the data can deviate from the model in a systematic loop, most obviously for the ANF showing peak splitting (row E; Fig. 1F). Systematic deviations are reduced or eliminated when distortion components are included. Also note that the steepness of the functions (and the scaling of the horizontal axes) varies with SPL.

The exponential transfer function changes with stimulus level

Figure 7 shows that the model with two distortion components describes the data well and that the exponential transfer function changes with stimulus level (except at low SPLs in Fig. 7D).

The slope factor

Figure 8 shows the dependence of the slope factor B on stimulus level, stimulus frequency, spontaneous rate, vector strength, and relative amplitude of the f₂ component. Figure 8D shows the combinations of stimulus frequency and ANF spontaneous rate for which parameter estimates were obtained. Horizontal dashed lines demarcate spontaneous-rate (SR) groups of <1 spike/s, between 1 and 10 spikes/s, and >10 spikes/s, whereas vertical dashed lines demarcate stimulus-frequency groups <500 Hz, between 500 and 2000 Hz, and >2000 Hz. For each point, the corresponding B-versus-level function is shown in the panel indicated in parentheses. Gray lines connect adjacent stimulus levels. Data points for one B-versus-level function are marked by squares. The example marked in Figure 8A is that shown in Figures 1C,D and 7C,D. The example marked in Figure 8H is that shown in Figures 1A,B and 7A,B. The color (except in Fig. 8G) indicates the vector strength for the stimulus frequency and shows that vector strength tends to increase with increasing stimulus level, decreasing stimulus frequency, and decreasing spontaneous rate, which is consistent with findings in other studies (for review, see Heil and Peterson, 2017).

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

A–C, E, F, H, I, Dependence of the slope factor B on stimulus level, stimulus frequency, and spontaneous rate. Gray lines connect adjacent stimulus levels. The color indicates the vector strength for the stimulus frequency (f₁). One example function in each panel is marked by squares. The dotted negative diagonal in each panel represents BP₁ = 1. The arrows in B show how sensitivity and vector strength increase within this coordinate system in the absence of harmonic distortions in the mechanical drive. D, Combinations of stimulus frequency and ANF spontaneous rate for which parameter estimates were obtained, with corresponding B-versus-level functions shown in the panel indicated in parentheses. G, Vector strength as a function of BP₁. Color indicates the relative amplitude of the f₂ component. The white line represents the theoretical relationship between vector strength and BP₁ in the absence of distortions.

For all ANFs and stimulus frequencies, B decreases with increasing stimulus level. For ANFs having high SRs and stimulated with tones of low or intermediate frequencies (Fig. 8A,B), B changes little at low stimulus levels and then decreases with a slope of approximately −1 in these log–log plots (i.e., B is approximately proportional to the inverse of the stimulus amplitude; we refer to this as a “linear decrease”). For ANFs having high SRs and stimulated with tones of relatively high frequencies (Fig. 8C) and for ANFs having intermediate (Fig. 8E,F) and low (Fig. 8H,I) SRs, B decreases linearly with stimulus amplitude over the entire range of stimulus levels for which parameter estimates were obtained. Nonetheless, these B-versus-level functions probably also have flat portions, which fail to show up here due to our criteria of significant vector strength and ≥125 spikes, a conclusion supported by the data of Horst et al. (2018). When stimulus frequencies are high (Fig. 8C,F,I), low-pass filtering (e.g., by the IHC membrane capacitance) would attenuate the AC component and prevent vector strengths from reaching significance when stimulus levels are low and spikes are too few. In other cases, vector strength is already near its maximum at stimulus levels at which the criterion of ≥125 spikes is reached (Fig. 1B).

The flat portion of a B-versus-level function corresponds to the range of stimulus levels over which the slope of the transfer function changes little or not at all, but the dimensionless product BP₁ of the slope factor and the stimulus amplitude (i.e., the concentration parameter) increases with increasing stimulus level (Fig. 9A). Over this range of stimulus levels, the vector strength increases with level and there is an overlap of the functions relating the recovered instantaneous event rate to the instantaneous mechanical drive for different stimulus levels (as in Fig. 7C,D at low stimulus levels; Horst et al., 2018). The linearly decreasing portion of each B-versus-level function corresponds to the range of stimulus levels over which BP₁ is approximately constant (Fig. 9A). This range probably extends to higher stimulus levels than were used here, because none of the B-versus-level functions shows any sign of asymptoting at high stimulus levels. These results are consistent with the proposals of Colburn (1973) and Heinz et al. (2001), who assert that the concentration parameter of the von Mises distribution increases with stimulus level before saturating. In the absence of distortion components, a constant BP₁ means that the shape of the period histogram and the vector strength are constant across stimulus levels due to the relationship between vector strength and BP₁ (Eq. 9). Figure 8G shows this relationship (white line) and also shows the vector strengths obtained from the period histograms of recovered events (symbols, where the color indicates P₂/P₁). Weak distortions in the mechanical drive yield vector strengths near the function described by Equation 9, whereas strong distortions yield predominantly lower vector strengths.

• Download figure
• Open in new tab
• Download powerpoint

Figure 9.

A, Example BP₁-versus-level functions corresponding to the B-versus-level functions in Figure 8B, with color also indicating vector strength. B, Maximum BP₁ obtained from each BP₁-versus-level function, plotted as a function of stimulus frequency.

The positions of the B-versus-level functions in Figure 8 along the negative diagonal (e.g., the positions of the knee points) might relate to the sensitivities of the ANFs to the stimulus frequency (Fig. 8B, left arrow). Such sensitivity is a function of CF and a function of the distance of the stimulus frequency to the CF. The linearly decreasing portions of the B-versus-level functions also differ in their positions along the positive diagonal. The higher the function along this diagonal, the larger BP₁ and therefore the higher the maximum vector strength (Figs. 8B, right arrow; 9A, vertical arrow). The positions of the functions along this diagonal, and therefore the maximum vector strength, depend little on spontaneous rate, which is consistent with the small effect of spontaneous rate on maximum vector strength reported in other studies (Johnson, 1980; Joris et al., 1994b; Louage et al., 2004; Dreyer and Delgutte, 2006; Temchin and Ruggero, 2010; for review, see Heil and Peterson, 2017). The position of the functions along the positive diagonal, and therefore the maximum vector strength, is primarily a function of stimulus frequency (Fig. 1E; see Fig. 11, compare columns; Palmer and Russell, 1986; Javel and Mott, 1988; Weiss and Rose, 1988; Joris et al., 1994a; Köppl, 1997; Temchin and Ruggero, 2010; Versteegh et al., 2011; Heil and Peterson, 2017) because low-pass filtering reduces B (and therefore BP₁) as frequency increases.

Figure 9A shows example BP₁-versus-level functions corresponding to the B-versus-level functions in Figure 8B. Figure 9B shows the maximum BP₁ obtained from each BP₁-versus-level function, plotted as a function of stimulus frequency. This plot is analogous to Figure 1E, which shows the maximum vector strength as a function of stimulus frequency. The maximum BP₁ decreases with increasing frequency over the entire frequency range and more clearly than does the vector strength. This is a consequence of the fact that the vector strength saturates but BP₁ does not. Using BP₁ as a measure of phase locking avoids the need to rescale vector strength near its upper limit when higher resolution is desired (Johnson, 1978, 1980; Joris et al., 1994a,b).

The event rate at the operating point

Figure 10 shows the dependence of the scale factor A on stimulus level, stimulus frequency, spontaneous rate, vector strength, and P₂/P₁, in the same format as Figure 8 and with the same examples marked. For ANFs having high SRs and stimulated with tones of low frequencies, A tends to decrease slightly with increasing stimulus level (Fig. 10A). Several of the functions are rather erratic, particularly at high stimulus levels. For ANFs having high SRs and stimulated with intermediate (Fig. 10B) and high (Fig. 10C) frequencies, A respectively changes little or increases slightly with increasing stimulus level. For ANFs having intermediate SRs, A increases more clearly with increasing stimulus level before saturating (Fig. 10E,F). For ANFs having low SRs, this increase is even more pronounced (Fig. 10H,I). For all ANFs and stimulus frequencies, A would have to approximately equal the spontaneous rate at stimulus levels near the rate threshold. This is seen in the data of Horst et al. (2018), but is not shown in Figure 10 because period histograms for such low levels did not satisfy our criteria of significant vector strength and ≥125 spikes. At high stimulus levels, A saturates at values considerably higher than the spontaneous rate. The ratio of the value of A at saturation to the spontaneous rate is largest for low-SR ANFs and smallest for high-SR ANFs. Vector strength is a function of A/R_mean (Eq. 11), which is shown in Figure 10G (white line). Weak distortions in the mechanical drive yield vector strengths near the function described by Equation 11, whereas strong distortions yield predominantly lower vector strengths. When both A and BP₁ are constant (as is often the case at moderate and high stimulus levels; compare Figs. 8, 10), and provided the relative amplitudes and phases of the distortion components remain constant, the mean rate, the shape of the period histogram, and the vector strength are constant across stimulus levels. This is equivalent to complete compression of the response (i.e., a growth of 0 dB/dB).

• Download figure
• Open in new tab
• Download powerpoint

Figure 10.

A–C,E,F,H,I, Dependence of the scale factor A on stimulus level, stimulus frequency, and spontaneous rate. Gray lines connect adjacent stimulus levels. The color indicates the vector strength for the stimulus frequency (f₁). One example function in each panel is marked by squares. D, Combinations of stimulus frequency and ANF spontaneous rate for which parameter estimates were obtained, with corresponding A-versus-level functions shown in the panel indicated in parentheses. G, Vector strength as a function of A/R_mean. Color indicates the relative amplitude of the f₂ component. The white line represents the theoretical relationship between vector strength and A/R_mean in the absence of distortions.

Exemplary changes in the transfer function with stimulus level

Figure 11 shows, for the examples marked in the middle and right columns of Figures 8 and 10, how the model exponential transfer function changes with stimulus level. The instantaneous event rate is plotted on a logarithmic axis, such that the transfer functions form straight lines. Color indicates stimulus level. The extent of each line shows the used range of the transfer function. This range is not necessarily symmetrical around P(t) = 0 Pa, because the minimum and maximum of the mechanical drive depend on the distortion components present (Fig. 6H). In the absence of distortions, however, the drive would be symmetrical. Furthermore, when plotted on a logarithmic axis, the instantaneous event rate in the period histogram would vary symmetrically around the value at the operating point. In each example in Figure 11, the transfer functions are steepest (i.e., the sensitivity is highest) at low stimulus levels, such that small changes in the pressure lead to large changes in the response. With increasing stimulus level, the transfer functions become shallower (i.e., the sensitivity decreases). This gain control is observed well below CF (Fig. 11A), at CF (Fig. 11B–E), and above CF (Fig. 11F). The ratio of the maximum to minimum instantaneous rates becomes attenuated when the stimulus frequency increases (e.g., across columns in Fig. 11), likely due to low-pass filtering. Low-SR ANFs show pronounced growth of the DC component with increasing stimulus level, and therefore their ratio of the mean rate to the spontaneous rate (Fig. 11, dashed horizontal lines) becomes higher.

• Download figure
• Open in new tab
• Download powerpoint

Figure 11.

A–F, Changes of the model transfer function with stimulus level, shown for the examples marked in the middle and right columns of Figures 8 and 10. Note that the instantaneous event rate is plotted on a logarithmic axis, such that the exponential transfer functions form straight lines. Color indicates stimulus level, with the lowest in blue and the highest in red. The extent of each line shows the range of the transfer function used. Dashed horizontal lines mark the estimated spontaneous rate of release events.

Rate adaptation does not account for the change in the slope of the transfer function

Here, we examine whether the change in the slope factor B can be accounted for by spike-rate adaptation occurring on time scales longer than the stimulus period. For each cycle of the tone stimulus, starting from the onset, we computed the vector strength and the mean number of spikes per cycle across all repetitions (proportional to the rate). This analysis was restricted to spike trains in which at least half of the cycles contained ≥10 spikes and had significant vector strength (622 spike trains). Figure 12A–E shows the results from five example ANFs, with only significant vector strengths shown. There is pronounced spike-rate adaptation in each example, which is consistent with previous reports (for review, see Heil and Peterson, 2015). Figure 12A–C shows data from ANFs for which the vector strength remains relatively constant across cycles. Vector strength was similarly constant for the other stimulus levels not shown here. Delgutte (1980, his Fig. 13) and Joris et al. (1994a,b, their Fig. 7 in each) also show examples of approximately constant vector strength across cycles despite spike-rate adaptation. Assuming distortions are constant (or absent) over time, the finding that vector strength is constant across cycles means that the slope factor B is also constant. This is because vector strength is a function of BP₁ (Eq. 9) and P₁ is constant (except during the brief stimulus onset and offset). Figure 12D shows data from an ANF for which the vector strength decreases over time, which could correspond to a gradual decrease in B or a gradual change in distortions. Figure 12E shows data from an ANF for which the vector strength increases over time, which could correspond to a gradual increase in B or a gradual change in distortions. More repetitions than were collected would be needed to investigate these possibilities. In any case, the changes in B over time that would potentially account for the changes in vector strength are small compared with the orders-of-magnitude changes in B occurring with changes in stimulus level (Fig. 8).

• Download figure
• Open in new tab
• Download powerpoint

Figure 12.

Vector strength and mean number of spikes per cycle as functions of the cycle number. A–C, Data from three spike trains for which the vector strength remains relatively constant across cycles while the spike rate adapts. D, Data from a spike train for which the vector strength decreases over time. E, Data from a spike train for which the vector strength increases over time. In A–E, vector strength is plotted only when significant. F, Values of Spearman's rank correlation coefficient ρ computed for each function relating significant vector strength to cycle number (as in A–E), plotted as a function of stimulus frequency.

To examine whether monotonic relationships between vector strength and cycle number are typical at the population level, we computed Spearman's rank-order correlation coefficient for each spike train. Figure 12F shows that most (532 of 622) coefficients scatter around zero and are not significant (p > 0.01), indicating that vector strength might typically be constant across cycles. Thirty of the 34 significant negative coefficients were obtained from a single ANF (A2-U3, same ANF as in Fig. 1F). Forty-nine of the 56 significant positive coefficients were observed for frequencies >1 kHz and were obtained from seven ANFs. These data show that the slope factor B is largely set by mechanisms independent of those responsible for spike-rate adaptation occurring on time scales longer than the stimulus period. The mechanism must also be very fast, because in most cases B is constant starting from the beginning of the response.

The power–law relationship between the synaptic Ca²⁺ signal and the rate of release events

Because period histograms are well described by an exponential transfer function, it is impossible to use fits of the model to the period histograms of recovered release events or spikes to estimate the power α in the model's relationship between the synaptic Ca²⁺ signal and the rate of release events (Eq. 16). To estimate the power, it would be necessary to know the Ca²⁺ signals at individual synapses in response to the acoustic stimuli. Such data are not available, but the assumed linearity between the membrane potential and the Ca²⁺ signal, and the experimental support for this assumption (Li et al., 2014), means that the membrane potential is a reasonable proxy for the Ca²⁺ signal. Only few in vivo recordings of the membrane potential of IHCs in response to low-frequency tone stimuli are available. These come from the pioneering work of Russell and Sellick (1983, their Figs. 1, 2) and of Dallos (1985, his Fig. 5; Dallos, 1986, his Fig. 6), and Dallos and Cheatham (1989, their Figs. 1, 10), all in guinea pig. No recordings have been published for cat. Sharp electrodes such as those used for these recordings likely alter the total conductance of the IHC. Such alterations in turn alter the magnitude of the changes in the membrane potential (Zeddies and Siegel, 2004; Altoè et al., 2018), but are unlikely to change the shape of the waveform. The waveforms shown in the literature are similar to one another, regardless of whether recorded from basal or apical IHCs. We therefore chose to compare our model results to the receptor potential recordings shown in Figure 1 of Russell and Sellick (1983), because this figure offers the highest resolution. Figure 13B shows these basal IHC recordings in response to 300 Hz stimuli presented at 90 (top) and 80 (bottom) dB SPL. Figure 13A shows, for our data, vector strength as a function of the rate asymmetry, which is the ratio of the difference between the maximum rate reached and the rate at zero pressure to the difference between the rate at zero pressure and the minimum rate reached. The color of the symbols indicates P₂/P₁. The white line shows an analogous function constructed using Equations 9 and 13. Because the receptor potential contains negative values, it is impossible to calculate vector strengths for the recordings from Russell and Sellick (1983) and to determine their corresponding positions along this function. Nonetheless, we can compare the shapes of these recorded waveforms to the shapes of the waveforms derived from our model fits. For this comparison, we chose to use period histograms having high vector strengths and relatively low distortions, located in the region marked by the small rectangle in Figure 13A. Figure 13C (top) shows, for each histogram in this region, two cycles of the model waveform corresponding to the membrane potential if the power were α = 1 (gray lines). An approximately linear relationship between the Ca²⁺ signal and the rate of release events has been proposed and interpreted to indicate a nanodomain control of release (Brandt et al., 2005; Kim et al., 2013). Each waveform is normalized to its peak value. The mean waveform (red line) and its value at 0 Pa (red dotted line) are also shown. Figure 13D (top) shows the normalized transfer functions (gray lines) that give rise to the waveforms in Figure 13C. The mean transfer function (red line) and its operating point (white circle) are also shown. Figure 13C,D (bottom) show the corresponding waveforms and transfer functions if the power were α = 3 (obtained by taking the cube root of the waveforms and transfer functions for α = 1). A power of ∼3, derived previously from the dependence of ANF spike rate on stimulus level (Heil et al., 2011) and of ANF first-spike latency on stimulus level and onset duration (Heil et al., 2008; Neubauer and Heil, 2008; Heil and Neubauer, 2010), has been interpreted to indicate a microdomain control of release. Note that for α = 3, the transfer function is shallower and the waveform more sinusoidal than for α = 1. Figure 13E shows six cycles of the mean waveforms for α = 1 (red line) and α = 3 (blue line) overlaying the center six cycles of the membrane-potential recordings of Russell and Sellick (gray lines, replotted from the boxed region in Fig. 13B). The model waveforms have been scaled to approximately match the peaks and troughs of the recordings. Note that the function for α = 3 matches the waveform of the recordings well, whereas the function for α = 1 has a much sharper peak and a much broader trough. The red and blue dotted lines mark the outputs of the respective transfer functions at 0 Pa. We also obtained rough estimates of the output at 0 Pa from the data of Russell and Sellick (1983). Precise estimates are impossible to derive due to the presence of distortions that affect the locations of the peaks and troughs (Fig. 6H). The rough estimates were obtained by identifying the voltages (white dots) at the time points halfway between the maxima and minima. The estimates for the leading and trailing edges differ due to the distortions, but the true value is likely within the range identified. These estimates are also better matched by the model when α = 3 than when α = 1. These comparisons show that, given our model, the IHC membrane-potential recordings and our data would be inconsistent if there were a linear relationship between the Ca²⁺ signal and the rate of release events, but would be consistent if there were a power-of-three relationship.

• Download figure
• Open in new tab
• Download powerpoint

Figure 13.

Recordings of the IHC receptor potential are consistent with a power of three, but not with a power of one, in the model relationship between the synaptic Ca²⁺ signal and the rate of release events. A, Vector strength as a function of the asymmetry of the output of each model transfer function. Color indicates the relative amplitude of the f₂ component. The white line represents the theoretical relationship in the absence of distortions. B, IHC receptor potential recordings from Russell and Sellick (1983, their Fig. 1) in response to 300 Hz tone stimuli presented at 90 (top) and 80 (bottom) dB SPL. C, Top, Two cycles of the waveform corresponding to the membrane potential of the model if the power were α = 1. Gray lines show the waveforms obtained from the histograms within the region marked by the rectangle in A, each normalized to its peak value. The red line shows the mean waveform, and the dotted red line shows its value at 0 Pa. Bottom, Corresponding waveforms if the power were α = 3. D, Transfer functions (gray, red and blue lines) and their operating points (white circles), which give rise to the waveforms in C. For α = 3, the transfer functions are shallower and the waveforms more sinusoidal than for α = 1. E, Six cycles of the mean model waveforms for α = 1 (red line) and α = 3 (blue line) overlaying the center six cycles of the membrane-potential recordings (B, boxed region) of Russell and Sellick. The model waveforms have been shifted and scaled to match the amplitude of the recordings. The red and blue dotted lines mark the outputs of the respective model transfer functions at 0 Pa. The white dots mark rough estimates of the output at 0 Pa from the data of Russell and Sellick. The function and operating point for α = 3 closely matches the waveform and operating point of the recordings of Russell and Sellick, whereas the function for α = 1 has a much sharper peak and a much broader trough, and its operating point is too low.

Relationship between the exponential transfer functions and pseudotransducer functions

Figure 14A,C shows model IHC receptor potentials estimated from model fits to the responses of two ANFs to tones having a range of stimulus levels. For clarity, the estimated distortions and phase shifts are omitted. Receptor potentials were calculated assuming a power-of-three relationship between the Ca²⁺ signal and the rate of release events. The model receptor potentials are similar to the receptor potentials measured in vivo in response to acoustic stimulation (Russell and Sellick, 1983; Dallos, 1985, 1986; Russell et al., 1986a; Cody and Russell, 1987; Dallos and Cheatham, 1989, 1990) or in vitro in response to sinusoidal deflections of the stereocilia (Russell et al., 1986a,b; He et al., 2004). The AC component increases with stimulus level before saturating. Despite this saturation, no peak clipping occurs. Figure 14B,D shows the corresponding exponential transfer functions relating the receptor potential to the instantaneous pressure. Their slopes decrease with increasing stimulus level. Several investigators have plotted the maxima and minima of the voltage responses of IHCs and outer hair cells (OHCs) against the maxima and minima of the stimulus amplitudes (Russell and Sellick, 1983; Russell et al., 1986a,b; Dallos and Cheatham, 1989; Cheatham and Dallos, 2000). The resulting functions have sometimes been referred to as “pseudotransducer” functions (Dallos and Cheatham, 1989, 1990; Cheatham and Dallos, 2000) and, aside from effects such as low-pass filtering by the hair-cell membrane, are assumed to approximate the true transducer functions (Russell, 2008, his Figs. 4, 18). In Figure 14B,D, the white dots connected by black lines represent pseudotransducer functions analogous and similar to those in the literature. These pseudotransducer functions fail to describe the transfer functions obtained for individual stimulus levels (colored lines), except those for low levels that operate far away from the saturating regions of the pseudotransducer functions. As in Figure 14D, many model pseudotransducer functions are nonmonotonic, with a sharp dip at low negative instantaneous pressures. Such dips occur because the transfer function must operate near the resting value when stimulus levels are low, whereas a rapid growth of the DC component can result in transfer functions that operate entirely above the resting value when stimulus levels are higher (Fig. 11C–F). Although the sharp dips are not apparent in published pseudotransducer functions, they should also occur in such cases.

• Download figure
• Open in new tab
• Download powerpoint

Figure 14.

Relationship between the exponential transfer functions and pseudotransducer functions. A, C, Model IHC receptor potentials estimated from fits to the responses of two ANFs. For clarity, the estimated distortions and phase shifts are omitted. B, D, Corresponding exponential transfer functions relating the instantaneous model receptor potential to the instantaneous pressure. White dots connected by black lines represent pseudotransducer functions analogous to those in the literature.