The influence of distributed source regions in the formation of the nonlinear distortion component of cubic distortion-product otoacoustic emissions

Otoacoustic emissions (OAEs) are acoustic signals generated in the cochlea either spontaneously or evoked by an external stimulus (Shera, 2004). They are a by-product of the electromechanical force produced by the soma of the outer hair cells (OHCs), which in turn is the basis of cochlear amplification (Dong and Olson, 2013). These signals are transmitted backward through the middle ear and are recorded as a pressure signal in the ear canal. Evoked OAEs are generated by two different mechanisms: (1) distortion due to the nonlinearity of the mechanoelectrical transducer in the OHC stereocilia and (2) coherent reflection of mechanical waves traveling along the basilar membrane (BM) (Zweig and Shera, 1995; Shera, 2004; Shera and Guinan, 2008). The level of the OAEs generated by both mechanisms is affected—usually decreasing—if cochlear amplification is reduced (Gorga et al., 1997; Kemp, 1978). Therefore, OAEs are routinely used clinically to estimate the functional status of the cochlear amplifier (Blankenship et al., 2018).

OAEs evoked by two simultaneously presented pure tones, with frequencies usually denoted as $f1$ and $f2$⁠, are called distortion product OAEs (DPOAEs). The largest DPOAE is the lower sideband, cubic component at $2f1−f2$⁠. The evoked emission, generated due to the nonlinearity of the mechanoelectrical transducers close to the place of maximum overlap of the two traveling waves produced by $f1$ and $f2$⁠, is not only transmitted retrograde to the ear canal, but also stimulates a wave traveling anterograde towards the tonotopic place of $2f1−f2$⁠. In the case of local perturbations of mechanical properties near the tonotopic place of $2f1−f2$⁠, this anterograde traveling wave can be coherently reflected back along the BM and add constructively or destructively to the initial retrograde traveling wave (Zweig and Shera, 1995; Shera and Guinan, 1999). Therefore, the overall DPOAE can be thought of as being composed of both the nonlinear-distortion and the coherent-reflection components. Conventionally, the generators of these two components are called the primary and the secondary sources, respectively.

The nonlinear distortion component of DPOAEs does not originate from a discrete point-like source. Instead, the primary source is spread along a region of the BM. Specifically, the nonlinear-distortion component is a superposition of distributed primary sources whose phase and amplitude may be different (Shera and Guinan, 2008). The primary-source region is located at and basal to the tonotopic place of $f2$⁠, which is where the traveling-wave envelopes evoked by both tones overlap maximally (Kim et al., 1980; Young et al., 2012). Cancelation between these sources due to phase differences may cause a decrease of DPOAE amplitude (Martin et al., 2010, 2013; Mills, 1997; Young et al., 2012). Martin et al. (2013) measured DPOAEs in rabbits and presented complexities (amplitude overshoots) in the onset and offset of the time-domain signal. Since the complexities disappeared and the steady-state amplitude of DPOAEs increased after adding a third tone with frequency above $f2$⁠, they proposed that the more basal sources of DPOAEs add in opposite phase to the sources near $f2$⁠. They showed this phenomenon for a small ratio of $f2/f1$ (1.07), and also for a larger ratio of 1.25, but using relatively high stimulus levels [$L1=L2=65–75$ dB sound pressure level (SPL)]. It was suggested that mutual cancelation of wavelets leads to a decrease of the amplitude of DPOAEs for $f2/f1$ ratios approaching 1. This suggestion agrees with modeling studies using nonlinear cochlear models (van Hengel, 1996; Talmadge et al., 1998; Shera, 2003; Liu and Neely, 2010; Sisto et al., 2018). On the other hand, the decrease of DPOAE amplitude for frequency ratios less than about 1.2 could be also explained by filtering by the tectorial membrane (TM) (Allen and Fahey, 1993).

Another phenomenon which may play a role in the reduction of DPOAE amplitude is two-tone suppression (Robles and Ruggero, 2001). Kanis and de Boer (1997) suggested that the DPOAE amplitude may decrease because the stimulus tones mutually suppress if their frequency separation is small $(f2/f1≈1.05)$⁠; also see Fahey et al. (2006). Two-tone suppression can partially explain the complexities in the onset and offset of DPOAE signals, as demonstrated with cochlear models (Gummer et al., 2018; Vencovský and Vetešník, 2018b) combined with experimental measurements (Zelle et al., 2018) from human subjects. Kummer et al. (2000) suggested that two-tone suppression can cause the non-monotonic growth of DPOAE input/output functions observed if one of the stimulus levels is held constant and the other is increased (Johnson et al., 2006). In contrast, Sisto et al. (2018) hypothesized that the non-monotonic DPOAE level growth may be due to mutual cancelation of wavelets from the primary DPOAE source.

