## Abstract

This paper reviews current understanding and presents new results on some of the nonlinear processes that underlie the function of the mammalian cochlea. These processes occur within mechano-sensory hair cells that form part of the organ of Corti. After a general overview of cochlear physiology, mathematical modelling results are presented in three parts. First, the dynamic interplay between ion channels within the sensory inner hair cells is used to explain some new electrophysiological recordings from early development. Next, the state of the art is reviewed in modelling the electro-motility present within the outer hair cells (OHCs), including the current debate concerning the role of cell body motility versus active hair bundle dynamics. A simplified model is introduced that combines both effects in order to explain observed amplification and compression in experiments. Finally, new modelling evidence is presented that structural longitudinal coupling between OHCs may be necessary in order to capture all features of the observed mechanical responses.

## 1. Cochlear physiology

Hearing is a complex biological process, which relies on many highly developed subsystems. The main concern of this paper is the mammalian cochlea, which houses the first active and nonlinear components of the auditory system. The cochlea's most basic function is to convert sound-evoked motion into an ensemble of neural signals which underlies an organism's ability to hear.

A schematic of the uncoiled mammalian cochlea can be seen in figure 1*a*,*b*. It is a fluid-filled tube that consists of three chambers: the *scala tympani*, *scala media* and *scala vestibuli*. At the basal end, there are two openings: the oval window, which is connected to the middle ear, and the round window, which is covered by an elastic membrane. The scala tympani and scala media are separated by the *organ of Corti*, which is built on the flexible *basilar membrane*. The scala media and scala vestibuli are separated only by the thin *Reissner's membrane* and are usually modelled as a single chamber.

A typical transverse cross section through the organ of Corti is depicted schematically in figure 1*b*. Below the basilar membrane there is the scala tympani, and above the organ of Corti is the scala media. The organ of Corti can also be seen from above in figure 1*c* with the tectorial membrane removed. The three rows of outer hair cells (OHCs) are connected to the basilar membrane by the *Deiters' cells*, and at the other end to the *reticular lamina*. The hair bundles of the OHCs are stimulated by the shearing motion of the *tectorial membrane* which slides above the reticular lamina as the basilar membrane is deflected by pressure differences between the fluid filled chambers. Sensing of sound occurs through the inner hair cells (IHCs), which are positioned on the opposite side of the *tunnel of Corti* (the triangular-shaped space between the IHC and OHCs depicted in figure 1*b*). The IHCs signal the characteristics of a sound to the brain via synaptic contacts with individual fibres of the auditory nerve.

Attempts to understand the function of the human hearing organ can be traced to the work of Aristotle, who supposed that acoustic resonance occurs in an air-filled cavity located behind the eardrum (see Gitter [1]). The discovery of the anatomy of the basilar membrane by Hensen [2] inspired Helmholtz [3] to formulate the theory of tonotopy, also known as the place theory of hearing. Specifically, every transverse cross section of the cochlea is most sensitive to a particular frequency, the so-called characteristic frequency (CF) of that location. The tonotopy map is such that low frequencies excite the basilar membrane at apical positions along the cochlear tube, while high frequencies stimulate the basal end. Helmholtz also recognized that merely having a resonant basilar membrane does not explain the observed sharp frequency discrimination. Therefore, he hypothesized that there are additional oscillators connecting the basilar membrane and the sensory nerves.

The experiments of Békésy [4] indicated that the fluid–structure interaction of the cochlear tube and membrane supports a travelling wave. This wave reaches a maximum amplitude and rapidly attenuates at a frequency-dependent distance from the stapes, thus explaining the tonotopic map (figure 1*d*,*e*). Such behaviour can be reproduced theoretically via a dynamic interplay between the fluid in the cochlear duct and the basilar membrane, with stiffness that varies with longitudinal position (e.g. [5–7]).

Békésy's experiments were conducted using high-amplitude stimulation of post-mortem cochleae, and revealed only relatively broadly tuned mechanics. Later experiments by Johnstone & Boyle [8], Rhode [9] and Sellick *et al.* [10] found much sharper tuning and active mechanical responses in living animals. Another piece of evidence pointing to an active process in the cochlea's mechanics was the discovery by Kemp [11] that sound is emitted from a healthy ear after a short stimulus, a phenomenon known as *otoacoustic emission*.

As experimental technology developed, it became possible to characterize the nonlinear vibrations of the active basilar membrane more fully (see Robles & Ruggero [12] for a review). Exemplar basilar membrane responses obtained by Rhode [13] are shown in figure 4*b*,*c* for the 9.1 kHz CF position of the chinchilla cochlea when subjected to a range of different pure tone frequencies and amplitudes. These data represent the experimental results by which we shall benchmark modelling results in §§3 and 4 below.

The active component of the cochlea's mechanics is widely accepted to be located within the OHCs of the organ of Corti. These particular cells have two types of active processes (e.g. [14,15]). The hair bundle on top of the cell is motile, and can hence exert force against the vibrating structures of the cochlea. For example, the experimental results of Kennedy *et al.* [16] show that direct mechanical stimulation of the hair bundles can result in a period of negative stiffness. On the other hand, it was observed by Brownell *et al.* [17] that if OHCs are depolarized, they shorten. Later Zheng *et al.* [18] discovered that a piezoelectric protein called *prestin*, which is embedded within the lateral membrane of the OHC body is responsible for this effect. Prestin can expand and contract in response to changes in the cell potential, and makes the cell body change its shape. Both bundle motility and somatic motility can be stimulated by cochlear vibrations through the mechanically sensitive ion channels (MET channels) found at the tips of the cilia that make up the hair bundle.

