An elemental approach to modelling the mechanics of the cochlea

&NA; The motion along the basilar membrane in the cochlea is due to the interaction between the micromechanical behaviour of the organ of Corti and the fluid movement in the scalae. By dividing the length of the cochlea into a finite number of elements and assuming a given radial distribution of the basilar membrane motion for each element, a set of equations can be separately derived for the micromechanics and for the fluid coupling. These equations can then be combined, using matrix methods, to give the fully coupled response. This elemental approach reduces to the classical transmission line model if the micromechanics are assumed to be locally‐reacting and the fluid coupling is assumed to be entirely one‐dimensional, but is also valid without these assumptions. The elemental model is most easily formulated in the frequency domain, assuming quasi‐linear behaviour, but a time domain formulation, using state space method, can readily incorporate local nonlinearities in the micromechanics. Examples of programs are included for the elemental model of a human cochlea that can be readily modified for other species. HighlightsGeneral formulation of an elemental model for cochlear mechanics.Reduce to the transmission line model for locally‐reacting micromechanical and 1D fluid coupling.Incorporation of non‐uniform areas, 3D fluid coupling and non locally‐reacting micromechanics.MATLAB programs for the elemental model in the frequency domain and time domain.


Introduction: modelling the mechanics of the cochlea
It is important to develop mathematical models of the mechanics of the cochlea in order to test our understanding of the physical processes involved in hearing. These models can have various levels of abstraction, depending on what features of the cochlear mechanics they are attempting to describe. Wave models, for example, are good at describing the global behaviour of the motion along the cochlea without getting too involved in the detailed physical processes that give rise to the wave motion. The most widely used wave model uses the WKB method (de Boer et al., 1982;Steele et al., 1979a;Steele et al., 1979b;Wang et al., 2016), in which the response along the cochlea can be predicted from an assumed distribution of complex wavenumbers. There are, however, several assumptions inherent in using the WKB method. Although it can account for both forward and backward travelling waves in the cochlea, its main assumption is that there is only one type of wave propagating. This is generally called the "slow wave", to distinguish it from the so-called "fast wave" that travels at the compressional speed of sound in the cochlear fluids. The fast wave is not thought to play an important role in normal hearing, since the pressure is the same in the two fluid chambers and so there is no pressure difference acting on the basilar membrane, BM, which is then not driven into motion. There are, however, several other kinds of wave that might play a significant role in determining the mechanical response of the cochlea, including those due to higher-order fluid modes (Watts, 2000;Elliott et al., 2013) and those due to mechanical coupling along various parts of the organ of Corti or the fluid within it (Zwislocki et al., 1979;Ghaffari et al., 2007;Ramamoorthy et al., 2007;Lamb et al., 2011;Meaud et al., 2010). The way that these and other modes couple into the slow wave is still a matter of some debate.
An alternative approach to modelling the mechanics of the cochlea is to divide the length of the cochlea into a finite number of elements and describe the individual physical processes involved in cochlear mechanics. These equations can be coupled together in a numerical model that makes no explicit assumptions about the type of wave propagation. The earliest example of this approach is the "transmission line" model (Shera et al., 1991;Zweig, 1991;Zwislocki, 1950), in which the inertance of the fluid in the cochlear chambers is represented as a series inductance and the response of the BM is represented as a shunt impedance. At the other extreme, the detailed behaviour of the fluid in the chambers and the motion of the organ of Corti could be represented by a finite element model, where both the fluid in the chambers and the different parts of the organ of Corti are meshed into many small elements and the M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 5 coupled set of equations of motion are, typically, solved using commercial software such as ANSYS (Ni et al., 2016). It is important to note that such finite element methods simultaneously solve for both the pressure in the fluid and the motion of the organ of Corti.
This can significantly increase the computational time compared to solving for the fluid pressure and the organ of Corti motion individually, since the computational time typically rises in approximate proportion to the square of the total number of degree-of-freedom in a finite element model. It also means that both the fluid coupling and the organ of Corti motion must be analysed numerically, even if there may be a simple analytic formulation in one case or the other.
An elemental approach to the formulation of cochlear mechanics can be viewed as a generalisation of the transmission line model, in that the fluid motion and the organ of Corti dynamics are analysed separately and then coupled together. In the elemental model, however, the organ of Corti dynamics are not restricted to being locally-reacting and the fluid flow is not restricted to being proportional to the pressure difference between adjacent elements. The elemental model reduces to a transmission line model if these restrictions are imposed, but can also account for the more general case. This paper first introduces the elemental model, for the case of a locally-reacting basilar membrane and 1D fluid coupling in a uniform box model. If the BM response is linear, the coupled solution can be solved efficiently in the frequency domain and this formulation is considered first. A time domain formulation is considered at the end of this paper which leads to a state space implementation of the elemental model that can be used to simulate the nonlinear response of the cochlea but is also important in assessing the stability of linear models.
It is shown how the fluid coupling part of the formulation can be adapted to model nonuniform distributions of fluid chamber area and the near field component of the pressure that is present when the fluid coupling is analysed in 3D. More general forms of the BM response are also then considered, where both longitudinal mechanical coupling along the BM and non-symmetrical feedforward behaviour along the organ of Corti can be accounted for.
The theoretical formulation behind the fluid coupling within the cochlea and its passive micromechanics are reasonably well understood. The intention of the present review paper is to demonstrate how these individual responses may be integrated into a numerical formulation that can be readily used to predict the coupled response of the cochlea under different conditions. Although there have been a number of micromechanical models proposed for the active cochlea, there is still considerable debate about the form of such a model and how it should couple into the fluid. Since this is still unclear, and also for reasons of brevity, the elemental model will be developed and illustrated here only using passive micromechanical models, on the understanding that active models can be incorporated with a more general form of BM admittance.