The aim of this study is to elucidate mechanisms responsible for the nonlinear-distortion component of DPOAEs up to moderate stimulus levels for $f2/f1=1.2$⁠. Simulations of a two-dimensional nonlinear cochlear box model are compared with experimental results measured in normal-hearing human subjects. Previous modeling studies revealed how the primary and secondary DPOAE sources affect both the steady-state DPOAE amplitude and also the envelope of the DPOAE signal (Talmadge et al., 1998; Talmadge et al., 1999; Tubis et al., 2000). Here, we focus on how the distributed nature of the primary source affects the steady-state DPOAE amplitude as well as the onset and offset of the DPOAE signal.

This section describes a hydrodynamical box model of the cochlea (Sec. II A) and an approximate analytical solution of the model at the angular frequency of the lower cubic distortion product (Sec. II B). The analytical solution is proportional to the DPOAE signal in the ear canal. This solution provides a tool for elucidating how DP wavelets generated along the BM might contribute to the nonlinear-distortion component of the DPOAEs. Two simulation techniques are used: (1) numerical solution of the box-model equations (Sec. II A) which provides pressure in the ear canal from which the DPOAE signal (the nonlinear-distortion component in the case of a “smooth” model) can be extracted, and (2) analytical approximation of the solution of the box-model equations (Sec. II B), which allows specification of the primary DPOAE source and hence analysis of the distributed nature of this source. The simulation techniques are described in Sec. II C and the experimental paradigm is presented in Sec. II D. If only interested in the results, the reader may skip most of the modeling details. However, to comprehend the results, the reader is advised to consult Eq. (16), which describes how the nonlinear-distortion component of the DPOAE results from the summation of distributed primary sources of unique phases. In this description, the summation is achieved by a cumulative integration.

The 2D cochlear model used in this study was formulated in Sondhi (1978), Allen and Sondhi (1979), and further developed in Vetešník and Nobili (2006), Vetešník and Gummer (2012). The model is supplemented with a 1D middle-ear model adapted from Talmadge et al. (1998); the equations are given in the  Appendix. The cochlear model assumes that the cochlea has the shape of a box. The box is divided by the BM with the top and bottom chambers filled with incompressible fluid. The BM is modeled as a continuous array of harmonic oscillators with mass $m(x)$⁠, damping $h(x)$⁠, and stiffness $k(x)$⁠, per unit BM length, where $x$ is measured from the basal end of the BM. Parameter values of the model are given in Table I and the most commonly used symbols in the model are listed in Table II. The main assumptions underlying the model are described in Mammano and Nobili (1993) and Nobili et al. (1998). The BM oscillators are coupled longitudinally by hydrodynamical forces and by shearing resistance between the adjacent BM oscillators. The input signal is fed into the model via vibrations of the stapes, which are also coupled to the BM oscillators by hydrodynamical forces. In the model, cochlear amplification derives from the OHC electromechanical force acting in phase with BM velocity and, therefore, against viscous forces acting on the BM, as found experimentally by Dong and Olson (2013). This phase relation is achieved by a second array of oscillators simulating the shearing displacement between the reticular lamina (RL) and the TM. The model is, therefore, based on the amplification mechanism postulated by Gummer et al. (1996).

The transversal displacement $ξ(x,t)$ (positive for displacement towards scala vestibuli) of the oscillators is calculated in the model using the integro-differential equation

where $∂t$ and $∂x$ represent partial derivatives with respect to time, $t$⁠, and $x$⁠, respectively, $GS(x)$ is the Green's function representing the hydrodynamic coupling between the stapes and the BM at $x$⁠, $σ(t)$ is the stapes displacement (positive for displacement into the cochlea), $∂xs(x)∂x$ accounts for the longitudinal coupling of the oscillators due to the shearing viscosity, $s(x)$⁠, and $G(x,x¯)$ is the Green's function representing longitudinal coupling between the BM segments caused by hydrodynamic forces from the surrounding fluid. The acceleration of the oscillator at $x¯$ with a length $dx¯$ produces force per unit BM length at $x$⁠. This force is given by the term $G(x,x¯)∂t2ξ(x¯,t)dx¯$⁠. The Green's function $G(x,x¯)$ can be expressed as a sum of two terms: the long-range term and the short-range term; analytical expressions for these two terms for the box model can be found in Vetešník and Nobili (2006).

