Lumped parameter model and Monte Carlo simulation to study middle ear
uncertainties
Lucas LOBATO(1) , Stephan PAUL(2) , Júlio A. CORDIOLI(3)(1) Acoustic and Vibration Laboratory, Federal University of Santa Catarina, Brazil. lucas.lobato@lva.ufsc.br(1) Acoustic and Vibration Laboratory, Federal University of Santa Catarina, Brazil. stephan.paul@ufsc.br(1) Acoustic and Vibration Laboratory, Federal University of Santa Catarina, Brazil. julio.cordioli@ufsc.br

Abstract
Several mathematical models of the human middle ear have been developed throughout the last 70 years fordifferent purposes. Although the experimental data of the dynamical behavior of the middle ear show a large
dispersion, only deterministic models have been proposed in the literature. This work aims to develop a probabilistic model of the human middle ear to study its dynamic behavior uncertainties. Then, a coupled lumped
parameter and analytical model is combined to Monte Carlo simulations (MCS). The parameters of the deterministic model was fitted to match the stapes velocity to one experimental results. Then, the input parameters
were modeled with Gaussian probability density functions in the MCS. It was observed the uncertainty of thenatural frequencies increases with the mode order. Also, the dispersion seen experimentally could not be reproduced defining the standard deviation (StD) of the input parameters as 20% of the mean value. Therefore, theStD of the input parameters was increased up to 30%, 40% and 50% of the mean values. Results suggest the
StD should be defined separately for each input parameter in order to reproduced the experimental dispersion.
At last, next steps are set for improvements in the proposed probabilistic approach. Keywords: Middle ear,
Biomechanics, Uncertainties, Probabilistic model1

INTRODUCTION

The middle ear is a biomechanical system that composes the peripheral auditory system of mammals together
with the outer and inner ears. The middle ear is mainly composed of a complex fibrous membrane called tympanic membrane (TM), an ossicular chain comprising three tiny bones (malleus, incus and stapes), two synovial
joints and various soft tissues such as ligaments, tendons and muscles [1]. Two functions are assigned to thephysiological working of the middle ear. At first, the middle ear matches the several different impedances ofthe outer and inner ears, improving the transmission of the sound energy between these two systems. Secondly,
the middle ear protects the inner ear against high static pressures and high sound stimuli through passive andactive working [2].
For a proper understand of the middle ear mechanics under healthy and pathological conditions as well as toevaluate the performance of prosthesis and hearing devices, several mathematical models of the middle ear havebeen proposed in the literature [3]. Since Esser [4] in 1947, who studied the deformation of the TM throughan analytical formulation, all mathematical models of the middle ear were developed using deterministic approaches. In contrast, experimental data of the human middle ear dynamics generally shows large discrepancies
among samples. Consequently, it is still debatable the applicability of deterministic models of the middle ear
once the real structure is a fuzzy system.
In this paper we use a lumped parameter model combined to Monte Carlo simulations in order to develop a
probabilistic model of the human middle ear for studying its dynamic behavior uncertainties. In addition, a
simplified analytical model of the outer ear is coupled to the middle ear model to expand the analysis. In
the end, we aims to provide a preliminary investigation on the statistical properties of the human middle earmechanics through a probabilistic approach.

5426

2

METHODS

Our analysis of the middle ear’s uncertainties is based on four different measures: the first three natural frequencies fn , the frequency response of the stapes velocity by the acoustic pressure at the TM Hfp , the input
impedance of the middle ear ZME and the wideband Energy Reflectance ER (also called as Power Reflectance).
The definition of the last three measure are given as
Hfp = Vfp /PTM ,(1)

ZME = PTM /VTM

(2)

ER = |RME |2 ,

(3)

