Explosive synchronization enhances selectivity: Example of the cochlea

Abstract

Acoustical signal transduction in the cochlea is an active process that involves nonlinear amplification and spontaneous otoacoustic emissions. Signal transduction involves individual subunits composed of globally coupled hair cells, which can be modeled as oscillators close to a Hopf bifurcation. The coupling may induce a transition toward synchronization, which in turn leads to a strong nonlinear response. In the model studied here, the synchronization transition of the subunit is discontinuous (explosive) in the absence of an external stimulus. We show that, in the presence of an external stimulus and for a coupling strength slightly lower than the critical value leading to explosive synchronization, the response of the subunit has better frequency selectivity and a larger signal-to-noise ratio. From physiological observations that subunits are themselves coupled together, we further propose a model of the complete cochlea, accounting for the ensemble of frequencies that the organ is able to detect.

Graphical abstract

Keywords

cochlea / frequency selectivity / periodical forcing / explosive synchronization

1 Introduction

Sound transduction, the process whereby information from incoming sound waves in the ear is processed and transmitted to the brain, occurs in the cochlea. This snail-shaped organ, located in the internal ear, i.e., the innermost part of the vertebrate ear [1], is wound into a coil with approximately two and a half turns. The cochlea is characterized by a remarkable frequency-selective sensitivity to mechanical stimulation; in humans, it can detect sound intensities over 6 orders of magnitude with a frequency selectivity as high as \text{Δ}\omega /\omega \approx 0.2\% [2]. The central part of the cochlea is the organ of Corti, which consists of three rows of outer hair cells and a single row of inner hair cells [3]; the former are responsible for enhancing the basilar membrane motion, and the latter translates sound-induced mechanical stimuli to the auditory nerve [4-6]. In contrast to passive mechanical systems, the cochlea can actively generate weak acoustic oscillations in a broad range of frequencies [7]; consequently, acoustical signal detection is an active process [8].

Understanding how the transduction system can achieve such remarkable performance is the theoretical challenge addressed here. We aim to elucidate the extreme sensitivity of the cochlea to incoming sound waves. Earlier models have shown that the generation of spontaneous otoacoustic emissions, as well as the self-tuning properties of the cochlea, can be understood in terms of nonlinear oscillators close to a Hopf bifurcation [9-11]. Moreover, several physiological models have been reduced exactly to a quantitative Hopf normal form [12]. At the same time, spatially extended models have been proposed, either involving the interaction of an ensemble of uncoupled oscillators with a deformation wave of the basilar membrane [13, 14, 16, 17] or more parsimoniously introducing viscoelastic coupling between adjacent oscillators [18-21]. One of us recently proposed that mammalian hearing can be enhanced by a balance between spontaneous otoacoustic emissions and spatial coupling [22].

To better understand the remarkable frequency sensitivity of the cochlea, we consider here the hierarchy of levels in the cochlea. As depicted in Fig. 1, hair cells are coupled into subunits via the sallets and further coupled to the basilar membrane. We consider here the collective behavior of oscillators coupled hierarchically, as shown schematically in Fig. 1. In this scenario, a hair cell or hair bundle is modeled as a single oscillator that is coupled viscously and elastically to its neighboring hair cells underneath the tectorial membrane (the sallet; see Fig. 1) [23]. Such a set of coupled cells defines a subunit. We further couple each hair cell to the neighboring subunits to account for the connection to the basilar membrane. For example, the basal segment of the bobtail lizard papilla reveals a tectorial membrane with an ensemble of about 80 sallets, where each sallet covers about 10 to 25 hair cells [24].

We therefore consider successively two levels of description. First, we focus on a subunit consisting of N globally coupled oscillators representing the hair cells coupled by the tectorial membrane or a sallet. At this intermediate level, and in the absence of an external signal, we observe that the subunit can undergo an explosive synchronization transition when the coupling strength is varied. We also study the response of the subunit for some fixed coupling strength slightly below the critical value by applying an external forcing. Second, we study a higher-level model of the cochlea in which hair cells are coupled to neighboring sallets.

A key factor in our model is the complex coupling strength, which includes both real and imaginary parts representing both viscous and elastic coupling separately [19], in contrast to previous models with only a real coupling strength [20, 22]. In this framework, the real and imaginary parts of the coupling strength reflect the interactions of the amplitude and phase, respectively. We find, surprisingly, that the competition between the real and imaginary parts of the coupling may lead to an explosive synchronization transition, which is currently a very hot topic [25-34]. Moreover, we find that this abrupt synchronization transition is very good for the implementation of sensitive frequency-selectivity, thus providing a new mechanism for signal amplification and selectivity in the cochlea.

The paper is organized as follows. In Section 2, we present a model of an isolated subunit consisting of N globally coupled cells and show that it can undergo explosive synchronization. In Section 3, we apply an external stimulus to the subunit, and we observe that synchronization is concomitant with frequency locking to the external stimulus. In Section 4, we present a numerical exploration of a complete model of coupled subunits. Finally, discussions and conclusions are presented in Section 5.

2 Isolated subunit

