A fully-coupled subwavelength resonance approach to modelling the passive cochlea

TL;DR

This paper models the cochlea using a fully-coupled subwavelength resonance approach, employing layer potential techniques and numerical methods to analyze acoustic scattering, resonant modes, and signal decomposition, providing mathematical insights into cochlear behavior.

Contribution

It introduces a novel mathematical model of the cochlea based on coupled subwavelength resonators, explaining key biological phenomena.

Findings

Mathematical explanation for the cochlea's 'travelling wave' behavior

Demonstration of the model's ability to decompose incoming signals

Computation of resonant modes matching biological observations

Abstract

The aim of this paper is to understand the behaviour of a large number of coupled subwavelength resonators. We use layer potential techniques in combination with numerical computations to study the acoustic pressure field due to scattering by a graded array of subwavelength resonators. Using this method, we study a graded-resonance model for the cochlea. We compute the resonant modes of the system and explore the model's ability to decompose incoming signals. We are able to offer mathematical explanations for the cochlea's so-called "travelling wave" behaviour and tonotopic frequency map.

Tables

Table 1. Table 1: Approximate values for the material parameters of a biological cochlea. Unless specified, the values are taken from [ 10 ] .

Quantity	Approximate Value
$L$: length of uncoiled cochlea	3.5cm
$w$: width of basilar membrane	0.015cm at base to 0.056cm at apex
$h$: average thickness of basilar membrane	0.002cm
$r$: average radius of scalae	0.1cm
$E$: Young’s modulus of basilar membrane	${10}^8$$N{m}^{-2}$ at base to ${10}^7$$N{m}^{-2}$ at apex
$\kappa$: bulk modulus of water	$2\times {10}^9$$\mathrm{Pa}$ [24]

Equations

$$\begin{cases}\left(\nabla\cdot\frac{1}{\rho}\nabla-\frac{1}{\kappa}\frac{\partial^{2}}{\partial t^{2}}\right)p=0,&\text{for }(x,t)\in\mathbb{R}^{2}\setminus\overline{D}\times\mathbb{R},\\ \left(\nabla\cdot\frac{1}{\rho_{b}}\nabla-\frac{1}{\kappa_{b}}\frac{\partial^{2}}{\partial t^{2}}\right)p=0,&\text{for }(x,t)\in D\times\mathbb{R},\\ p_{+}-p_{-}=0,&\text{for }(x,t)\in\partial D\times\mathbb{R},\\ \frac{1}{\rho}\frac{\partial p}{\partial\nu_{x}}\big{|}_{+}-\frac{1}{\rho_{b}}\frac{\partial p}{\partial\nu_{x}}\big{|}_{-}=0,&\text{for }(x,t)\in\partial D\times\mathbb{R},\\ p^{s}:=p-p^{in}=0,&\text{for }x\in\mathbb{R}^{2},t\ll 0,\end{cases}$$

$$v=\sqrt{\frac{\kappa }{\rho }},v_b=\sqrt{\frac{{\kappa }_b}{{\rho }_b}},k=\frac{\omega }{v},k_b=\frac{\omega }{v_b},$$

$$\mu :=\frac{{\kappa }_b}{\kappa },\delta :=\frac{{\rho }_b}{\rho },\tau :=\frac{v_b}{v}=\sqrt{\frac{\rho {\kappa }_b}{{\rho }_b\kappa }}.$$

$$\mu \ll 1,$$

$$\begin{cases}\left(\Delta+k^{2}\right)u(x,\omega)=0,&\text{in }\mathbb{R}^{2}\setminus\overline{D},\\ \left(\Delta+k_{b}^{2}\right)u(x,\omega)=0,&\text{in }D,\\ u_{+}-u_{-}=0,&\text{on }\partial D,\\ \delta\frac{\partial u}{\partial\nu}\big{|}_{+}-\frac{\partial u}{\partial\nu}\big{|}_{-}=0,&\text{on }\partial D,\\ u^{s}:=u-u^{in}\text{ satisfies the SRC},&\text{as }|x|\to\infty.\end{cases}$$

$$\underset{|x|\to \infty }{\lim }|x{|}^{1/2}\left(\frac{\partial }{\partial |x|}-ik\right)u(x,\omega )=0.$$

$${\mathcal{S}}_D^k[\phi ](x):={\int }_{\partial D}{\Gamma }^k(x-y)\phi (y)d\sigma (y),x\in \partial D,\phi \in L^2(\partial D),$$

$${\mathcal{K}}_D^{k,\ast }[\phi ](x)={\int }_{\partial D}\frac{\partial {\Gamma }^k(x-y)}{\partial {\nu }_x}\phi (y)d\sigma (y),x\in \partial D,\phi \in L^2(\partial D).$$

