Abstract

A review on mathematical and mechanical models of the vocal cords is given. The basic model is a two-mass nonlinear oscillator system which is accepted to be the basic one for mechanical description in voice production. The model is not only extended into three, five, and more mass systems, systems with time variable parameters and three-dimensional systems, but also simplified into one-mass system with coupled two-direction deflection and damping functions. The corresponding mathematical models are the systems of coupled second-order differential equations which describe the vibrations of the symmetric and asymmetric vocal folds. The models give the conditions for the regular and irregular motions like bifurcation and deterministic chaos in vocal folds. The obtained results are of special interest for detecting the pathology of vocal cords, when there are no visual effects of disease. Based on the results given in the paper, the objectives for future investigation in this matter are given.

1. Introduction

For a long time, the researchers are trying to simulate the human voice production. Various mechanical and mathematical models are developed not only for describing of the human organs which are connected with voice production, but also of the process of phonation. In essence, voice is produced by movement of the two lateral opposing vocal folds located in the larynx (see textbooks in Otolaryngology, e.g., [1]) caused by the air flow through the trachea and generated by lung.

The vocal cords, also known as vocal folds, represent the part of the mechanism for sound generation. The vocal cords are composed of twin infoldings of mucous membrane stretched horizontally across the larynx. Open during inhalation, closed when holding one’s breath, and vibrating for speech or singing; the folds are controlled via the vagus nerve. The vocal cords are brought near enough together such that air pressure builds up beneath the larynx. The cords are pushed apart by this increased subglottal pressure with the inferior part of each cord leading the superior part. Under the correct conditions, this oscillation pattern will sustain itself. Modulating the flow of air being expelled from the lungs during phonation, the vibration of the vocal cords appears. Rhythmic opening and closing of the vocal folds, that is, their oscillation in the larynx, causes chopping up of a steady flow of air around the glottal into little puffs of sound waves. The studies of vocal fold biomechanics give an insight into voice production and also provide important information about laryngeal pathology development. Detecting of pathology of voice production system when there is no visual evidence for morphological laryngeal abnormalities is of special interest.

Complexity of vocal folds, their histology, shape, position, and so forth, give us a possibility to treat the problem in quite different manner. Same ability is evident for modeling of the process of phonation. It is the reason that a numerous aspects of the problem are investigated, and a great number of results are published (more than 1000).

This paper gives a review only of the models which are suitable for analyzing of the vocal fold vibration separated from the vocal tract dynamics. We are aware that neglecting the influence of the vocal tract and subglotal resonances is a strong simplification of the voice source, but it gives us the possibility to analyze extremely reduced equations describing solely the vocal fold dynamics.

The intention of the paper is to give the reviews of the vocal fold models but also to give objectives for the future investigation in this complex matter by including the nonlinear properties of the system.

2. One-Degree-of-Freedom One-Mass Model of the Vocal Fold

One of the most simple models of the vocal cord is an one-mass system with one degree of freedom which was introduced by Flanagan [2] and after that investigated by a significant number of authors (see Mermelstein [3], Titze [46], Cronjaeger [7], etc.).

Flanagan [2] used the electroacoustics analogy to describe the dynamics of the vocal cord. He gave the model of the cord as an one-mass oscillatory system with damping, where the self-excited vibrations are caused by the flow of the air with the variable pressure. The model is shown in Figure 1. The excitation force is obtained experimentally and included into the model where is the velocity of air flow. The model is a linear one but is useful for explanation of the basic oscillatory properties of the vocal cord. Unfortunately, the model was not able to explain the mechanism of the vocal fold excitation.

Fulcher et al. [8] exceed this problem by including of the negative Coulomb damping into the model. They stated that the action of the aerodynamic forces on the vocal folds is captured by negative Coulomb damping which causes the vibration of the vocal cords. Their model is based on the idea that the oscillation of vocal folds maintains their motion by deriving energy from the flow of the air through the glottis. The aerodynamic forces which act on the vocal cords energize them by negative Coulomb damping. This force adds energy to the oscillator instead of removing it. A viscous force is added to include the effects of the tissue damping. Adding a viscous damping term makes steady-state motion possible. In the long-time limit the analytical solutions approach a limit cycle, and the amplitude and velocity lose their dependence on the history of the motion. Due to combination of these two forces, the limit cycle is reached. The mathematical model considered in the paper is where is the negative Coulomb force, is the lung pressure, and is the medial surface of the glottis. An elevated lung pressure gives rise to a flow of air through the glottis and produces a series of alternating converging and diverging shapes of the vertical dimensions of the vocal folds. The pressure distributions in the glottis resulting from the series of shapes are alternately higher and lower than the pressures in the vocal tract. These pressure variations are in phase with the motion of the vocal folds and add energy to the oscillator in the same way as negative Coulomb damping does. Limit cycle of the oscillator with negative Coulomb damping provides a natural explanation of the self-oscillation property of the model.