There is currently a debate in the literature on the relative importance of these two forms of motility in the active amplification process within mammalian hearing (e.g. [19]). While somatic motility can exert larger forces, and may hence seem more potent from a mechanical point of view, amplification is also known to occur in lower organisms, such as birds and reptiles, whose hair cells do not contain prestin but have rather more pronounced stereocilia. On the other hand, transmission to the basilar membrane motion is more straightforward with a somatic motility hypothesis, and the experimental results of Lagarde *et al.* [20], in mutant mice, suggest that somatic motility alone is sufficient to match wild-type behaviour. It has also been suggested that somatic motility may be too slow. However, Frank *et al.* [21] have recently shown that somatic motility can operate up to very high frequencies, for example much higher than the highest that are perceived by humans. The main limitation on the speed of somatic motility seems to be the time needed to charge up the capacitance of the OHC's cell body—a process which involves low-pass filtering with a corner frequency of about 1 kHz. However, taking into account ion current adaptation in the ion channels of the hair bundle, even this filter can be compensated for (see Ó Maoiléidigh & Jülicher [22]). Indeed, the model we present in §3 below suggests that, from a modelling point of view, the observed behaviour can be explained by a combination of current adaptation and somatic motility.

There have been many other explanations of the source of the active process in hearing (see Hudspeth *et al.* [23], Ashmore *et al.* [24] and references therein). Of particular note is the work of Zweig [25], who added a complex, but linear, feedback term to a simple travelling-wave model of the basilar membrane; the results were found to reproduce the active response of the cochlea well. See also the works of Steele & Lim [26] and Geisler [27], who found similar results using a spatial feed-forward term within the travelling-wave equation. A discussion and extension of these results forms the subject of §4 below.

The measurements in figures 1 and 4 all concern the motion of the basilar membrane, rather than of hearing itself. The transduction into neural signals occurs within the IHCs. It is widely held that increased basilar membrane motion is transferred into increased stimulation of the IHCs through the fluid flow in the subtectorial space. This flow causes displacement of the cilia of the IHCs, which then causes ion channels to open, which in turn causes depolarization of the IHCs, allowing the connected nerve fibres to fire.

This paper uses the theory and techniques of nonlinear dynamical systems. Bifurcation analysis refers to the study of qualitative changes to the state of a system as a parameter is varied. This question is the basis of bifurcation theory and it underlies a qualitative understanding of many biological processes and transitions, such as the onset of oscillation, switching, multi-stability, emergence and localization. Nonlinear systems can have multiple stable steady states or attractors to which solutions with nearby initial conditions converge. They can also support periodic oscillations. Moreover, as a parameter is varied, a new attractor can emerge out of ‘thin air’, typically when an unstable state gains local stability. For this reason usual numerical methods and computer simulation for ordinary differential equations (ODEs) can be highly unreliable in understanding the true dynamics of the system, if the method is based on simulation from fixed initial conditions. Bifurcation theory can be especially useful in this context as it enables the tracing of paths of unstable states as parameters vary, and hence determine precisely the transitions (or ‘bifurcation points’) at which qualitatively distinct stable behaviour emerges. There are dedicated numerical implementations of bifurcation theory using path-following, or numerical continuation [28]. A codimension-one bifurcation can often be represented in a *bifurcation diagram* that depicts a measure, or norm, of a system state against a single parameter (e.g. figure 2*c*). Codimension-one bifurcations can also be used to divide regions in a parameter plane in which qualitatively distinct bifurcations can occur.

The rest of this paper is organized as follows. First, we present recent progress on mathematical modelling and experimental investigation of hair cell dynamics. Section 2 introduces some recent progress by Tsaneva-Atanasova and Kennedy on mathematical models of the function of developing IHCs. The model is validated using electrophysiological recordings reported here for the first time and simulation results are explained with the aid of numerical bifurcation analysis. Section 3 reviews local nonlinear models of active OHC behaviour, within a fixed transverse section of the organ of Corti. We describe a simple model, first proposed in Szalai *et al.* [29], which is based on a comprehensive model by Ó Maoiléidigh & Jülicher [22] with the emphasis on the role of somatic motility, transduction current adaptation and nonlinearity that is consistent with the experimental results in Kennedy *et al.* [16]. We also argue that it is not necessary for local models of the OHC dynamics to take the form of the popular ‘Hopf oscillator’ [30] to provide results in agreement with experimental data. We then go on (in §4) to discuss how local models (specifically, of the organ of Corti) can be incorporated into a longitudinal cochlear model, which includes fluid coupling as well as coupling through a feed-forward mechanism. Preliminary results are presented that show how a longitudinally coupled version of the local OHC model can qualitatively reproduce experimentally observed vibration characteristics of the basilar membrane. Finally, §5 draws conclusions and suggests avenues for further work.

## 2. Modelling the development of inner hair cell function