$$u=\{\begin{array}{ll}{u}^{in}(x)+{\mathcal{S}}_D^k[\psi ](x),& x\in {\mathbb{R}}^2\backslash \overline{D},\\ {\mathcal{S}}_D^{k_b}[\varphi ](x),& x\in D,\end{array}$$

$$\mathcal{A}(\omega ,\delta )\left(\begin{array}{c}\varphi \\ \psi \end{array}\right)=\left(\begin{array}{c}{u}^{in}\\ \delta \frac{\partial u^{in}}{\partial {\nu }_x}\end{array}\right),$$

$$\mathcal{A}(\omega,\delta):=\left[{\begin{array}[]{cc}\mathcal{S}_{D}^{k_{b}}&-\mathcal{S}_{D}^{k}\\ -\frac{1}{2}Id+\mathcal{K}_{D}^{k_{b},*}&-\delta(\frac{1}{2}Id+\mathcal{K}_{D}^{k,*})\\ \end{array}}\right].$$

$$\mathcal{A}(\omega ,\delta )={\mathcal{A}}_0+{\omega }^2\ln \omega {\mathcal{A}}_{1,1,0}+{\omega }^2{\mathcal{A}}_{1,2,0}+\delta {\mathcal{A}}_{0,1}+O(\delta {\omega }^2\ln \omega )+O({\omega }^4\ln \omega ),$$

$$\mathcal{A}_{0}:=\left[{\begin{array}[]{cc}\hat{\mathcal{S}}_{D}^{k_{b}}&-\hat{\mathcal{S}}_{D}^{k}\\ -\frac{1}{2}Id+\mathcal{K}_{D}^{*}&0\\ \end{array}}\right],\mathcal{A}_{1,1,0}:=\left[{\begin{array}[]{cc}v_{b}^{-2}\mathcal{S}_{D,1}^{(1)}&-v^{-2}\mathcal{S}_{D,1}^{(1)}\\ v_{b}^{-2}\mathcal{K}_{D,1}^{(1)}&0\\ \end{array}}\right],$$

$$\mathcal{A}_{1,2,0}:=\left[{\begin{array}[]{cc}v_{b}^{-2}(-\ln v_{b}\mathcal{S}_{D,1}^{(1)}+\mathcal{S}_{D,1}^{(2)})&-v^{-2}(-\ln v\mathcal{S}_{D,1}^{(1)}+\mathcal{S}_{D,1}^{(2)})\\ v_{b}^{-2}(-\ln v_{b}\mathcal{K}_{D,1}^{(1)}+\mathcal{K}_{D,1}^{(2)})&0\\ \end{array}}\right],$$

$$\mathcal{A}_{0,1}:=\left[{\begin{array}[]{cc}0&0\\ 0&-(\frac{1}{2}Id+\mathcal{K}_{D}^{*})\\ \end{array}}\right].$$

$${\mathcal{S}}_D[\varphi ](x)$$

$${\hat{\mathcal{S}}}_D^k[\varphi ](x)$$

$${\mathcal{S}}_{D,1}^{(1)}[\varphi ](x)$$

$${\mathcal{S}}_{D,1}^{(2)}[\varphi ](x)$$

$${\mathcal{K}}_{D,1}^{(1)}[\varphi ](x)$$

$${\mathcal{K}}_{D,1}^{(2)}[\varphi ](x)$$

$$\dim \ker {\mathcal{S}}_D\le 1.$$

$$f=\frac{\psi }{{\int }_{\partial D}\psi }-\frac{\tilde{\psi }}{{\int }_{\partial D}\tilde{\psi }}.$$

$${\hat{\mathcal{S}}}_D^k[y]={\mathcal{S}}_D[y]+{\eta }_k{\int }_{\partial D}y=0.$$

$${\eta }_k\left({\int }_{\partial D}y\right)\left({\int }_{\partial D}{\psi }_0\right)=0.$$

$${\mathcal{S}}_D[y]=-{\eta }_k{\int }_{\partial D}y,$$

$$c\left(1+{\eta }_k{\int }_{\partial D}{\psi }_0\right)=0.$$

$$\mathcal{A}(\omega ,\delta )\left(\begin{array}{c}\varphi \\ \psi \end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),$$

$$u=\{\begin{array}{ll}{\mathcal{S}}_D^k[\psi ](x),& x\in {\mathbb{R}}^2\backslash D,\\ {\mathcal{S}}_D^{k_b}[\varphi ](x),& x\in D.\end{array}$$

$${\hat{\mathcal{S}}}_D^{k_b}[\varphi ]-{\hat{\mathcal{S}}}_D^k[\psi ]={\hat{\mathcal{S}}}_D^k[\varphi -\psi ]+\frac{1}{2\pi }\ln \frac{v}{v_b}{\int }_{\partial D}\varphi ,$$

$$\psi =\varphi +\frac{1}{2\pi }\ln \frac{v}{v_b}\left({\int }_{\partial D}\varphi \right)({\hat{\mathcal{S}}}_D^k{)}^{-1}[{\chi }_{\partial D}]+O({\omega }^2),$$

$${\int }_{\partial D_j}{\mathcal{K}}_{D,1}^{(1)}[\varphi ](x)d\sigma (x)$$

$$=4b_1|D_j|{\int }_{\partial D}\varphi (y)d\sigma (y).$$

Full text

Mathematics subject classification: 35R30, 35C20

Keywords: subwavelength resonance, cochlear mechanics, coupled resonators, hybridisation, passive cochlea, signal processing

1 Introduction

The development of the understanding of the cochlea has largely been a dichotomy between two classes of models [12]. The first, proposed by Hermann von Helmholtz in the 1850s, is based on resonators tuned to different audible frequencies being distributed along the length of the cochlea [19]. Later, Georg von Békésy demonstrated that when the cochlea is stimulated a wave travels from the base to the apex along the basilar membrane [11]. This discovery won him a Nobel Prize in 1961 and lead to the creation of models based on each receptor cell being excited in sequence as the signal travels through the cochlea.

The cochlea is, at its simplest, a long tube filled with fluid, into which sound waves enter through the oval window. An elastic membrane, known as the basilar membrane, is suspended across the centre and upon its surface sit bundles of cylindrical cells, known as hair cells. These bundles of hair cells are the receptor cells of the ear, which produce electrical signals when deflected laterally [21, 20]. The tips of the hair cells are attached to a membrane known as the tectorial membrane and, as a result, motion of the basilar membrane displaces the hair cell and a signal is produced [27]. The cochlea’s ability to filter sounds by pitch is based on the fact that the basilar membrane is graded in both size and stiffness. Thus, cochlear mechanics is, at its heart, a question of studying the motion of a graded elastic membrane in response to a sound wave.

Given what is now know about the structure of the cochlea, and that its function can be reduced to studying membrane motion, one might be inclined to think that Helmholtz’ resonance ideas are no longer relevant. However, the basilar membrane is much stiffer across the cochlea’s width (perpendicular to the page in Figure 1) than along its length. As a result, it has been shown, e.g. by Charles Babbs in [10], that its motion can be modelled by an array of harmonic oscillators. Comparing the cochlea’s length (approximately 3cm) to the wavelength of audible sound (a few centimetres to several metres) and it is clear that these resonators operate in a subwavelength regime.

Babbs shows that if the basilar membrane is split into segments (and considered as an array of oscillators) then the interactions between each part due to elastic tension can be neglected. However, the oscillators will also be coupled by variations in the pressure of the cochlear fluid, by which they are surrounded. The mathematical complexity of modelling these interactions has been one of the main impediments facing the development of this class of cochlear models.

In this paper, we apply boundary integral techniques to understand the complex interactions between coupled subwavelength resonators [4, 9]. In order to simulate the basilar membrane, we consider the problem of acoustic wave scattering by compressible elements in two-dimensional space. Similar layer-potential techniques have previously been applied to other materials that exhibit subwavelength resonance, the classical example being the Minnaert resonance of air bubbles in water [2, 4]. This analysis (in Sections 2.2 & 2.3) relies on the use of layer potential techniques [3, 8, 6].

It is found that a graded array of hybridised resonators has a set of resonant frequencies that becomes increasingly dense (within a finite range) as the number of resonators is increased. We study the eigenmodes and present a scheme (in Section 2.4) for how the model processes incoming signals, filtering them into the system’s resonant frequencies. Finally, in Sections 2.5 & 2.6 we present the important observations that our graded-resonance model predicts the existence of a travelling wave in the pressure field and a basis for the tonotopic map. This acoustic pressure wave is complementary to the wave seen in the motion of the basilar membrane by Békésy and has itself been observed experimentally [26].

2 Response of the coupled resonators

2.1 Preliminaries

We consider a domain $D$ in $\mathbb{R}^{2}$ which is the disjoint union of $N\in\mathbb{N}$ bounded and simply connected subdomains $\{D_{1},\ldots,D_{N}\}$ such that, for each $n=1,\dots,N$, there is $0<s<1$ so that $\partial D_{n}\in C^{1,s}$ (that is, each $\partial D_{n}$ is locally the graph of a differentiable function whose derivatives are Hölder continuous with exponent $s$). We will consider the resonators arranged in a straight line since the curvature of the cochlea does not contribute to its mechanical behaviour [17]. Figure 2 shows an example of such an arrangement (in the special case of circular subdomains, which we will consider for the numerical simulations in Section 2.3 onwards).

We denote by $\rho_{b}$ and $\kappa_{b}$ the density and bulk modulus of the interior of the resonators, respectively, and denote by $\rho$ and $\kappa$ the corresponding parameters for the auditory fluid (which we assume occupies $\mathbb{R}^{2}\setminus\overline{D}$).

We consider an incident acoustic pressure wave $p^{in}(x,t)$ (where $x=(x_{1},x_{2})\in\mathbb{R}^{2}$ and $t\in\mathbb{R}$) that is scattered by $D$. This problem is given by

[TABLE]

where $\frac{\partial}{\partial\nu_{x}}$ denotes the outward normal derivative in $x$ and the subscripts + and - are used to denote evaluation from outside and inside $\partial D$ respectively.

We then introduce the auxiliary parameters

[TABLE]

which are the wave speeds and wavenumbers in $\mathbb{R}^{2}\setminus\overline{D}$ and in $D$ respectively.

We introduce the dimensionless contrast parameters

[TABLE]

By rescaling the dimensions of the physical problem, we may assume that $v=O(1)$ and that the resonators $\{D_{1},\ldots,D_{N}\}$ have widths that are $O(1)$. We further assume that $\tau=O(1)$. On the other hand, since each resonator $D_{n}$ behaves as a harmonic oscillator [15], we can use the calculations of Babbs [10] to show that in order for the system of resonators to replicate the elastic properties of the basilar membrane it must be the case that

[TABLE]

the details of which are given in Appendix A. Since $\tau=\sqrt{\mu/\delta}$, these assumptions give that $\delta\ll 1$. It is important to note that this material contrast condition is an essential prerequisite for the structure $D$ to exhibit resonant behaviours at subwavelength frequencies [3].

We transform problem (1) into the complex frequency domain by making the transformation $u(x,\omega):=\int_{-\infty}^{\infty}p(x,t)e^{i\omega t}\,dt$, $\omega\in\mathbb{C}$ to reach

[TABLE]

‘SRC’ is used to denote the Sommerfeld radiation condition

[TABLE]

The SRC is the condition required to ensure that we select the solution that is outgoing (rather than incoming from infinity) and gives the well-posedness of problem (4).

We wish to use integral operators known as layer potentials to represent the solution to the scattering problem (4).

Definition 2.1**.**

We define the Helmholtz single layer potential associated with the domain $D$ and wavenumber $k$ as

[TABLE]

where $\Gamma^{k}$ is the outgoing (i.e. satisfying the SRC) fundamental solution to the Helmholtz operator $\Delta+k^{2}$ in $\mathbb{R}^{2}$. We similarly define the Neumann-Poincaré operator associated with $D$ and $k$ as

[TABLE]

We can then represent the solution to (4) as

[TABLE]

for some surface potentials $(\phi,\psi)\in L^{2}(\partial D)\times L^{2}(\partial D)$.

We define the space $H^{1}(\partial D):=\{u\in L^{2}(\partial D):\nabla u\in L^{2}(\partial D)\}$ in the usual way and use $Id$ to denote the identity on $L^{2}(\partial D)$. Then, using the representation (8), problem (4) is equivalent [6, 8] to finding $(\phi,\psi)\in L^{2}(\partial D)\times L^{2}(\partial D)$ such that

[TABLE]

where

[TABLE]

We now recall from e.g. [2, 3] the main result that will allow us to understand the leading order behaviour of $\mathcal{A}$ in (9).

Lemma 2.2**.**

In the space $\mathcal{L}(L^{2}(\partial D)\times L^{2}(\partial D),H^{1}(\partial D)\times L^{2}(\partial D))$ we have

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The above operators are defined as

[TABLE]

where $b_{1}:=-\frac{1}{8\pi}$, $c_{1}:=-\frac{1}{8\pi}(\gamma-\ln 2-1-\frac{i\pi}{2})$ and $\gamma=0.5772\ldots$ is the Euler constant.

The operator $\mathcal{S}_{D}$ is the Laplace single layer potential associated with $D$. Since we are working in two dimensions this is not generally invertible however the following two lemmas help us understand the extent of its degeneracy.

Lemma 2.3**.**

If for some $\phi\in L^{2}(\partial D)$ with $\int_{\partial D}\phi=0$ it holds that $\mathcal{S}_{D}[\phi](x)=0$ for all $x\in\partial D$, then $\phi=0$ on $\partial D$.

Proof.

The arguments given in [7, Lemma 2.25] can be easily generalised to the case where $D$ is the disjoint union of a finite number of bounded Lipschitz domains in $\mathbb{R}^{2}$. ∎

Proposition 2.4**.**

Independent of the number $N\in\mathbb{N}$ of connected components making up $D$, we have that

[TABLE]

Proof.

Let $\psi\in\ker\mathcal{S}_{D}$. Then by Lemma 2.3 if $\int_{\partial D}\psi=0$ then $\psi=0$. Suppose that $\int_{\partial D}\psi\neq 0$ then take $\tilde{\psi}\in\ker\mathcal{S}_{D}$ with $\int_{\partial D}\tilde{\psi}\neq 0$ and then consider the function

[TABLE]

Then $f$ satisfies $\mathcal{S}_{D}[f]=0$ and $\int_{\partial D}f=0$ so by Lemma 2.3 we have that $f=0$. Therefore $\psi=(\int_{\partial D}\psi/\int_{\partial D}\tilde{\psi})\tilde{\psi}$. ∎

There are two cases to consider, in light of 2.4:

• •

Case I: $\dim\ker\mathcal{S}_{D}=1$,

• •

Case II: $\dim\ker\mathcal{S}_{D}=0$.

By the Fredholm Alternative Theorem, an equivalent formulation is:

• •

Case I: $\mathcal{S}_{D}$ is not invertible,

• •

Case II: $\mathcal{S}_{D}$ is invertible,

as an operator in $\mathcal{L}(L^{2}(\partial D),H^{1}(\partial D))$. We are now in a position to prove an important property of the operator $\hat{\mathcal{S}}_{D}^{k}$ that was defined in Lemma 2.2 and is the leading order approximation to $\mathcal{S}_{D}^{k}$ as $k\to 0$.

Lemma 2.5**.**

For any fixed $k\in\mathbb{C}\setminus\{z\in\mathbb{C}:\Re(z)=0,\Im(z)\geq 0\}$, $\hat{\mathcal{S}}_{D}^{k}$ is invertible in $\mathcal{L}(L^{2}(\partial D),H^{1}(\partial D))$.

Proof.

Since $\hat{\mathcal{S}}_{D}^{k}$ is Fredholm with index 0 we need only to show that it is injective. To this end, assume that $y\in L^{2}(\partial D)$ is such that

[TABLE]

Case I: Let $\psi_{0}$ be the unique element of $\ker\mathcal{S}_{D}$ with $\int_{\partial D}\psi_{0}=1$ (which exists as a result of Lemma 2.3). We then find that $\mathcal{S}_{D}[y]\perp\psi_{0}$ in $L^{2}(\partial D)$ and hence (11) becomes

[TABLE]

Thus $\int_{\partial D}y=0$. It follows from (11) that $\mathcal{S}_{D}[y]=0$ and further by Lemma 2.3 we have that $y=0$.

Case II: Define $\psi_{0}=\mathcal{S}_{D}^{-1}(1)$ then (11) gives us that

[TABLE]

is constant so, since $\mathcal{S}_{D}$ is injective, we find that $y=c\psi_{0}$ for some $c$. Substituting back into (11) gives

[TABLE]

Everything within the brackets is real with the one exception of $\eta_{k}$ (which has nonzero imaginary part, thanks to the choice of $k$) so we must have that $c=0$. ∎

2.2 Resonant modes

Definition 2.6**.**

For a fixed $\delta$ we define a resonant frequency to be $\omega\in\mathbb{C}$ with positive real part and negative imaginary part such that there exists a nontrivial solution to

[TABLE]

where $\mathcal{A}(\omega,\delta)$ is defined in (10). For each resonant frequency $\omega$ we define the corresponding eigenmode (or resonant mode or normal mode) as

[TABLE]

Remark 2.7**.**

The reason for the choices of sign in 2.6 is to give a physical meaning to a complex resonant frequency. The real part represents the frequency of oscillation and the imaginary part describes the rate of attenuation (hence it should be negative, to give a solution that decays over time).

Remark 2.8**.**

We will see from Figure 5 that 2.6 is equivalent to the notion that resonant frequencies are those at which the system will oscillate at much greater amplitude than is generally the case.

We wish to now compute the resonant frequencies and associated eigenmodes for our system. Manipulating the first entry of (12) we find that

[TABLE]

hence

[TABLE]

since an application of $(\hat{\mathcal{S}}_{D}^{k})^{-1}$ rescales like $O(1/\ln\omega)$. Here, $\chi_{\partial D}$ is used to denote the characteristic function of $\partial D$.

To deal with the second component of (12) we first prove some technical lemmas.

Lemma 2.9**.**

For any $\phi\in L^{2}(\partial D)$ and $j=1,\ldots,N$, we have that

(i) $\int_{\partial D_{j}}(\frac{1}{2}I-\mathcal{K}_{D}^{*})[\phi]=0$,

(ii) $\int_{\partial D_{j}}(\frac{1}{2}I+\mathcal{K}_{D}^{*})[\phi]=\int_{\partial D_{j}}\phi$.

Proof.

(i) follows from the jump relations for single layer potentials and the fact $\mathcal{S}_{D}[\phi]$ is harmonic in $D$ [8, 7]. Then (ii) is immediate. ∎

Lemma 2.10**.**

For any $\phi\in L^{2}(\partial D)$ and $j=1,\ldots,N$, we have that

(i) $\int_{\partial D_{j}}\mathcal{K}_{D,1}^{(1)}[\phi]=4b_{1}|D_{j}|\int_{\partial D}\phi$,

(ii) $\int_{\partial D_{j}}\mathcal{K}_{D,1}^{(2)}[\phi]=-\int_{D_{j}}\mathcal{S}_{D}[\phi]+(4b_{1}+4c_{1})|D_{j}|\int_{\partial D}\phi$,

*where $|D_{i}|$ is the area of $D_{i}$. *

Proof.

(i) follows from the divergence theorem

[TABLE]

Similarly for (ii) we can show that

[TABLE]

making use of the fact that $b_{1}=-1/8\pi$. ∎

Turning now to the second component of (12) we see that

[TABLE]

We substitute expression (14) for $\psi$ to see that $\phi$ satisfies the equation

[TABLE]

At leading order (15) is just $(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})[\phi]=0$ so it would be useful to understand this kernel, which we achieve with the following two lemmas.

Lemma 2.11**.**

If $\phi\in L^{2}(\partial D)$ is such that $\phi\in\ker(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})$ then there exist constants $b_{j}$ such that $\mathcal{S}_{D}[\phi]=\sum_{j=1}^{N}b_{j}\mathcal{X}_{\partial D_{j}}$.