Each oscillator is locally amplified by the undamping force, $U(x,t)$⁠, which represents the electromechanical force produced by the OHC soma (Nobili et al., 1998). Similar to Nobili and Mammano (1996), we assume that $U(x,t)$ depends nonlinearly on the hair-bundle displacement $η(x,t)$⁠:

where $a=1 (cm−1)$ and $u(x)$ is a suitable function controlling the degree of undamping along the BM and is given by

$g1$⁠, $g2$⁠, and $g3$ were set empirically so that the desired cochlear amplifier gain at CF was achieved; the chosen values are $g1=1.44$ (⁠$g/s2$⁠), $g2=5.88$⁠, and $g3=0.15$⁠. $S(y)$ is a sigmoidal function (second-order Boltzmann function) given in Nobili and Mammano (1996) as

where $y$ is equivalent to $aη(x,t)$ in Eq. (2); the parameter values of the function are $y1=0.1139$⁠, $y2=0.3736$⁠, $c1=0.7293,$$c2=1.4974$⁠, $b=0.309 91$⁠. This function simulates the mechanoelectrical transducer which is the dominant nonlinearity in the system (Santos-Sacchi, 1993).

The hair-bundle displacement, $η(x,t)$⁠, caused by the BM motion, is defined as the shearing displacement of the RL relative to the TM, and is calculated using the differential equation

where $γTM(x)$ is the damping caused by the viscosity of the subtectorial space and viscoelasticity of the TM, and $ωTM(x)$ is the TM resonance frequency. The sigmoidal shape of the nonlinear function, $S(y)$⁠, limits cochlear amplification at higher intensities and causes the input/output (I/O) functions of BM motion to be compressive.

This section provides an analytical approximation to the solution of the model equations at the radial frequency of the lower cubic DPOAE. The solution, Eq. (16), determines how adjacent parts of the distributed DPOAE primary sources contribute to the resulting DPOAE. The equations presented below are based on the equations in Vetešník and Gummer (2012) and, therefore, the reader is referred to that publication for a comprehensive presentation of the details.

The model is analyzed in the frequency domain. Since the stimulus is composed of two pure tones and the model is nonlinear, the solution becomes tractable if we assume that the input is composed of tones which are multiples of some fundamental angular frequency, $ω0$⁠; that is, $ω1=n1ω0$⁠, $ω2=n2ω0$⁠. The steady-state solution of Eqs. (1) and (A1), $ξstat(x,t)$⁠, can then be approximated by a truncated Fourier series (Nobili and Mammano, 1996)

because the amplitudes of the Fourier components, $ξ̂n(x)$⁠, are negligible for $|n|>N$⁠, where $N$ is limited by the frequency range of the cochlear model. (The circumflex designates a phasor.) Each term of the Fourier series satisfies the integral equation

where $ĤnTM(x)$ is a dimensionless transfer function relating $η̂n(x)$ to $ξ̂n(x)$⁠; namely,

and $ÛnNL(x)$ is the nth Fourier component of the nonlinear undamping force defined as

Finally, $Ω̂n$ is the nth Fourier component of the stapes displacement, $σ(t)$⁠, which is calculated from Eq. (A1) as

The approximate analytical solution at the radial frequency, $ωDP=2n1ω0−n2ω0$⁠, (⁠$nDP=2n1−n2$⁠), was derived in Vetešník and Gummer (2012). It assumes (1) the DP generation site is restricted to a narrow region which is well-separated from the BM characteristic frequency (CF) site of $ωDP$ and (2) the longitudinal coupling due to shearing viscosity and the short-range term of the Green's function do not qualitatively affect the steady-state solution. The effects of both couplings can be then included into local properties of the BM oscillators; i.e., into the BM transfer function $ĤnDP(x)$ defined by Eq. (22) in Vetešník and Gummer (2012). This transfer function relates the local hydrodynamic force per unit length across the cochlear partition to BM displacement at the DP frequency and includes both passive and active components. Under these assumptions, Eq. (7) is equivalent to a second-order differential equation with boundary conditions at $x=0$ and $x=L$⁠, and allows derivation of the approximate analytical solution.

This solution is the sum of two terms (Vetešník and Gummer, 2012); namely,