It is worth to say that the simple, one-mass vocal fold model with one degree of freedom is useful for explanation of the connection between the kinematic features of the vocal cord movement and the muscle stiffness [9] and in connection with vocal tract dynamics [10], too.

3. Two-Degrees-of-Freedom One-Mass Model of the Vocal Fold

Adachi and Yu [11] extended the one-mass-model of the vocal cord assuming that the vocal folds can vibrate both parallel and perpendicular to the air flow. The model is with one mass but with motions in two orthogonal directions and (Figure 2).

As the mass moves, the parallelogram, which models the vocal cord, is deformed and the vocal fold simultaneously executes both swinging and elastic motions. The motions are coupled and mathematically described with a complex function with displacements in and direction, and is the imaginary unit.

If the left and the right vocal folds vibrate symmetrically, the mathematical model of the oscillating of each system is as follows: where is the resistance coefficient with the quality factor . The authors introduced the factor (1/2) in the differential (3.1) implying that the acceleration and velocity of the center of the vocal fold are half as much as those of the edge of the vocal fold.

The forces which act on the vocal cords and are considered in the paper are as follows: where is the pressure in the glottis, is the subglottal pressure, is the pressure at the entrance of the vocal tract, is the complex conjugate function to , is the complex function of the rest position with displacements and in two orthogonal directions, and are the coordinates originating from the rest position (,), and is a coefficient representing the nonlinear elastic property. The contact force is generated during closure of the glottis and is a linear deflection function with constant . The elastic force is assumed to be the cubic displacement function. This type of nonlinear Duffing model (3.1) with (3.2) is widely discussed in the papers of Cveticanin [1217].

In the papers [18, 19] the model (3.1) with (3.2) is modified by assumption that the elastic property of the vocal cord differs in and direction (Figure 3). For the horizontal and vertical displacements and , the system of differential equations is where the change of the air pressure satisfies the relation is the damping coefficient of the vocal cord, is the coefficient of the additional damping, and are a horizontal and vertical stiffness coefficient of the vocal cord, is a stiffness of the coupling between the two directions of motion, is the coefficient of the cubic stiffness term, is a hyperbolic type stiffness, is the unloaded equilibrium position (), is the average pressure, and is air flow and . The system of (3.3) describes the self-excited vibrations of the vocal cords. For the interval of damping considered, it has been shown that this system possesses two Hopf bifurcation points and three period-doubling bifurcation points. There is an interval of damping for which two stable periodic and periodic solutions exist. It is possible to jump from one solution to another. Harmonic and subharmonic unstable solutions also exist for this interval of damping.

Both of the two-dimensional one-mass motion models of the vocal fold can successfully simulate the self-excited oscillation in a wide frequency region.

4. Basic Two-Mass Model of the Vocal Fold

Another type of the vocal cord motion model is introduced by Ishizaka and Flanagan [20]. The model is based on the assumption that the energy transfer from the air flow to the vocal cord is generated by phase difference between the oscillations , the lower, and upper edge of the vocal cord, where the vibrations are driven by lung pressure. This effect can be modeled by representing each fold by two coupled oscillators (Figure 4). For normal voice, the oscillators are bilaterally symmetric. The displacing tissue of each cord is considered to be approximated by two mass-damper-spring systems ( and ), coupled through a spring . The two vocal folds are assumed identical, and they move symmetrically with respect to the glottal midline. The mathematical model of the system is given with two coupled second-order differential equations where is the horizontal displacement measured from a rest (neutral) position , is the driving vibration force produced by the air flow, and . The damping and elastic properties of the vocal cord tissue are assumed to be the linear functions of the velocity and displacement, respectively.

The mass is permitted displacement in both the lateral () direction and the longitudinal () direction. Motion in lateral direction is opposed by restoring stiffness and viscous damping. Longitudinal displacement of masses and is conditioned by the internal coupling stiffness which permits realistic phase differences in the lateral displacements.