Proof.

Let $u:=\mathcal{S}_{D}[\phi]$. Then $\Delta u=0$ in $D$ and \frac{\partial u}{\partial\nu}\big{|}_{-}=(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})[\phi]=0 on $\partial D$ (known as a “jump condition”) [6, 8] so $u$ satisfies a homogeneous interior Neumann problem on each of the $N$ connected components $D_{1},\dots,D_{N}$ of $D$. It is known that such problems are uniquely solvable up to the addition of a constant. ∎

Lemma 2.12**.**

Fix some ${k_{0}}\in\mathbb{C}\setminus\{0\}$. The set of vectors $\{\psi_{1},\ldots,\psi_{N}\}$ defined as

[TABLE]

forms a basis for the space $\ker(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})$.

Proof.

The linear independence of $\{\psi_{1},\ldots,\psi_{N}\}$ follows from the linearity and injectivity of $\hat{\mathcal{S}}_{D}^{k_{0}}$, plus the independence of $\{\mathcal{X}_{\partial D_{1}},\ldots,\mathcal{X}_{\partial D_{N}}\}$.

For $\phi\in L^{2}(\partial D)$ the difference between $\hat{\mathcal{S}}_{D}^{k_{0}}[\phi](x)$ and $\mathcal{S}_{D}[\phi](x)$ is a constant (in $x$) so they will have the same derivatives. In particular, they are both harmonic and satisfy the same jump conditions across $\partial D$. Therefore, using arguments as in Lemma 2.11, we see that if $\phi\in\ker(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})$ then $\hat{\mathcal{S}}_{D}^{k_{0}}[\phi]\in\text{span}\{\mathcal{X}_{\partial D_{1}},\ldots,\mathcal{X}_{\partial D_{N}}\}$. Thus $\phi\in\text{span}\{\psi_{1},\ldots,\psi_{N}\}$. ∎