The hair bundles in the cochlea can be modeled as oscillators undergoing a Hopf bifurcation [9]. On the level of an individual subunit, hair bundles are both elastically and viscously coupled [19-21]. Considering that hair cells in one subunit are very similar and the coupling among them is homogeneous, we may treat them as all-to-all (globally) coupled and consider a mean field coupling. The following equation for the reduced variable z_j describing the displacement and velocity of the tip of the j-th hair bundle has been derived [19]:

\begin{array}{}\text{(1)}& {\dot{z}}_j=({\mu }_j+\text{i}{\omega }_j)z_j-B_j|z_j{|}^2{z}_j+k_j(d_{\text{R}}+\text{i}{d}_I)(\overline{z}-z_j)\mathrm{,}\end{array}

where N represents the number of oscillators; j=\mathrm{1,2,}\cdots \mathrm{,}N; {\omega }_j is the natural angular frequency of the j-th oscillator; {\mu }_j is a parameter representing the nondimensional distance from the Hopf bifurcation onset; B_j describes the intrinsic nonlinearity of the oscillator; d_{\text{I}} and d_{\text{R}} are proportional to the elastic and viscous coupling constant [19], respectively; k_j describes the relative coupling strength of the j-th oscillator; and \overline{z}=\frac{1}{N}{\displaystyle {\sum }_{j=1}^N}{z}_j is the centroid, defined as the average of z_j over all the oscillators. For an isolated oscillator (k_j=0), the stationary solution z_j=0 is stable for ; in contrast, when {\mu }_j>0, z_j oscillates in a limit cycle corresponding to a spontaneous otoacoustic emission [7]. Unless stated otherwise, we use N=101 oscillators throughout this paper. We have checked that the results do not qualitatively depend on N unless N is too small.

A key feature of Eq. (1) is that its coupling consists of both a real part d_{\text{R}} and an imaginary part d_{\text{I}}. To study the collective behavior described by Eq. (1), we introduce the real, positive order parameter R=|\overline{z}|, together with the phase ψ of the centroid, as R{\text{e}}^{\text{i}\psi }=\overline{z} [35]. A high (low) value of R implies that the oscillators are well (poorly) synchronized. Moreover, R represents the amplitude of the oscillation of the centroid, i.e., for the synchronized state, it gives the amplitude of the oscillation.

For the sake of simplicity, all the oscillators are assumed to be described by the same parameters, {\mu }_j=\mu, B_j=B, and k_j=k. Oscillators in the group differ only in their natural angular frequencies \left\{{\omega }_j\right\}\mathrm{,1}\le j\le N, for which we choose a symmetrical distribution around a central frequency {\omega }_0. Because Eq. (1) is invariant under any phase shift, including z_j\mapsto z_j{\text{e}}^{-\text{i}{\omega }_{0}t}, we further choose {\omega }_0=0 without loss of generality. The results reported in this article were obtained by choosing a uniform distribution of width 2\pi \text{Δ} with \text{Δ}=0.1 Hz, or, when further indicated, a deterministic set of angular frequencies equispaced by 2\pi \text{Δ}/(N-1) in the interval -\pi \text{Δ} and +\pi \text{Δ}, which we call a linear distribution. In the rest of this article, we choose \mu =\mathrm{1,}B=1, and d_{\text{R}}=0.15, as in Refs. [19, 21]. Eq. (1) is numerically integrated using the fourth-order Runge–Kutta algorithm. We typically allowed a relaxation time of 40 000 time units to avoid transient behaviors. We define the time average, {\langle X\rangle }_T, of any quantity X as

\begin{array}{}\text{(1)}& {\langle X\rangle }_T=\frac{1}{T}{\displaystyle {\int }_t^{t+T}}X(\tau )\text{d}\tau \mathrm{,}\end{array}

and we choose an integration time T=20000; this is sufficiently large compared to the largest period in the system, which we approximate as the smallest angular frequency in the distribution of the natural angular frequencies \left\{{\omega }_j\right\}, i.e., T\gg 2(N-1)/\text{Δ}.

As the coupling k is increased, we expect the oscillators to synchronize; i.e., we expect R to take large values. We approximate R by its time-averaged value {\langle |\overline{z}|\rangle }_T. The initial conditions are chosen such that all the oscillators have an initial radius \sqrt{\mu }=1 and initial phases randomly distributed in the interval \left[\mathrm{0,}\frac{\pi }{2}\right]. Figure 2 shows the dependence of R on the coupling strength k for several values of the imaginary part of the coupling, d_{\text{I}}, which reveals a change in the order of the transition. The results shown in the main part of Fig. 2 were obtained by starting a simulation of the model (1) with random initial conditions for each set of values of d_{\text{I}} and k. For small |d_{\text{I}}|, the transition looks continuous, whereas large values of |d_{\text{I}}| — corresponding to a large elastic coupling — yield a discontinuous transition to synchronization. Note that the sign of d_{\text{I}} does not matter here, as Eq. (1) can be complex conjugate provided the distribution of frequencies \left\{{\omega }_j\right\} is symmetrical around {\omega }_0=0. The way that the system switches from a second- to a first-order transition when the parameters are varied is the subject of a separate publication [36]. In Ref. [36], the discontinuous transition is shown to result from the coexistence of a particular synchronized state together with the trivial incoherent state in a range of parameters (k and d_{\text{I}}) in which the latter is linearly stable. In contrast to the paradigmatic Kuramoto model [37], this particular state observed at the synchronization transition contains a finite, nonvanishing number of synchronized oscillators; hence, it results in a first-order transition. In the following, we choose d_{\text{I}}=-1.2, which corresponds to a first-order transition to synchronization, also called explosive synchronization [25-34]. Adiabatically increasing the coupling k by small steps, \delta k=0.01, and adiabatically decreasing k by the same steps, we find a hysteresis region where bistability is observed (inset of Fig. 2). In this region, the system can be either synchronized or nonsynchronized depending on the initial conditions.

