Introduction

Non-linear dynamic systems are mathematical models for several problems in physics, chemistry, biology, economy, engineering, and correlated areas of knowledge (Chavarette, Balthazar, Rafikov, & Hermini, 2009;Zhang, Liu, & Wang, 2009;Gupta, Dalal, & Mishra, 2014;2015; Ferreira, Chavarette, & Peruzzi, 2017;Mishra, Sen, & Mohapatra, 2017). These non-linear dynamic systems are usually associated with non-linear ordinary differential equations. When studying the behavior of dynamic systems that are modeled by ordinary differential equations, some typical behaviors are found, such as periodic, chaotic, and synchronized (Rafikov, Balthazar, & Tusset, 2008;Zhang et al., 2009).

For some systems originated from a physical formulation, a specific behavior is desired, and synchronization is mostly the case. Synchronization is the process by which the elements of a dynamic system start exhibiting a collective behavior. Natural systems such as fireflies blinking, cardiac pacemakers, and neurons firing are prone to synchronization (Pikovsky, Rosenblum, & Kurths, 2003).

The synchronization phenomenon was discovered in the 17th century by the scientist Christiaan Huygens (Willms, Kitanov, & Langford, 2017) when he studied the movement of two oscillating pendulums, which were hung over the same wooden support. Huygens, inventor of the pendulum clock, was observing two clocks, which he had recently manufactured, hung side by side. He noticed that, regardless of the initial position, the two pendulums oscillated in synchrony after some time; these oscillations coincided and the pendulums moved in opposite directions. He suspected, therefore, that the clocks might be influenced by each other somehow.

At first, the synchronization of periodic dynamic systems, which present a limit cycle attractor, was studied. However, in the past few decades, this concept has been extended in order to be able to understand other rhythmic adjustments processes in chaotic dynamic systems, which has proven to be even better for understanding what rhythm truly means in the context of dynamic systems (Pecora & Carrol, 1990;Heagy, Carroll, & Pecora, 1994;Rosenblum, Pikovsky, & Kurths, 1996;Rosenblum, Pikovsky, Kurths, Schafer, & Tass, 2001;Boccaletti, Kurths, Osipov, Valladares, & Zhou, 2002).

Dynamic systems do not always present a synchronized behavior. Nevertheless, some naturally non-synchronized systems may present a synchronized movement when certain techniques that force such synchronization are applied. In many problems related to engineering, ecology, among other areas, the goal is to choose a control rule that moves the system from a non-controlled or chaotic regime into a fixed equilibrium point or periodic or non-periodic (chaotic) desired orbit. To that extent, this paper aims at synchronizing the movement of vocal folds through optimal linear control (LQR), which presented satisfactory results in many cases analyzed (Berry, Herzel, Titze, Krischer, & Story, 1994;Berry, Herzel, & Titze, 1996;Giovanni et al., 1999;Giovanni, Ouakinine, & Triglia, 1999;Tusset, Baltthazar, Chavarette, & Felix, 2012;Chavarette, 2011).

Material ands methods

This section describes the material of this research and a brief physiological and mathematic description of the requirements for understanding the applied methods and achieved results.

Physiology

The human vocal apparatus is composed of the lungs, larynx, pharynx, mouth, and nasal cavity. On average, the larynx is 5 to 6 cm long in adults. Its anatomy is divided into:

- Subglottic region: the bottom portion of the larynx, which extends from the lower edge of the vocal folds up to the beginning of the larynx;

- Glottic region or glottis: the gap between the vocal folds;

- Supraglottic region: region that comprises the upper edge of the vocal folds up the end of the larynx.