Formulation of the elemental model in the frequency domain
In this section, we assume a linear micromechanical response, and an excitation of the cochlea at a single frequency, so that the amplitude and phase of all the dynamic variables can be considered complex and proportional to e iωt . The essence of the elemental model is that the length of the cochlea is divided into a finite number of sections and that the interaction between the fluid in the chambers and the mechanical response of the BM is described in terms of the transverse modes of the element in each section. The general form of the assumed model is illustrated in Figure 1, which also shows the coordinate system that will be used below. The cochlea is shown uncoiled since the coiling is not thought to play an important role in its dynamics (Steele et al., 1985). The origin of the coordinates system is on the basilar membrane at the base of the cochlea and it consists of a longitudinal variable, x, a radial variable, y, and a transverse variable, z. The aim of the elemental model is to reduce the complicated interaction that occurs between the fluid and the motion of the BM in all three directions to a single set of discrete equations, as a function of only the longitudinal variable, x.

Frequency domain formulation
In general, there may be multiple transverse modes across the basilar membrane, but in practice the fluid coupling is relatively insensitive to the exact form of the radial BM velocity distribution  and so only a single transverse mode is assumed. If the normalised shape of this transverse mode is ψ(y) and v(x, y) is the complex BM velocity distribution in the longitudinal and radial directions at a frequency ω, the dependence of which is suppressed for notational convenience, then the single variable can be used to express the complex BM velocity at the longitudinal position x, as defined as where W is the width of the cochlear partition and the normalisation of the mode shape is assumed to be such that (2) In a similar way, a single variable describing the difference in pressures in the two main fluid chambers that act on the element is defined to be where p 1 (x, y, z) and p 2 (x, y, z) are the three-dimensional pressure distributions in the upper and lower fluid chambers.
The continuous distribution of BM velocity and pressure, v(x) and p(x), are now spatially sampled to give a finite set of BM velocities and pressures at a number of discrete points, v n and p n , where the integer variable, n, runs from 2 to N. The first element in the set of velocities and pressures, n=1, in these arrays is reserved to describe the excitation at the oval window from the middle ear, and the final element, n=N, is used to define the conditions at the helicotrema. Although it is not essential to the formulation, we assume here that the spatial sampling is uniform, so that if the length of the cochlear is L then the length of the elements in the longitudinal direction is always equal to L/N, which is defined as ∆. The number of elements and hence their length is determined by the shortest wavelength that needs to be reproduced in the longitudinal direction. A rule of thumb commonly used in finite element analysis is that they should be at least six elements within the shortest wavelength, so in practice it is common to use about 500 elements to describe the wave motion, which may have a wavelength as small as 0.5 mm (Shera, 2007), in a human cochlea of length 35 mm.
Vectors of these discrete spatial samples of complex pressure and BM velocity are now defined as which are coupled both because of the fluid in the chambers and also because of the micromechanical response of the BM. The fluid coupling in the chambers is now assumed to generate a set of pressures, p, equal to a fluid coupling matrix, Z FC , multiplied by the vector of elemental velocities v, The vector of BM velocities, v, is also assumed to be equal to the sum of an excitation vector, which contains the imposed stapes velocity, v s , minus the matrix of BM admittance, Y BM , multiplied by the pressure vector, p, The negative sign in equation (6) accounts for the fact that in our convention a positive pressure difference drives the basilar membrane in a downward direction, whereas the elements of the vector v are defined to be those in an upward direction. Normally all the elements of the excitation vector v s are taken to be zero except the first element, which is equal to the normalised stapes velocity. More generally, it is also possible to represent internal velocity excitation of the elements along the basilar membrane by setting other elements of v s to be equal to this internal excitation. This can be used, for example, to model the generation of DPOAEs (Kanis et al., 1997;Young et al., 2012).
By combining equations (5) and (6) This formula is the same as that used to solve similar coupled fluid-structural interaction problems in engineering, when they are formulated in elemental terms (see, for example, Fahy and Gardonio, 2007).
Substituting equation (7) into equation (5), the expression for the vector of complex pressures is also obtained as which will be used below.