While OHCs serve to mechanically pre-amplify the mechanical vibrations (i.e. sound) that reach the inner ear, IHCs are the proper receptor cells of hearing and are connected to the afferent nerves. Sound transduction by IHCs generates a receptor potential which triggers the coordinated release of neurotransmitters onto the ensemble of auditory nerve fibres. The membrane filter properties that define the input–output function of IHCs are dependent on membrane conductance and capacitance [31]. In normal hearing, deflection of the cilia of an IHC causes the opening of mechano-electrical transducer ion (MET) channels. This allows a flow of positively charged ions, which depolarizes the cell's plasma membrane, and causes the nerve connections at the basal end of the cell to respond. These signals are then sent through nerve pathways to the brain, which makes sense of them as sounds. This process must be performed at exceptional speed in order to encode the complex characteristics of sound such as speech and music. It is therefore important that the excitability of the membrane is carefully controlled to be ready to respond rapidly.

The main control of the membrane potential in an IHC occurs through the action of calcium on the excitability of the hair cells. Calcium is involved in both the mechanotransduction at the cell's apex and neurotransmitter release at the base (e.g. [32]). Thus, regulation of intracellular Ca^{2+} is important for local dynamic fine-tuning of the IHC membrane filter [31]. A delicate calcium balance must be maintained in order for hearing to occur, but more insight is required into how this is achieved. An interesting approach to such a problem is to study how calcium also controls the developmental processes within IHCs. Before we are able to hear, the ear must develop the unique sensory properties of the hair cells and form the intricate connections to the brain. It is believed that a key factor in this process is an increase in calcium ions (Ca^{2+}), that is thought to cause neurotransmitters to be released, which in turn helps guide neural connections to form at the basal end of the cell.

During development, IHCs do not respond to sound but instead are capable of firing spontaneous and evoked Ca^{2+}-based spikes that disappear close to the onset of hearing (e.g. [33]). These spikes generate sufficient Ca^{2+} influx to trigger neurotransmitter release (e.g. [34]), and activate small conductance Ca^{2+}-activated currents (so-called SK2 currents) that are known to be expressed in the cochlea (e.g. [35]). Under basal conditions, the membrane potential of these cells is constant near its equilibrium value and the intracellular calcium concentration, [Ca^{2+}]_{i}, is low, so the K_{Ca} channels are closed. The negative membrane potential keeps the voltage-dependent Ca^{2+} channels and [Ca^{2+}]_{i}, low, and no significant neurotransmitter release takes place. Upon stimulation, which generally leads to depolarization, the membrane potential increases, opening voltage-gated Ca^{2+} channels. This results in Ca^{2+} influx into the cell cytosol, which stimulates the release of neurotransmitters (e.g. [36]). The increase of [Ca^{2+}]_{i}, opens K_{Ca} channels, generating outflow of K^{+} ions, which repolarizes the membrane potential. This repolarization leads in turn to closure of voltage-gated Ca^{2+} channels and subsequent decrease in the levels of intracellular calcium concentration. In immature IHCs, the combined sequence of events described above triggers repetitive rises (oscillations) in the intracellular calcium concentration that are accompanied by action potential firing in the form of single spikes, as well as plateau-bursting electrical activity (see new data in the right-hand panel of figure 2*b*, more details of which will be presented elsewhere). Plateau-bursting electrical activity is characterized by periodic switches between an active (depolarized) phase accompanied by increase in [Ca^{2+}]_{i}, and a silent (repolarization) phase during which [Ca^{2+}]_{i}, decreases.

We have developed a simple preliminary model based on realistic interaction between ion channels to describe how Ca^{2+} influx and buffering generate a [Ca^{2+}]_{i}, signal within developing IHCs. Our model is of Hodgkin & Huxley [37] type, and is based on the voltage-clamp experimental data obtained from immature IHCs (e.g. [32,33,38]; further details, including the parameter fitting technique, will be presented elsewhere). The model equations take the form of four coupled ODEs for the membrane potential *V*_{m}, the concentration of intracellular calcium *c*=[Ca^{2+}]_{i}, an activation variable *n* and an inactivation variable *h* associated with the voltage-dependent K^{+} channels:
2.1
2.2
2.3
and
2.4
Here, *C*_{m} is the membrane capacitance, *τ*_{n}(*V*_{m}) is the activation time constant for the voltage dependent K^{+} channels, *τ*_{h} is the inactivation time constant for the voltage-dependent K^{+} channels, and are the steady-state functions for *n* and *h*, respectively. Equation (2.1) describes the overall dynamics in which the rate of change of membrane potential *V*_{m} is balanced by the various ionic currents. These currents include the voltage-gated calcium current *I*_{Ca}, the voltage-gated K^{+} current *I*_{K} and the Ca^{2+}-sensitive potassium current *I*_{KCa}. Specific forms of the currents included in the model are
and the steady-state activation–inactivation functions are
where is the steady-state Ca^{2+}-dependent inactivation of the voltage-gated calcium current [39] and is the steady-state Ca^{2+}-dependent activation of the K_{Ca} channels.