From Lemma 2.12 we know that $\ker(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})$ has dimension equal to the number of connected components of $D$ (a wider discussion can be found in e.g. [1]). Thus we can take a basis

[TABLE]

of the null space $\ker(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})$. Then, in light of the fact that at leading order (15) is just $(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})[\phi]=0$, it is natural to seek a solution of the form

[TABLE]

for some non-trivial constants $a_{j}$ with $\sum_{j}|a_{j}|=O(1)$. The solutions $(\phi,\psi)$ to (12) are determined only up to multiplication by a constant (and hence so are $a_{1},\dots,a_{N}$). We fix the scaling to be such that the eigenmodes are normalised in the $L^{2}(D)$-norm

[TABLE]

We now integrate (15) over each $\partial D_{i},i=1\ldots N$ and use the results of Lemmas 2.9 and 2.10 to find that, up to an error of $O(\delta\omega^{2}\ln\omega)+O(\omega^{4}\ln\omega)$,

[TABLE]

When we substitute the expression (17) for $\phi$ in (19) we find the system of equations, up to an error of order $O(\delta\omega^{2}\ln\omega)+O(\omega^{4}\ln\omega)$,

[TABLE]

Remark 2.13**.**

Thanks to the linearity of the operators $B_{\delta}^{(i)}$, the solutions $\omega(\delta)$ to (20) (as well as the associated eigenmodes) are independent of the choice of basis $\{\phi_{1},\dots,\phi_{N}\}$.