1D fluid dynamics
The cross-sectional area of the two fluid chambers is initially assumed to be constant along the cochlea, as in the uniform box model (de Boer, 1996). If the pressure in each of the fluid chambers is assumed to be uniform over each cross-section, then the fluid dynamics can be described as being 1D. Physically, the condition under which it is valid to use 1D fluid dynamics is when the shortest wavelength of the slow wave in the cochlea is greater than the cross-sectional dimensions of the fluid chamber. This condition is not a bad approximation in the passive cochlea, but becomes less valid when the cochlea is active and the shortest wavelength becomes smaller, in which case a 3D description of the fluid dynamics, as discussed below, should be used.
Assuming 1D fluid coupling, the equation for the fluid mass and momentum conservation can be combined to provide a differential equation relating the longitudinal distribution of the complex pressure difference across the BM, p(x), to that of the complex BM velocity distribution, v(x), both proportional to e iωt , as where h is the effective height of the fluid chambers, given by π 2 A/8B, where A is the area of each chamber and B is the width of the BM, as discussed below. Assuming that the wavelength of the slow wave is small compared with the length of the elements, then the second spatial derivative in equation (10) can be approximated using a finite difference approach as The boundary condition at the base of the cochlea can be written in terms of the momentum equation for the complex pressure as where v 1 is the stapes velocity when p 1 is equal to 0, i.e. when the middle ear in unloaded, as described below. Using the finite difference approximation for the spatial derivative in equation (12), we obtain The boundary condition at the apex of the cochlea, at the helicotrema, is generally taken to be pressure release, i.e. the pressure difference across the BM is equalised in the two fluid chambers so that p N =0.
These finite difference equations can be written in matrix form as or more compactly in terms of the vectors of velocity and pressure as where F is the matrix of finite difference terms above. Assuming that this matrix is not singular then we can write the vector of pressures as which is of the form assumed above for the fluid dynamics with The need to be able to invert this matrix requires that a finite value be entered into its bottom right element, taken here as ∆ 2 /h 2 but discussed more fully below.
which is of the form originally used by Neely (1981), so that s BM 1 iω The

