Clustering in an array of nonlinear and active oscillators as a model of
spontaneous otoacoustic emissions
Liv Moretto Sørensen(1) , Peter Leer Bysted(2) , Bastian Epp(3)(1) Niels Bohr Institute, University of Copenhagen, bvf365@alumni.ku.dk(2) Department of Health Technology, Technical University of Denmark, s144045@student.dtu.dk(3) Department of Health Technology, Technical University of Denmark, bepp@dtu.dk

Abstract
Spontanoeus otoacoustic emissions (SOAEs) are sound signals emitted by the inner ear of many species in theabsence of an external stimulus. The precise mechanism underlying the generation of SOAEs is still unknown.
This paper investigates the phenomenon of frequency clustering in a system of coupled, nonlinear oscillators
as a possible mechanism underlying the generation of SOAEs. A generic system consisting of 100 Van der
Pol oscillators was implemented with dissipative and reactive coupling elements of varying strength. Frequency
clusters were investigated in terms of size, frequency and phase coherence as a function of coupling strength andscaling of the nonlinearity. One conﬁguration of the modelled system was synchronised by an external driving
force. Phase relations between oscillators were analysed and compared both before and after synchronisation.
Frequency clusters appeared for some conﬁgurations of the coupled system. Dissipative coupling resulted inclusters of equal frequency spacing and similar size, while reactive coupling resulted in fewer clusters and much
more complex dynamics of the resulting system. Even though the phenomenon of frequency clustering is welldescribed in the literature, less is known about both the phase relation and coherence within and across clusters.
These two factors play an important role in the generation of the observed SOAEs and might help to identify thecriticial parameters of the underlying mechanism in inner ear mechanics.
Keywords: Cochlea, SOAE, nonlinear dynamics, otoacoustic emission, oscillators

1

INTRODUCTION

The human ear’s ability to detect extremely small differences in sound levels and to discriminate frequencies
with high precision has yet to be explained fully. It was argued in [6] that it would not be possible to achieve
the healthy ear’s sensitivity with passive mechanisms. In this connection the presence of an active mechanicalfeedback mechanism in the ear was postulated, often referred to as the "active process" [13, 5]. This was
supported by the discovery of the spontaneous otoacoustic emissions (SOAEs) [7, 15]. SOAEs are sound signalsthat can be measured in the absence external stimuli, pointing towards self-sustained activity in the inner ear,possibly linked to outer hair cell motility [1]. An interesting characteristic of these emissions are the narrow
distributions of energy in the spectra of the emitted signal. SOAEs are vulnerable to a damage of outer haircells and tend to disappear even for mild hearing impairment [10]. Although the lizard basilar papilla shows a
simpler structure than the human cochlea, the SOAEs share some interesting characteristics [2]. Hence lizard
inner ear models might provide important insights for mammalian SOAE generation [13, 5]. It was reported [8]that SOAEs in bobtail lizards contain roughly 10 evenly spaced peaks in the power spectrum. This number ofpeaks is much lower than the number of hair cells in the bobtail lizard. Under the assumption that outer haircells play a key role, it is therefore believed that numerous cells must be synchronised within the basilar papillato produce sound with the characteristics of SOAEs. A model that has been of main interest in the investigationof the active process in the inner ear is a system of nonlinear, active, coupled limit cycle oscillators [11]. Thiswas ﬁrst suggested by [6] and modeled in different varieties [4, 3, 13, 5, 14]. These models can account for
many key aspects of the phenomena observed in SOAEs, including effects of entrainment. Entrainment is thephenomenon that causes oscillators with different natural frequencies to group in clusters with neighbouring688