For each oscillator, we write z_j=r_j{\text{e}}^{\text{i}{\theta }_j} and introduce the effective angular frequency of the oscillator as the time-averaged time derivative of its phase {\theta }_j [25]:

\begin{array}{}\text{(2)}& {\omega }_j^{\text{e}\text{f}\text{f}}={\langle {\dot{\theta }}_j\rangle }_T=\frac{{\theta }_j(t+T)-{\theta }_j(t)}{T}\mathrm{.}\end{array}

The distribution of {\omega }_j^{\text{e}\text{f}\text{f}} can be characterized by its mean, \overline{{\omega }^{\text{e}\text{f}\text{f}}}=\frac{1}{N}{\displaystyle {\sum }_{j=1}^N}{\omega }_j^{\text{e}\text{f}\text{f}}, and its standard deviation σ. Figure 3 shows the evolution of the average and standard deviation of the set of effective frequencies \left\{{\omega }_j^{\text{e}\text{f}\text{f}}\right\} as a function of the decreasing coupling k. For large coupling, k>k_c^{″}\approx 0.71, the standard deviation of the effective frequency vanishes, and R is large, which implies that all the oscillators are fully frequency-synchronized. For small coupling, , there is no synchronization, and the oscillators all have different effective frequencies, which are nothing but their natural angular frequencies shifted by the imaginary part of the coupling: {\omega }_j^{\text{e}\text{f}\text{f}}={\omega }_j-k{d}_{\text{I}}. Thus, the standard deviation has the same constant value, 2\pi \text{Δ}/(2\sqrt{3}), as the distribution of natural angular frequencies, and the average evolves linearly with the coupling strength k.

We also computed the angular frequency \Omega of the centroid, defined as the time derivative of the phase ψ of the centroid, by using the same averaging procedure as in Eq. (2). The phase of the order parameter is a nonlinear function of the phases of the individual oscillators, so \Omega behaves differently from \overline{{\omega }^{\text{e}\text{f}\text{f}}}, unless the oscillators are fully frequency-synchronized. As seen in Fig. 3, for large coupling, k>k_c^{″}, \Omega is equal to the average effective angular frequency. For intermediate values of the coupling strength k, the order of magnitude of \Omega is unchanged at around {10}^{-1}.

Let us now describe qualitatively the existence of the bistable region in Fig. 2. We begin by rewriting Eq. (1) for each oscillator, z_j, in terms of its amplitude, r_j, and phase, {\theta }_j:

\begin{array}{}\text{(3)}& \begin{array}{l}{\dot{r}}_j=(\overline{u}-r_j^2)r_j+\overline{k}R\cos (\psi -{\theta }_j+\alpha )\mathrm{,}\\ {\dot{\theta }}_j={\omega }_j-\overline{k}\sin \alpha +\overline{k}\frac{R}{r_j}\sin (\psi -{\theta }_j+\alpha )\end{array}\end{array}

for j=\mathrm{1,2,}\cdots \mathrm{,}N, where we have introduced the notation \overline{k}=k\sqrt{d_{\text{R}}^2+d_{\text{I}}^2}, \alpha =\arctan (d_{\text{I}}/d_{\text{R}}), and \overline{\mu }=\mu -\overline{k}\cos \alpha [recall that \overline{z}=R\exp (\text{i}\psi )].

The difference between the instantaneous angular frequencies of two distinct oscillators i and j thus reads

\begin{array}{}\text{(4)}& {\dot{\theta }}_i-{\dot{\theta }}_j=({\omega }_i-{\omega }_j)-\overline{k}R\left[\frac{\sin (\psi -{\theta }_j+\alpha )}{r_j}-\frac{\sin (\psi -{\theta }_i+\alpha )}{r_i}\right]\text{}\mathrm{,}\end{array}

which allows us to discuss the two regimes.

First, we can easily see that R=0 will give a trivial incoherent state of the system. In such a state, each oscillator rotates freely on a circle of radius \overline{\mu } with its natural angular frequency shifted by k{d}_{\text{I}}, so not a single pair of oscillators will synchronize. For small coupling strength, k\ll 1, the second term on the right-hand side of Eq. (4) can be neglected, and {\dot{\theta }}_i-{\dot{\theta }}_j\approx {\omega }_i-{\omega }_j\ne 0, so the phase difference {\theta }_i-{\theta }_j increases with time for any set of initial conditions. This is the case for nearly all (i\mathrm{,}j), which implies that R is close to zero in the steady-state regime.