Locally-reacting micromechanics
The mechanical coupling along the length of the BM is relatively weak and longitudinal coupling in the cochlea is dominated by that due to the fluid in the chamber, so that a reasonable first approximation is to assume that the BM is locally reacting. The BM is said to be locally reacting if the velocity of an element at one position depends only on the pressure difference at that point, so that where Y BM (n) is the mechanical admittance of the BM and the negative sign has been discussed above. If the passive BM dynamics are approximated by a single degree of freedom system, then where ω is the excitation frequency and m n , s n and r n are the mass, stiffness and damping, per unit area, of the BM at the position of the n-th element, i.e. at a distance of x n =n∆ along the cochlea.
It is assumed in the example below that BM mass, m n , is constant along the cochlea, and is denoted, m 0 . The BM stiffness is chosen to match the distribution of natural frequencies along the cochlea, which is assumed to take the form where f B is the natural frequency at the base and l is a characteristic length for the frequency variation, so that The BM damping is chosen to ensure that the Q factor of the micromechanical elements is constant along the cochlea (de Boer, 1996), and equal to Q 0 , so that 0 0 n n s m r Q = . (24) To account for the mechanical loading of the middle ear by the cochlear input impedance we can write the velocity of the first element as where v st is the unloaded stapes velocity and Y ME is the admittance looking into the middle ear from the cochlear. This admittance corresponds to the reciprocal of the parameter M 3 in the analysis of Puria (2003), who shows that this admittance is well approximated by that of a single degree-of-freedom system and who gives appropriate values for the mass, stiffness and damping per unit area in this case.

Transmission line model
For the special case of 1D fluid coupling and locally-reacting micromechanics, we show that the elemental model can be interpreted in terms of a transmission line, as shown in Figure 2.
The BM velocity is the analogue of current in this network, and pressure difference is the analogue of the voltage. The velocity into the n-th shunt element, shown in Figure 3, where n is greater than one, is given by where L is the inertance of the longitudinal coupling elements, p n is the pressure at the n-th elements and v n denotes the BM velocity in the upward direction. This equation can also be written as By comparing this with equation (11) we find that the series inertance due to the 1D fluid coupling is generally equal to The exception to this is the first series element, which can be seen from equation (12) to have the value The velocity of the n-th elements, for n greater than one, is equal to the BM admittance multiplied by the pressure at this point, as expressed in equation (20). The velocity into the first element is given by equation (25), so that the source term at the stapes is indicated as a current source in Figure 2 and Y ME is the first shunt element. The transmission line is terminated by short-circuit at the apical end to impose the condition at the helicotrema: that p N is equal to zero. Although the elemental approach reduces to the transmission line model, which provides considerable insight from an electrical engineering perspective under these restrictive assumptions, it should again be emphasised that it is not necessary to make these assumptions for the elemental approach to be valid, as discussed below. The elemental model should thus be seen as a generalisation of the transmission line approach.

Example of an elemental model of the human cochlea
The human cochlea is chosen as an example to illustrate the elemental approach with the parameters listed in Table 1. The choice of the number of elements used in the model, N, is a ... ...  Figure 3. This elemental approach, however, is capable of modelling more complicated cases, for example, taking cross-sectional area variation and 3D fluid dynamics into account, as discussed later in Section 3. The coupled BM responses are also shown in Figure 3, for comparison, when both area non-uniformity, assumed to be linearly tapered (Ni et al., 2017), and 3D fluid dynamics are considered. The linear variation in the physical dimensions of the fluid chambers along the length of the cochlea is based on Thorne et al. (1999), and the average values of those at the two ends are used for the uniform model. The modelling parameters are listed in Table 1 and example code can be found in supplementary material.   (Wever, 1949) BM width at the apex (mm) B A 0.50 (Wever, 1949) Frequency at the base (kHz)

Non-uniform scala areas
It has been assumed above that the two fluid chambers have a uniform area along the length.
This restriction can be removed by modifying equation (10)  and assuming the pressure difference cross the BM is zero at x = L, account for the helicotrema, the 1D, far-field, fluid dynamics can be given as where ‫′ݔ‬ is dummy integration variable, A e (x) is the effective area, given by the harmonic mean of the upper and lower areas of the fluid chambers. It should be noted that equation (30) is for the pressure distribution basal to the excitation point and equation (31)  x to L to give such a variation in pressure.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
If the areas of the fluid chambers in the cochlear models are divided up into N discrete sections, as for the BM, the integrals in equations (30) and (31) can be approximated by summations to give the pressure at the n-th element as where n 0 = x 0 /∆.