where the first term represents the direct response of the BM to $ÛnDPNL(x)$ and the second term represents responses due to fluid coupling. $QnDP$ is a constant defined as (Vetešník and Gummer, 2012)

where $ρ$ is the fluid density given together with the remaining parameters in Table I, and $K̂nDP(x,x¯)$ is the Green's function of the boundary-value problem, and is given by (Vetešník and Gummer, 2012)

where $Ŵ$ is the Wronskian determinant of the homogeneous differential equation and is a constant for this application (Vetešník and Gummer, 2012). The functions $ξ̂nDP(1)(x)$ and $ξ̂nDP(2)(x)$ are linearly independent homogeneous solutions (base functions) which satisfy the boundary conditions at $x=0$ and $x=L$⁠, respectively. $ξ̂nDP(1)(x)$ can be expressed as a superposition of two waves with opposite phase slope (a backward wave and a forward wave due to the reflection at the stapes),

where $R̂$ represents the stapes reflection coefficient viewed from within the cochlea, $â$ is the complex constant determined by the boundary condition at $x=0$⁠, and $AnDPB,F(x)$ and $φnDPB,F(x)$ represent the amplitudes and phases of the backward $(B)$ and forward $(F)$ waves, respectively. $ξ̂nDP(2)(x)$ is approximated by a single wave traveling from the base to the apex of the cochlea,

where $b̂$ is the complex constant determined by the boundary condition at $x=L$⁠.

For studying the generation of the DPOAE nonlinear component, the solution of Eq. (11) at $x=0$ is the relevant BM displacement because it is this spatial value that drives the stapes in reverse. Provided the tonotopic place of $f2$ is located sufficiently distant from the basal end of the BM, the nonlinear force at the DP frequency at the basal end, $ÛnDPNL(0)$⁠, is negligibly small. Therefore, $ξ̂nDP(0)$ is given (approximately) by the second term in Eq. (11). Furthermore, the integration interval in Eq. (11) can be narrowed by assuming that $ÛnDPNL(x)$ due to interaction between the $f1$ and $f2$ traveling waves is negligible except for $x1<x<x2$⁠.¹ That is, the fluid pressure difference at the BM base, $P̂nDPP(0),$ which is proportional to the DPOAE pressure in the ear canal, is given by

In this interval, the nonlinear force, $ÛnDPNL(x),$ produces DPs which may contribute to the DPOAE signal. In Sec. III, we determine this interval by the method described in footnote ¹ and use Eq. (16) to analyze the dependence of DPOAE amplitude on stimulus intensities and on the spatial extent of the nonlinear force. We emphasize that the intensity dependence of the DPOAE signal does not depend on $ξ̂nDP(2)(x)$ because being a solution of a homogeneous equation, $ξ̂nDP(2)(x)$ is independent of a stimulus. Because the first three right-hand factors of Eq. (16), $ξ̂nDP(1)(0)(QnDP/Ŵ)$⁠, are independent of the stimulus amplitudes, we treat those terms as a constant for a given $nDP$ and focus the analysis only on the integral. The theoretical analysis and the experimental measurements use a frequency ratio of $f2/f1=1.2$⁠, the value commonly used in routine measurements in the laboratory and the clinic. This choice provided a relatively localized primary source up to stimulus levels of at least 50 dB SPL.

According to Vetešník and Gummer (2012), the nonlinear force $ÛnDPNL(x)$ can be written as

where $φn1F(x)$ and $φn2F(x)$ are the spatial phases of the steady-state solution for the stapedial input at the stimulus frequencies, $C(x)$ is a slowly changing function of $x$⁠, and $UANL(x)$ is the spatial amplitude of $ÛnDPNL(x)$⁠. If the generation region is restricted to a sufficiently narrow interval, then the phases of $ÛnDPNL(x)$ and $ξ̂nDP(2)(x)$ do not change significantly and the DPOAE signal depends on the constructive contribution from distributed DPOAE sources.

If sufficiently slow onset and offset of the DPOAE eliciting tone is used, we can assume instantaneous coupling between the DP generation sites and the DPOAE signal at the BM base. Numerical simulations suggest that for an elicitor duration longer than about 10 ms (Vencovský and Vetešník, 2018b), Eq. (16) can be extended to the time domain to approximate the time courses of the onset and offset responses of the DPOAE fluid pressure difference at the BM base; that is, $P̂nDPP(0)$ is now extended to $PnDPP(0,t)$ by the approximation

where