and
being Vfp the velocity of the stapes, PTM the sound pressure at the TM, VTM the velocity of the TM and RME
the reflection coefficient of the middle ear seen by the ear canal. These quantities were chosen to reach asmechanical measures ( fn , Hfp and ZME ) as well as vibroacoustical measures (ER). While a pure mechanicalmodel can be used to obtain the mechanical parameters, a vibroacoustical model is necessary to calculate ER
taking account the geometry of the ear canal. Thus, the two needed approaches are described in Section 2.1
and 2.2. Lastly, Section 2.3 details the Monte Carlo simulation used in this work.
2.1 Lumped model of the middle ear
The lumped parameter model used here is based on the model proposed by Feng and Gan [5]. Our model
has four degrees of freedom (DOF) related to the normal displacement of the four presented masses, neglecting
the rotational movements. The masses represent the TM (m1 ), malleus (m2 ), incus (m3 ) and the stapes (m4 ).
The stiffness (ki ) and dampers (ci ) with index 2, 4 and 6 represent the connection of TM and malleus andthe incudomallear joint and incudostapedial joint, respectively. The stiffness and dampers with index 1, 3 and5 represent the ligaments and tendons that hold the TM and ossicular chain in the middle ear cavity. Finally,
stiffness k7 and dampers c7 model the cochlear load over the stapes. Figure 1(a) shows the lumped parametermodel scheme.
Middle ear

k1

f

k3

k5

Ear canal

k7

m1 c1

k2

c2

m2 c3 k4 m3 c5k6 m4
c7
c4

c6

j�t

ZME

u(0,t)= e

0

a)

b)

l

x

Figure 1. (a) the lumped parameter mass-spring model of the middle ear and (b) the analytical model of ear
canal.
Assembling the matrices of mass [M], stiffness [K] and damping [C], the displacement of each DOF is obtained
through the solution of the equation of motion[M]{ẍ} + [C]{ẋ} + [K]{x} = { f },(4)

being {x}T = {x1 , x2 , x3 , x4 } the displacement vector of all the DOFs, {ẋ} and {ẍ} the time derivative of firstand second order, and { f } the excitation vector, defined as{ f }T = {(ATM PTM ) 0

0

0},

(5)

5427

being ATM the TM area defined as π(0.01/2)2 [m2 ]. Note that the excitation acts only at the TM defining zerofor the remaining DOFs. It should to be note that√ defining PTM , the quantities Hfp and ZME can be calculatedusing VTM = i2π f x1 and Vfp = i2π f x4 , being i = −1.
To define the baseline lumped parameter model we fitted the mechanical parameters in order to match the Hfp
from our model to one experimental data of a cadaveric temporal bone. The experimental Hfp was measured ina fresh temporal bone, whereas Vfp and PTM were measured using Laser Doppler Vibrometry and a microphonewith a probe tube. We carried out the experiment reproducing the procedure described by Aibara et al. [6].
Then, the resultant mechanical parameters fitting the theoretical-experimental Hfp are exposed in Table 1. The
comparison of theoretical-experimental Hfp is shown in Figure 2(a). Also, both theoretical-experimental results
are compared to the experimental mean ± standard deviation (StD) range over 11 samples experimented by
Aibara et al. [6].
Following the fitting process of Hfp , we then calculated ZME . The result is shown in Figure 2(b), being thetheoretical result compared to the mean ± StD range over 13 samples of cadaveric temporal bones experimented
by Voss et al. [7]. Our result falls in the StD range up to 1 [kHz], being above the experimental range between
1 [kHz] and 4 [kHz]. Note that although we display ZME up to 10 [kHz], Voss et al. [7] measured only up to
4 [kHz].
Table 1. Mechanical parameters of the baseline model.
Stiffness [kN/m]
k1 = 0.2
k2 = 94.7
k3 = 0.96
k4 = 10k5 = 0.2
k6 = 10k7 = 1.16

Mass [mg]m1 = 2,7m2 = 35.3
m3 = 35.8
m4 = 4.4-

Damping [mNs/m]
c1 = 105c2 = 20c3 = 120
c4 = 66
c5 = 30
c6 = 0.28
c7 = 124