oscillators, all oscillating at the same frequency [9]. This phenomenon can either be triggered by intrinsic
oscillations within the system or by an external driving force. Even though both systems of coupled oscillators
and SOAEs are well investigated, the methods to quantify the observations differ. In order to transfer insightsfrom one ﬁeld to the other, the implications, differences, and similarities of the various analysis methods need
to be linked in order to draw conclusions about the underlying mechanisms. The only way to analyse SOAEs inthe ear canal is by means of Fourier analysis while it is common to describe nonlinear dymamics with the helpof, e.g., phase diagrams. Because only one state variable is available for SOAEs (the pressure in the ear canal)
the assumption of a periodic oscillation is commonly made, even though nonlinear and active oscillators often
shown non-harmonic dynamics. In this paper, a model- and analysis framework inspired by previous modelsof coupled, nonlinear, active oscillators is described. A system of coupled Van der Pol oscillators [12, 9] iscreated to investigate the changes in behaviour both globally and locally while the parameters that deﬁne thedynamics and the coupling between oscillators are varied. Speciﬁcally, the effects of dissipative and reactive
coupling between nearest neighbours are investigated with focus on the formation of clusters. These clusters areanalysed with respect to size, stability and phase coherence. Synchronisation effects are quantiﬁed with respect
to cluster size and cluster frequency.

2

METHODS

2.1 Model
The modeled system consisted of N = 100 Van der Pol oscillators, each described by the equation of motion:ẍ = −ω 2 x + μ(1 − x2 )ẋ(1)where x is the displacement from equilibrium, ẋ and ẍ are the velocity and acceleration,
 respectively. The factor
μ scales the nonlinearity of the system. In the case μ = 0, the eigenfrequency ω = K/M of the oscillator isdeﬁned by the spring constant K and the mass M (for this study, all masses were set to M = 1). To accountfor tonotopy, a linear gradient in the eigenfrequency was introduced by(n − 1)ωmax 2
Kn = ωmin +
· εn
N −1

n ∈ {1...N}εn = N (1, ν)(2)The value of Kn was varied stochastically for each oscillator by a random variable drawn from a normal distribution N with mean one and standard deviation ν to aid generation of clusters [9]. Because the phase plane
trajectory of a Van der Pol oscillator varies both with μ and ω, a scaling of the time differential was introduced
by a scaling of μ to keep the ratio ω/μ constant and to result in identical phase trajectories of all N oscillators.
This yields the equation of motion for oscillator nẍn = −xn +μ(1 − x2 )
ẋn
ωn

(3)

The force between two coupled neighbouring oscillators is given by
Fn = dn+1 (ẋn+1 − ẋn ) + kn+1 (xn+1 − xn )+dn (ẋn−1 − ẋn ) + kn (xn−1 − xn )(4)

where the dissipative and reactive coupling are described by d and k, respectively. Positions x0 and xN+1 were
held constant to zero. For the systems analysed in this paper, the boundary conditions were ﬁxed, which meansthat oscillators n = 1 and n = N are coupled to the ﬁxed boundary points. Extending (3) with (4) and adding
an external driving force, Fext , yields a second order differential equation. This can be rewritten as a set of ﬁrst

689

order differential equations for numerical treatment:
dxn
= ẋn
dt
d ẋn
μ(1 − x2 )= −xn +
ẋn + Fn + Fext
dt
Kn

(5)

2.2 Analysis
2.2.1 Steady state criteria
The dynamics of the system was simulated for time intervals of length Tint and stopped as soon as a steady stateof the system was reached. This state was deﬁned by a “nearly constant” energy over time and the requirementthat the minimum energy in the system must exceed a value of 10−2 · N. The energy of each oscillator wasapproximated by:
En (m) =(m+1)Tint∑

xn (t)2 , m = 0, 1, 2...

(6)

t=mTint