Vocal cords generate a sound rich in harmonics [21]. One of the generators of the harmonics is the collision of the vocal cords between themselves. It required the inclusion of the collision effect of the vocal cords into the mathematical model (4.1). The mathematical model is [22] where is the rigidity coefficient during the collision between two sides of the vocal cords, and denote the lower and upper glottal areas, respectively. For the symmetric model, the left side of vocal folds is the same as the right side, so that the glottal areas can be given as where is the glottal rest area. The collision function is assumed as The driving forces are expressed as with where is the pressure inside the glottis acting on the upper and lower masses of vocal cords, is the steady pressure component in Bernoulli region, and is the pressure in the turbulent jet region. Analysing the relation (4.5), it is evident that the tissue collision during phonation produces a very large impact pressure which correlates with the lung pressure and glottal width . Larger lung pressure and a narrower glottal width increase the impact pressure. Namely, it is known that by the narrowest glottal gap, the effects of airflow are separated into two regimes: Bernoulli regime and jet regime (Figure 5). Any point upstream, that is, the region from the subglottal region to the minimum glottal diameter satisfies the steady Bernoulli equation. However, downstream the region from the minimum glottal diameter to the glottal exit fully developed turbulence generates numerous vortices at the glottis exit. The masses in the turbulent region are perturbed by extremely irregular fluctuations. The turbulence in the unsteady jet region acts as an external random force to drive vocal folds, and the additional pressure perturbation in this region is a random force. The pressure in the turbulent jet region is different from zero. Dependently on the type of air flow, in (4.6) is 0 for Bernoulli flow and 1 for turbulent jet flow of the air. The results given in [23] demonstrate that this excitation causes irregular vocal fold vibrations and inhibit fine period-doubling bifurcations. Finally, it is concluded that the bifurcations are closely related to observations in voice pathology.

During the last forty years, the two-mass model evolved to a standard for exploring the voice producing system. The model of vocal cord represents a self-excited oscillatory system which includes the basic requirement of the phase difference of the lower and the upper mass caused by varying of the area of the vocal cords due to air flow (see, e.g., [24]). The induced phase difference of the upper and lower mass enables the energy transfer from the air stream to the vocal folds [25]. Due to its simplicity, the model was convenient for application and explanation of the phenomena which appears in voice production [26, 27], for connection with dynamics of vocal tract [28] and for detection of various anomalies of the vocal cord.

Remarks
Suggested mechanical model is convenient not only for voice folds vibration, but also for their posturing. Namely, vocal fold dynamics for phonation can be treated in two parts [29]: (i)large and relatively slow deformation occurring when the vocal folds are positioned for voicing, coughing, and breathing—this part is referred as vocal fold posturing which is subdivided into:(1)adducting or abducting the medial surfaces of the vocal folds,(2)elongating or shortening of the vocal folds;(ii)small and relatively fast deformations occurring when the tissue is driven into self-sustained oscillation by aerodynamic and acoustic pressures—this part is referred as fold vibration.
 The posturing occurs in a nonperiodic but ultimately always cyclic fashion at small frequencies, but the vocal fold vibration occurs at higher frequencies. Although vocal fold posturing and vocal fold vibration are often thought to be separate mechanical processes, many parameters of vibration (e.g., fundamental frequency, amplitude of vibration, and voice onset time) are dependent on posturing. Analyzing the model simulating the adduction of the medial surface of the vocal cord during posturation, it is concluded that it affects the intensity of vocal folds vibration and involves into fundamental frequency regulation [30]. This result is previously obtained by clinical observation [31], too.
The characteristics of all mechanical elements are based upon dates obtained by clinical observations (see [3234]) and measures [35, 36].

5. Nonlinear Two-Mass Model of the Symmetrical Vocal Folds

Lucero and Koenig [37] extended the basic two-mass model of the vocal fold (4.1) by introducing the nonlinear properties of the vocal fold tissue.

The two vocal folds are assumed identical and to move symmetrically with respect to the glottal midline. When the glottis is open, the equations of motion for a vocal cord may be written as where is the horizontal displacement measured from a rest (neutral) position and . In the work, they introduce the cubic characteristic for the tissue elastic forces and for the damping forces, instead of the usual linear term , the nonlinear characteristic of the form The reason was the need to limit the amplitude of the vocal fold oscillation as the glottal width increased.

For the case when there is a contact between the opposite vocal folds, at the rest position, the both masses are at a distance from the glottal midline. Then each, mass collides with its opposite counterpart at a displacement . During contact, the stiffness is increased and the damping coefficient is increased by adding 1 to the damping ratio: if is the value of the damping coefficient when the glottis is open, for a given damping ratio , then, during contact of mass , it assumes the value.

This low-dimensional vocal fold model can reproduce the vocal fold vibration onset-offset patterns observed experimentally during speech.

6. Two-Mass Models of the Asymmetric Vocal Folds