Remark 2.14**.**

One can think of the step where we integrated (15) over each $\partial D_{i}$, for $i=1\dots,N$, to give (19) as the point where the hybridisation (between the $N$ resonators) was performed (see also e.g. [4]).

2.3 Numerical computations of resonant modes

In order to improve computational efficiency, we will assume from here onward that the resonators are circular. This means that we can use the so-called multipole expansion method, an explanation of which is provided in e.g. [5, Appendix C]. The method relies on the idea that functions in $L^{2}(\partial D)$ are, on each circular $\partial D_{i}$, $2\pi$-periodic so we may approximate by the leading order terms of a Fourier series representation. We found that as few as seven terms was sufficient to give satisfactory results.

Using such an approach we can find, for each fixed $\delta>0$, the $N$ values of $\omega\in\mathbb{C}$ such that there exists a nontrivial solution to (20). For the case where $N=50$ the results are shown in Figure 3.

We see that there is a range of frequencies where the (the real part of the) resonances occur most commonly. As $N$ is increased, the resonances become increasingly dense in this region. In fact, with the current arrangement, this range of frequencies does not change as $N$ increases. Instead, the region becomes increasingly densely filled.

It is also seen from Figure 3 that the imaginary parts of the resonances is smallest in the region where they are most dense. This means that these frequencies experience the least significant attenuation, suggesting that tones in this range will be most easily audible. The reason $\omega_{1}=0.0002284-0.0000526i$ has been omitted from Figure 3 is due to its $O(10^{-4})$ imaginary part. This is not only inconvenient for plotting but also means that this resonant mode will suffer much greater attenuation and thus will be a less significant part of the motion over time.