Moreover, if one then increases the coupling strength k but still sets the initial state totally randomly to ensure that R is close to zero, from Eq. (4), we can see that it has the same effect as a small k, which means that the second term can still be neglected initially. As long as the incoherent state is linearly stable under this k value, we can be sure that R does not increase significantly while Eq. (1) evolves, and hence we have a nearly zero R. We define the linear stability boundary as k_c^{\prime } [36].

However, for the same k value, if we set the initial state in such a way that R is quite large, we can no longer neglect the second term, which can be large enough to balance the frequency difference {\omega }_i-{\omega }_j, in order to ensure a constant phase difference {\theta }_i-{\theta }_j for some pairs of oscillators that have similar natural frequencies. This synchronization state can exist as long as there is a large enough synchronized cluster for such k. This can give another critical coupling strength, k_c [36].

In cases such as the classical Kuramoto model [37], k_c>k_c^{\prime }, so a k value slightly larger than k_c^{\prime } will first cause a small synchronization cluster, and then the cluster size will increase up to the system size when k reaches k_c. In this case, the transition to synchronization is continuous. On the other hand, leads to a bistable regime. In this case, when k is decreased around k_c, the system will suddenly jump from a high R state down to the steady incoherent state with R\approx 0, giving a discontinuous or first-order transition, or explosive synchronization. Previous theoretical results for a similar system [36] showed that k_c can be smaller than k_c^{\prime } if |d_{\text{I}}/d_{\text{R}}| is large enough, as seen in Fig. 2.

3 Forced subunit

An interesting corollary of the discontinuous transition is the remarkable frequency selectivity of the cochlea. To address this problem, we now consider a slightly subcritical regime by choosing a coupling strength slightly below the critical value [16], , and study the effect of external monochromatic forcing with amplitude f and frequency {\omega }_{\text{f}}. This forcing, f{\text{e}}^{\text{i}{\omega }_{\text{f}}t}, which represents the effect of an incoming sound wave, is assumed to act additively and with the same strength on all the oscillators:

\begin{array}{}\text{(5)}& {\dot{z}}_j=(\mu +\text{i}{\omega }_j)z_j-B|z_j{|}^2{z}_j+k(d_{\text{R}}+\text{i}{d}_{\text{I}})(\overline{z}-z_j)+f{\text{e}}^{\text{i}{\omega }_{\text{f}}t}\mathrm{.}\end{array}

Our key observation is that synchronization can be achieved if , provided that the angular frequency of the external stimulus lies within some range, which depends on its amplitude. For fixed values of the external stimulus amplitude f and coupling coefficient k, varying the angular frequency {\omega }_{\text{f}} reveals two distinct solutions with either a large or a small value of the order parameter R. Figure 4 shows the boundary delimiting the region of ({\omega }_{\text{f}}\mathrm{,}f) where a large R solution is observed. This region shrinks when k is increased below k_c.

For a given value of the coupling strength k, the oscillators are synchronized inside the tongue in Fig. 4 for a sufficiently large forcing amplitude, whereas they are not synchronized outside this tongue for a lower forcing amplitude. Thus, the minimum forcing f_{\text{m}\text{i}\text{n}} at the tip of the tongue represents the sensitivity to a weak signal. We find that the sensitivity depends on the coupling strength k and that the system can detect weaker signals when k is closer to the critical value k_c; i.e., f_{\text{m}\text{i}\text{n}} decreases as the coupling strength k increases toward k_c. Focusing on the dependence of f_{\text{m}\text{i}\text{n}} on k, we observe that when the coupling strength k is increased toward k_c, the minimum amplitude f_{\text{m}\text{i}\text{n}} required to synchronize the oscillators decreases linearly toward 0 for k=k_c, as shown in the inset of Fig. 4. This property can be used to obtain a more accurate estimate of k_c than that obtained from Fig. 2 by performing a linear fit of f_{\text{m}\text{i}\text{n}} as a function of k.

Is the synchronized solution phase-locked with the external signal? To answer this question, we measure the amplitude Q of the Fourier mode of \overline{z} at the forcing frequency {\omega }_{\text{f}} as follows:

\begin{array}{}\text{(6)}& Q=|\frac{1}{T}{\displaystyle {\int }_t^{t+T}}\overline{z}(t){\text{e}}^{-\text{i}{\omega }_{\text{f}}t}\text{d}t|\mathrm{.}\end{array}

