Recent advances in periodic vibrations of the middle ear with a floating mass transducer

The paper investigates periodic solutions of a nonlinear model of the middle ear with a floating mass transducer. A multi degree of freedom model is used to obtain a solution near the first resonance. The model is solved analytically by means of the multiple time scales method. Next, the stability of obtained periodic solutions is analysed in order to identify the parameters of the floating mass transducer that affect the middle ear dynamics. Moreover, some parameters of the middle ear structure are investigated with respect to their impact on obtained periodic solutions.


Introduction
The middle ear is the smallest and one of the most complicated biomechanical structures in the human body, therefore its modelling is especially demanding and difficult. The human middle ear (HME) is usually modeled by the finite element method (FEM) or as a multi-degree of freedom system of lumped masses, the latter is more difficult to perform but gives deeper inside into the system dynamics than FEM. Lumped mass models of 3 [18], 4 [15] and even 6 [7] degrees of freedom (dof) are usually built to analyse sound transfer through the human ear. Most lumped masses models are solved numerically, but some papers offer analytical solutions by the multiple time scale method (MTSM) can be found, e.g. [21]. It is practically impossible to find in the literature an analytical solution for a system having more than 3dof. However, numerically analysed multi-degree-of-freedom systems are used in various applications. For instance, to predict the stability of an autonomous vehicle platoon system, the authors base on a nonlinear 5-degree-offreedom vehicle model [23]. A 5dof spatial medical manipulator is analyzed for path planning [16]. Owning to unusual joint configurations of the considered manipulator it is very difficult to analytically obtain closed-form solutions. To overcome these shortcomings of the analytical method, a new approach called an adaptive neuro-fuzzy inference system is proposed. In the article [8], a physical prototype and a specially designed control system for a new 5dof hybrid robot manipulator are developed. The mechanical structure, kinematics, dynamics, and control system of this robot manipulator are presented. A new set of closed-form solutions is proposed to address the inverse kinematics problem of the robot manipulator. A 5dof model is also presented in [2] to describe an automated system designed to study the complete flexibility functions of the knee in vitro. A rigid shaft supported by a pair of angular contact ball bearings is modeled as a 5dof system as well [1]. The model simulates an existing precision grinding machine tool spindle. Another application of a 5dof model is presented in [5,6], where a mathematical model is analysed to determine equilibrium and associated load distribution in rolling element bearings. Other applications of multi-degreeof-freedom model can be found in [24,25].
A multi-degrees of freedom model must be used for both the intact ear and the middle ear with passive and active hearing aid devices. An active device is implanted to the human ear to cure sensorineural hearing loss. Then, an implantable middle ear hearing device (IMEHD) is placed into the intact ossicular chain of the middle ear. Therefore, at least a 3dof model of the middle ear and a 2dof system of a transducer as an active element of the implant are required to describe the implanted middle ear (IME). The human middle ear is composed of three bones that are the smallest in the human body, i.e. the malleus, the incus and the stapes. The bones are connected to each other and to the temporal bone by means of ligaments and tendons. In contrast, a typical IMEHD consists of three parts: a microphone, a signal processor and an vibrating output transducer called the floating mass transducer (FMT). The transducer is the main active part of the IMEHD and provides a direct mechanical stimulation to the human middle ear [9,12] or directly to the cochlea by driving the round window [3,12]. Clinical studies on the IMEHD's application have been reported by several authors [13,14,27]. FEM analysis is also applied to investigate the IME [11,12,26]. Undoubtedly, there is no mechanical lumped mass model that could explain the implanted human middle ear behaviour. Only several papers that treat the middle ear as a lumped mass structure have been published [20,22]. The paper [22] shows a simple linear model without deep investigation, while the model proposed in [20] is more complex and nonlinear. The paper [20] focuses on irregular vibrations of the middle ear bones and the transducer elements. In light of the above, the present paper uses the model taken from [20] to explore periodic solutions under various excitation conditions and system parameters.
The paper is organized as follows: Sect. 2 presents a 5dof nonlinear model of the human middle ear with the FMT. In Sect. 3 the multiple time scales method is employed to obtained approximate periodic solutions for the IME model, near to the main resonance. Next, the stability of periodic solutions is determined. In Sect. 4 shows resonance curves obtained by means of the multiple time scales method and numerical verification of the results are compared. Then, nonlinear parameters and damping coefficients are investigated with respect to solution stability. Finally, Sect. 5 provides discussion of the results and conclusions.