$x1(t)$ and $x2(t)$ specify the generation region of the primary sources, but in contrast to the foregoing steady-state analysis, this region now changes with time during the onset and offset of the DPOAE-eliciting stimulus tone, as does the nonlinear force, $ÛnDPNL(x,t)$⁠. To determine the generation region at a given time instant, we used the same approach as for the steady-state analysis as given in footnote ¹. Equation (9) calculates the nonlinear force in the frequency domain and assumes a stationary signal. However, the calculation can be also done numerically in the time domain. The method used to evaluate Eq. (18) is given in Sec. II C 2.

It has been shown in the study by Vencovský and Vetešník (2018b) that the suppression of BM responses, $ξ̂n1,2(x)$⁠, at one stimulus frequency during the onset of a pulsed tone causes complexities in the onset and offset of $PnDPP(0,t)$⁠. In the present study, we used Eq. (18) to investigate the temporal development of the complexities within the onset and offset regions of analytically simulated $PnDPP(0,t)$ signals. The results are presented in Sec. III A 2.

This section describes the methods used for the simulations. Subsections have the same names as corresponding subsections in Sec. III, to improve orientation and clarity. In the text, we use the terms “numerical simulations” and “analytical simulations.” By “numerical simulations” we mean simulations obtained by numerical solution of the cochlear box-model equations in the time domain [Eqs. (1), (5), and (6)] which, together with the middle-ear transmission model ( Appendix), yield the acoustical pressure of the nonlinear-distortion component of the DPOAE signal in the external auditory canal. The term “analytical simulations” denotes calculation of Eqs. (16) and (18), using $ξ̂nDP(2)(x)$ and $ÛnDPNL(x)$ derived from the numerical simulations, to yield analytical approximations of the nonlinear-distortion component of the DPOAE fluid pressure across the organ of Corti at the BM base.

The cochlear box model was discretized in the spatial dimension with 800 equal-length segments. We repeated some of the simulations with a larger number of segments (1200) and did not find any significant differences. The time-domain solution of the model was performed with a matlab implementation of the Dormand–Prince method, a member of the Runge-Kutta family of ODE solvers. The sampling frequency was 600 kHz. Source codes of the model will be provided on request. The box model was not calibrated to yield absolute values of BM displacement and DPOAE pressure expected for the human cochlea. In addition, only those parts of Eqs. (16) and (18) are calculated which depend on the examined stimulus parameter (stimulus levels). Therefore, arbitrary units (a.u.) are used to represent simulated responses.

Because of the assumed near-scale invariance of the box model² (Vetešník and Nobili, 2006), we mainly focused our analysis on one setting of the stimulus frequencies: $f1=2$ kHz, $f2=2.4$ kHz $(f2/f1=1.2)$⁠, which gives the cubic distortion product at $fDP=2f1−f2=1.6$ kHz. We kept the level of one of the stimulus tones fixed (50 dB SPL) and changed the other tone level from 30 to 70 dB SPL with a step of 10 dB SPL. In addition, we simulated DPOAEs for the same stimulus tones as in the experiment described below. Those simulations were conducted for two $f2$ values: 5.5 and 5.9 kHz, and a frequency ratio $f2/f1$ of 1.2. Stimulus levels were changed in the same way as described above for $f1=2$ kHz and $f2=2.4$ kHz stimuli.

In the analytical solution, only the integral over the product of the base wave $ξ̂nDP(2)(x)$ and the nonlinear force $ÛnDPNL(x)$ in Eq. (16) depends on the stimulus levels. Therefore, an estimate of the variation of the steady-state DPOAE signal as a function of $L1$ or $L2$ can be obtained by calculating only this integral. The base wave $ξ̂nDP(2)(x)$ was calculated from the model equations solved in the frequency domain [Eqs. (1), (5), and (6)]. The nonlinear force $ÛnDPNL(x)$ was calculated numerically in the time domain by first subtracting the linear approximation of $U(x,t)$ from $U(x,t)$⁠; that is, referring to Eq. (3), by calculating the difference $u(x)S[aη(x,t)]−u(x)aη(x,t)$⁠. To extract the cubic distortion component, this difference function was filtered with a zero-phase bandpass filter centered at $fDP$⁠, which was designed as a window-based finite impulse response (FIR) filter of order 800 with filter coefficients computed using a Hamming window. The 6-dB bandwidth was chosen to be dependent on $f2$ and calculated according to

where the unit of frequency in this equation is Hz. Before filtering, the time-domain signal was resampled to the sampling frequency used in the experiments; namely, 102.4 kHz.