Figure 5 shows that the dependences of Q and R on the forcing frequency are very similar. We further checked that the ratio Q^2/\langle R^2(t)\rangle, which measures the fraction of energy located at the forcing frequency, always equals 1 in the synchronized region within a 1\% error. Additionally, we can check whether a single oscillator is frequency-locked with the external forcing. We define the fraction {\rho }^{\text{l}\text{o}\text{c}\text{k}}=N^{\text{l}\text{o}\text{c}\text{k}}/N, where N^{\text{l}\text{o}\text{c}\text{k}} is the number of oscillators that is frequency-locked with the external stimulus. As seen in Fig. 5, this fraction is close to zero in the nonsynchronized state, whereas it is large in the synchronized solution. Alternatively, we can use {\rho }^{\text{l}\text{o}\text{c}\text{k}} to define the synchronized region in Fig. 4, which gives the same boundary as that obtained by using R or Q. Thus, the synchronized solution observed inside the tongue in Fig. 4 is frequency-locked to the forcing stimulus. We will further refer to the synchronized region as an Arnol’d tongue [15], by analogy with the case of a single forced nonlinear oscillator. We should emphasize that the phase-locked area of the tongue in Fig. 4 is induced by the external driving. We denote a point on the boundary line in Fig. 4 as ({\omega }_{\text{f}}^{\text{b}}\mathrm{,}{f}^{\text{b}}). It has two states, synchronized and nonsynchronized states. When f>f^{\text{b}}, the synchronized state is stable, whereas the nonsynchronized state is unstable. When , the nonsynchronized state is stable, whereas the synchronized state is unstable. This is how signal detection is guaranteed. On the other hand, from Fig. 5, we see that even in the area of the tongue, the fraction of synchronized oscillators is about 75% but not 100\%, implying that 25\% of the oscillators are still in the nonsynchronized state. In this sense, we can understand the signal frequency selectivity.

Another observation is that the frequency {\omega }_{\text{m}\text{i}\text{n}} at the tip of the tongue, which is the frequency that is the easiest to force, is in the range [8-9]\times {10}^{-2}. This is much larger than the smallest natural angular frequency in the distribution \left\{{\omega }_j\right\}, which is 6.3\times {10}^{-3} for N=101. In fact, {\omega }_{\text{m}\text{i}\text{n}} is on the order of the angular frequency \Omega that would be obtained in the absence of external forcing.

We now ask how the number of oscillators in the system, N, influences our results. Figure 6(a) shows that in the synchronized state, the order parameter R_{\text{s}} shows little variation and depends linearly on 1/N. In contrast, in the nonsynchronized state, R_{\text{n}\text{s}}, as well as the minimum amplitude f_{\text{m}\text{i}\text{n}}, decreases significantly in proportion to 1/N as N increases; see Figs. 6(b) and (c). Thus, for a realistic subunit in the cochlea system, which is typically composed of dozens to hundreds of hair cells, finite size effects may be relevant.

From Eq. (5), we can derive the dynamical equation of the mean field \overline{z}(t):

\begin{array}{}\text{(7)}& \dot{\overline{z}}=(\mu +\text{i}\overline{\Omega })\overline{z}-\frac{1}{N}{\displaystyle \sum _{j=1}^N}|z_j{|}^2{z}_j+f{\text{e}}^{\text{i}{\omega }_{\text{f}}}\mathrm{,}\end{array}

where we have written \overline{\Omega }\equiv \frac{1}{N\overline{z}}{\displaystyle {\sum }_{j=1}^N}{\omega }_j{z}_j. In the synchronized state, oscillators are populated around the mean field; hence, the amplitude |z_j|\approx R. We can write |z_j{|}^2=R^2+{\delta }_j with {\delta }_j\ll 1, so we have \frac{1}{N}{\displaystyle {\sum }_{j=1}^N}|z_j{|}^2{z}_j=R^2\overline{z}+\xi \overline{z} with \xi \overline{z}\equiv \frac{1}{N}{\displaystyle {\sum }_{j=1}^N}{\delta }_j{z}_j and \xi \ll 1. Because \left\{{\omega }_j\right\} is symmetrical around 0, \overline{\Omega } has a small modulus and can be approximated to the leading order by the small mean field frequency \Omega. Eq. (7) then reads

\begin{array}{}\text{(8)}& \dot{\overline{z}}=(\mu +\text{i}\Omega -R^2)\overline{z}-\xi \overline{z}+f{\text{e}}^{\text{i}{\omega }_{\text{f}}t}\mathrm{.}\end{array}

To find a 1:1 frequency-locked solution [22, 39], we further impose the condition that the mean field frequency \Omega equals the forcing frequency {\omega }_{\text{f}} up to a slowly varying phase {\psi }^{\prime }, \Omega t={\omega }_{\text{f}}t+{\psi }^{\prime }, so Eq. (8) can be further rewritten as

\begin{array}{}\text{(9)}& \begin{array}{l}\dot{R}=R(\mu -R^2)-\xi R+f\cos {\psi }^{\prime },\hfill \\ {\dot{\psi }}^{\prime }=\Omega -{\omega }_{\text{f}}-\frac{f}{R}\sin {\psi }^{\prime }.\hfill \end{array}\end{array}

For weak external forcing, i.e., when f is small compared to {\mu }^{\mathrm{3/2}}, and neglecting ξ, the synchronized solution has an amplitude R\simeq \sqrt{\mu } [39], and it is frequency-locked to the external forcing if {\dot{\psi }}^{\prime }=0, with a constant phase difference {\psi }^{\prime } given by f\sin {\psi }^{\prime }=({\omega }_0-\omega )\sqrt{\mu }, if and only if the forcing amplitude f is larger than \sqrt{\mu }|{\omega }_0-\omega |, which defines the tongue in Fig. 4.