The vocal folds, formerly called vocal cords, are 5 to 17 mm long. It is possible to highlight their role in the process of speech as an essential structure. They have an anatomy that allows the mucosal layers to vibrate freely over smoother structure rigid bottom layers. As Figure 1 shows, the vocal fold can be divided into three levels: mucosa, vocal ligament and vocal muscle. The mucosa is mainly responsible for vibrating in the phonation process, and can be divided into squamous epithelium and superficial lamina propria, with gelatinous texture, which move freely upon other structures below. Under the superficial lamina propria, there exist two other laminas that compose the vocal ligament. These three layers of lamina are particularly different because as they approach the vocal muscle (the stiffest layer of the vocal folds), they get stiffer and stiffer (Hirano, 1974). Therefore, this layer and its structure determine the stiffness coefficient for each vocal fold.

The vocal folds are responsible for emitting sounds in human beings and other species of animals, such as dogs. These structures, located in the larynx, vibrate due to the air pressure coming from the lungs, generating the emission of sounds, which are in turn modified according to the articulation performed by the mouth, and then amplified by the resonance chamber, composed of the larynx, pharynx, mouth and nose.

Each person has two vocal folds, and each one is composed of body and cover. The body is a muscle structure composed of the vocal muscle, whereas the cover is composed of an epithelium and the superficial lamina propria. In between them there exists a transition layer composed of both intermediate and deep lamina propria. The body stiffens up at the moment of speech, whereas the cover is flexible and thus able to vibrate.

A person's voice can be modified if the vocal folds are under the action of a disease, such as tumors, polyps, Reinke's edema, nodules, and partial or full folds paralysis. An example of unilateral paralysis can be given when one of the vocal folds is sick, that is, when there is a difference in stiffness compared to the other healthy fold (flaccid), which causes an asymmetry between them and a non-synchronized movement in most of the cases (varying the asymmetry parameter).

Based on this context, the technique developed on this work yields a better understanding of the vibrating process in vocal folds and, despite their stiffness asymmetry, synchronization can be attained, thus producing a harmonic and synchronized movement between the folds and providing a good functioning, which renders possible sound emission.

Mathematical model

Modeling the vibration of vocal folds considering all their peculiarities would be very tricky since the glottis is tridimensional, the glottal surface is asymmetrically deformable, the vocal folds themselves are not identical (due to possible asymmetries, such as difference in length, width and stiffness, for instance), among other reasons. In other words, it is a truly complex matter. The geometry to be considered is not an easy task and would make the model confusing, with several parameters.

Therefore, the movement of each vocal fold (right/left) is assumed to be the movement of wave propagation towards the air flow. Besides, a simpler structure will be considered, such that the system's non-linearities (geometry, viscosity, abrupt collision of vocal folds) do not present any substantial losses during the analysis. All properties of the tissue are assumed to be concentrated at the mid-point of the vocal folds, that is, the model considered is to be grouped. By performing a balance of forces in the direction of the movement (Titze, 1988), the motion equations for both right and left vocal folds, respectively (see Figure 2), can be written as follows Equation 1 and 2.

where:

ξr,l are the mid-point displacements of the glottis tissue, Mr,l, br,l and Kr,l are the mass, linear damping and stifiness by mean unit of surface area, respectively, nr,l are the non-linear damping coefficients, and Fg is the external force acting on the glottis, which will be replaced here by the mean pressure on the glottis (Laje, Gardner, & Mindlin, 2001) indicated as Pg considered a more appropriate term.

In Figure 3, L is the vocal fold length and xr and xl are, respectively, the displacement of the right and left fold; x0 is half the glottis' width when the folds are at rest, and T is the glottis' height.

According to (Titze, 1988), for small displacements of xr and xl and as well as for small time lag τ, where τ is the time lag for the surface wave to cover half of the glottis' height T it follows Equation 3.

where:

kt is the transglottal pressure and Ps is the subglottic pressure. As only stiffness asymmetry is considered , .

Considering , and assuming P_s high enough so that , changing the variable and performing a time scale such that a van der pol system coupled with 2 equations is achieved, given as Equation 4 and 5.

where:

is the asymmetry parameter and is the coupling parameter.