The plateau spikes during a burst are caused by the interaction between the voltage-dependent calcium and potassium currents *I*_{Ca} and *I*_{K}. There is also a slow process included, which is responsible for the shifts between the active and the silent phases. The slow modulation is provided by the intracellular calcium concentration *c*, which is a negative feedback variable. The rate of change of *c* is also given by a balance equation (2.4) for the fluxes involved in its dynamics. The Ca^{2+} flux across the plasma membrane is given by the difference between the Ca^{2+} influx, represented by the voltage-gated calcium current, and Ca^{2+} efflux through the plasma membrane Ca^{2+} pump. Here, *C*_{m}=10^{−5}×*A*_{cell} is the membrane capacitance, and *f*_{c} is the fraction of free to total cytosolic Ca^{2+}, *α*=10^{5}/(2×9.65×*A*_{cell}) is a factor that converts current to flux, is the area of the cell, *k*_{PMCA} is the maximum rate and *K*_{p} is the half-maximal activation of the plasma membrane Ca^{2+} ATPase pumps, *k*_{SERCA} is the endoplasmic reticulum (ER) membrane Ca^{2+} ATPase pump rate, *p*_{ER} is the permeability of the ER membrane and *c*_{ER} represents the free Ca^{2+} concentration in the ER. Since *c* represents the free Ca^{2+} concentration in the cytosol, the corresponding fluxes are multiplied by the fraction of free to total cytosolic calcium, *f*_{c}.

The values of all parameters used in the model simulations are as follows: *k*_{PMCA}=3.6 s^{−1}, *I*_{app}=0 pA, *k*_{SERCA}=1.2 s^{−1}, *K*_{p}=0.08 μM, *d*_{cell}=15 μm, *k*_{s}=1.25 nS, *V*_{K}=−60 mV, *g*_{Ca}=2.4 nS, *V*_{Ca}=60 mV, *g*_{K}=2.85 nS, *V*_{leak}=−20 mV, *g*_{leak}=0.12 nS, *V*_{mL}=−26.7 mV, *V*_{n}=−16 mV, *s*_{m}=11.5 mV , *s*_{n}=10 mV, *V*_{h1}=−60.5 mV, *V*_{h2}=−17.8 mV, *s*_{h1}=6.8 mV , *s*_{h2}=7.1 mV, *k*_{q}=0.6 μM, *τ*_{h}=0.55 s^{−1}, *c*_{ER}=500 μM and *f*_{c}=0.01.

It has been established experimentally that the activation of K_{Ca} channels controls the magnitude and duration of IHC response to stimulation [38]. Figure 2*a* presents a validation of our IHC model, illustrated via simulations of a gradual reduction in the amount of *I*_{KCa}, modelled by a time-dependent exponential decay in the maximum channel conductance *g*_{KCa}. These simulations accurately reproduce previously published experimental data of continuous recordings of voltage responses in the presence of K_{Ca} channel blocker (e.g. fig. 3C in Marcotti *et al.* [38]). As K_{Ca} channel activity is decreased in the model, spike repolarization gradually slows and eventually begins to fail leading to prolonged depolarizations. Moreover, the model predicts that, in response to spikes, signals of intracellular calcium *c* initially rise sharply and are brief. However, as repolarization begins to fail the resting *c*-value rises and the *c*-increases are less sharp, rising throughout depolarization. Interestingly, our experimental recordings of spikes (top traces in the right-hand panel of figure 2*a*), and their corresponding near membrane *c*-signals (bottom), confirm that during prolonged action potentials (bursts) *c*-signals appear larger in magnitude and duration. This would presumably result in increased neurotransmitter release probability [36]. In the left-hand panel of figure 2*b*, we show typical IHC model solutions for a fixed value of *g*_{KCa}=4 nS. These results show that our model supports periodic solutions of mixed type, which are composed of a single spike and a burst with several small amplitude spikes. We have also observed such patterns of electrical activity in immature IHCs as shown in the right-hand panel of figure 2*b*.

In order to gain an insight into the effect of changing maximum K_{Ca} channel conductance, *g*_{KCa}, we have performed a bifurcation analysis of the full model system equations (2.1)–(2.4) using *g*_{KCa} as a bifurcation parameter. Note that *g*_{KCa} also indirectly affects the time scale at which *c* changes in the model (equation (2.4)). The corresponding bifurcation diagram of the full IHC model, presented in three-dimensional (*g*_{KCa} *c*, *V*_{m})-space is shown in figure 2*c*. The bifurcation diagram of our model shows the emergence of complex bursting from a primordial, large amplitude spiking solution that undergoes a cascade of bifurcations in the full system as the parameter *g*_{KCa} decreases. For *g*_{KCa}≈9.25 the branch of periodic orbits that corresponds to stable large amplitude spiking loses stability in a period-doubling (PD) bifurcation. As *g*_{KCa} decreases a new stable period-2 orbit emanates from PD and corresponds to a 2-large-spikes bursting solution. This period-2 bursting attractor loses stability in another PD bifurcation that gives rise to a complex period-4 bursting solution (not shown) that is stable for considerably smaller range of values of *g*_{KCa} and, in turn, leads to a classical PD cascade to chaos (not shown). However, the range of values for *g*_{KCa} for which this behaviour is found is rather narrow, and we will not focus on it further. Instead, we concentrate on the stable complex bursting behaviour of the IHC model that occurs over a wider range of *g*_{KCa} values.

Direct numerical simulations of equations (2.1)–(2.4) show that decreasing *g*_{KCa} results in complex bursting solutions with more small spikes. Thus, we continue such periodic orbits, namely with 1 large spike and 2, 3, 4, 5, 6 or 7 small spikes, in parameter *g*_{KCa}. We find that they lie on isolated solution branches, or isolas, as illustrated in figure 2*c*. As *g*_{KCa} decreases each of these isolas that we compute loses stability in PD bifurcations (not shown) that presumably give rise to PD solutions and generate further complex behaviour that is beyond the scope of this paper. Examples of complex bursting trajectories with 1 large spike and 2, 3, 4, 5, 6 or 7 small spikes are superimposed on the bifurcation diagram in figure 2*c*. Note also that as *g*_{KCa} decreases, the stable portions of the isolas become smaller, resulting in smaller regions in parameter space where bursts with large number of small spikes exist.