If En (m + 1) did not differ from En (m) by more than ten percent for all n, the system was deﬁned to be insteady state.
2.2.2 Frequency and cluster analysis
Three methods were used to analyze the frequency of each oscillator: 1) “Fourier mean method” (Ωf,mean ),computed as the weighted mean of all frequencies Ωi in the Fourier domain, 2) “Fourier peak method” (Ωf,peak ),computed by localising the frequency with the highest normalized peak m in the Fourier domain, and 3) the“Zero crossing mean method” (Ωzc,mean ) where the frequency was determined from the inverse of the mean
distance between zero crossings Δtzc in the time domain:Ωf,mean = ∑mi Ωi;Ωf,peak = Ω(find(m = 1))

;

i


−1
1
∑ ΔtzcΩzc,mean = = 2
P
nzc

(7)

Clusters of minimum four oscillators were analyzed with the system in steady state. In order to qualify as
a cluster, the frequency differences between oscillators within one cluster must have been smaller than half
the frequency difference between to neighboured oscillators assuming a linear gradient between the natural
frequencies of the two oscillators with lowest and highest index n in the cluster:
stepcluster <(Ωhigh − Ωlow )2N

where frequencies Ωhigh and Ωlow can be found using any of the frequency analysis methods mentioned above.
The phase coherence of the oscillators was computed using the estimated oscillation period combined with theamplitude and velocity of each oscillator. The normalized [displacement, velocity]-vectors after each oscillation
period were added in the phase plane and divided by the number of periods the oscillator has gone through.
2.3 Numerical evaluation
The coupled equations of motion were solved using MATLAB’s built in ODE45 solver. For all simulations,values in (2) were set to
N = 100 , ωmin = 2π , ωmax = 5 · 2π
(8)The sampling frequency was set to fs = 1000Hz and the maximum step size used by ODE45 to 0.01s. The
simulation was run in intervals of length Tint = 20s until the system met the steady state criteria or until the
maximum simulation time of 1000s was reached. All oscillators had initial displacement and velocity equal to

690

B

0

20
40
60
80
Oscillator number

Magnitude (normalised)

Frequency [Hz]Time [s]

20
60
100

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

A

1
2
3
4
5
6

100

Figure 1. System Pd . A) Spectral representation of the system dynamics. A 20 s steady state time interval
was windowed using a Hann window prior to Fourier transformation. The spectrum for each oscillator was
normalised to its maximum value to improve the graphical representation. B) State of each oscillator as functionof time (in intervals of length Tint = 20s). Steady state criteria are met for values of 1 and are not met for valuesof 0. C) Clusters found when analysed using Ωf,peak . The system exhibited a regular cluster structure for theseven equally spaced clusters. The cluster size was biggest at the ends of the row of oscillators. D) Clusters
for analysis with Ωzc,mean . Here some oscillators were located between clusters, showing intermediate cluster
structure.
0 and were excited by a short rectangular pulse at time t = 0s with width w = 0.05s and unity height. The
systems were analysed for parameter combinations P = {d, k, μ, ν}. To get a more explicit understanding of
the impact of dissipative and reactive coupling, systems were coupled by either dissipative or reactive couplingelements of varying strength. The following parameters were used for the simulations:d ∈ {0, 1, 2, 5, 10, 20, 50, 100},k ∈ {0, 1, 2, 5, 10, 20, 50, 100},μ ∈ {0.5, 1, 1.5 ... 4, 4.5, 5},ν ∈ {0, 0.1} (9)In the absence of stochastic variations of Kn (i.e., ν = 0) the system was simulated once, otherwise 25 times.
The system with dissipative coupling Pd = {100, 0, 1, 0} and the system with reactive coupling Pk = {0, 20, 1, 0}
will be described in more detail because they showed the most comprehensive clustering behaviour.
To analyse synchronisation effects by an external driving force, a force Fext = Fsin was applied to all oscillators:
Fsin (t) = α sin(2π f t)
with α the amplitude and f the frequency of the sinusoidal force. The cluster size and -frequencies wereanalysed as a function of both α and f :fd ∈ {0, 0.25, 0.5 ... 6.5, 6.75, 7},

fk ∈ {0, 0.25, 0.5 ... 9.5, 9.75, 10},α ∈ {0, 1, 10, 50, 100, 150}