Equation 4 and 5 in turn, can be written as Equations 6 to 9 of states.

Rewriting the system in a more compact form, Equation 6 to 9 become, according Equation 10.

where:

The variable ‘X’ stands for the derivative of X with respect to t.

Note that the only equilibrium point is E_q = (0 0 0 0)^T.

Results and discussion

Below are the results found through various simulations where different situations were considered.

Numerical simulation and scenario analysis

Three scenarios are considered in this section, assuming in all simulations that the right vocal fold is healthy, whereas the left vocal fold is flaccid, which is classified as a disease.

Stability

In order to begin an investigation on the asymptotic behavior of the system modeled by equation (10), the stability theory for non-linear systems is applied (Hale, 2009).

Note that G(X) = o(|x|), that is, G(x) → 0 when |x| → 0. Besides, G(0) = 0, where O = (0000)^T. Therefore, based on the stability theory for disturbed linear systems, it is possible to study the qualitative behavior of the system above in the neighborhood of the equilibrium point Eq simply by addressing the linear system's behavior (Hale, 2009). To determine the stability of the linear system X = AX the characteristic polynomial associated with matrix A must be found. In the present case, it is as follows Equation 11.

By varying the parameters of coupling (α) of asymmetry (Q) and setting the parameter (µ), the stability diagram (Figure 4) for the only equilibrium point - which is the origin Eq = (0 0 0 0)^T - can be found. The variable Q is considered fixed in the interval using a 0.05 step size, and then α varies in the same interval, for the same step size. Varies in the same interval, for the same step size 4 shows, where the circles indicate the equilibrium point's stability, and the asterisks indicate the equilibrium point's instability.

Scenarios considered

To obtain a more realistic picture of the problem, the parameters considered are based on an adult man according to Titze (1988), where M = 0.5 g cm^-2, B = 50 dyne s cm^-3, K = 200000 dyne cm^-3, τ = 1 ms, kt = 1.1, ζ0 = 0.1 cm, Ps = 800 Pa, which yield μ = 0.3 and α = 0.23. The parameter α can also vary in the interval [0 1] (which corresponds to real situations). All simulations are based on the same initial conditions (x_rx_rx_lx_l)^T = (0.10 - 0.10)^T. Three different situations are presented as follows:

Situation 1: As illustrated in Figure 4, as shown in Figure 5 the vocal fold movement's amplitude increases until it becomes constant (equilibrium point's instability, see the instability diagram in Figure 4), and then remains constant after a certain time. It can be seen that the vocal folds are practically in phase.

Situation 2: In this case, as shown in Figure 6, the oscillating movement decreases in amplitude until it stops eventually. It is based on these parameters in the stability diagram that the equilibrium point, the origin, is stable (see Figure 4).

Situation 3: As shown in Figure 7 it can be seen that there exists no relationship between the phases from the left and right vocal folds, that is, the movement is not synchronized. After a long time, the amplitudes of each fold start decreasing. Based on the stability diagram, Figure 7, one notices that there exists stability for the chosen parameters. In this figure, a close-up has been applied to the second graph in order to better understand the analysis.

Lyapunov exponent

Lyapunov exponent's negative values λ characterize a convergent behavior, either periodic or of fixed-point, in which the system's trajectories over time that initially started in states close to each other tend to get closer. Points whose values are zero indicate qualitative behavior changes that are called bifurcation points. On the other hand, if the values are positive, they indicate the occurrence of chaos, thus describing the distancing rate of the system's trajectories that initially started in states close to each other.

The more negative the exponent is, the faster the series converges. Therefore, the Lyapunov exponent is a good chaos indicator and also measures how fast two arbitrarily close orbits move closer or away from each other over time (Vakakis, 2001).