In order to better understand the dependence of the steady-state DPOAE signal on stimulus level for the box model, which contains a longitudinal array of coupled nonlinearities, we investigated the two-tone steady-state response inherent to the sigmoidal nonlinearity used in the box model [Eq. (4)]. A purely analytical approach for quantifying the two-tone response is complicated because of the saturating character of the sigmoidal function. The Taylor series expansion of the sigmoidal function converges too slowly in the saturation regions. Therefore, the two-tone response was simulated in matlab. A systematic analysis of the behavior of the nonlinearity can be found in Duifhuis (1989) and Lukashkin and Russell (1998, 1999, 2001).

The stimulus-level dependence of the DP amplitude inherent to the (isolated) sigmoidal function will be shown to be qualitatively similar to that for a single nonlinear BM oscillator generating a DP in the box model (Sec. III A 1). However, the similarity cannot be anything but qualitative because not being embedded in the cochlear box model, the output and, therefore, the input of the isolated nonlinearity is not coupled to neighbouring nonlinearities. Therefore, DPOAEs generated by a single oscillator within the cochlea are, in this publication, determined analytically by using only one model-segment in the cumulative integral given in Eq. (16). For this purpose, we used the model segment at the place of maximum nonlinear force. This approach allows us to determine how the distributed nature of the primary source contributes to the observed level dependence of the DPOAE in the box model.

DPOAE signals were simulated by numerical solution of the model equations in the time domain. In order to remove the responses at the stimulus frequencies, the stimulus tones were presented with different phases (multiples of 90° and 180°) and the signals averaged in the time domain (Whitehead et al., 1996). For the analytical simulations, the nonlinear force $ÛnDPNL(x,t)$ and the base wave $ξ̂nDP(2)(x)$ were calculated in the same way as described in Sec. II C 1 for analytical simulation of the steady-state DPOAE response. Then, the DPOAE signals for fluid pressure at the basal end of the BM were derived analytically from Eq. (18). In all simulations, onset/offset of the DPOAE elicitors were shaped with 10 ms raised-cosine ramps.

To investigate how the coherent-reflection source may affect the DPOAE steady-state signal for the stimulus parameters used in this publication, random inhomogeneities (roughness) were introduced into the model used for numerical simulations. For this purpose, the function $u(x)$ given by Eq. (3), which determines the undamping force in the model, was modified to $ũ(x)$ by adding Gaussian noise,

where $ϵ=0.05$ and $N(0,1)$ denotes a Gaussian distribution with zero mean and variance of unity. In order not to affect the nonlinear-distortion component generated near the tonotopic place of the $f2$ traveling wave, we added the inhomogeneities only to the model sections between $x=1.50$ cm and $1.67$ cm (CF of 1.8 kHz and 1.4 kHz) for $f2=2.4$ kHz, and $x=0.90$ and $1.10$ cm (CF of 4.0 kHz and 3.0 kHz) for $f2=5.5$ kHz. The onset and offset of the inhomogeneities were shaped with a Tukey window defined as

where $r=0.3$ and $z$ ranges between 0 and 1. The window covered all the model sections into which the inhomogeneities were added; e.g., for $f2=2.4$ kHz the model section for $x=1.50$ cm corresponds to $z=0$⁠, and the model section for $x=1.67$ cm corresponds to $z=1$⁠; the positions of the remaining model sections were linearly scaled between 0 and 1.

DPOAEs were measured in two normal-hearing (pure-tone hearing threshold <20 dB hearing level) men of age 29 and 46 years, denoted as s01 and s02, respectively. The experiments were approved by the Ethics Committee of the University of Tübingen in accordance with the Declaration of Helsinki for human experiments. All measurements were performed in a sound-proof chamber (Industrial Acoustics Company, Niederkrüchten, Germany) using an ER-10C DPOAE probe-microphone system (Etymotic Research, Elk Grove Village, IL) connected to a PC-based setup (Zelle et al., 2015, 2017b).

The DPOAE measurements were conducted for two different $f2$ values: 5.5 and 5.9 kHz. These frequencies were chosen so as to find a frequency region which was free of spontaneous otoacoustic emissions (SOAEs) and which had a relatively small secondary-source amplitude (≤0.27 of the primary-source amplitude) at low stimulus levels (≤45 dB SPL). SOAEs can interfere with small-amplitude DPOAEs (van Dijk and Wit, 1990; Zelle et al., 2017b). In both subjects, the secondary source component was small at high frequencies: in s01 at 5.5 kHz and in s02 at 5.9 kHz. DPOAE measurements were conducted for a frequency ratio $f2/f1=1.2$⁠.