Although we have focused on the mathematical effects of varying the parameter *g*_{KCa}, which represents the maximum K_{Ca} channel conductance in the model, perhaps most significant from a physiological point of view is the fact that it indirectly affects the calcium component of the bursting orbits in the system (2.1)–(2.4). Smaller values of *g*_{KCa} decrease the activation of the K_{Ca} channels and in turn slow down the rise in Ca^{2+}. As shown here, when *g*_{KCa} is small the cells exhibit bursts with more spikes and, hence, longer depolarized plateaus. Since the trajectory spends more and more time in the direction in the region of phase space with large *c* (that is represented by the third dimension in figure 2*c*) as a result it generates more spikes within a burst. Consequently, the longer bursts are accompanied by larger Ca^{2+} signals.

These results on the dynamics of IHCs, while not describing actual hearing, are likely to be important in development since the levels of intracellular calcium are instrumental for the regulation of neurotransmitter release. More widely, by correctly fitting mathematical models to the more experimentally amenable developing IHCs, we can in principle change parameters to represent mature IHCs in order to understand how they operate within hearing. For example, the transition between various oscillatory regimes could be regarded as a form of plasticity of the intrinsic membrane properties of IHCs, and thus could have a profound effect on their function. The mathematical model described here, carefully validated with new experimental results, will play a key role in developing this understanding.

## 3. Micromechanical models of outer hair cell amplification

While the IHCs are believed to be primarily responsible for encoding acoustic stimuli, the primary role of the OHCs appears to be cochlear amplification [24]. As described in §1, OHCs provide two forms of motility: hair bundle motility [16] and somatic motility [17].

A schematic of an OHC can be seen in figure 3*c*. Like the IHC, the OHC has MET (ion) channels that are connected to adjacent cilia in the bundle through elastic fibres, i.e. the tip links. The tip links are presumed, as in non-mammalian hair cells, to be coupled to myosin motors that pre-tense the tip link, perhaps to the most sensitive state of the MET channels. The MET channels are not only mechanically sensitive, but also Ca^{2+} sensitive. This provides a fast adaptation mechanism to tune the operating point of the entire hair bundle [40]. The myosin motors are also Ca^{2+} sensitive: their force decreases with increasing Ca^{2+} concentration.

It has been suggested that the negative stiffness seen in mammalian OHCs, together with spontaneous oscillations of the hair bundle that has been observed in lower vertebrates, e.g. in bullfrog, mean that the nonlinearity of hair cells can be essentially modelled as that of a nonlinear oscillator undergoing a Hopf bifurcation [30]. See Hudspeth [15] for a detailed review, including appealing evidence that a forced Hopf bifurcation normal form can give rise to a compressive nonlinearity and a possible explanation for otoacoustic emissions and two-tone interference in hearing [41].

In order to understand the possible physiological basis for the Hopf bifurcation hypothesis, Tinevez *et al.* [42] adapted a model of the myosin motor mechanism constructed by Assad & Corey [43] that can account for gating currents observed on a voltage clamped cell, as well as the mechanical properties of the hair bundle in bullfrogs. The current through the ion channels is qualitatively similar to that illustrated in figure 3*b*. In particular, Tinevez *et al.* [42] adapted the model to describe mammalian hair bundle dynamics, based on the experimental data of Kennedy *et al.* [16]. It was found that this model can indeed undergo Hopf bifurcation.

Ó Maoiléidigh & Jülicher [22] incorporated this myosin motor dynamics in their unified ODE model of an OHC which also includes the charge dynamics and linear coupling to all the structural elements in the cochlear cross section. Recently, Szalai *et al.* [29] analysed a version of this model and various other simplified ODE models using bifurcation analysis. They find that the dynamics of the system contains many local bifurcations including saddle-node, cusp and Hopf bifurcations of the equilibrium state. They argue that the presence of a bifurcation *per se* is not the most important feature needed to reproduce the observed compressive nonlinear response; rather it is the form of the nonlinearity that is key.

We present here a simplified model of the dynamics of the organ of Corti that captures what we believe to be the simplest ingredients that are consistent with physiological observations. For more details of how this model may be seen as a reduction of that of Ó Maoiléidigh & Jülicher [22], see Szalai *et al.* [29]. The model is depicted schematically in figure 3*d*. The model includes adaptation owing to hair bundle displacement and nonlinearity of the transduction current. The actual force to move the basilar membrane is exerted via the OHC's somatic motility mechanism. We assume an elastic and underdamped basilar membrane. The hair cells connect to the basilar membrane through the elastic Deiters' cells on the one hand, and on the other hand, are rigidly fixed to the reticular lamina. The passive organ of Corti mechanism is assumed to stimulate the hair bundle, which in our model does not exert any force.