In the present work, the Wolf method was applied (Wolf, Swift, Swinney, & Vastano, 1985). In this case, the Lyapunov exponents are obtained through comparison between the trajectory obtained from the non-linear differential equation and the trajectory obtained through the original equation's linearized solution. In this way, successive approximations are made and, after some time evolution, an orthonormalization process of the linearized solutions is carried out.

The Lypuanov exponent is applied to the same values as the simulated parameters in the previous section (see Table 1).

It can be seen that in Figure 8 the system is chaotic for the first and third situations, that is, there is no synchronization between the systems, which does not cause sound emission or causes disordered voice emission (white noise). For the second situation, the system is non-chaotic.

Optimal linear control

Given the previous results, that is, for situations 1 and 3, where chaos occurs, the vibration is not synchronized, which entails very noisy sound emission or, in the worst case scenario, total absence of sound. The implementation of a method to fix this behavior, or disorder, is thus necessary, and to do so a new optimal control design shall be addressed in the next section.

For the second situation, in which chaos does not occur, the method is also applied to explore the possibility of making the vocal folds vibrate in a new desired frequency, thus enhancing the quality of sound emitted: high frequencies present high pitch sounds, and low frequencies present low pitch sounds. In this case, vibration control is also important, which justifies the control method applied.

A controller design usually consists of determining a control law that makes the system comply with certain performance requirements. In the case of optimal control systems, the control law can be found through minimization of a cost functional (or performance index) J= F(x, u, t).

The control law must be such that it ensures the system's stability, as well as minimizes the system's efforts on keeping this stability. In compliance with these two requirements, the Linear Quadratic Regulator (LQR) may be used.

For the vocal folds, it is interesting to note that in all situations, even for the case in which chaos occurs, control can be attained such that the system exhibits synchronized solutions, leading to the good functioning of vocal folds, thus causing a synchronous vibration, and ultimately producing sound.

Linear control formulation

Consider the following controlled non-linear system (Rafikov et al., 2008), according Equation 12.

where:

is the vector of control. Let x be the vector defined as the desired trajectory of U composed of two parts, according Equation 13.

where:

u is the feedback control given by Equation 14.

and u is the linear response that can be rewritten as Equation 15.

where:

is a constant matrix. Defining Equation 16.

Considering the desired trajectory and both Equation 14 and 15, it is obtained, through derivation Equation 17.

Considering now Equation 18.

the system (Equation 17) may be written as Equation 19.

According to Bryson (1975) there exists a positive definite symmetric matrix Q(t) and a positive definite R(t) such that the matrix, according Equation 20.

is positive definite. So, the linear control with feedback, according Equation 21.

is optimal and moves the system (Equation 19) from the initial state to the final state, according Equatin 22.

Minimizing the functional, according Equation 23.

where:

the real symmetric matrix P(t) is calculated by solving the Ricatti Equation 24.

in agreement with the final condition, according Equation 25.

The case where t_f →∞ and the matrices A, B, Q and R are matrices with constant components, matrix P is obtained from the Ricatti Equation's 26 solution.

LQR control application

In this section, the optimal linear control is applied with the aim of obtaining a stable periodic solution for the system, where the synchronized movement represents the good functioning of the vocal folds, considering that this movement starts and, after some short time, it synchronizes itself in order to make the most of the force related to the pressure due to air expelled from the lungs.

From Equation 10 it follows that Equation 27.

The desired trajectory is Equation 28 and 29.

where:

B is a fixed and arbitrarily chosen positive definite matrix, according Equation 30.

where:

Q is a positive semi-definite matrix, according Equation 31.

Matrix Q and the scalar R weigh the effects of minimization over the states and the input, respectively.

Situation 1: By applying the control method, a controllability value of 0.1224, is attained, which means that the system is controllable. For the gain factor matrix, K = [2.50701.11481.75610.5111]. From Figure 9, one notices that the vocal folds quickly reach synchronization. Therefore, the movement that previously showed to be almost synchronized, now becomes synchronized in phase.