In order to reduce the amplitude responses at $f1$ and $f2$⁠, the stimulus tones were presented with different phases (multiples of 90° and 180°) and the signals averaged in the time domain (Whitehead et al., 1996). Two different acquisition schemes were used. In the first scheme, the focus was on the onset, offset, and steady-state parts of the DPOAE signal. The longer of the two stimulus tones was turned on shortly after the onset of the acquisition block; the other stimulus tone was turned on after a further time delay and turned off before the offset of the first stimulus tone. The shorter stimulus tone is called the “elicitor” since its presentation evokes DPOAEs. Both $f1$ and $f2$ tones were tested as elicitors, similarly to Zelle et al. (2018). Each acquisition block had duration of 90 ms. The longer (80-ms duration) stimulus tone was turned on 5 ms after the onset of the acquisition block; the onset and offset envelopes were shaped with 4-ms raised-cosine ramps. The elicitor tone was turned on 20 ms after the onset of the acquisition block and had duration of 30 ms. The onset and offset of the elicitor envelopes were shaped with 10-ms raised-cosine ramps. The measurements were performed for two level paradigms of the stimulus tones: either for $L1$ kept fixed at 50 dB SPL and $L2$ varied from 30 to 70 dB SPL with a 10-dB SPL step, or $L2$ kept fixed at 50 dB SPL and $L1$ varied between 30 and 70 dB SPL with a 10-dB SPL step. The second acquisition scheme used shorter elicitor pulses (duration of 9 ms instead of 30 ms) to allow separation and acquisition of both the nonlinear-distortion and coherent-reflection components (Zelle et al., 2017a). The $f2$ tone was the elicitor and the relative levels of the stimuli were optimized for each $f2$ and $L2$⁠, by choosing $L1$ such that the generated DPOAE amplitude would be as large as possible (Zelle et al., 2015). The duration of the acquisition block was 60 ms. The $f1$ tone was turned on 5 ms after the onset of the block and had duration of 50 ms. The onset and offset envelopes of the $f1$ tone were shaped with 4-ms raised-cosine ramps. The elicitor was turned on 20 ms after the onset of the acquisition block; its onset and offset envelopes were shaped with 2.5-ms raised-cosine ramps. For both acquisition schemes, the frequency ratio was $f2/f1=1.2$⁠.

The same averaging procedure was used for both acquisition schemes. Averaging was performed for a maximum of 100 acquisition blocks; acquisition was terminated if the signal-to-noise ratio (SNR) between the DPOAE amplitude and the noise, estimated in the time domain, was at least 10 dB; see Zelle et al. (2017b) for details. The average included only those acquisition blocks which enhanced the SNR.

The short-pulse DPOAEs, acquired with the second scheme described in Sec. II D 1, were used to extract information about the coherent-reflection source. Due to the different latencies of the two sources, the onset of the DPOAE signal is dominated by the nonlinear-distortion source and the offset by the coherent-reflection source (Zelle et al., 2013). Therefore, the DPOAE signal can be described as a superposition of the nonlinear-distortion and coherent-reflection components, each of which is estimated in amplitude and phase by fitting a pulse-basis function (PBF) to the recorded DPOAE signal (Zelle et al., 2013, 2017a). Details of the PBF method are given in Zelle et al. (2013). The PBF method was applied to the DPOAE signals averaged in the time domain (sampling frequency 102.4 kHz) and band-pass filtered using a zero-phase finite impulse response (FIR) filter with an order of 800 and filter coefficients computed using a Hamming window.

Section III A presents a theoretical analysis of the steady-state (Sec. III A 1) and time course (Sec. III A 2) of the DPOAE signal. To this end, the DPOAE generated by a single oscillator is compared with the DPOAE generated by the box-model containing an array of oscillators, operating in a smooth cochlea (without mechanical inhomogeneities) and producing mutually interfering wavelets. Section III B compares the simulations (without mechanical inhomogeneities) with experiments (Sec. III B 1) and then examines a possible contribution from the secondary source to the DPOAE onset and offset responses (Sec. III B 2). It is concluded that, for the stimulus intensities and frequencies used in the present experiments, the secondary source had negligible effect on the DPOAE signal and, therefore, that complexities in the onset and offset responses of the DPOAE were due to two-tone suppression between stimuli.