2.2 Analytical model of ear canal
As shown in the previous section, Hfp and ZME can be calculated using the outcomes of the lumped parametermodel. However, to calculate the ER taking account the length of the ear canal, an acoustical model of theear canal is needed. Thus, we assumed the ear canal as an one dimension duct of constant cross section with
an impedance ZME at the end and a velocity condition at the entrance, as shown in Figure 1(b). The acousticpressure field in an one dimension duct can be described byp(x,t) = Aeikx + Be−ikx ei2π f t(6)being k the wave number given as 2π f /c0 , c0 the sound velocity in air and A and B the complex constants
obtained from the boundary conditions. Then, the two needed boundary conditions arep(l,t)= ZME
u(l,t)(7)u(0,t) = ei2π f t ,(8)being l the length of the ear canal. We assumed l = 32 [mm], c0 = 343 [m/s] and the air density ρ0 = 1.21
[kg/m3 ]. Note that the impedance ZME used as boundary condition is previously calculated with the lumped
parameter model described in Section 2.1. Finally, we calculate the reflection coefficient at the end of the ear
canal as RME = B/A and the ER using Equation 3. The computed ER is shown in Figure 2(c), where it is5428

compared to the confidence interval (CI) of 80% measured by Shahnaz et al. [8] in healthy subjects. Our resultof ER is in the CI range in almost all frequency range, being slightly out around 2 [kHz].

10−1

10−2

Exp. LVA
Exp. mean ± std
Model

100

1k

10k