Situation 2: In this case controllability has a value of 0.5112, the gain is K = [0.41721.05721.03251.0814] and the system, which was previously stable and non-synchronized, now becomes synchronized in phase, as shown in Figure 10.

Situation 3: For the last simulation, a controllability value of 0.1251, is attained, and the gain is K = [0.75221.60861.15901.0044]. For this system, which previously was chaotic, now presents a synchronized behavior, as shown in Figure 11.

Therefore, Figure 9, 10 and 11 show that the optimal linear control method is effective for the situations analyzed, which ensures the good functioning of the vocal folds, rendering possible sound emission again.

Conclusion

A new application of an optimal linear control design has been proposed in this work concerning the synchronization of asymmetric vocal fold's movement.

The simulations showed that the system under optimal linear control has proven to be robust and effective on controlling the movement of vocal folds which present different stiffness. Thus, a periodic and in phase movement was produced with the aim of achieving the good functioning of the vocal folds, rendering possible sound emission even in diseased folds (one stiff, or healthy, and the other one flaccid, or diseased), such as the case of flaccid dysphonia.

Acknowledgements

The authors are grateful to the Laboratory of Complex Systems (Sisplexos), where the project has been developed, as well as to the Post-Graduate Program of Mechanical Engineering (PPGEM), and the Federal University of Mato Grosso (UFMT).

References

Berry, D. A., Herzel, H., Titze, I. R., & Krischer, K. (1994). Interpretation of biomechanical simulations of normal and chaotic vocal fold oscillations with empirical eigenfunctions. Journal of the Acoustical Society of America, 95(6), 3595-3604. doi: 10.1121/1.409875

Berry, D. A., Herzel, H., Titze, I. R., & Story, B. H. (1996). Bifurcations in excised larynx experiments. Journal of Voice, 10(2), 129-138. doi: 10.1016/S0892-1997(96)80039-7

Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. L., & Zhou, C. S. (2002). The synchronization of chaotic systems. Physics Reports, 366(1-2), 1-101. doi: 10.1016/S0370-1573(02)00137-0

Bryson, A. E., & Ho, Y.-C. (1975). Applied optimal control. Washington, DC: Hemisphere Publishing Corporation.

Chavarette, F. R. (2011). On an optimal linear control of a chaotic non-ideal duffing system. Applied Mechanics and Materials, 138-139, 50-55. doi: 10.4028/www.scientific.net/AMM.138-139.50

Chavarette, F. R., Balthazar, J. M., Rafikov, M., & Hermini, H. A. (2009). On non-linear dynamics and an optimal control synthesis of the action potential of membranes (ideal and non-ideal cases) of the Hodgkin-Huxley (HH) mathematical model. Chaos, Solitons & Fractals, 39(4), 1651-1666. doi: 10.1016/chaos

Ferreira, D. C., Chavarette, F. R., & Peruzzi, N. J. (2017). Linear matrix inequalities control driven for non-ideal power source energy harvesting. Journal of Theoretical and Applied Mechanics, 53(3), 605-616. doi: 10.15632/jtam-pl.53.3.605

Giovanni, A., Ouakinine, M., & Triglia, J. M. (1999). Determination of largest Lyapunov exponents of vocal signal: application to unilateral laryngeal paralysis. Journal of Voice, 13(3), 341-354. doi: 10.1016/S0892-1997(99)80040-X

Giovanni, A., Ouakinine, M., Guelfucci, B., Yu, T., Zanaret, M., & Triglia, J. M. (1999). Nonlinear behavior of vocal fold vibration: the role of coupling between the vocal folds. Journal of Voice, 13(4), 465-476. doi: 10.1016/s0892-1997(99)80002-2

Gupta, S., Dalal, U., & Mishra, V. N. (2014). Novel analytical approach of non conventional mapping scheme with discrete hartley transform in OFDM system. American Journal of Operations Research, 4(5), 281-292. doi: 10.4236/ajor.2014.45027