3D fluid dynamics
The vector of pressures in the elemental model is related to the vector of velocities by the fluid coupling matrix, equation (5). The physical significance of the columns of the matrix Z FC in this equation are that they give the pressure distribution along the cochlear due to a single element of the BM vibrating, with all the other elements being fixed. If Z FC is calculated assuming 1D fluid coupling in a uniform cochlear, so that Z FC is given by equation (17), then this pressure distribution is given by the grey dashed line in of excitation. The near-field pressure is the difference between these two cases and its magnitude and longitudinal distribution will depend on the local geometry around the BM.
Since the near-field does not generally extend to the boundaries of the fluid chambers, however, it does not depend significantly on their cross-sectional shape, which considerably simplifies the modelling of the fluid coupling.
The factors that have the greatest effect on the near-field pressure are the width and the position of the BM across the cochlear partition, which separates the main two fluid chambers. It is shown in Elliott et al. (2011) and Ni and Elliott (2015) that to a good approximation the near-field pressure distribution due to the vibration of a single element is a decaying exponential function where x 0 is the position of vibrating element, ω is angular frequency, v 0 is the linear velocity of the BM element, ∆ is its width, B is the width of the BM, W is the width of the cochlear partition, and a and b are fitting coefficients (Ni and Elliott, 2015). Assuming that the BM begins from one side of the cochlear partition, the spatial ligament side, and that B/W is about 0.3, the peak near field pressure can be expressed as approximately, where the factor of 2 is a parameter calculated based on fitted values of a and b (Ni and Elliott, 2015). Similarly, the decay distance, ݀ ≈ √ܾܹ, is approximately equal to In fact, the 3D fluid coupling analysis has already been used in the definition of the effective height of the fluid chambers in the 1D case (Elliott et al., 2011). There is, however, only a small difference between the distribution of BM velocity along the cochlear predicted for a passive model with either 1D or 3D fluid coupling, as shown for example in Figure 8 of Elliott et al. (2011). This comparison assumes that the mass of the BM in the 1D case is increased to account for the fluid loading. This increases the BM mass from about 0.05 kg/m 2 in the 3D case, corresponding to a physically reasonable thickness for the organ of Corti of 50 µm, to a value of about 0.3 kg/m 2 in the 1D case, since there is effectively a layer of fluid of total thickness about 250 µm that moves with the BM in this case and has to be accounted for in the micromechanics rather than the fluid coupling (Neely, 1981;Elliott et al., 2011). Figure 4 The pressure distribution along a uniform (grey lines) and non-uniform (black lines) cochlea due to excitation by the vibration of a single BM element, at x 0 in the case of 1D fluid coupling (dashed) and 3D fluid coupling (solid). In this example, the human cochlear geometrical data for the non-uniform model is the same as that used in Ni et al. (2017). MATLAB Codes for generating the figure can be found in supplementary material.

Non locally-reacting micromechanics
The assumption made in Section 2.2 above is that the BM velocity at one position along the cochlea depends only on the pressure difference at that same point, i.e. it is locally reacting.
There are, however, a number of mechanisms for longitudinal coupling along the BM, which will give rise to non locally-reacting behaviour, and this may be readily incorporated into the element model. The effect of these complexities in the model on its predicted response will depend on the extent of the assumed longitudinal coupling and the magnitude of the BM admittance. If reasonable values are assumed for the longitudinal coupling it is found that it again makes little difference to the predicted response in the passive cochlea, but would make a greater difference in a fully active one.
The first case of longitudinal coupling along the BM that we consider is due to its orthotropic behaviour, i.e. although its stiffness in the longitudinal direction is less than that in the radial direction, this longitudinal stiffness is not negligible (Emadi et al., 2004;Meaud et al., 2010; Naidu et al., 2001). Other authors have also emphasised the role of longitudinal damping along the BM (Chan et al., 2015;Fisher et al., 2012;Reichenbach et al., 2014). The simplest method of accounting for such behaviour in the elemental model is assuming that adjacent BM elements are coupled by a longitudinal stiffness, k L , and damping c L , as shown in Figure   5.
In order to incorporate these elements into the BM admittance matrix, we first write down the equation for the pressure difference on the n-th element, from Figure 5 where m n , k n and c n are the locally-reacting mas, stiffness and damping of the n-th BM element, assumed to be passive here for convenience, and k L and c L are the longitudinal coupling stiffness and damping terms. This equation can be written as A tri-diagonal matrix of BM impedances can thus be assembled from these elements as illustrated in Figure 6(a), so that in which case, assuming Z BM is not singular, the BM admittance matrix required in the elemental model can be obtained by inverting Z BM . The resulting form of the BM admittance matrix is shown in Figure 6(c).