4 Coupled subunits

In the previous section, we studied synchronization and frequency locking within a subunit that represented a set of hair cells under a sallet (green arrow in Fig. 1). We now add another level to the hierarchy and study an assembly of M different subunits (red arrow in Fig. 1). Within each subunit, we consider, as before, mean-field coupling between oscillators, whereas the subunits are themselves coupled to their nearest neighbors by the subunits mean fields. Each subunit contains N oscillators and is characterized by its natural frequency distribution. We use the same uniform distribution inside a subunit as in the previous sections but shift the average natural frequency from one subunit to another. These average natural frequencies are linearly distributed so that the complete set of frequencies \left\{{\omega }_{lj}\right\} reads

\begin{array}{}\text{(10)}& {\omega }_{lj}={\overline{\omega }}_l+\frac{2j-N-1}{N-1}\pi \text{Δ}\mathrm{,}{\overline{\omega }}_l=l\delta \omega \mathrm{.}\end{array}

The indices j=\mathrm{1,2,}\cdots \mathrm{,}N and l=\mathrm{1,2,}\cdots \mathrm{,}M indicate the j-th oscillator in the l-th subunit. \delta \omega is the difference between the average natural frequencies of two neighboring subunits [38]. We choose to distribute the average frequencies {\omega }_l linearly, but a more realistic logarithmic spacing can be considered as well and does not change our results.

The dynamics of the j-th oscillator in the l-th subunit is given by

\begin{array}{}\text{(11)}& \begin{array}{ll}{\dot{z}}_{lj}=\hfill & (1+\text{i}{\omega }_{lj})z_{lj}-|z_{lj}{|}^2{z}_{lj}+k_{\text{i}\text{n}}(d_R+\text{i}{d}_I)[({\overline{z}}_l-z_{lj})\hfill \\ \hfill & +k_{\text{o}\text{u}\text{t}}({\overline{z}}_{l+1}+{\overline{z}}_{l-1}-2z_{lj})]+f{\text{e}}^{\text{i}{\omega }_{\text{f}}t}\mathrm{.}\hfill \end{array}\end{array}

The mean field {\overline{z}}_l is defined locally in each subunit l as {\overline{z}}_l\equiv \frac{1}{N}{\displaystyle {\sum }_{j=1}^N}{z}_{lj}. The coupling strength within a subunit is denoted as k_{\text{i}\text{n}}, whereas the coupling between neighboring subunits is denoted as k_{\text{o}\text{u}\text{t}}. In our numerical simulations, we used M=10 subunits, \delta \omega =0.5, and Dirichlet-type boundary conditions for the first (l=0) and last (l=M) subunits.

We first study the unforced case (f=0) and ask how the synchronization transition of a subunit is altered by the coupling to neighboring subunits. We find that a nonzero positive k_{\text{o}\text{u}\text{t}} does not affect the nature of the transition observed in a single subunit, but it reduces the critical value of k_{\text{i}\text{n}} at which it occurs. Figure 7 shows the (k_{\text{i}\text{n}}\mathrm{,}{k}_{\text{o}\text{u}\text{t}}) parameter plane, where explosive synchronization occurs in the subunit l=4. Similar behavior is observed in the other subunits. This relationship between the critical k_{\text{o}\text{u}\text{t}} and k_{\text{i}\text{n}} can be briefly explained as follows. k_{\text{i}\text{n}} quantifies the strength of the restoring force that brings oscillator z_{lj} to the value {\overline{z}}_l, i.e., that synchronizes the oscillator in subunit l with the average frequency of the subunit. The role of k_{\text{o}\text{u}\text{t}} is to synchronize the oscillator with its two neighboring subunits. In other words, because each individual oscillator is coupled to its subunit (by k_{\text{i}\text{n}}) and to its two neighboring subunits (by k_{\text{o}\text{u}\text{t}}), all the oscillators are synchronized to a common frequency, i.e., their average frequency. Therefore, k_{\text{i}\text{n}} and k_{\text{o}\text{u}\text{t}} have the same role in the synchronization process, and increasing either one or the other has the same effect, which is exactly what we observe in Fig. 7.

We then consider an incoming sound wave (f>0) and focus on how the number of oscillators in a subunit, N, affects signal detection. We choose k_{\text{i}\text{n}}=0.6 and k_{\text{o}\text{u}\text{t}}=0.2 in the nonsynchronized region of the coupling parameter space (Fig. 7) and measure the local order parameter R_l=|{\overline{z}}_l| in the subunit l as a function of the frequency {\omega }_{\text{f}} of the incoming sound wave. The upper panels of Fig. 8 show the results for l=\mathrm{4,5,6}, and 7 (corresponding to {\overline{\omega }}_l values of 2.0 to 3.5) with N=11 [Fig. 8(a)] and N=101 [Fig. 8(b)]. As in the case of a single subunit, the resonance of the subunit occurs at a forcing frequency close to the average of the natural frequency distribution of the subunit, although it is shifted slightly upward. The resonance peak is sharp, which indicates good frequency selectivity. We also observe that the signal-to-noise ratio is enhanced when N is increased, mainly because the value of R away from the resonance is lower when N is larger, as seen in Fig. 6. The behavior of the system is easily deduced from the behavior of the M subunits composing it, and our model predicts a discrete set of resonances at frequencies given by the distribution of the average frequencies in each subunit.