with fd and fk denoting the frequencies of the forces driving Pd and Pk , respectively.

3

RESULTS

3.1 Dissipative coupling
Figure 1A shows the spectral analysis of system Pd in steady state (reached after 80s of simulation, see Figure
1B). Clusters appeared at frequencies between 1 Hz to 5.5 Hz with equal spacing and size. This representationshows harmonics of each oscillator at frequencies of neighbouring clusters. The clusters overlap at the cluster
edges, indicating two dominant frequencies for these oscillators. When extracting clusters using Ωf,peak (Figure1C), all oscillators were assigned to a cluster. Analysis of the same data with Ωzc,mean (Figure 1D), resulted
in a collection of oscillators that was not discrete but more smooth with some oscillators in between clusters,
consistent with [9].

691Mean clustersize
Synchronisation (%)

20

A

B

C

D

E

15
10
5

1
0.5
0

0

1
d=100
d=50

2

3

4

5
d=20
d=10

0

1

2
d=5
d=2

3

4

5
d=1

Figure 2. A-D) Mean cluster size, the number of oscillators being part of a cluster (expressed as synchronisation) for systems with varied combinations of coupling strength d and nonlinearity μ and two methods offrequency analysis (Fourier mean method, A,C; zeros crossing mean method, B,D). E) Cluster analysis usingΩzc,mean with estimates of the standard deviation of the estimated frequency (grey), the estimated phase coherence (red) and the temporal dynamics of two selected oscillators (insets) of system Pd in steady state.
The number and size of the found clusters depended on the combination of coupling strength (dn ), the value ofthe of nonlinearity parameter μ and the method used to analyze the dominant oscillation period of each oscillator. Figure 2 shows the size and stability of the oscillation period (expressed as synchronization) of identiﬁed
clusters for different combinations of dn and μ for Ωf,peak (A,C) and Ωzc,mean (B,D). Larger clusters were foundfor stronger coupling and for higher degrees of nonlinearity. The mean cluster size and the synchronization
were slightly larger when using Ωf,peak compared to Ωzc,mean . While the maximum synchronization for Ωf,peak
was close to 1, the maximum value for Ωzc,mean was around 0.8.
The dynamics of the oscillators within a cluster or at the corners of a cluster are shown in Figure 2E forthe zero crossings mean method. The phase coherence in the middle of clusters is much higher than at theborders, where the uncertainty on the frequency analysis is correspondingly higher. The temporal dynamics
showed a clear beating of one oscillator located between two clusters, consistent with two dominant frequencies
as indicated in Figure 1.
3.2 Synchronization with harmonic driving force
For the phase analysis the central cluster of system Pd was used. This cluster consisted of 11 oscillators witha mean frequency of 3.3 Hz. The phase analysis of each oscillator in each cluster is shown in the absence(upper panel) and presence (lower panel) of an external driving force in Figure 3. The phase is shown foreach time period p = 1/3.3s during the time interval Tint for each oscillator. In both cases, the ﬁrst and last
oscillators showed a much wider spread of phases compared to the oscillator in the middle of the cluster. The
spread of phase is consistent with the estimate of the phase coherence (red lines). Furthermore, the mean of
the phase distribution changed along the oscillators in the analysed cluster, from a mean phase of around 5/4π
for oscillator 42 and a mean phase of around 1/3π for oscillator 52.
In the presence of an external driving force with frequency 3.25 Hz and amplitude α = 10, the cluster grows insize and the phase distributions of all oscillators become more narrow. Especially in the middle of the cluster,
the phases collapse into a very narrow region around the mean, indicating a very stable oscillation period.
Neighboured oscillators show though a slight shift in mean phase, similar to the result in the absence of the

692

A

E