It is also important to understand the eigenmodes $u_{n}$ associated with each resonant frequency $\omega_{n}$. The six resonant modes for the case of six resonators are shown in Figure 4. They take the form of increasingly oscillating patterns that inherit the asymmetry of the resonator array.

It is also notable that the solution is approximately constant on each resonator. This is because the solution, taking the form (8), is given by $\hat{\mathcal{S}}_{D}^{k_{b}}[\phi]$ at leading order which by Lemma 2.11 is constant for $\phi\in\ker(-\frac{1}{2}Id+\mathcal{K}_{D}^{*})$.

2.4 Signal processing

We wish to offer an explanation of how, given an incident wave $p^{in}(x,t)$, our system of coupled resonators is able to classify (and hence identify) the sound. The system of resonators $D$ is able to decompose the signal over its resonant modes. It is clear that the $N$ eigenmodes are linearly independent so we may define the relevant $N$-dimensional solution spaces.

Definition 2.15**.**

We define the $N$-dimensional spaces $X$ and $Y$ as

[TABLE]

[TABLE]

We will approximate the solution by a decomposition in the frequency domain. The fact that, for $n=1,\dots,N$, the Fourier transform of $e^{-i\omega_{n}t}$ for $t>0$ is given by $i/(\omega-\omega_{n})$ motivates us to employ the ansatz