Model of the implanted human middle ear
The middle ear with the FMT, shown in Fig. 1a is for brevity called the implanted middle ear (IME). A detailed description of the model can be found in [20]. However, in this paper more important features of the IME are presented for clarity. The middle ear is represented by three lumped masses: the malleus (m M ), the incus (m I ) and the stapes (m S ), that are connected to each other by means of joints (IMJ, ISJ) and to the temporal bone by means of ligaments (AML, PIL, AL). Damping and stiffness properties of the elements are denoted as c and k, respectively. Moreover, the middle ear is attached to the tympanic membrane on the left which is modelled as k TM and c TM . On the right, the ossicular chain is connected to the cochlea, its properties being denoted as k c and c c ). It should be observed that the AL has nonlinear stiffness characteristics, as reported in [10,20]. The FMT is the main executive element of the IMEHD and is usually connected to the incus (Fig. 1).
The FMT is composed of a floating mass (magnet, M m ) and a metal case (M c ). The floating mass, suspended with dampers (c m ) and silicon springs (k m ), and it is excited by electromagnetic field. As a result, the mass (M m ) is moved by an external force with an amplitude P and a frequency x, as shown in Fig. 1. The FMT is fixed to the incus long process by means of a coupler (clip). Its damping and stiffness coefficients are denoted as c CLIP , k CLIP and k CLIP3 . The coupler's stiffness is assumed to be nonlinear because its design can easily be adopted to individual patient's needs. A typical structure of the FMT consists of a silicon rubber suspension of the magnet (M m ) which is also nonlinear as reported in [4,28]. Therefore, to describe elasticity of the magnet suspension, a third order polynomial is used with the coefficients (k m and k m3 ). Thus, the governing differential equations of the IME system in the dimensional form are defined as: where: To simplify a forthcoming analysis, the dimensionless time s, the frequency X and the coordinates x 1 À x 5 are introduced according to following expressions: Finally, the dimensionless equations of motion take the form: where the dimensionless parameters are defined as follows: Pcos t FMT Fig. 1 The 5dof model of the human middle ear with FMT The presented model has already been used in some papers [19,20], but only a numerical analysis of the clip fixation [19] and an investigation into irregular behaviours of the ear and the implant [20] have been done.

Multiple time scales method analysis
In order to find a primary resonance curve, the MTSM is used. Since the nonlinear 5dof system needs long and time-consuming calculations, the second order approximation is performed here. To obtain an analytical solution of Eq. 4, the following procedure is applied.

Transformation to quasi-normal coordinates
First, a linear eigenvalue problem is solved. Considering only linear terms and neglecting damping and excitation in (4) we get Next, the system is transformed from physical X to normal g coordinates where X ¼ ½x 1 ; x 2 ; x 3 ; x 4 ; x 5 T is the column matrix of physical coordinates and g ¼ g 1 ; g 2 ; g 3 ; g 4 ; g 5 ½ T the column matrix of normal coordinates, while u is a modal matrix defined by Eq. 19. Some of the extensive transformations of the equations are moved to the appendix. Then, the modal vectors are normalised to be M-orthonormal, thus where d ij is the Kronecker delta, i; j ¼ 1; 2; 3; 4; 5. Transforming the nonlinear system to normal coordinates by substituting Eq. (7) into Eq. (4) and taking into account the orthogonality in condition (8) we get: where the matrixes M; C; K; F NL ; F e are defined by Eq. 20. Now, Eqs.(9) have the form where P i are the modal amplitudes of external excitation. Expanding Eqs. (10) we get an ODE of motion in the normal coordinates (Eq. 23) Damping can be decoupled to be independent for every single mode only, if In our case, After transforming to normal coordinates and taking M-orthonormal eigenvectors, we get the modal damping matrix (Eq. 26) where the coefficients of the modal matrix f 1 . . .f 5 are selected on the basis of the original system response in order to represent its naturally damped oscillations. Assuming that the system is weakly nonlinear, we introduce a formal small parameter, e, keeping all small terms.