We use a simplified approach to describe the transduction current *I* as a function of the hair bundle displacement. To keep a consistent notation with earlier models, in particular, the work of Ó Maoiléidigh & Jülicher [22], we refer to the deflection of the tip link from its equilibrium position via an adaptation variable *y*_{a}. Specifically, we assume
3.1
and
3.2
where *Z* is the deflection of the hair bundle, and the open probability function *P*_{o} is given by the Boltzmann function
where the induced current has been scaled to unity if the ion channels are fully open. In contrast to previous work, we do not include separate dynamics for the hair bundle displacement *y* and the tip link *y*_{a}, as this is only a phenomenological description of the data in Kennedy *et al.* [16]. The behaviour of the model is illustrated in figure 3*a*,*b*. As in the experiments, we stimulate the hair bundle in the model with step displacements of different magnitude (shown in figure 3*a*). This results in peaks of the transduction current whose maximum value depends on the height of the stimulating step with the function *P*_{o}. When the stimulation is held, the current falls to almost zero; that is, the adaptation is nearly complete. This means that the transduction current is more of a function of the hair bundle velocity than its position. We found that the values *α*=0.1 and *β*=1 form a good fit with the experimental observations (see Szalai *et al.* [29] for more details). The parameter Δ is an overall scale factor that controls the slope of *P*_{o} at equilibrium, and can be determined from energetic considerations of the ion channels [44].

To include equations (3.1) and (3.2) as a feedback mechanism, we assume that the basilar membrane is modelled as a linear oscillator with natural frequency *ω*_{0}, the CF of the longitudinal position in question, and damping ratio *ζ*. This basilar membrane motion is supposed to be coupled to the OHC through the Deiters' cells that are again represented by a spring and a damper. We suppose that the spring and damping constants of this coupling are chosen so that the extension *l*(*t*) of the OHC causes an additional force on the basilar membrane equivalent to a reduced restoring force, , and reduced damping, .

We suppose that the extension *l*(*t*) of the OHC is controlled linearly by somatic motility. Specifically, the excess charge *q* in the hair cell is assumed to leak at rate *γ* and that the hair cell expands and contracts from its equilibrium length *l*_{0} in direct proportion to the excess charge. For simplicity, the units of the charge are chosen so that *l*(*t*)=*l*_{0}+*q*(*t*). The hair cell is charged by the transduction current. The hair bundle excitation is driven by the basilar membrane motion *z*(*t*) and indirectly by the assumed sinusoidal pressure difference . Even though the indirect forcing by the pressure is small, its effect can be rather large owing to the sensitivity of the hair bundle. A parameter *ϵ* is included to represent the relative importance of *p*(*t*) on the motion of the basilar membrane *z* and on the deflection of the hair bundles *y*_{a}.

Under the above assumptions, the governing set of equations of the organ of Corti model is
3.3
Note that the only nonlinearity is the open probability function *P*_{o}. The time constant of the OHC low-pass filter is compensated for by the current adaptation, because the open probability of the ion channels is mostly dependent on the speed of adaptation variable *y*_{a}.

The bifurcation diagram of equation (3.3) can be seen in figure 3*e*. The white region is stable and bounded by a pitchfork bifurcation (horizontal line) and a Hopf bifurcation curve (diagonal line). Both bifurcations are supercritical. This means that at the pitchfork bifurcation while the actual equilibrium loses its stability two other stable equilibria arise as the parameters depart from the stable region. In case of the Hopf bifurcation, the stability loss comes with a periodic oscillation which increases in amplitude outside the stable region. Therefore, the unforced system in the stable region has no other attractor than the equilibrium. We have found that within the white stability region, the system is rather robust to parameter changes. At our chosen set of parameters (denoted by a cross in figure 3*e*) the response of the model can be seen in figure 3*f*. The response has the desired qualitative feature of being linear for small forcing amplitude, compressive for medium amplitudes and linear again for higher amplitudes. One might notice that its tuning is rather sharp, sharper than is possible by a 1-degree of freedom (d.f.) oscillator, but less sharp than a 2-d.f. oscillator. This is due to the additional dynamics from the adaptation and charging of the hair cell.

One conclusion from these results is that the Hopf bifurcation normal form is not the only model that produces qualitatively similar results to experimental data. Even though a Hopf bifurcation is present within the bifurcation diagram of figure 3*e*, the model cannot be reduced to the Hopf oscillator as other features play a significant role. It can be argued that this Hopf bifurcation is the possible source of otoacoustic emissions owing to the consequent limit cycle oscillation that might result from a perturbation to the system's parameters away from its stable operating point. It has been suggested that spatial inhomogeneity in the cochlea might cause such a perturbation [45]. However, it may be misleading to focus solely on local dynamics because, as we outline in the next section, the combination of longitudinal coupling together with local nonlinearity seems to be the most important feature of a model to explain experimental data.

## 4. Longitudinal coupling

The organ of Corti model in §3, even though it has some desirable features, is far from complete. The cochlea cannot simply be described as an ensemble of disconnected oscillators each subject to a simple pure-tone input. In truth, there is significant longitudinal coupling between individual organ of Corti cross sections. Effects that cause this coupling include the longitudinal stiffness of the basilar membrane (which is significantly weaker than its transverse stiffness), the motion of the tectorial membrane, the fluid motion, the fact that the hair cells themselves tilt in the longitudinal direction, and that the Deiters' cells are coupled longitudinally through their so-called phalangeal processes to the reticular lamina (see figure 4*a* for a schematic).