Isshiki et al. [38] discussed the clinical significance of asymmetrical vocal cord tension. Asymmetry may signify the pathology of the vocal cord. Due to pathological changes of the mechanical properties, and asymmetry of the vocal folds instabilities may appear which lead to qualitative different sounds being produced in the larynx. The dynamic model has to include the effect of a tension imbalance between the two vocal folds (left and right).

Steinecke and Herzel [39] formed the model of the asymmetric vocal folds (Figure 6) based on the well-known two-mass one (4.2). Differential equations of the motion of the upper and lower mass remain the same as previously given (4.2), but the left and the right vocal folds are not symmetric and have various motion. The following four second-order differential equations are suitable for consideration where For the derivation of the driving forces , it is usually assumed that the air jet separates at the point of minimum area , inducing an immediate pressure drop to zero. It causes that the force is zero for all glottal configurations. On the lower part of the vocal cord, the force (4.5) acts where is the pressure which can be obtained according to the Bernoulli law and is the glottal volume flow. For the pressure is Using the asymmetric vocal fold model (6.1), it could be seen that for certain parameter values, irregular motion and bifurcation occur. The different stiffness values of the left and right vocal fold and the overcritical tension imbalance result in diverse oscillatory patterns of the vocal folds like: subharmonic vibrations, two-frequency oscillations (biphonation), and chaos for which a distributed or hoarse voice arises. These irregularities are often caused by anatomical asymmetries between left and right vocal fold. The correctness of the model is proved experimentally for the unilateral laryngeal paralysis [40, 41]. For a unilateral polyp, the vibrations of the vocal folds are also irregular, that is, chaotic (see [42, 43]).

Using the model (6.1), Zhang [44] investigated the phonation in asymmetric system, recently. He analyzed the vibration in a self-oscillating vocal fold model with left-right asymmetry in body-layer stiffness. In spite of the fact that the frequencies of the twofolds have a ratio close to 1 : 3, the subharmonic synchronization was not observed in the asymmetric model. Instead, the vibratory behavior was dominated by the dynamics of onefold only, and the other fold was enslaved to vibrate at the same frequency. Increasing subglottal pressure induced a shift in relative dominance between the twofolds, leading to abrupt changes in both vibratory pattern and frequency.

Comparing the properties of vibration asymmetry obtained by the model (6.1) with the measured vocal fold vibratory asymmetry [45], it can be seen that in spite of the fact that the model is a two-mass one, the results are comparable.

7. Nonstationary Two-Mass Model of the Asymmetric Vocal Folds

The conventional clinic examination of vocal fold vibrations is done during stationary, sustained phonation. However, the conclusions drawn from a stationary phonation are restricted to the observed steady-state vocal fold vibrations and cannot be generalized to voice mechanisms during running speech. Classification of vocal fold vibrations is an essential task of the objective assessment of voice disorders. The two-mass model-based approach is suitable to be used for the analysis and the quantitative interpretation of nonstationary vocal fold vibrations. In order to simulate nonstationary vibrations, the original constant parameters are modified to be time dependent [46]. In Figure 7, the two-mass model of asymmetric vocal folds with time variable parameters is plotted. The corresponding mathematical model is a system of four coupled differential equations Comparing (7.1) with (6.1), it is obvious that the form of the equations is the same, but in (7.1), the parameters are time variable, and also the reactive force is generated. Namely, for the system where mass is a time variable function, the reactive force acts which is defined as [4749] The damping coefficient is The model is suitable for study of the nonstationary phonation. The obtained solutions allow the intuitive assessment of vocal fold instabilities and are applicable to an objective quantification of asymmetries.

Nonstationary vocal fold oscillations are typical for many voice disorders. Namely, the fine adjustment of the model requires the determination of the time dependence of the asymmetry coefficient and of the parameter deviations from the standard set.

8. Multimass Model of the Asymmetric Vocal Cord

The generalization of the previously mentioned model is recently given by Yang et al. [50]. Combining the idea of one-mass system with two degrees of freedom and of the two-mass system where each of masses have one degree of freedom is extended to 5-mass system where each of the masses has three degrees of freedom: not only lateral and longitudinal (as it was the case in the two-dimensional model), but also vertical vibrations (Figure 8). The model is the most appropriate one as it describes the human voice originates from the three-dimensional (3D) oscillations of the vocal folds. The mathematical model of the 3D system with five masses on the left-hand side () and five on the right-hand side () is where denotes the position of each mass element in the Cartesian coordinate system, is the vertical coupling force, which serves as internal force of the vocal cord tissue in vertical direction, is the anchor force, which acts as the function of the thyroarytenoid muscle in lateral direction, is the longitudinal coupling force which describes actions of the thyroarytenoid muscle and the vocal ligament in longitudinal direction, is the collision action between the two vocal cords during phonation, which causes the elastic structure deformation of the vocal folds, and is the fluid driving force.