[TABLE]

where $\alpha_{1},\dots,\alpha_{N}$ are complex-valued functions of a real variable.

It is important to understand whether knowing the value of the solution on each resonator (which is the information that a cochlea is able to capture) means that one can recover the weight functions $\alpha_{1},\dots,\alpha_{N}$ in (23).

Remark 2.16**.**

The eigenmodes $u_{1},\dots,u_{N}$ are not orthogonal in $L^{2}(D)$. It turns out, however, that they are nearly orthogonal. For example, the normalised eigenmodes shown in Figure 4 satisfy $(u_{n},u_{m})_{L^{2}(D)}=O(10^{-3})$ for $n\neq m$.

Proposition 2.17**.**

Let $\{\omega_{1},\ldots,\omega_{N}\}$ be the resonances of the system $D=D_{1}\cup\ldots\cup D_{N}$ and denote by $u_{1},\ldots,u_{N}$ the corresponding eigenmodes. Then the matrix $\gamma\in\mathbb{C}^{N\times N}$ defined by

[TABLE]

is invertible.

Proof.

We can apply the Gram-Schmidt procedure to produce a basis $\{v_{1},\ldots,v_{N}\}$ for $X$ that is orthonormal with respect to $(\cdot,\cdot)_{L^{2}(D)}$. This procedure produces a nonsingular lower triangular matrix $P\in\mathbb{C}^{N\times N}$ such that $(v_{1},\dots,v_{N})^{T}=P(u_{1},\dots,u_{N})^{T}$ (superscript $T$ denotes the matrix transpose). If we define $Q\in\mathbb{C}^{N\times N}$ as $Q:=P^{-1}$ then $Q$ is also nonsingular and lower triangular. We can then calculate that

[TABLE]

Integrating (25) componentwise gives that, for $i,j=1,\dots,N$, it holds that

[TABLE]

and thus

[TABLE]

∎

In order to find the weight functions $\alpha_{1},\dots,\alpha_{N}$ in Equation 23 we must take the $L^{2}(D)$-product with $u_{n}(x)$ for $n=1,\dots,N$ and then invert $\gamma$. This gives that

[TABLE]

Thanks to its representation (8) in terms of single layer potentials, $u(\cdot,\omega)$ is an analytic function of $\omega\in\mathbb{C}$. Thus, from (28) we can see that $\alpha_{1},\dots,\alpha_{N}$ are analytic and hence we can recover a similar decomposition for $p(x,t)$ using the Laplace inversion theorem

[TABLE]

Example 2.18**.**

$p^{in}(x,t)$* is a plane wave*

We take as an example the case where $p^{in}(x,t)$ is a pulse of a plane wave with frequency $\omega_{in}\in\mathbb{R}$ travelling in the $x_{1}$ direction. This is given by

[TABLE]

This has Fourier transform

[TABLE]

We can then compute $\alpha_{1}(\omega),\dots,\alpha_{N}(\omega)$ as in (28).

In Figure 5 we show firstly how the $L^{2}(D)$-norm of the solution to the scattering problem (4) varies as a function of $\omega_{in}$. As is expected, the response is (locally) much greater when $\omega_{in}$ is close to $\Re(\omega_{n})$ for some $n=1,\dots,N$. We also show how the weights $\alpha_{1}(\omega_{1}),\dots,\alpha_{N}(\omega_{N})$ in (29) vary as a function of $\omega_{in}$. Each constant is small except in a region of the associated resonant frequency when the corresponding eigenmode is excited most strongly.

It should also be noted that $|\alpha_{n}(\omega_{n})|$ decreases in $n$. If we considered higher order resonances the corresponding constants would be significantly smaller. This justifies our choice to approximate $p$ as an element of $Y$ in (29) (i.e. to only consider the $N$ subwavelength modes).

2.5 Travelling waves

In trying to resolve the differences between the two main classes of cochlear model a crucial realisation is that our (resonance) model for the acoustic pressure exhibits the travelling wave behaviour. This is easy to see in models based on graded arrays of uncoupled resonators, since a resonator’s response time increases with decreasing characteristic frequency [18, 12], but is also true of our hybridised model.