Early models of the cochlea assumed no coupling between basilar membrane cross sections, other than the fluid motion in the ducts of the cochlea. These models could not account for every detail of the basilar membrane vibration, whether or not they included higher spatial dimensional fluid motion. On the other hand, Zweig [25] determined the impedance of the basilar membrane necessary to produce observed basilar membrane response. He found that the effective basilar membrane oscillator must have a negative damping coefficient and a linear time-delayed feedback term. Such a model is rather difficult to underpin physiologically. The negative damping could be thought of as a crude model of OHC dynamics, but there is no obvious physiological evidence for a local time delay of the appropriate duration in the organ of Corti.

However, the puzzle of Zweig can be resolved by other means. Steele & Lim [26] and Geisler [27] suggested that a longitudinal forward coupling mechanism in the cochlea can produce the same result. Recent work by Szalai *et al.* [46] provides further evidence in this direction, and an explanation of the similarities and differences between temporal delay and spatial feed-forward. Other forms of longitudinal coupling could arise from a second travelling wave, as suggested by Hubbard [47], that could be related to either the tunnel of Corti flow [48], or the waves of the tectorial membrane [49]. Recently, Meaud & Grosh [50] considered simple plate models for the basilar membrane and the tectorial membrane. They showed that the combined effects of these yield longitudinal coupling that is bidirectional. This bidirectional coupling was found to lead to a cochlear response that has broader, more realistic, peaks as well as to lead to increased stability.

In this paper, we use a combination of passive bidirectional coupling and active forward coupling. The bidirectional coupling is inspired by Meaud & Grosh [50] while the spatial feed-forward coupling is used as more realistic surrogate of the time delay in Zweig's model [46]. We then assume that the micromechanics is described by equation (3.3) with some simplifications.

Specifically, the model consists of two parts: the fluid dynamics in the chambers of the cochlea, and the micromechanical model of the organ of Corti. As is common in modelling the fluid mechanics, we assume an incompressible and inviscid fluid description of the perilymph, which yields a wave equation for the pressure difference *p* between the scala tympani and scala vestibuli, and the basilar membrane displacement *z*,
4.1
Here, the lumped parameter *ϵ*≈0.025 is a function of the geometry of the cochlear chambers and the density of the perilymph fluid, while *m* is the mass surface density of the basilar membrane. We scale the cochlea length to unity, therefore .

To describe the basilar membrane motion, we use a simplified version of our local model (3.3). First it is no longer necessary to introduce the applied forcing term into the stimulation of the hair cells. Instead, we suppose that the basilar membrane (the main elastic member within the organ of Corti) drives all mechanical components within a cross section, including the motion between the tectorial membrane and reticular lamina that stimulates the OHC hair bundles. Hence, the variable *Z* can be replaced with something strictly proportional to the basilar membrane displacement *z*. Also, we assume that the current adaptation is rather fast compared with the other dynamics, that is, the adaptation variable *y*_{a} and the basilar membrane displacement *z* are also strictly proportional. Without loss of generality, we can assume *y*_{a}=*z*. The governing equations are, therefore,
4.2

Another simplification to the model (3.3) is that we use only the charge *q* and not its derivative as the feed-forward term as a modelling convenience. However, we introduce feed-forward and feed-backward terms from the basilar membrane motion. These two terms are equal in magnitude, as they approximate the longitudinal stiffness of basilar membrane, meaning that the membrane is more like an inhomogeneous elastic plate rather than an ensemble of strings (cf. [50]). For computational simplicity here, we assume that the feed-forward distance *h*(*x*) of the action of the hair cells is the same as the longitudinal coupling distance of the basilar membrane (although this assumption can easily be relaxed [51]). A realistic estimate for the magnitude of *h*(*x*) can be computed from the length of the hair cells based on the assumption that the geometry of the organ of Corti is self-similar at different positions along the length of the cochlea. Specifically, we take *h*(*x*)=1.2*l*(*x*), where *l* is the length of the OHC.

In order to match the experimental data of Rhode [13], shown in figure 4*b*,*c* we assume the following mechanics for the basilar membrane: *l*(*x*)=(0.0027+0.0027*x*), natural frequency *ω*_{0}(*x*)=2*π*20.8(e^{−4.8354x}−0.1455) and damping coefficient *ζ*(*x*)=0.03(1−e^{−40(1−x)})+0.1e^{−40(1−x)}. We also include spatial inhomogeneity by assuming that the feedback coefficient *f*_{b} has a normal random distribution around with standard deviation . Other parameters are the same as in equation (3.3), with the exception of *γ*=1.

We have used a boundary-value problem formulation to solve equations (4.1) and (4.2). We also assume that the stapes is moved by a harmonic velocity, that is, , where *u* is the fluid velocity of the scala vestibuli. From the one-dimensional fluid flow equation, we find that , where *ρ* is the fluid density. Hence, assuming that the pressure difference at the cochlear apex vanishes, we find pressure boundary conditions
4.3
For the basilar membrane displacement *z* and charge *q*, since the response of the cochlea for a single-tone excitation is periodic, we prescribe a periodic boundary condition in time
4.4

The numerical method used to solve equations (4.1) and (4.2) subject to equations (4.3) and (4.4) is a finite differences discretization in space and a piecewise-polynomial approximation, so-called orthogonal collocation, in time [52]. The solution of the nonlinear equation is propagated with increasing forcing amplitude *A* using pseudo-arclength continuation [28]. Preliminary numerical results are shown in figure 4*d*,*e*. More details will be presented elsewhere.