Figure 3. A) Phase analysis of the central cluster in the system Pd in steady state. The cluster consisted ofoscillators n = 42 to n = 52. The mean frequency of the cluster was 3.3 Hz. The phase of each oscillator wasinvestigated by marking points in phase space for each time period p = 1/3.3 s. Same legend as in 2. B) Phase
analysis of central cluster in the system Pd when in steady state while being synchronised by external force
with frequency f = 3.25 Hz and amplitude α = 10. The cluster consisted of oscillators n = 33 to n = 62.
external driving force.
3.3 Reactive coupling
Figure 4 shows the spectral analysis of a system using reactive coupling (A,B) and an analysis of cluster sizeand synchronization analogue to Figure 2. Two big clusters with frequencies around fl = 3.5Hz and fh = 9Hz
are apparent in Figure 4A. The cluster frequency around 9 Hz was higher than the maximum eigenfrequency
in the uncoupled system, indicating a contribution of the coupling elements to the dominant frequency of theoscillators. Inbetween the clusters, the oscillators showed a complex distribution of energy and no consistentclustering behaviour. Consistent with the results using dissipative coupling, clusters tended to be larger for
higher amount of coupling and higher degree of nonlinearity, but the behaviour was much more complex (Figure 4C-F). The high variability in the synchrony of the single oscillators (panels D,F) are consistent with thecomplex distribution of energy found in the spectral analysis of the system (panel A). These results indicate
a more complex dynamics and impact on the resulting frequencies when using reactive instead of dissipative
coupling.

4

DISCUSSION

For systems with dissipative coupling, similar cluster patterns compared to both [13, 5, 4] appeared in the systems. Although no frequency analysis method has been deﬁned by [13], their results show a high similarity to
the results in the current study using Ωf,peak . The result that stronger coupling results in larger clusters (Figures2,4) can be explained by assuming that the interaction across oscillators happens exclusively through this parameter. With increasing magnitude, the eigenfrequency of each oscillator and the local mechanical parameters
play a decreasing role. Close by oscillators will start to synchronize and, as a compound effect, increase the693Frequency [Hz]Time [s]

3
5
7
9
B

20
60
0

20
40
60
80
Oscillator number

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

Magnitude (normalised)

A

1

C

D

E

F

100

Figure 4. A-B) Same as Figure 1, but for system Pk with reactive coupling. C-F) Same as in 2, but forvariations of reactive coupling constant k.
effect on neighboring oscillators. This will be the case until the balance between across-oscillator and local
oscillator effects will balance and move towards the next oscillator cluster. Given that clustering is an effect ofactive, nonlinear oscillators, it is also intuitive that increasing the amount of nonlinearity increases the effect ofclustering, leading to larger clusters.
Regarding the reactive coupling, similarities to the data presented by [4] can be observed in the spectral representation in terms of some clearly deﬁned clusters and some more blurry areas between the clusters. Consistent
with this, [5] found that the number of clusters decreases for bigger values of the reactive coupling coefﬁcient,which might be connected to the increasing mean cluster size documented in this paper. One argument why
reactive coupling shows a more complex behaviour compared to dissipative coupling is that the coupling spring
constants will affect the local eigenfrequency and hence the local tuning and by that the frequency gradientacross oscillators.
The often neglected information of oscillator clusters is the phase relation across oscillators within one cluster.
With the results in Figure 3 it becomes clear that even though oscillators might be entrained to a speciﬁcfrequency, the local properties of the oscillators will have an effect on the phase relation within a cluster. For
phenomena that represent the compound result of a whole cluster like SOAEs measured in the ear canal, thisphase shift will have important consequences on the measured magnitude. The measured magnitude can hence
not be used to estimate the size of the cluster because the addition of phase-delayed components will reduce theresulting compound amplitude, in the extreme case to zero for components in anti phase. Phase distributions are,in general, more narrow for oscillators located centrally in a cluster that for oscillators remote from the cluster
center. This might be connected to frequency difference between the cluster frequency and the eigenfrequencyof the local oscillator. An external driving force heavily reduces the width of the phase distributions, making
them very narrow for oscillators in the center of the cluster. It is, however, interesting to note that the phaseshift of oscillators within a cluster remains the same as in the absence of the external driving force.