Simultaneously exciting a graded array that is initially at rest produces the evolution shown in Figure 6. The existence of a wave travelling from the small high-frequency resonators at the base of the cochlea to the larger low-frequency resonators at the apex is clear. This wave is the movement of the position of maximum acoustic pressure along the array of resonators. It is a consequence of the asymmetric eigenmodes (shown in Figures 7a-d) growing from rest at different rates.

The parallels between the travelling wave in Figure 6 and that observed (e.g by Békésy) on the basilar membrane are clear. While it is true that acoustic waves enter the cochlea at the base and travel through the fluid to the apex, the wave observed by Békésy moves much more slowly than this. The speed of sound in cochlear fluid is approximately 1500$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ whereas the travelling wave is observed at speeds close to 10$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$[12, 16]. This justifies the choice to assume that all the resonators are excited simultaneously by an incoming signal [12]. Since pressure changes in the fluid and the motion of the membrane will be physically linked, it is not surprising that the wave in Figure 6 shares a number of characteristics with Békésy’s observations. For instance, the amplitude initially grows before quickly diminishing and the wave is seen to slow as it moves through the array [18, 11, 16]. It should also be noted that travelling waves have been observed in the cochlear fluid (as predicted here) as well as on the membrane [26].

2.6 Tonotopic map

Békésy’s famous experiments further revealed the existence of a relationship between signal frequency and the position in the cochlea where the sound is most strongly detected. His results showed that the frequency $f(x)$ giving rise to maximum excitation at a distance $x$ from the base of the cochlea satisfies a tonotopic map of the form

[TABLE]

for some $a,d,c\in\mathbb{R}$ [11]. In Figure 7 we show the relationship between the position of maximum amplitude of each eigenmode and the associated resonant frequency. We see that, if some of the lowest frequency modes are ignored, the pattern follows a relationship that is approximately of the form (32) (with $a=0.0126,d=-0.0117,c=0.0060$). The eigenmodes shown in Figures 7b-d demonstrate the basis for the tonotopic map. Each features oscillations with a clear peak followed by a rapid decrease in amplitude (which explains the growth and then rapid decay of the travelling wave that was observed in Figure 6).

It is not clear why the lowest frequency modes (e.g. Figure 7a) do not fit the pattern that is established by the majority of the eigenmodes, or what the implications of this could be. However the relatively large negative imaginary parts of the associated resonant frequencies mean this phenomenon has a less significant impact on the evolution of the acoustic pressure field.

3 Concluding remarks

In this paper, we have modelled the passive cochlea by approximating its response by a graded array of subwavelength resonators. We have used boundary integral methods to compute leading order approximations to the resonant frequencies and associated eigenmodes this fully-coupled system. This model has the ability to decompose incoming signals into these resonant modes. As the number of resonators is increased, the resonant frequencies densely fill a finite range meaning that a large system can capture signals with a high frequency resolution.

It is a significant observation that a simple graded-resonance model, with appropriate coupling between the resonators, predicts travelling wave behaviour in the acoustic pressure field, and that this has similar characteristics to that observed in the membrane motion. In some sense, these models represent the unification of Helmholtz’ and Békésy’s ideas [12, 13].

It is well known that the cochlea is an active organ and even emits sounds (known as otoacoustic emissions) as part of its response to a signal [22, 23, 14, 25, 28]. For instance, a key feature that our current model lacks is the ability to amplify quiet sounds more greatly than louder ones. Such non-linear amplification is needed in order to account for the ear’s remarkable ability to hear sounds over a large range of amplitudes. In this work, we have only considered a passive system of resonators but have presented a model which, by introducing appropriate non-linear forcing terms in (1), can be modified to include active elements in future work.

The code developed for this study is available online at

https://github.com/davies-b/cochlea_passive

Acknowledgements

The authors are grateful to Fabrice Lemoult for their participation in valuable discussions and to Andrew Bell for insightful comments made on an early version of this manuscript.

Appendix A Material parameters

Here, we estimate the appropriate material parameters for the system of subwavelength resonators used in our version of the model from [10]. In particular, using physical values for the cochlea and representations of subwavelength acoustic resonators as harmonic oscillators, we estimate the appropriate value for the bulk modulus contrast $\mu$, defined in (2).

In [10] it is shown that a suitable approximation to basilar membrane motion can be achieved by considering an array of masses on springs. For a piece of membrane with area $A$, thickness $h$, width $w$ and Young’s modulus $E$, the spring constant of the equivalent resonator is shown to be given by

[TABLE]

where $C\approx 30$ is a dimensionless constant.

In [15] it is shown that a subwavelength acoustic resonator behaves as a one-degree-of-freedom harmonic oscillator. In the case where the resonator $D_{i}$ is spherical with radius $R$, it is shown that its stiffness is given by

[TABLE]

Combining (33) and (34), we see that the contrast $\mu$ is given by

[TABLE]

In Table 1, we give values for the relevant material parameters, derived by experimentalists working on biological cochleas. Using the orders of magnitude of these values we find that, if the membrane is approximated by $N=O(10^{2})$ subwavelength resonators, then it should hold that

[TABLE]