Noise is always unavoidable in realistic situations. To probe the robustness of our model (11) to noise, we added a noise term to the external force and replaced f{\text{e}}^{\text{i}{\omega }_{\text{f}}t} with f{\text{e}}^{\text{i}{\omega }_{\text{f}}t}+\xi (t), where \xi (t) is a Gaussian white noise with complex values and strength σ: \xi (t)={\xi }_1(t)+\text{i}{\xi }_2(t), where \langle {\xi }_d(t)\rangle =0, and \langle {\xi }_d(t){\xi }_{d^{\prime }}(t^{\prime })\rangle ={\sigma }^2{\delta }_{d\mathrm{,}{d}^{\prime }}\delta (t-t^{\prime }) \forall d\mathrm{,}{d}^{\prime }\in \left\{\mathrm{1,2}\right\}. We find that when σ is reasonably small, the results are the same, indicating that our model (11) is robust against noise. Figures 8(c) and (d) show the results for \sigma =0.1; the main features of Figs. 8(a) and (b) still appear, although the signal-to-noise ratio is reduced. Increasing the noise amplitude σ increases the baseline and reduces the signal-to-noise ratio.

5 Discussion and conclusions

In summary, we investigated a model of globally coupled oscillators having two different layers of organization inspired by the known structure of the cochlea. Our results show that the synchronization of single oscillators can be concomitant with phase-locking to an external stimulus, provided that the coupling is strong enough and the external stimulus is powerful enough. In this situation, the selectivity is high and increases with the number of oscillators at the level of the subunit, whereas the number of discrete frequencies that are detectable is directly related to the number of subunits.

The results obtained with our simple model may give interesting insight about the properties of the transduction process. In particular, we found that, as the number of oscillators in a subunit increases, both f_{\text{m}\text{i}\text{n}} — which is related to the minimum amplitude of the incoming sound to be detected — and R_{\text{n}\text{s}} — which is related to the noise level that should be overcome by the output signal — decrease as 1/N. This indicates that larger organs, or more precisely larger subunits, are not only sensitive to weaker acoustical signals, but also have a better signal-to-noise ratio. This property is conserved considering that the organ is composed of coupled subunits.

References

[1] L. Robles and M. A. Ruggero, Mechanics of the mammalian cochlea, Physiol. Rev. 81, 1305 (2001)

[2] A. Hudspeth, Hearing, in: Principles of Neural Science, 4th Ed., McGraw-Hill, 2000

[3] See, for example, J. O. Pickles, An Introduction to the Physiology of Hearing, 2nd Ed., Academic Press, 1988

[4] S. S. Narayan, A. N. Temchin, A. Recio, and M. A. Ruggero, Frequency tuning of basilar membrane and auditory nerve fibers in the same cochleae, Science 282(5395), 1882 (1998)

[5] J. F. Ashmore, G. S. Géléoc, and L. Harbott, Molecular mechanisms of sound amplification in the mammalian cochlea, Proc. Natl. Acad. Sci. USA 97(22), 11759 (2000)

[6] K. E. Nilsen and I. J. Russell, The spatial and temporal representation of a tone on the guinea pig basilar membrane, Proc. Natl. Acad. Sci. USA 97(22), 11751 (2000)

[7] D. T. Kemp, Stimulated acoustic emissions from within the human auditory system, J. Acoust. Soc. Am. 64(5), 1386 (1978)

[8] R. Probst, B. L. Lonsbury-Martin, G. K. Martin, B. L. Lonsburymartin, and G. K. Martin, A review of otoacoustic emissions, J. Acoust. Soc. Am. 89(5), 2027 (1991)

[9] V. M. Eguíluz, M. Ospeck, Y. Choe, A. J. Hudspeth, and M. O. Magnasco, Essential nonlinearities in hearing, Phys. Rev. Lett. 84(22), 5232 (2000)

[10] S. Camalet, T. Duke, F. Jüicher, and J. Prost, Auditory sensitivity provided by self-tuned critical oscillations of hair cells, Proc. Natl. Acad. Sci. USA 97(7), 3183 (2000)

[11] J. Cartwright, D. González, and O. Piro, Nonlinear dynamics of the perceived pitch of complex sounds, Phys. Rev. Lett. 82(26), 5389 (1999)

[12] K. A. Montgomery, M. Silber, and S. A. Solla, Amplification in the auditory periphery: The effect of coupling tuning mechanisms, Phys. Rev. E 75(5), 051924 (2007)

[13] M. O. Magnasco, A wave traveling over a Hopf instability shapes the cochlear tuning curve, Phys. Rev. Lett. 90(5), 058101 (2003)