5

CONCLUSIONS

A system of coupled Van der Pol oscillators with a frequency gradient showed clustering effects for both dissipative and reactive coupling. The number of and distance between clusters is dependent on the magnitude andtype of the coupling between oscillators. Hence, the distance observed between SOAEs might provide information about the type and strength of coupling in the inner ear. The size and frequency of clusters depends,

694however, on the metric used to quantify the oscillation period of each oscillator. Finally, oscillators within one
cluster showed a linear phase gradient. Hence, a better understanding of the parameters deﬁning this phasegradient and the interaction with the size and frequency distance of clusters is required to be able to link theseproperties to analogue characteristics observed in SOAEs.

ACKNOWLEDGEMENTS
This work was supported by the William Demant foundation through a travel grant (LMS). We would like tothank Mogens Høgh Jensen for constructive discussions about this project.

References[1] J. Ashmore. Cochlear outer hair cell motility. Physiological Reviews, 88(1):173–210, 2008.
[2] C. Bergevin and A. Salerno. Dynamics of spontaneous otoacoustic emissions. Aip Conference Proceedings,
1703(1):090024, 2015.
[3] B. Epp, J. L. Verhey, and M. Mauermann. Modeling cochlear dynamics: Interrelation between cochlea
mechanics and psychoacoustics. Journal of the Acoustical Society of America, 128(4):1870–1883, 2010.
[4] B. Epp, H. Wit, and P. van Dijk. Clustering of cochlear oscillations in frequency plateaus as a tool toinvestigate soae generation. Aip Conference Proceedings, 1703(1):090025, 2016.
[5] M. Gelfand, O. Piro, M. O. Magnasco, and A. J. Hudspeth. Interactions between hair cells shape spontaneous otoacoustic emissions in a model of the tokay gecko’s cochlea. Plos One, 5(6):e11116, 2010.
[6] T. Gold. Hearing. ii. the physical basis of the action of the cochlea. Proceedings of the Royal Society of
London. Series B - Biological Sciences, 135(881):492–498, 1948.
[7] D. Kemp. Stimulated acoustic emissions from within human auditory-system. Journal of the AcousticalSociety of America, 64(5):1386–1391, 1978.
[8] C. Köppl and G. Manley. Spontaneous otoacoustic emissions in the bobtail lizard .1. general-characteristics.
Hearing Research, 71(1-2):157–169, 1993.
[9] G. V. Osipov, J. Kurths, and C. Zhou. Synchronization in oscillatory networks. 2007.
[10] R. Probst, B. Lonsbuty-Martin, and G. Martin. A review of otoacoustic emissions. Journal of the AcousticalSociety of America, 89(5):2027–2067, 1991.
[11] C. Talmadge, A. Tubis, H. Wit, and G. Long. Are spontaneous otoacoustic emissions generated by selfsustained cochlear oscillators. Journal of the Acoustical Society of America, 89(5):2391–2399, 1991.
[12] B. van der Pol. Lxxxviii.on relaxation-oscillations. London, Edinburgh, and Dublin Philosophical Magazine
and Journal of Science, 2(11):978–992, 1926.
[13] A. Vilfan and T. Duke. Frequency clustering in spontaneous otoacoustic emissions from a lizard’s ear.
Biophysical Journal, 95(10):4622–4630, 2008.
[14] H. P. Wit and P. van Dijk. Are human spontaneous otoacoustic emissions generated by a chain of coupled
nonlinear oscillators? Journal of the Acoustical Society of America, 132(2):918–926, 2012.
[15] P. Zurek. Spontaneous narrowband acoustic-signals emitted by human ears. Journal of the Acoustical Societyof America, 69(2):514–523, 1981.

695