Figure 4*b*,*c* shows, for comparison, experimental results from Rhode [13] at 9.1 kHz CF. Figure 4*b* shows the input–output functions for different stimulus frequencies. For frequencies substantially lower than CF, the response is linear. For frequencies closer to CF, the response becomes compressive for mid-range amplitudes and tends towards linear for high amplitudes. The compressive nature remains for very high frequencies but the linear region for higher amplitudes starts earlier. In figure 4*c*, the cochlear response is shown as a function of frequency at constant amplitudes. For low amplitudes, the response has a distinct broad peak at the CF. With increasing amplitude the peak gets less pronounced and slightly moves towards lower frequencies. The amplitude independent variation at lower frequencies observed in experiments could be due to the middle ear, while wiggles at higher than CF frequencies might be explained by interference of vibrations from the travelling wave and a faster wave [53]. Another observation is that in figure 4*b* the response at 13 kHz is nearly linear, even a little expansive. This might also be explained by a second sound source perhaps from the compressive wave or the basilar membrane longitudinally coupled dynamics.

In comparison, our models (4.1) and (4.2) yield similar results in figure 4*d*,*e*. The significant differences are that the numerical results show sharper responses, which yield somewhat thinner peaks in figure 4*e*, and curves that are somewhat further apart in figure 4*d*. Because we do not model the middle ear, amplitude disturbances due to it (at low frequencies) are absent in figure 4*e*. However, surprisingly the disturbances at high frequencies are reproduced in figure 4*e*. This is especially striking since the model does not include fast compression waves or any other source that the travelling wave could interfere with. A control simulation with longitudinally smooth feedback shows identical dips in the tuning curves, hence the inhomogeneity of parameters is similarly ruled out as a possible cause. Also, surprisingly, in figure 4*d* for the 13 kHz curve we get a similar response to the experimental result. Our speculation is that the early linearity of high frequency input–output functions is because of the longitudinal coupling, without which the curve would be flat up until high-amplitude forcing. The exact mechanism, however, is still under investigation.

## 5. Conclusions

This paper has reviewed the dynamical properties of the mammalian cochlea, and presented some preliminary results on modelling its specific nonlinear components.

First, we presented a new model, supported by new data, of how IHC calcium currents might play a role in the development of hearing. IHC Ca^{2+} signalling pathways are thought to be critical to both the intrinsic development of IHCs and the auditory system as a whole. Understanding the basic elements of these signalling pathways in the healthy ear is essential in order to approach future strategies for hair cell regeneration or replacement therapies that may be used to treat conditions such as age-related hearing loss. Such understanding could be greatly facilitated by using complementary approaches of experimentation and mathematical modelling.

Next, we reviewed the dynamics of the OHCs, which are believed to be the source of cochlear amplification. The nonlinearity of these cells appears to play a significant role in cochlear function, leading to compressive nonlinearity in the basilar membrane response and in sound perception more generally. Based on our synthesis of earlier micromechanical models, we introduced a longitudinally coupled cochlear model that includes many aspects of the cochlear mechanism already identified by others. In particular, we argue that it is the combination of spatial feed-forward, longitudinal coupling, somatic motility and nonlinear MET channel open probabilities that enables our model to agree qualitatively and quantitatively with data. As a deficiency, our model shows extreme compression near saturation. This can be remedied by an additional linear coupling between the charge of the OHC and the basilar membrane displacement that represents the reverse piezoelectric mechanism of prestin.

Note that all of this agreement is found by physiological modelling without assuming *a priori* that the hair bundles are tuned to undergo a local bifurcation. This finding seems to run counter to some current opinion in the literature which suggests the normal form for a Hopf bifurcation is a suitable local model for the mammalian cochlea (e.g. Hudspeth *et al.* [23]). Indeed, proponents of a Hopf bifurcation theory point to a supposed 1/3-power law in the amplification data (which does not seem to be universal). In fact, it can be shown by harmonic balance [29] that a 1/3-power law comes from the presence of the lowest order nonlinear term in any nonlinear function.

One key question that remains, then, is what is the true source of the nonlinearity in the OHC dynamics? Here, we have assumed an open probability function that depends on a mixture of the displacement and rate of displacement of the hair bundles. This provides a good agreement with the experimental data, but does not shed light on the biochemical processes that might cause such a response. This is where future work that incorporates ion-channel modelling as we described for the IHCs might be helpful. In particular, Verpy *et al.* [54] have found that stereocilin-deficient mice lack signatures of nonlinearity, but preserve amplification and sharp tuning that progressively vanishes. Determining the effect of stereocilia links on MET is therefore another avenue to pursue.

One key mathematical technique we have used in this paper and seems especially important when investigating nonlinear models is the use of numerical bifurcation analysis. This method can be useful to test for robustness of mechanisms, to point to possible multi-stability or plasticity and the proximity of instability. Bifurcation analysis works best on simplified models which capture the essential features believed to be at work.

More generally though, as the nature of the cochlear amplifier is still debated, we need to use more comprehensive simulation models that incorporate the true geometry of the cochlea in order to decide what modelling assumptions play essential roles. Even though linear, the model of Meaude and Grosh described in Ashmore *et al.* [24] comes closest to this requirement. Careful parametrization of such models remains a major issue, as does the correct treatment of noise, inhomogeneity and within-species heterogeneity.

## Footnotes

One contribution of 11 to a Theme Issue ‘Towards the virtual physiological human: mathematical and computational case studies’.

- This journal is © 2011 The Royal Society