[14] T. Duke and F. Jülicher, Active traveling wave in the cochlea, Phys. Rev. Lett. 90(15), 158101 (2003)

[15] P. L. Boyland, Bifurcations of circle maps: Arnol’d tongues, bistability and rotation intervals, Commun. Math. Phys. 106(3), 353 (1986)

[16] A. Kern and R. Stoop, Essential role of couplings between hearing nonlinearities, Phys. Rev. Lett. 91(12), 128101 (2003)

[17] R. Stoop and A. Kern, Two-tone suppression and combination tone generation as computations performed by the hopf cochlea, Phys. Rev. Lett. 93(26), 8103 (2004)

[18] K. Dierkes, B. Lindner, and F. Jüicher, Enhancement of sensitivity gain and frequency tuning by coupling of active hair bundles, Proc. Natl. Acad. Sci. USA 105(48), 18669 (2008)

[19] A. Vilfan, and T. Duke, Frequency clustering in spontaneous otoacoustic emissions from a Lizard’s ear, Biophys. J. 95(10), 4622 (2008)

[20] M. Gelfand, O. Piro, M. O. Magnasco, and A. J. Hudspeth, Hudspeth a J. Interactions between hair cells shape spontaneous otoacoustic emissions in a model of the Tokay Gecko’s cochlea, PLoS One 5(6), e11116 (2010)

[21] H. P. Wit and P. van Dijk, Are human spontaneous otoacoustic emissions generated by a chain of coupled nonlinear oscillators? J. Acoust. Soc. Am. 132(2), 918 (2012)

[22] Z. Liu, B. Li, and Y.-C. Lai, Enhancing mammalian hearing by a balancing between spontaneous otoacoustic emissions and spatial coupling, Europhys. Lett. 98(2), 20005 (2012)

[23] G. A. Manley, Cochlear mechanisms from a phylogenetic viewpoint, Proc. Natl. Acad. Sci. USA 97(22), 11736 (2000)

[24] C. Köppl, Morphology of the basilar papilla of the bobtail lizard Tiliqua rugosa, Hear. Res. 35(2–3), 209 (1988)

[25] J. Gómez-Gardenes, S. Gómez, A. Arenas, and Y. Moreno, Explosive synchronization transitions in scale-free networks, Phys. Rev. Lett. 106(12), 128701 (2011)

[26] I. Leyva, R. Sevilla-Escoboza, J. M. Buldú, I. Sendina-Nadal, J. Gomez-Gardenes, A. Arenas, Y. Moreno, S. Gómez, R. Jaimes-Reátegui, and S. Boccaletti, Explosive first-order transition to synchrony in networked chaotic oscillators, Phys. Rev. Lett. 108(16), 168702 (2012)

[27] P. Ji, T. K. D. Peron, P. J. Menck, F. A. Rodrigues, and J. Kurths, Cluster explosive synchronization in complex networks, Phys. Rev. Lett. 110(21), 218701 (2013)

[28] X. Zhang, X. Hu, J. Kurths, and Z. Liu, Explosive synchronization in a general complex network, Phys. Rev. E 88(1), 010802 (2013)

[29] I. Leyva, A. Navas, I. Sendina-Nadal, J. A. Almendral, J. M. Buldú, M. Zanin, D. Papo, and S. Boccaletti, Explosive transitions to synchronization in networks of phase oscillators, Sci. Rep. 3, 1281 (2013)

[30] Y. Zou, T. Pereira, M. Small, Z. Liu, and J. Kurths, Basin of attraction determines hysteresis in explosive synchronization, Phys. Rev. Lett. 112(11), 114102 (2014)

[31] X. Zhang, Y. Zou, S. Boccaletti, and Z. Liu, Explosive synchronization as a process of explosive percolation in dynamical phase space, Sci. Rep. 4, 5200 (2014)

[32] X. Hu, S. Boccaletti, W. Huang, X. Zhang, Z. Liu, S. Guan, and C. H. Lai, Exact solution for first-order synchronization transition in a generalized Kuramoto model, Sci. Rep. 4, 7262 (2014)

[33] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Explosive synchronization in adaptive and multilayer networks, Phys. Rev. Lett. 114(3), 038701 (2015)

[34] T. Qiu, Y. Zhang, J. Liu, H. Bi, S. Boccaletti, Z. Liu, and S. Guan, Landau damping effects in the synchronization of conformist and contrarian oscillators, Sci. Rep. 5, 18235 (2015)

[35] P. C. Matthews, R. E. Mirollo, and S. H. Strogatz, Dynamics of a large system of coupled nonlinear oscillators, Physica D 52(2–3), 293 (1991)

[36] C. Wang and N. Garnier, Continuous and discontinuous transitions to synchronization, arXiv: 1609.05584

[37] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77(1), 137 (2005)

[38] P. D. Welsby, The 12, 24, or is it 26 cranial nerves? Postgrad. Med. J. 80(948), 602 (2004)

[39] L. Fredrickson-Hemsing, S. Ji, R. Bruinsma, and D. Bozovic, Mode-locking dynamics of hair cells of the inner ear, Phys. Rev. E 86(2), 21915 (2012)

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2017