Approximate solution of the first resonance
Approximate analytical solutions of Eg. 23 are determined by the MTSM [17] presented in details in Appendix 2. Finally, the steady state solution takes the form where the first column of the modal matrix (Eq. 19 is renumbered as follows: u Eqs. 13 can be used to perform a bifurcation analysis around the first natural frequency and resonance curves. It is worth stressing that the solution obtained by the modulation Eqs.(13) depends on selected system parameters, i.e., the excitation amplitude p and frequency X, the modal damping f 1 , the linear stiffness of the system (modal vector u) and the nonlinear parameters c 3 and c 45 .

Stability analysis
The stability of periodic solutions is determined based on the Jacobian matrix (J) of the modulation equations (13) in the form Generally, stability of the solution is determined by the sign of real part of eigenvalues of the Jacobian matrix (J). When the real part of eigenvalues is negative, the solution is stable. However, the eigenvalues of the Jacobian are long and mathematically complicated, therefore the stability is defined with the help of the trace (Tr(J)) and determinant (Det(J)) of the Jacobian are defined as The solution is stable when DetðJÞ [ 0 and TrðJÞ\0. Substituting Eq. 17 into Eq. 15, the trace finally takes the form In the present case the trace is always negative (TrðJÞ\0), then the stability depends only on Det(J).
To check the solution stability, first the amplitude (a) must be calculated and next Det(J). This analysis is performed numerically in subsequent sections of this paper.

Resonance curves
To obtain a resonance curve for the stapes (a 3 ), the amplitude (a) calculated from Eq. 13 should be multiplied by u from the modal matrix Eq. 19. Then, the resonance curves near the first (main) resonance are analysed for the system parameters listed in Tab. 1. It should be noted, that two sets of damping coefficients are presented: normal (c Ã ) and pathological (c 1Ã ). The former is typical of the intact (healthy) ossicluar chain while the latter describes the pathological ear with decreased damping properties. However, damping does not affect the modal matrix (Eq. 19), the value of which are calculated on the basis of system stiffness. Since, it is only the first resonance that is analysed, only the first modal vector is necessary to obtain the solution. Obtained values of this vector are as follows u 11 ¼ 0:591082; u 21 ¼ 0:5 91289; u 31 ¼ 0:376607; u 41 ¼ 0:709169; u 51 ¼ 0: 8 6 1285. Resonance curves of the normal and pathological IME obtained for three variants of the external excitation p (1p, 2p, 5p), are shown in Figs. 2 and 3, respectively. The black line denotes the stable periodic solution (SPS) while the unstable periodic solution (UPS) is marked with a red line. To verify the analytical results, numerical simulations were performed using the MATLAB Simulink package. The Runge-Kutta fourth-order integration method (ode45) with a relative tolerance of 1e À10 and a variable step size was applied to ensure appropriate accuracy. The numerical results are marked as blue points and stars. For the normal ear IME (Fig. 2) the stapes vibrates like in a linear system, where only the stable periodic solutions exist regardless of the excitation level, while for the pathological IME classical nonlinear resonance curves were obtained. The numerical results (marked as blue points) of the pathological IME are consistent with the approximate analytical solutions, in contrast to results of the normal IME where the compatibility is significantly worse. The presence of unstable solutions strongly depends on the excitation amplitude (p) and the excitation frequency (X). Regions of the stable (SPS) and unstable periodic solutions (UPS) as well as the number of these solutions are shown in Fig. 4. One UPS occurs when X [ 0:98 and the excitation amplitude p is strong enough, however for X [ 1:1 and a very small p the region of unstable solutions vanishes and the system exhibits only one periodic solution (1SPS).  An analysis of the resonance curves for different damping coefficients shows the importance of modal damping (f 1 ). Therefore, the effect of f 1 is plotted in Figs. 5 and 6 as a function of X and p, respectively. When the excitation amplitude is 1p this leads to the critical value of f crX 1 ¼ 13:5e À3 (Fig. 5). When f 1 [ f cr 1 , the system has only 1SPS regardless of the excitation frequency . The lower the damping (f 1 ) is, the larger the gray regions of 1UPS and 2SPS become. Interestingly, the strong excitation (p [ 4:5e À4 ) prevents the occurrence of the grey region of three solutions (Fig. 6). On the other hand, in order to eliminate the three-solution (grey) zone, near the primary resonance, the damping should be greater than f crp 1 ¼ 0:034. Thus, the final critical damping coefficient ensuring 1SPS, regardless of the excitation amplitude (p) and the excitation frequency (X), is f cr 1 ¼ 0:034. The presence of the 2SPS zone can be dangerous because this means that the system is not clearly defined. The solution amplitude depends on initial conditions.

Effect of system nonlinearity on the periodic solutions
System nonlinearities are equally important for the periodic solutions as excitation. This is true as far as c 3 is concerned (Fig. 7). It can be observed that near resonance the gray region depends on c 3 . However, when c 3 \0:1e 4 , the system has only 1SPS. It should be noted that c 45 has no effect on the stability of the solutions (Fig. 8). Nevertheless, one can distinguish a borderline marked by c 3 , where an increase in c 45 causes a transition from the white to gray region (Fig. 9). To sum up, c 3 is crucial for solution stability whereas the influence of c 45 is insignificant. Thus, the nonlinearity of the silicon rubber characteristics can be omitted in the analysed range of parameters.

Discussion and conclusions
From a medical and practical point of view, the IME system should have only 1SPS. With more solutions, Fig. 4 Regions of stable and unstable periodic solutions for the pathological IME versus excitation amplitude p and excitation frequency X  Fig. 6 Influence of damping f 1 and excitation amplitude p on stable and unstable periodic solutions for the pathological IME and X ¼ 1 the implant control becomes more difficult. For this reason, it is crucial to estimate the system parameters when there is only 1SPS. From the point of view of mechanical engineering however, finding more then one periodic solutions of the 5dof system is a challenging task. Due to strong nonlinearity, the analytical solution shows agreement with the numerical results only for relatively small damping. The nonlinear properties of the silicon rubber (c 45 ) have an insignificant effect on the stability of periodic solutions. The nonlinearities of AL (c 3 ) have a significantly bigger influence on stability. Therefore, various kinds of silicon rubber or another viscoelastic element can be applied in the transducer to connect the can and the magnet. Both strongly nonlinear and pure linear element can be used to build the floating mass transducer. That gives various possibilities for engineers to use the materials having wide range of stiffness without a loss of periodic solution stability. Nonlinearity of the FMT coupler (c 24 ) is an another key parameter of the implant but the performed analysis shows that it does not influence the periodic solution. However, it can have an effect on irregular vibrations and in consequence on the coupler construction.
System damping is also an important parameter of the system. The normal middle ear is free from nonlinear effects even under strong excitation due to relatively big damping. Thus, the use of the IMEHD is a good solution in this case. A problem arises when the ear is pathological with decreased damping. In a situation like this, there can be three periodic solutions, stable or unstable. From a scientific point of view, this poses an interesting case, but from a practical viewpoint it can be dangerous.
The multiple time scales method (MTSM) that was employed in this study to find periodic solutions of the 5dof system is characterized by a relatively complicated and time-consuming procedure, but the results are positive because they yield a closed-form solution that is especially required for to system dynamics identification. Thanks to the use of the MTSM, it was possible to determine the stability of periodic solutions, which is marked as the line constraining the grey regions in Figs. 4, 5, 6, 7 and 9. ð Appendix 2 According to the MTSM, the solution of Eg. 23 is expressed by a series of the small parameter e. g j ðs; eÞ ¼g j0 ðT 0 ; Time is also expressed as series of the small parameter in the fast and the slow scales in the following way where T 0 is fast and T 1 , T 2 slow scales of time. Thus new definitions of time derivatives are obtained Excluding an internal resonance case, the solutions are sought near the first natural frequency resonance zones. Thus, the excitation frequency can be expressed as where r i is the frequency detuning parameter. Assuming that the excitation frequency is close to the first natural frequency the following dependence can written Substituting the analytical solution (Eq. 27) into Eq. 23 and taking into account the derivatives definitions (Eqs. 29 and 30) the natural frequency according to Eq. 32, after grouping terms with respect to the order, a set of differential equations in the successive perturbation orders is defined as follows • e 0 -order • e 1 -order D 2 0 g 11 þ x 2 g 11 ¼ À2D 0 D 1 g 10 À f 1 D 0 g 10 þ ::: þ pu 51 cos xs • e 2 -order D 2 0 g 12 þ x 2 g 12 ¼ À2D 0 D 1 g 11 À 2D 0 D 2 g 10 À D 2 1 g 10 À f 1 ðD 0 g 11 þ D 1 g 10 Þ þ ::: where i=2,...,5.
Because the internal resonance, near the first natural frequency are excluded, the coordinate plays dominant role. Therefore, the response of the coordinate is taken in higher perturbation orders. Then, the solution of Eq. 33 can be expressed in the form where i ¼ ffiffiffiffiffiffi ffi À1 p is the imaginary unit, A 1 is complex amplitude and A 1 is its conjugate. The solution Eq. 36 is substituted into Eq. 34 and, after grouping terms in proper exponential functions, one can get Equalling ST 1 terms to zero, particular solutions of Eq. 37 are determined. Substituting the solutions into equations in the second order approximation (Eq. 35) one can get D 2 0 g 12 þ x 2 g 12 ¼ ST 2 e ix T 0 þ NST 1 þ cc ð39Þ NST 1 and NST 2 represent nonsecular generating terms. To avoid secular terms in the solution of Eq. 39 vanishing secular generating terms ST 2 is required, therefore the term with e ixT 0 should vanish thus the following expression is obtained where, AT 1 means the additional terms that are too long to show them.
By means of so called the reconstitution method, Eq. 38 is reformulated into the first order differential equation for complex amplitude A 1 . The amplitude derivative with respect to the dimensionless time is define as Derivatives À e 2 pu 51 8x where, AT 2 means the additional terms that are too long to show them. Expressing complex amplitude A 1 in the polar form then substituting it into Eq. 43 and separating real and imaginary parts, modulation equations for the amplitude ( _ a) and the phase ( _ /) of the stapes motion are obtained. In a steady state _ a ¼ 0 and _ / ¼ 0, then the amplitude (a) and the phase (/) can be calculated.