Energy Reflectance [-Exp. mean ± std
Model|ZME | [Ns/m5 ]|Hfp | [mm/s/Pa]

1010

109

108

107
100

1k

10k

Exp. (CI 80%)
Model

1

0.5

0
102

103

Frequency [Hz]Frequency [Hz]Frequency [Hz]

(a)

(b)

(c)

104

Figure 2. Magnitude of Hfp and ZME obtained from the lumped parameter model described in Section 2.1 andEnergy Reflectance obtained from analytical model of the ear canal described in Section 2.2. Model outcomesare compared to the experimental mean ± StD from [6] in (a), [7] in (b) and [8] in (c).
Although our model provides a good results compared to the experimental references, some limitations of themodel should be pointed out. Firstly, it is assumed that the TM acts as a rigid body, what is not true as
observed by Cheng et al. [9]. Furthermore, the model is able to described only an one dimensional movementof the middle ear, also called as “piston-like movement”, despite the ossicular chain presents a more complex
movement in high frequencies [7]. On the other hand, the analytical model of the ear canal does not represent
the ear canal curvature and variation of the cross section area.
2.3 Monte Carlo Simulation description
Monte Carlo simulation (MCS) is a probabilistic parametric method widely employed to analyze the statisticalresponse of structures with uncertainties [10]. The procedure consists in generating n pseudo-random samplesof the random variables of the system, with a given probability density functions (PDF). In what follows, thedeterministic models described in Section 2.1 and 2.2 are computed for each pseudo-random sample of theensemble, resulting in n outcomes. Finally, statistical analyzes can be made over the n results, computing thefirst and second statistics moments, for example.
Due the lack of experimental data of the mechanical parameters of the middle ear, some assumption on thestatistical properties of them have been taken. We assumed all parameters described in Table 1 and the area of
the TM (ATM ) as random variables with a gaussian PDF. The mean value of each parameter is the same usedin the baseline deterministic model. On the other hand, the StDs were defined as a percentage of the mean.
Initially, we assumed a StD of 20% of the mean for all random variables, excepted for the ATM , whereas a StD
of 12% was taken based on the De Greef et al. [1] observations.
In order to define a reasonable n number of pseudo-random samples, a convergence analysis was carried out.
The mean and StD of the first three natural frequencies of the lumped parameter model were calculated varying
the n number of samples. Figure 3 shows the convergence of the mean and StD of the first three natural
frequencies for 10 up to 3000 n samples, with a step sample of 10, resulting in 299 cases computed. In
addition, it is shown the mean of the means and StDs over the all computed cases. The convergence criteria
were defined as ±1% of the mean of the means and ±5% for the mean of the StDs. It was observed a safe
convergence for a number of samples n ≥ 1500, where the variation of the mean for each natural frequency
did not exceeded 1% of the overall mean of the means and 5% of the overall mean of the StDs.

5429

First fn
Frequency [Hz]Third fn

Second fn
3,900

960
Convergence
Mean of means
Mean ± 1%

940
920

8,800
3,800
8,600

900

3,700

880
860

8,400
0

1,000

2,000

3,000

3,600

0

1,000

2,000

3,000

Frequency [Hz]550
80
Convergence
Mean of StD
Mean ± 5%

40

0

1,000

2,000

Sample size n [-]

3,000

2,000

3,000

0

1,000

2,000

3,000

1,200

450

1,000

400

800

350

1,000

1,400

500

60

0

0

1,000

2,000

Sample size n [-]

3,000

600

Sample size n [-]Figure 3. Convergence analysis of mean (upper panel) and standard deviation (lower panel) of the three firstmodes over the number of samples used in the MCS.

3RESULTS

3.1 Dispersion of the natural frequencies fn
To obtain the undamped natural frequencies of the system the following eigenvalue problem was solved(λ [M] − [K])Ψ = 0,

(9)begin√λ the eigenvalues and Ψ the eigenvectors of the system. The natural frequencies are then calculated as
fn = λ /2π. Then, the Equation 9 was solved for the n samples in the MCS giving n outcomes of each natural
frequencies.
Figure 4 groups the results presenting the distribution of the first three natural frequencies of the system interms of histograms. The first three natural frequencies of the system present a distribution similar to gaussian
distribution. It is expected since the input parameters were statistically modeled by gaussian PDFs. For the
first mode, the mean over 1500 samples is 890 [Hz] with StD of 78 [Hz], what represents 8% of the mean.
On the other hand, for the second and third modes the StDs increased to 11% and 13% of the mean values,
respectively. The fourth mode falls out of the frequency range of interest, being the mean equal to 32 [kHz]with a StD of almost 15% of the mean. These results suggest that the uncertainty of natural frequencies of thehuman middle ear increases along the frequency.
At last, it should to be noted that the natural frequencies found here are related to the one dimensional movement of ossicular chain, named in the literature as “piston-like movement”. Due to the ossicular chain complexity, other resonances can be found in the in the same frequency range we analyzed. Homma et al. [11] observed
through a 3D Finite Element model the first mode close to 1 [kHz] (similar to our result), but the second mode
at 1.7 [kHz], not observed in our results.
3.2 Dispersion of the Hfp , ZME and ER
The measures computed with the deterministic model and exposed in Figure 2 were then computed using the
MCS. Figure 5 shows the results of Hfp , ZME and ER by means of minimum (min) and maximum (max) values
at each frequency and the mean ± StD. Also, these results are compared to experimental mean ± StD from
Aibara et al. [6], Voss et al. [7] and Shahnaz et al. [8].

5430

100
50

150
Occurrence [-]150
Occurrence [-]Occurrence [-]

150

100
50

0
800

1,000 1,200

2,000

Frequency [Hz]

50
0

0
600

100

4,000

6,000

1

0.5

Frequency [Hz](a) Mode 1: µ=890 [Hz] and σ =78 [Hz]Frequency [Hz](b) Mode 2: µ=3755 [Hz] and σ =437 [Hz]

1.5
·104

(c) Mode 3: µ=8568 [Hz] and σ =1106 [Hz]Figure 4. Histograms (with 40 bins each) of the first three natural frequencies of the system with n = 1500,
being µ the mean and σ the StD. These results consider a StD of 20% of the input parameters in the MCS.
In general, the max and min results have a range larger than the experimental mean ± StD reference for the
Hfp and ZME . However, the theoretical mean ± StD of Hfp is much smaller than the experimental reference.
On the other hand, the theoretical mean ± StD of the ZME is comparable to the experimental reference below
1 [kHz], but has diverged from the reference above this frequency. In addition, the theoretical statistical resultsof the ER show the variation of the middle ear parameters has a minor influence on the ER dispersion.
Lastly, the experimental data of Hfp , ZME and ER come from different samples of temporal bones/subjects.
Therefore, it is not possible to conclude what is the greater dispersion among Hfp , ZME and ER based on thoseexperimental data.

10−2

100

1k

10k

Energy Reflectance [-]10−1|ZME | [Ns/m5 ]|Hfp | [mm/s/Pa]

1010

109

108

107
100

1k

10k

1

0.5

0
100

1k

Frequency [Hz]Frequency [Hz]Frequency [Hz]

(a)

(b)

(c)

10k

Figure 5. The dispersion of the |Hfp |, |ZME | and ER by means of max and min values () and mean ±), being () the experimental mean ± StD from [6] in (a), [7] in (b) and [8] in (c). These results
StD (consider a StD of 20% of the input parameters and n = 1500 in the MCS.
Once the outcomes in terms of StD of our probabilistic model showed to be considerably smaller than theexperimental data, we increase the dispersion of the input parameters. Thus, the StD of the parameters which
compound [M], [C] and [K] were increased from 20% of the mean to 30%, 40% and 50%. Figure 6 shows themean ± StD of the |Hfp |, |ZME | and ER varying the StD of the input parameters in the MCS. The theoreticalresults are once more compared to the experimental data from Aibara et al. [6], Voss et al. [7] and Shahnaz
et al. [8]. The increment of input parameters’ StD from 20% up to 50% made the dispersion range of |Hfp |5431

similar to the experimental reference below 300 [Hz] and above 2 [kHz]. On the other hand, the range of mean
± StD of |ZME | diverged at high frequencies while the StD of input parameters increased. In the end, the ER
dispersion has not presented relevant increment above 500 [Hz] increasing the StD of the input parameters.

10−2

100

1k

10k

Energy Reflectance [-]10−1|ZME | [Ns/m5 ]|Hfp | [mm/s/Pa]

1010

109

108

107
100

1k

10k

1

0.5

0
100

1k

Frequency [Hz]Frequency [Hz]Frequency [Hz]

(a)

(b)

(c)

10k

Figure 6. Mean ± StD of the |Hfp |, |ZME | and ER varying the StD of the input parameters in the MCS, being() StD of 50%, () 40%, () 30%, () 20% and () the experimental reference ([6] in (a), [7]in (b) and [8] in (c)).

4FINAL REMARKS AND FUTURE INVESTIGATIONSIn this work we develop a model of the human middle ear using a probabilistic approach for studying middle
ear uncertainties. A coupled lumped parameter and analytical model of the middle and outer ear was combined
to Monte Carlo simulation and four different measures were evaluated. As expected, the first three natural
frequencies of the middle ear are distributed with gaussian shape for gaussian ensembles of input parameters.
Moreover, the uncertainties of the natural frequencies increase with the mode order.
In Section 3.2 we showed the model did not reproduce the dispersion of Hfp , ZME and ER seen in the experimental data assuming the StD of the input parameters as 20% of the mean with Gaussian distribution. Twohypotheses can be stated: 1. the dispersion of the input parameters was underestimated; and 2. the chosen
Gaussian PDF can not be representative for one or more input parameter. The first hypothesis was tested here
rising different remarks for the three measures analyzed. For the Hfp , the increment of input parameters’ StD
made the model outcomes closer to the experimental data in low and high frequencies, but not around 1 [kHz].
On the other hand, the dispersion of ZME became higher than the experimental reference above 1 [kHz] fora input parameters’ StD of 40% and 50%. In other words, the increase of the input parameters’ StD resultin a non proportional increase of the model output’s StD along the frequency. These result suggest the input
parameters should hold different StD in the MCS, once we assumed the StD percentage equal for the all input
parameters. Furthermore, the ER showed a low sensitiveness to the middle ear parameters variation. Also, ER
grater than 1 could be seen in some cases, what is not expected. At last, our probabilistic model was based ina deterministic model fitted to one specific experimental result. This approach lead to a biased statistical result
once the shape of the dispersion is similar to the deterministic case responses.
Although our probabilistic model developed in this work to study middle ear uncertainties has presented some
limitations, it helped to rise some interesting questions for future investigations towards a better probabilistic
approach. Therefore, some topics are intended to be assessed in the next steps, such as: 1. to use other PDF
to model the distribution of the input parameters; 2. to determine the suitable StD of each input parameter
separately; 3. since the present probabilistic results seemed to be biased by the deterministic model, to basethe probabilistic model in more than one deterministic model; 4. to estimate experimentally the distribution of5432

the natural frequencies; and 5. to develop a more realistic model of the ear canal considering its curvature anddissipation in order to explain the values of ER grater than 1.

Acknowledgment
The authors would like to thank CNPq (National Council for Scientific and Technological Development) and
FINEP (Funding Authority for Studies and Projects) for the financial support.

REFERENCES[1] D. De Greef, J. A. Buytaert, J. R. Aerts, L. Van Hoorebeke, M. Dierick, and J. Dirckx, “Details of human
middle ear morphology based on micro-ct imaging of phosphotungstic acid stained samples,” Journal ofMorphology, vol. 276, no. 9, pp. 1025–1046, 2015.
[2] P. K. Gottlieb, Y. Vaisbuch, and S. Puria, “Human ossicular-joint flexibility transforms the peak amplitude
and width of impulsive acoustic stimuli,” The Journal of the Acoustical Society of America, vol. 143, no. 6,
pp. 3418–3433, 2018.
[3] A. D. Paolis, M. Bikson, J. T. Nelson, J. A. de Ru, M. Packer, and L. Cardoso, “Analytical and numericalmodeling of the hearing system: Advances towards the assessment of hearing damage,” Hearing Research,
vol. 349, pp. 111 – 128, 2017. Noise in the Military.
[4] M. H. M. Esser, “The mechanism of the middle ear: Part ii. the drum,” The bulletin of mathematicalbiophysics, vol. 9, pp. 75–91, Jun 1947.
[5] B. Feng and R. Z. Gan, “Lumped parametric model of the human ear for sound transmission,” Biomechan
Model Mechanobiol, no. 3, pp. 33–47, 2004.
[6] R. Aibara, J. T. Welsh, S. Puria, and R. L. Goode, “Human middle-ear sound transfer function and cochlear
input impedance,” Hearing Research, vol. 152, no. 1, pp. 100 – 109, 2001.
[7] S. E. Voss, J. J. Rosowski, S. N. Merchant, and W. T. Peake, “Acoustic responses of the human middle
ear,” Hearing Research, vol. 150, no. 1, pp. 43 – 69, 2000.
[8] N. Shahnaz, K. Bork, L. Polka, N. Longridge, D. Bell, and B. D. Westerberg, “Energy reflectance andtympanometry in normal and otosclerotic ears,” Ear and Hearing, vol. 30, no. 2, pp. 219–233, 2009.
[9] J. T. Cheng, M. Hamade, S. N. Merchant, J. J. Rosowski, E. Harrington, and C. Furlong, “Wave motion
on the surface of the human tympanic membrane: Holographic measurement and modeling analysis,” The
Journal of the Acoustical Society of America, vol. 133, no. 2, pp. 918–937, 2013.
[10] Daouk, Sami, Louf, François, Dorival, Olivier, Champaney, Laurent, and Audebert, Sylvie, “Uncertainties
in structural dynamics: overview and comparative analysis of methods,” Mechanics & Industry, vol. 16,
no. 4, p. 404, 2015.
[11] K. Homma, Y. Du, Y. Shimizu, and S. Puria, “Ossicular resonance modes of the human middle ear forbone and air conduction,” The Journal of the Acoustical Society of America, vol. 125, no. 2, pp. 968–979,
2009.

5433