Gupta, S., Dalal, U., & Mishra, V. N. (2015). Performance on ICI self-cancellation in FFT-OFDM and DCT-OFDM system. Journal of Function Spaces, 2015, 1-7. doi: 10.1155/2015/854753

Hale, J. K. (2009). Ordinary differential equations. New York, NY: Dover Publications.

Heagy, J. F., Carroll, T. L., & Pecora, L. M. (1994). Synchronous chaos in coupled oscillator systems. Physical Review E, 50(3), 1874-1885. doi: 10.1103/PhysRevE.50.1874

Hirano, M. (1974). Morphological structure of the vocal cord as a vibrator and its variations. Folia Phoniatrica, 26(2), 89-94. doi: 10.1159/000263771

Hirano, M., & McCormick, K. R. (1986). Clinical examination of the voice. The Journal of the Acoustical Society of America, 80(4), 1273. doi: 10.1121/1.393788

Laje, R., Gardner, T., & Mindlin, G. B. (2001). Continuous model for vocal fold oscillations to study the effect of feedback. Physical Review E, 64(5 Pt 2), 056201. doi: 10.1103/PhysRevE.64.056201

Mishra, L. N., Sen, M., & Mohapatra, R. N. (2017). On existence Theorems for some generalized nonlinear functional-integral equations with applications. Filomat, 31(7), 2081-2091. doi: 10.2298/FIL1707081N

Pecora, L. M., & Carrol, T. L. (1990). Synchronization in chaotic system. Physical Review Letters, 64(8), 821-824. doi: 10.1063/1.4917383

Pikovsky, A., Rosenblum, M., & Kurths, J. (2003). Synchronization: a universal concept in nonlinear sciences. New York, NY: Cambridge Nonlinear Science Series.

Rafikov, M., Balthazar, J. M., & Tusset, A. M. (2008). An optimal linear control design for nonlinear systems. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 30(4), 279-284. doi: 10.1590/S1678-58782008000400002

Rosenblum, M. G., Pikovsky, A. S., & Kurths, J. (1996). Phase synchronization of chaotic oscillators. Physical Review Letters, 76(11), 1804-1807. doi: 10.1103/PhysRevLett.76.1804

Rosenblum, M. G., Pikovsky, A. S., Kurths, J., Schafer, C., & Tass, P. A. (2001). Chapter 9 phase synchronization: From theory to data analysis. Handbook of Biological Physics, 4, 279-321. doi: 10.1016/S1383-8121(01)80012-9

Titze, I. R. (1988). The physics of small-amplitude oscillation of the vocal folds. Journal of the Acoustical Society of America, 83(4), 1536-1552. doi: 10.1121/1.395910

Tusset, A. M., Baltthazar, J. M., Chavarette, F. R., & Felix, J. L. P. (2012). On energy transfer phenomena, in a nonlinear ideal and nonideal essential vibrating systems, coupled to a (MR) Magneto-rheological damper. Nonlinear Dynamics, 69(4), 1859-1880. doi: 10.1007/s11071-012-0391-5

Vakakis, A. F. (2001). Inducing passive nonlinear energy sinks in vibrating systems. Journal of Vibration and Acoustics, 123(3), 324-332. doi: 10.1115/1.1368883

Willms, A. R., Kitanov, P. M., & Langford, W. F. (2017). Huygens’ clocks revisited. Royal Society Open Science, 4(9), 170777. doi: 10.1098/rsos.170777

Wolf, A., Swift, J. B., Swinney, H. L., & Vastano, J. A. (1985). Determining lyapunov exponents from a time series. Physica D: Nonlinear Phenomena, 16(3), 285-317. doi: 10.1016/0167-2789(85)90011-9

Zhang, H., Liu, D., & Wang, Z. (2009). Controlling Chaos: suppression, synchronization and chaotification. London, UK: Springer-Verlag.

Notas de autor

fabio.chavarette@unesp.br