M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT Figure 5 Longitudinal coupling along the cochlea represented by stiffnesses, k L , and dampings, c L , between adjacent elements (a) and the feedforward connection of the OHC (b) (Reprinted from Fig. 1 (a) in Yoon et al. (2011) with permission).
Another form of longitudinal coupling that exists within the organ of Corti is illustrated in Figure 5 (b) and is due to the way that the bottom of each outer hair cell, OHC, is attached to both to the BM by the Deiter's cell but also to more apical position along the reticular lamina by the phalangeal process (Fukazawa, 2002;Geisler et al., 1995;Steele et al., 1993;Yoon et al., 2011). Several models for this "feedforward" behaviour have been proposed, but one of the simplest was suggested by de Boer (2007), as a similar modification to the BM impedance model used for longitudinal coupling above, but acting only in one direction, so where ∆x 3 is the longitudinal length scale over which the phalangeal process acts which is equal to ∆x 2 -∆x 1 , in Figure 5 (b). Although this feedforward distance varies from species to species and in one species varies along the length of the cochlea (Yoon et al., 2011), it is typically about 20 µm. This is much less than a wavelength of the slow wave, even in an active cochlea and so significantly less than the spatial discretisation that would normally be used in an elemental model. If a larger number of elements were used, so as to incorporate this feedforward effect, so that L/N=∆x 3 , then equation (42) can be written in discrete form as The impedance matrix now takes the form shown in Figures 6(b) and the columns of the BM admittance matrix take the form shown in Figures 6(d). In order to provide an illustration of this formulation, however, simulation have again been performed with 512 elements and it has been assumed, rather arbitrarily, that k L and c L are 10% of k n and c n for longitudinal coupling in equation (36) Figure 6 (d). Figure 6 The structure of the BM impedance matrix in the case of symmetrical longitudinal coupling (a) as illustrated in Figure 5 (a), and feedforward action in the organ of Corti (b), as illustrated in Figure 5 (b), together with the longitudinal distribution of BM admittance corresponding to the velocity distribution along the cochlea when forced at a single point, given by a column of Y BM =Z BM -1 . In the longitudinal coupling case, Z L is assumed to be 10% of Z BM at each location, whereas a value of 50% is assumed for the feedforward case to more clearly show the effect.

Time domain formulation
The elemental model can also be formulated in the time domain, to give a state space description of the cochlear mechanics (Elliott et al., 2007). It is more convenient in this formulation to include the excitation vector at the stapes in the fluid coupling part of the Noting that the matrix F defined above has entirely real and frequencyindependent elements, we can thus write the fluid coupling as where F is exactly as above, where x n (t) is the vector of state variables for the n-th element, which has a pressure difference p n (t) acting on it, resulting in a velocity v n (t).
For example, if the passive micromechanics are described by a single degree-of-freedom system as above, then there are two state variables per element, x 11 (t) and x 12 (t), corresponding to the displacement and velocity of the element, and the equations above can be written as  x v x .
The dynamics of the middle ear can be similarly formulated to relate p 1 (t) to v 1 (t).
The uncoupled micromechanics of all the elements can thus be gathered together in a combined matrix equation where M A N U S C R I P T A C C E P T E D  T  T  T  T  1  2   T  T  T  T  T  1 Substituting the expression for p(t) above into this and noting that ( ) ( ) ɺ . Since all the elements of C E are real, the state space model for the coupled elemental model can be written as

Summary and conclusions
This paper describes an elemental approach to modelling the mechanics of the cochlea, particularly the interaction between the BM dynamics and the fluid coupling. The elemental model reduces to the classical transmission line formulation for a box model if the BM dynamics are locally-reacting and the fluid coupling is 1D, but is seen to be valid even without these restrictive assumptions. In particular, modifications are described for the model to incorporate non-uniform scala and BM dimensions, 3D fluid coupling and both symmetric and feedforward forms of longitudinal coupling along the BM.
Elemental formulations are described in the frequency domain, for linear models, and in the time domain, with the potential for extension to nonlinear models, and MATLAB programmes are included for these models in the supplementary material.