The anchor forces, which connect the masses to fixed body, and the vertical and longitudinal coupling internal forces between masses are supposed to be nonlinear deflection functions. The collision impact force is included into the model, too. The influence of the aerodynamic force, which causes the three-dimensional mass elements of vocal fold to oscillate, is investigated. The driving force is of Bernoulli type produced by the glottal flow originating from the lung and acting on the vocal folds from inferior to superior through the whole larynx. The driving force depends not only on the subglottal pressure, but also on the geometric dimensions (thickness and length of the vocal folds and on the rest positions). The resulting model enables visualization of the three-dimensional dynamics of the human vocal folds during phonation for both symmetric and asymmetric vibrations.

9. Nonstationary Multimass Model of the Vocal Fold

Dynamics of a multimass model with time-dependent parameters are matched to vocal fold vibrations extracted at dorsal, medial, and ventral positions by optimizing the time-dependent model parameters [51]. The parameters serve as an objective measure for quantitatively assessing left-right asymmetry measures of amplitude and phase shifts in the medial part of the vocal folds as well as anterior-posterior longitudinal asymmetries along a vocal fold side. The vibration behaviour at the dorsal, medial, and ventral vocal fold regions contributes information for characterizing the vibration pattern. In Figure 9, the sketch of the multimass time-dependent model showing the longitudinal and lateral vibrations for each mass in the system is plotted.

The equations of vibration of the -mass system with motion in two directions (2D) which models the asymmetric vocal fold is where the time variable mass element is denoted as with for the left-hand side and for the right-hand side. For the lower mass, it is , and for the upper mass . is the vertical coupling force which results from the connection between the upper and lower mass. As there is assumed that the motion in vertical direction is impossible, the force gives the connection between lateral and longitudinal vibrations. is the anchor force which is the function of coefficient of rigidity and damping, and is the longitudinal coupling force which is the function of rigidity and damping coefficients, while is the collision force between the two vocal cords. The reactive force which is connected with mass variation is .

The graphical representation of the time-dependent parameters gives an intuitive view of the actual degree of asymmetry. A rating value derived from the parameters provides a classification into normal and pathological vocal fold vibrations.

10. Conclusions and Directions to Future Investigation

Based on the published results, it is concluded that mathematical models of vocal cord vibrations give very good qualitative description of the phenomena, but in spite of the fact that the clinically observed and measured parameters are used for modeling, the obtained results quantitatively differ from real one. It requires the improvement of the accuracy of models. The following is suggested.(1)In the mathematical models, the nonlinear properties of the vocal cord have to be included. The differential equations would have nonlinear terms of integer or noninteger order. Not only numerical, but also approximate analytic methods have to be developed for solving these differential equations.(2)For all masses in the system, modeling the vocal fold, it is suggested to introduce the motion in three directions: lateral, longitudinal, and vertical, which is not the case at the moment. Such a model would give a better interpretation of the vocal fold vibration.(3)In the model of vocal cord, the negative Coulomb damping is assumed to be a linear one. To improve the model, we suggest the introduction of the nonlinear damping of the integer or noninteger order which is obtained empirically by clinical measuring. We believe that such models would give more accurate results. To solve the differential equations of motion, the known analytical methods have to be extended, and also new solving methods have to be developed.(4)Nonstationary vocal fold oscillations are typical for many voice disorders. Inclusion of the time variable parameters and the reactive force, which acts due to mass variation in time, gives an additional influence on the vibration of the fold. The stability of nonstationary vocal folds vibration has to be more intensively analyzed. (5)The special attention has to be paid to irregular motion of the model of the vocal cord and to qualitative and quantitative analysis of the differential equations which describe these motions. The instabilities and irregular motion of the vocal cord may signify the diseases or an anomaly.

Abbreviations

Tissue damping force
:Depth, that is, thickness of the vocal cord
:Driving force of the vocal cord
:Coefficient of rigidity of the vocal cord
:Coupling stiffness coefficient
:Length of the vocal cord
:Length of the glottis
:Mass of the vocal cord
:Pressure in the vocal tract
:Damping coefficient of the vocal cord
:Tissue elastic force
:Width of the vocal cord
:Specific mass of the air for the unit length.

Acknowledgment

This work has been supported by the Ministry of Science of Serbia through projects ON 174028 and IT 41007.