This section presents model simulations of DPOAEs, focusing on their steady-state properties and the effect of changing one of the stimulus levels. It combines numerical simulations in the time-domain with analytical simulations in order to reveal how stimulus levels affect the DPOAE primary source.

Before presenting the DPOAE simulations, we demonstrate the consistency of model BM responses with experimental data reported in the literature for single-frequency and impulse stimulation. To this end, the I/O function of the BM displacement amplitude, obtained from numerical solution of the box model, at the BM position tuned to 2.4 kHz (upper panel) and 5.5 kHz (lower panel) at 0 dB SPL is depicted in Fig. 1. At these BM positions, the gain of the cochlear amplifier³ was 54 and 58 dB, respectively. The simulations at 5.5 kHz were performed because this frequency was also used for experimental DPOAE data (Sec. III B 1).

Figure 2(A) shows displacement amplitude and phase of the iso-intensity model responses at the BM position tuned to 2.4 kHz (left panels) and 5.5 kHz (right panels) at 0 dB SPL. Notice the broadening of the amplitude responses and the basal shift of their maxima with increasing stimulus level. The phase responses are relatively independent of level. The level dependence of the iso-intensity responses is in qualitative agreement with experimental data (Robles and Ruggero, 2001). Figure 2(B) depicts impulse responses obtained at the same BM position as for the iso-intensity frequency responses. The zero crossing points of the impulse responses are independent of level, agreeing with experimental data (Recio and Rhode, 2000; Verhulst et al., 2012).

Taken together, this agreement between simulation and experiment provides confidence in the cochlear model.

To check the validity of the analytical approximation for the DPOAE amplitude given by Eq. (16), steady-state DPOAE amplitudes were compared for the numerical and analytical simulations of the box model. Figure 3 compares the intensity dependence of the steady-state DPOAE amplitude calculated by both methods for $f2=2.4$ kHz with $f2/f1=1.2$⁠. The analytical and numerical amplitudes well-nigh superimpose, thus supporting the use of Eq. (16) in our analysis of suppression phenomena. For each intensity condition, either constant $L1=50$ dB SPL or constant $L2=50$ dB SPL, the maximum DPOAE amplitude is reached for the optimal combinations of $L1$ and $L2$ given in Kummer et al. (1998) and Zelle et al. (2015). At the optimal levels, according to definition (Kummer et al., 1998), the amplitudes of both traveling waves at the $f2$ tonotopic place are approximately equal.

Now, using the analytical solution [Eq. (16)], we describe the spatial distribution of the DPOAE primary source, which is given by the nonlinear undamping force $ÛnDPNL(x)$⁠.

Figures 4(A) and 4(D) depict the steady-state amplitude and phase of the traveling-wave displacement of the BM when stimulated, respectively, with fixed $L1=50$ dB SPL or $L2=50$ dB SPL; the variable stimulus level ranged from 30 to 70 dB SPL with a step of 10 dB SPL. The traveling waves are depicted for the stimuli and the DP component. Figures 4(B), 4(E) show the corresponding nonlinear forces, $ÛnDPNL(x)$⁠, which are the source of the nonlinear-distortion component of DPOAEs. Figures 4(C) and 4(F) depict the corresponding cumulative integrals calculated from these forces as an integral of the product of the nonlinear force and the DP base wave $ξ̂nDP(2)(x)$ [Eq. (16)]. Evaluated on the interval $[0, x]$⁠, the result of the integration gives a value proportional to the DPOAE signal for primary source generators distributed along the BM between the base (⁠$x=0$⁠) and $x.$ In other words, the integral sums all the primary sources (wavelets) from the base to the place $x$⁠, with each generator having a unique phase. Notice that the amplitude of the cumulative integral does not change for $x$ greater than some apical value. This means that there are no sources of DPOAE apically from that place. The region in which the cumulative integral is drawn with a thickened line designates the region which should mainly contribute to the generation of the DP.¹ That region corresponds to the integration interval $[x1, x2]$ in Eq. (16); the procedure for estimating this interval is given in footnote ¹. The ringing in some of the cumulative integrals [Figs. 4(C) and 4(F)] around the $fDP$ tonotopic place, apical to $x2$⁠, is associated with the rapid roll-off of the phase of the nonlinear force in that region, which leads to interference between wavelets. Since the ringing is approximately symmetrical about the asymptotically constant value of the cumulative integral at large $x$⁠, the wavelets in that ringing region mutually cancel and, therefore, do not contribute to the DPOAE signal.