Experimental and Computational Modeling of the Effects of Voice Therapy Using Tubes

1 EXPERIMENTAL AND COMPUTATIONAL MODELING OF THE EFFECTS OF VOICE THERAPY USING TUBES J. Horáček1, V. Radolf1, A.-M. Laukkanen2 1Institute of Thermomechanics, Czech Academy of Sciences, Dolejškova 1402/5, 182 00 Prague, Czech Republic, [email protected], [email protected] 2Speech and Voice Research Laboratory, Faculty of Education, Åkerlundink. 5, University of Tampere, 33100 Tampere, Finland, [email protected] Corresponding Author: Jaromír Horáček, Institute of Thermomechanics, Czech Academy of Sciences, Dolejškova 1402/5, 182 00 Prague, Czech Republic, [email protected], Conflict of Interest: There are no relevant conflicts of interest. Funding Statement: This study was supported by a grant from the Czech Science Foundation: No. 16-01246S “Computational and experimental modelling of self-induced vibrations of vocal folds and influence of their impairments on human voice.” This is the accepted manuscript of the article, which has been published in Journal of Speech, Language, and Hearing Research. 2019, 62(7), 2227-2244. https://doi.org/10.1044/2019_JSLHR-S-17-0490 2 Abstract: Phonations into a tube with the distal end either in the air or submerged in water are used for voice therapy. This study explores the effective mechanisms of these therapy methods. The study applied a physical model complemented by calculations from a computational model, and the results were compared to those that have been reported for humans. The effects of tube phonation on vocal tract resonances and oral pressure variation were studied. The relationships of transglottic pressure Ptrans(t) variation in time vs. glottal area variation GA(t), were constructed. The physical model revealed that, for the phonation on [u:] vowel through a glass resonance tube ending in the air, the first formant frequency F1 decreased by 67%, from 315 Hz to 105 Hz, thus slightly above the fundamental frequency F0 that was set to 90–94 Hz. For phonation through the tube into water, F1 decreased by 91-92%, reaching 26–28 Hz, and the water bubbling frequency Fb≅19– 24 Hz was just below F1. The relationships of Ptrans(t) vs. GA(t) clearly differentiate vowel phonation from both therapy methods, and show a physical background for voice therapy with tubes. It is shown that comparable results have been measured in humans during tube therapy. For the tube in air, F1 descends closer to F0, while for the tube in water, the frequency Fb occurs close to the acoustic-mechanical resonance of the human vocal tract. In both therapy methods, part of the airflow energy required for phonation is substituted by the acoustic energy utilizing the first acoustic resonance. Thus, less flow energy is needed for vocal fold vibration, which results in improved vocal efficiency. The effect can be stronger in water resistance therapy if the frequency Fb approaches the acoustic-mechanical resonance of the vocal tract, while simultaneously F0 is voluntarily changed close to F1. Keywords: Biomechanics of voice; maximum glottal area declination rate (MADR), vocal tract acoustics; formant frequencies; phonation into tubes; water resistance voice therapy; water bubbling frequency, vocal efficiency. 3 INTRODUCTION Phonation into different kinds of tubes is widely used for voice training and therapy (Sovijärvi, 1969; Laukkanen, 1992; Titze, Finnegan, Laukkanen, & Jaiswal, 2002; Granqvist et al., 2015; Amarante et al., 2016; Guzman et al., 2017). According to practical experience, phonation into a tube may feel easier than ordinary vowel phonation, and speech immediately after tube phonation often sounds louder and less strained. Voice therapy has exploited phonation through a tube into water – i.e., water resistance therapy – for treating both hypofunctional and hyperfunctional voice disorders. Two variants of tubes have been increasingly widely used in water resistance therapy: a glass tube called a “resonance tube” (Sovijärvi, 1969; Laukkanen, 1992; Simberg & Laine, 2007) with a length of 24–28 cm and an inner diameter of ca. 8–9 mm, and a Lax Vox tube made of flexible silicon, 35 cm in length and 1 cm in inner diameter (Sihvo, 2006; Sihvo & Denizoglu, 2007; Sihvo, 2017). During the last two decades, many studies have been carried out to explain the basis of tube training and therapy. It is known that artificial prolongation and narrowing of the vocal tract (VT) increases the impedance. In particular, phonation through narrow tubes and straws and phonation through a tube into water increase supraglottic resistance (air pressure/flow) (Titze et al., 2002; Horáček, Radolf, Bula, & Laukkanen, 2014; Amarante et al., 2016). Increased supraglottic resistance offers training for breath support (Carroll & Sataloff, 1991) and decreases transglottic pressure – i.e., the difference between subglottic and supraglottic pressure – which provides the driving force of vocal fold vibration (Titze & Laukkanen, 2007). The phonation threshold and amplitude of the vocal folds’ (VFs) vibration decrease; thus, the biomechanical loading related to phonation decreases (Titze, 2006b). Simultaneously, especially if the first acoustic resonance of the VT plus the tube is lowered effectively close to the fundamental frequency (F0) of the VFs’ vibration, the positive inertance of the VT may improve the voice quality by assisting the VFs’ vibration and/or strengthening the amplitude of the harmonics (Story, Laukkanen, & Titze, 2000; Fant & Lin, 1987; Rothenberg, 1981). 4 Earlier modelling studies have shown that positive inertance of the vocal tract promotes vocal fold vibration as the oscillation of the air pressure above the vocal folds is in phase with the velocity of the vocal folds, and the phonation threshold is also reduced (Titze, 2001). There are also findings for a female singer showing that the vocal fold vibration is positively affected when F0 is close to F1 (Rothenberg, 1988). Additionally, the increased supraglottic resistance intensifies sensations of resonatory vibrations in the VT (Titze, 2006b). Thus, the exercise with tubes and other semiocclusions of the VT offers the trainee sensations of more economic and efficient voice production. These sensations may then be a guide to establish similar acoustic-mechanical conditions in the VT also after the exercise, since the trainee may learn to apply epilaryngeal narrowing as a source of impedance-matching between the glottis and the VT (Titze, 2006b). Moreover, water resistance therapy offers the element of water bubbling, which is reflected in the oscillation of supraglottic pressure and results in variation of the amplitude of the VFs’ vibration (Granqvist et al., 2015). This bubbling effect has been reported to feel like a massage of the VT and the larynx, and it has been suggested that it could offer similar beneficial effects as massage – e.g., improved blood circulation in the tissue (Mori et al., 2004). Although water bubbling has been regarded as relaxing for the muscles and potentially healing for the tissues, there are also some results showing that the relative glottal closing speed and the closed time of the glottis may increase (Guzman et al., 2017) when the tube immersion depth is high (e.g., 10–18 cm below the water surface) which suggests higher mechanical loading of the VFs. It is also possible for the water bubbling frequency to approach the acoustic–mechanical resonance of the laryngeal tissues, which may increase the effect of bubbling both beneficially and adversely (Horáček, Radolf, & Laukkanen, 2017a). This study aims to shed more light on the similarities and differences between tube training with the distal end of the tube in air and tube training with the distal end of the tube immersed in water. This study focuses on the following research questions: 5 1. Are there any essential differences between the glottal area variation in time GA(t) versus transglottic pressure Ptrans(t) variation measured for phonation on the vowel [u:] and phonation into the resonance tube with the distal end either in air or in water? An important reason for studying this relation is the fact that it gives an estimate of the average work done by airflow during vocal fold vibration thus allowing estimation of changes in vocal efficiency. Such a study is difficult to realize in humans, although investigations have been conducted with high speed filming of the VF using larger tubes (Laukkanen et al., 2007; Guzman et al., 2017) and filming through the nose (Granqvist et al., 2015). Both cases obviously result in unnatural VT conditions compared to ordinary phonation into a resonance tube. 2. Are there any important differences between the evaluated glottal area time derivatives for all three cases of phonation considered? Our focus is especially on the potential differences in the maximum area declination rate (MADR), which is considered to be a measure of the stress loading imposed on the VFs during collisions (Titze & Laukkanen, 2007). 3. What are the effects of the acoustic-mechanical resonance and the first acoustic resonance in phonation through the tube in air and in water? Are there any fundamental differences in the principles of these two therapy methods? This information would be useful to guide the proper choice and use of the tubes. In this study, a physical model will be applied that includes a lung model, artificial VFs, and a plexiglass model of the VT (Horáček, Bula, Košina, & Radolf, 2016a). The results obtained with this hard-walled VT model are further compared with (1) those obtained previously in humans and with (2) calculations obtained from a computer model, where yielding walls were implemented to better resemble the VT of humans.: Yielding walls affect the formant frequencies and also lead to appearance of a low-frequency acoustic-mechanical resonance (Radolf, Laukkanen, Horáček, & Liu, 2014; Radolf et al., 2016; Horáček et al., 2017a). 6 For simplicity, we use F1 to refer both to the first formant frequencies measured in the model and in humans and to the first computed acoustic resonance frequencies. METHODS Physical model A model of the human lungs, which includes the splitting of the airways up to fourth-order branching, was built in the subglottic part of the experimental facility (see Figure 1 and more in Horáček, Radolf, Bula, & Košina, 2017b). The air flowed through the model of the lungs to the trachea, which was modeled by a metal tube. The total length of the trachea was 23 cm, and the inner diameter was 18 mm. The experiments were performed with a 1:1 scaled three-layered vocal fold model with a total length of 20 mm, vertical thickness 10.3 mm and width 8 mm (Horáček et al., 2016b; Horáček et al., 2017b). The middle cross-section of the vocal fold model was based on CT measurements of a female subject (Hampala et al., 2015). The layers were formed by a silicon wedge modeling the vocal fold body inside the VFs, a 2-3 mm thick water layer modeling the lamina propria, and a silicon cover with the thickness slightly below 2 mm. Filling with water lowers the F0 of the VF model to ca. 80 Hz and the phonation threshold flow (PTF; i.e., the lowest flow rate to start sustained phonation) to ca. 0.02 l/s, covering also the lowest limits of human phonation (see, e.g., Baken & Orlikoff, 2000). The fundamental phonation frequency was preset by an initial pre-stressing of the VFs. This was accomplished by a slight variation of the distance between the frames enclosing the VFs, followed by changing the hydrostatic pressure inside the VFs. After filling the “lamina propria” layer in the VFs with water and completing the tuning procedure, the model was fixed in a plexiglass frame and connected to the model of the subglottal spaces on one side and to the model of the VT on the other side. The geometrical configuration of the plexiglass model of the vocal tract corresponded to the Czech vowel [u:]. It consisted of 46 circular cross-sectional areas, perpendicular to the midline along the VT from the VFs up to the lips. These areas were obtained from magnetic resonance images (MRI) 7 recorded during the phonation of a 35-year-old male subject (Vampola, Horáček, & Švec, 2008). The total length of the VT model was 18.7 cm. Set-up and measurements on the physical model Measurements were performed for: a) simulated phonation on the vowel [u:], b) phonation with [u:]- vowel-shaped VT attached to a glass resonance tube (length 27 cm, inner diameter 7.8 mm) ending in air, c) phonation with [u:]-vowel-shaped VT attached to the same tube ending in water, and d) blowing through the setup mentioned in c) but without the vocal folds vibrating (no phonation). The (d) condition was included to compare subglottic pressure, oral (back) pressure, and water bubbling frequency with and without phonation. Such a comparison may reveal more on the interplay between different types of oscillations of the air pressure occurring simultaneously in the VT during phonation into water. In the experiments with phonation, flow rates ranged between the PTF and the maximum rate beyond which the vibration amplitudes would endanger the material coherence of the artificial VFs. The glottal width between the VFs was preset to zero and not changed. In the experiment without phonation, the VFs were abducted by deflation, thus simulating the breathing position of the VFs. In this (d) condition the glottal width was preset to ca.1 mm. The following parameters were measured: 1) mean airflow, Q; 2) phonation threshold flow, PTF; 3) mean subglottic air pressure, Psub; 4) phonation threshold pressure, PTP (i.e. the lowest Psub to start and sustain phonation); 5) peak-to-peak variation in Psub(t); 6) frequency of the lowest subglottic resonance; 7)mean oral air pressure, Poral; 8) peak-to-peak variation in Poral; 9) frequency F1 of the lowest oral resonance (formant); 10) 3 dB bandwidth B1 of the lowest oral resonance, 11) fundamental frequency F0 of phonation (measured from Poral variation); 12) water bubbling frequency, Fb (measured from Poral variation); and 13) variation of glottal area, GA(t). The time derivative of glottal area variation dGA/dt (representing the speed of glottal opening and closing) was calculated by the method of central differences (Mathews, & Fink, 2004). Parameters 1–13 were measured because the VF vibrations and 8 thus glottal area variation are dependent on aerodynamic variables below and above the glottis. The time variation in the pressures, in turn, is related to supraand subglottic resonances. A schema of the measurement set-up is shown in Figure 1. The VFs were excited by airflow coming from a regulated central pressure supply. The mean airflow rate was measured by a float flowmeter (EMKO type DF3-09K5), and by an orifice flow meter with a pressure sensor (LP1000 DRUCK model LPX1812-C1S) using an analogue output. The mean subglottic and oral pressures were measured by integrated pressure semiconductor sensors (NXP Freescale MPXV5010GC6U) mounted in the walls of the trachea and oral cavity models. The mean subglottic pressure at the entrance to the VFs was registered by a digital manometer (Greisinger Electronic, GDH07AN) connected by a short compliant tube to the subglottic space of the glottal cavity. A second digital manometer, a Greisinger Electronic GMH3151, measured the mean pressure of the water inside the VFs, which was adjusted with a waterfilled syringe. The fluctuations of the subglottic and oral pressures were measured by a miniature microphone (Bruel & Kjaer 4138, range 6.5 Hz–140 kHz) and by a special microphone probe (B&K 4182, range 1 Hz– 20 kHz), respectively. The spectra of the pressure signals were calculated by an in-house developed program in Matlab, and the formant frequencies were estimated from the peaks of the averaged spectral data (Horáček et al., 2017a). All the measured signals were simultaneously sampled at a frequency of 16.4 kHz and registered by a B&K PULSE type 3560 C measurement system with type 7537A and 3109 Input/Output Controller Modules controlled by a personal computer (PC I) equipped with SW PULSE LabShop Version 10. A high-speed CCD NanoSense Mk. III camera (maximum resolution 1280 x 1024 pixels) with a Nikon AF micro Nikkor 60mm zoom lens was included in the measurement set-up for the analysis of the VFs’ vibration. The camera was positioned at a 90-degree bend of the trachea model where a glass window was installed; this enabled the viewing of the VFs’ vibration from the subglottal side. A personal computer (PC II) was used for recording the VFs’ vibration. During image recording at a frequency of 9 3000 frames/second, three intensive LED lights (2 x 13 W+1 x 35 W) were focused on the vibrating VFs. The images recorded by the camera were synchronized with the time records of the pressure signals. Figure 1. Somewhere here. Computational analysis of the measured relations between Ptrans(t) and GA(t) The relationship of Ptrans(t) vs. GA(t) measured during regular periodic self-sustained VF oscillations can be visualized as a closed, clockwise-oriented cyclic curve (loop). Similar graphs showing the relationship between the subglottic pressure and the glottal opening time variation were introduced by Mišun, Švancara, & Vašek (2011). They measured these relationships on artificial VFs for various selfoscillating regimes without any VT. In the present study, we use similar graphs to compare phonation on a vowel with phonation through the resonance tube with the distal end in air and in water. By computing the average area inside the loop, we get an estimate of a part of the average mechanical work done by the airflow per one cycle of the VFs’ vibration (i.e. part of the flow energy needed to excite the self-sustained VF oscillations). This can be used as a measure of vocal efficiency. We note that the acoustic energy is not included in our considerations, because it was not possible to measure the airflow volume velocity waveform Q(t) synchronously with the pressure signals. The normalized area A inside the graph (loop) was computed by an in-house developed program in Matlab that uses the following mathematical formula for the counter-integral along the boundary A of the area A :         00 00 11 d dt = dt nT nT eq eq eq eq trans A dGO t dGA t A F GO F t T P t n dt n dt        , (1) where   () eq trans F t P t T l is a roughly estimated force loading the VF surface approximated by the area Tl ; 10 mmT and l=20 mm are the VFs’ thickness and length, respectively;     / eq GO t GA t l is an equivalent glottal width;   GA t and Ptrans(t) are the measured glottal area and transglottic pressure, 16 Figure 8 presents an example of the relationships Ptrans(t) vs. GA(t) and Ptrans(t) vs. dGA(t)/dt for phonation through the resonance tube with the distal end 10 cm in water for the same mean flow rate Q=0.08 l/s as in the previous case for phonation through the tube into air. The measured loops are very similar as in the previous case. This suggests that the effect of the near coincidence of the formant frequency F1≅26 Hz and the bubbling frequency Fb≅23 Hz is similarly strong as in the case of the near coincidence of F0 and F1 for phonation through the tube into air, where the difference F1–F0 was 11 Hz. The mean flow rate is again much lower than the flow rate Q=0.25 l/s for the phonation on the vowel [u:]; furthermore, as in the previous case, the mean transglottic pressure Ptrans≅0.78 kPa is lower than the mean Ptrans≅1.95 kPa for phonation on the vowel [u:]. Figure 8. Somewhere here. Figure 9 demonstrates that phonation into water is more complicated due to the bubbling process, which creates waves on the water’s surface. When the variation of the glottal area, the derivative of the glottal area, transglottic pressure, and the constructed loops are observed for a longer time, the irregularities in amplitudes of all quantities caused by the water bubbling are evident. The areas inside the loops vary irregularly from one oscillation cycle to the next. Therefore, when studying phonation through a tube into water, the energy transfer between the flow and the VFs needs to be investigated for a longer time, taking into account more cycles in order to obtain an average value of the area ∆A. This procedure resulted in ∆A≅-0.041 mJ which is quantitatively comparable value with that obtained for phonation through the tube into air, where the behavior of the system is periodic with only small disturbances. The principle of supporting the self-sustained VF oscillations by the effect of the near coincidence of the formant frequency F1 with the fundamental frequency F0 or bubbling frequency Fb is valid in both cases of phonation into the tube. Figure 9. Somewhere here. This principle can be documented by the flow power ΔFP=ΔA F0 [W] that was computed from all the loops measured within the airflow rate interval Q=0.02–0.25 l/s considered in this study. Figure 10 17 shows the measured power ΔFP as a function of the input steady aerodynamic flow power defined as the mean subglottic pressure Psub multiplied by mean flow rate Q (Schutte, 1980). The power ΔFP measured for the vowel phonation approximates the part of the flow power consumed by the selfoscillating vocal folds. If ΔFP for phonation with the tube is negative, it implicates that the part of the flow power consumed by the self-oscillating vocal folds is smaller. The flow power ΔFP, evaluated from the area inside the loops, is the highest for phonation on the vowel [u:], for which ΔFP≅0.8–7 mW, while in the case of phonation through the tubes ΔFP is lower and mostly negative with a minimum ΔFP≅-8 mW for the tube phonation into air, where the difference F1-F0=14-15 Hz between the formant frequency F1=105 Hz and the frequency F0=91 Hz was minimal (see Figures 2b and 5, and Table 1). The flow power ΔFP for the tube phonation into water is negative having the minimum ΔFP≅-6 mW at the highest input flow power used in the experiments for tube phonation into water where the difference between the formant frequency F1=28 Hz and the bubbling frequency Fb was ca. 4-8 Hz (see Figures 2c and 5, and Table 1). The experiments on vowel phonation and on tube phonations were performed in only slightly overlapping intervals of steady input flow power, because of substantial lowering of the PTF caused by the tube. Thus in only one case shown in Figure 10, it was approximately possible to compare directly the power ΔFP for phonation on vowel with that for phonation into water. The steady input flow power for the vowel was 159 mW and for the tube in water it was 150 mW, which are close values. The power ΔFP decreased from 0.8 mW for vowel to -4.2 mW for tube in water. This shows that 5 mW of flow power needed for vocal folds’ self-sustained vibration was saved. Figure 10. Somewhere here. DISCUSSION An important finding in the present study is that for phonation through the resonance tube in air, not only F1 decreased, but the fundamental phonation frequency F0 also decreased slightly below F1 (see Table 1). Similarly, for phonation through the tube into water, F1 decreased further, slightly above the 18 water bubbling frequency Fb. When F1 is close to F0 or Fb, it assists the VFs’ vibration, which enables easier phonation utilizing the resonance effect of the acoustic system. In practice, it is possible to tune the system in order to gain maximum support for phonation when phonating through the tube either into air or into water. In the former case it can be accomplished by changing F0, while in the latter case by changing Fb. The relationships Ptrans(t) vs. GA(t) constructed for vowel [u:] resulted in clockwise-oriented cyclic graphs (loops), similar to those measured by Mišun et al. (2011) for artificial VFs without a VT, while the loops obtained for both therapy methods are more complicated. They show that less airflow energy is needed for phonation because the first acoustic resonance at the frequency F1 is closely above the excitation frequency F0 or Fb. This effect reduces the airflow energy needed for the VFs’ vibration and the excitation of acoustic waves. In principle the same situation happens, if the fundamental frequency F0 of ordinary vowel phonation gets close to the first formant frequency F1 of the vocal tract and increases the inertance of the supraglottic spaces (Titze, 2001). Comparison of results from physical modeling, computer modeling, and measurements in humans The experimental results presented in this paper for phonation through tubes into air and into water can be compared to measurements on humans published by Horáček et al. (2017a) for a male subject and by Radolf et al. (2014) for a female subject. The main difference between the results is caused by the existence of the low frequency acoustic-mechanical resonance in the human VT, whereas in the physical model with hard VT walls, such a resonance does not exist. The pressure and frequency data measured here on the physical model are compared with measurements on humans in Table 1, where the data from the model are carefully chosen according to the same magnitudes of the mean subglottic pressure Psub as was measured in humans. From the measurement using the model, the values of the mean oral pressure Poral, peak-to-peak oral pressure 19 Poral_p-t-p, the water bubbling frequency Fb, and the fundamental frequencies were interpolated from the data given by the graphs shown in Figures 3 and 5, considering the same subglottic pressure Psub as measured in humans. The acoustic resonance frequencies computed for the mathematical models with yielding and hard walls are also included in Table 1 for comparison with the measured data. Good agreement was found between the measurements on human subjects and the experiments with the physical model for the mean oral pressure Poral. However, the peak-to-peak fluctuations of the oral pressure Poral_p-t-p measured in the model are far higher than in humans. The reason for this is that the yielding walls of the human VT may cause a higher damping of the pressure fluctuations, and that the difference between the first resonance frequency and the pressure fluctuation frequency, was smaller in the model. Really, for phonation through the tube with the distal end in air, the differences F1-F0≅90 Hz for the male and F1-F0≅34 Hz for the female are evidently higher than the differences F1-F0≅14–15 Hz for the measurement with the VT model (see Table 1). Therefore, the acoustic-structural coupling was stronger in the VT model, because the excitation frequency F0 was nearer to the formant frequency F1, and consequently the amplitudes of oral pressure fluctuations were higher in the model than in humans. The acoustic-mechanical resonance at the computed frequency Fa-m≅27 Hz is not so important in vocal exercises using the resonance tube with the distal end in air, because the frequency Fa-m is considerably below the lower limit of F0 for the ordinary human voice (at least when vocal fry is not considered). However, at least in theory, the mechanical resonance associated with the yielding walls may raise the computed value of F1 for the hard VT from 97 Hz to 181 Hz, which comes close to F1=190 Hz measured in humans (see Table 1). This in turn would be beneficial, especially for the female speaking voice, since F0 could be easily set close to F1. In phonation into water the Fa-m becomes important. For the model, the VF vibrations are supported by water bubbling, because the difference F1-Fb ≅4–8 Hz between the first formant frequency and the bubbling frequency is very small (see Table 1), while the fundamental phonation frequency F0≅79–82 20 Hz is far from the acoustic resonance at F1≅28 Hz. We remark that the vibration amplitude of the VF in the vertical direction should be essential at the bubbling frequency. It should be proved in a next study, e.g. using laser vibrometry. In humans, the bubbling frequency Fb≅11.5–14 Hz was considerably lower than in the model, where Fb≅20–24 Hz (see Table 1). The computed frequency of the acoustic-mechanical resonance Fa-m≅9 Hz is very close to the bubbling frequency Fb measured in humans. This means that in humans, the acoustic-mechanical resonance can be excited by the water bubbling because the difference of only 2–5 Hz can be theoretically estimated between Fa-m and Fb (see Table 1). Moreover, in the case of a female subject, also the first acoustic resonance at the frequency F1≅150 Hz was excited by the very close fundamental frequency F0≅149 Hz, which is also close to the acoustic resonance frequency F1≅143 Hz computed for the VT with the yielding walls (see Table 1). This means that for the female subject, the effect of the water voice therapy can be doubled, because Fb can be close to Fa-m, and at the same time F0 is close to F1. It is worth noting that according to the measurements on a human subject the amplitudes of the bubbling induced vibrations of the larynx were higher than those caused by the fundamental frequency of phonation (Laukkanen et al., 2018). From a comparison of the results obtained from the model with hard VT walls and the measurements made in humans, it can be concluded that the effects of tube phonation in humans and in the model are the same, but the situation is more complicated in humans due to the existence of the low frequency acoustic-mechanical resonance. The results show that optimized conditions for a maximal utilization of the resonances for supporting phonation into the tubes could be found. It would be necessary to tune the system. Tuning can be accomplished by changing several parameters:  the inner diameter and length of the tube following the theoretical results (Horáček et al., 2017a), since the tube length and diameter affect the first resonance;  the fundamental frequency of phonation, within physiological limits, for males and females;  the airflow rate which changes the bubbling frequency, as shown in Figure 5. 21 The most problematic seems to be to change the bubbling frequency Fb to a wider frequency range covering the acoustic-mechanical resonance of the VT (ca. 9 Hz), because an exact theoretical modeling of water bubbling is excluded and only special experiments for chosen tubes can be carried out, as recently presented in Wistbacka et al. (2017), where for a hard-walled VT, the Fb increased from ca. 6 Hz at Q≅0.001 l/s to ca. 21 Hz at Q≅0.04 l/s. It is known that the bubbling frequency increases with the flow rate very quickly from zero for low flow rates to a limit value over which is not possible to cross (see Figure 5 and Davidson & Amick, 1956). General discussion Our modeling results are in line with the computations of Story et al. (2000) for humans when it comes to the lowering of F1 with an artificial prolongation of the VT by a tube in air. Our F1 values were lower when a hard-walled VT was considered. However, when the values were calculated considering a yielding-walled VT, as Story et al. (2000) did, our results (F1=181 Hz) for tube in air are roughly similar to the results by Story et al. for a 30 cm tube with a ca. 1 cm inner diameter (F1 about 210 Hz). Their results for the bilabial plosive [b] (F1 about 170 Hz,) in turn, resemble our results for the tube in water (F1=143 Hz; see Table 1). The differences can be caused by the different VT and slightly different tube length and inner diameter. As far as we know, the present study and the study by Horáček et al. (2017a) are the first to calculate F1 for a water resistance exercise showing the effect of the low frequency acoustic-mechanical resonance in the VT. The Psub and airflow rates considered in the physical modeling were within the range reported for humans (Baken & Orlikoff, 2000), and they remained in the same region as the one that has been found in humans during tube exercises (Tyrmi et al., 2017; Radolf et al., 2014). The bubbling frequency Fb=19–24 Hz in phonation through a tube into water was close to the 21 Hz reported by Wistbacka et al. (2017) for a hard-walled VT. Granqvist et al. (2015) measured Fb in humans in the range 10.5– 12.7 Hz and Tyrmi et al. (2017) in the range 15.8-16.5. The results agree with our finding that the yielding walls in the human VT lower the bubbling frequency (see Table 1). 22 The analysis of the measured relationships Ptrans(t) vs. GA(t) showed the principle of supporting the self-sustained VF oscillations by a near coincidence of the first resonance frequency with the fundamental or bubbling frequency in both cases of phonation into the tube. The measured relationships showed that when F1 was closer to F0 in the tube phonation, less airflow energy was needed for the VFs’ vibration. Thus, the results are in line with the computer modeling results presented by Titze (2006) and Titze & Laukkanen (2007), suggesting improved vocal efficiency with the tubes. The results of our study are also in line with the computer modeling results of Fant & Lin (1987) and Titze (1988) as well as the results of Rothenberg (1986) for a female singer, which showed that the condition where F0 and F1 are close improves the VFs’ vibration. A further study is warranted to quantify, how much the flow energy reduction and thus improvement of VF vibration depends on the difference between the fundamental phonation frequency and the first resonance frequency for phonation through the tube into air, and on the difference between the bubbling frequency and the first acoustic (or the acoustic-mechanical) resonance frequency of the vocal tract for phonation through the tube into water. Hypothetically, it is possible to obtain an efficiency maximum in tube therapy by varying (a) the fundamental phonation frequency, (b) the lowest acoustic resonance frequencies using different tube lengths and inner diameters, and (c) the bubbling frequency (through changing the airflow rate). A question of interest is whether the enhancement of VF vibration may increase the mechanical loading imposed on the VFs. According to the results of the present study, for phonation through the tube into air or into water, PTP, PFT and Ptrans decreased, while MADR did not increase, which suggests that the mechanical loading of the VFs did not increase. These results are in line with the previous modelling study that showed that the impact stress was even smaller during water resistance therapy compared to vowel phonation (Horacek et al. 2018). 23 On the other hand, Guzman et al. (2017) have found in a high-speed study with humans that in some cases of phonation through a tube into water (especially when the immersion depth is deep, 10-18 cm), there is an increased amplitude-to-VF length ratio (substitute of maximum glottal amplitude), closing quotient (substitute of MADR), and glottal spectral flatness; all of these changes suggesting increased impact stress in phonation. Similarly, the high-speed results reported by Laukkanen et al. (2007) showed a trend that in phonation with a longer tube compared to a shorter one, the Psub was higher and glottal open time shorter (i.e. closed time longer). These observations seem to imply increased phonatory effort as a compensation for increased supraglottic resistance. In the present study, the method used was modeling. It is possible that in some cases humans increase laryngeal adduction and Psub, so that the impact stress may rise somewhat; on the other hand, the human high-speed imaging data did not include data for Psub and airflow. The results of the present study suggest that maximum glottal area and the area derivative need a reference to Psub and airflow for drawing conclusions regarding impact stress. Furthermore, as our results show, maximum area declination rate is not a correct parameter for estimating the impact stress; instead, the area declination rate just prior to vocal fold collision should be studied in future. In any case, during tube voice therapy, the subject him-/herself should avoid phonation that may potentially be harmful due to exaggerated enhancement of VF vibration near the first resonance or in water resistance therapy when there is a double effect near the first two resonances (F1 and Fa-m). Such a situation is likely to be felt as unpleasantly strong vibrations in the vocal tract and larynx. This can be expected especially for loud phonation with high subglottic pressure. CONCLUSIONS The main findings from the modeling of the voice therapy based on phonation through a tube into air and into water can be summarized as follows. 1) Glottal area variation measured simultaneously with the variation of the transglottic pressure showed that the airflow energy needed for VF vibration and the excitation of acoustic waves can be 24 noticeably reduced by phonation through a resonance tube into air and into water compared to phonation on the vowel [u:]. In both voice therapy methods, part of the airflow energy required for phonation is substituted utilizing the first acoustic or acoustic mechanical resonance. 2) The basic principle in vocal exercises with a resonance tube with the distal end in air and in water (water resistance therapy) is the same. For phonation into air, the fundamental phonation frequency excites the acoustic resonance at the first formant frequency, and for phonation into water, the bubbling frequency excites a low frequency acoustic mechanical resonance in the human VT with yielding walls. 3) The maximum glottal area time derivative during the glottal closure (so-called maximum area declination rate; MADR) is considered a measure of the maximum impact stress between the colliding VFs. However, it was found in the present study that the glottal closing speed was noticeably reduced just before the VFs’ collision. Therefore, MADR does not seem to be an adequate parameter for impact stress estimation. The glottal closing speed should be investigated in detail in a future study. Acknowledgments This study was supported by a grant from the Czech Science Foundation: No. 16-01246S “Computational and experimental modelling of self-induced vibrations of vocal folds and influence of their impairments on human voice.” 25 References Amarante Andrade, P., Wistbacka, G., Larsson, H., Södersten, M., Hammarberg, B., Švec J. G., & Granqvist S. (2016). The Flow and Pressure Relationships in Different Tubes Commonly Used for Semi-occluded Vocal Tract Exercises. Journal of Voice, 30(1), 36–41. Baken, R. J., & Orlikoff, R. F. (2000). Clinical measurement of speech and voice (2nd ed.). San Diego, CA: Singular/Thomson Learning. Carroll, L. M., & Sataloff, R. T. (1991). The singing voice. In R. T. Sataloff (Ed.): Professional voice: The science and art of clinical care (pp. 381-383). New York, NY: Raven Press Ltd. Davidson, L., & Amick, E. H. JR. (1956). Formation of gas bubbles at horizontal orifices. A.I.CH.E. Journal, 2(3), 337-342. Fant, G., & Lin, Q. (1987) Glottal source –vocal tract acoustic interaction. Quarterly Progress and Status Report STL-QPSR, 1, 13-27. Granqvist, S., Simberg, S., Hertegård, S., Holmqvist, S., Larsson, H., Lindestad, P. A., Södersten, M., & Hammarberg, B. (2015). Resonance tube phonation in water: High-speed imaging electrographic and oral pressure observations of vocal fold vibrations – a pilot study. Logopedics Phoniatrics Vocology, 40, 113-121. Guzman, M., Laukkanen, A. M., Traser, L., Geneid, A., Richter, B., Muñoz, D.,& Echternach, M. (2017). Effect of water resistance therapy on vocal fold vibration: A high speed registration study. Logopedics Phoniatrics Vocology, 42(3), 99-107. Hampala, V., Laukkanen, A. M., Guzman, M., Horáček, J., & Švec, J. G. (2015). Vocal fold adjustment caused by phonation into a tube: A double-case study using computed tomography. Journal of Voice, 29(6), 733-742. 32 Figure 3. Measured values of: a) mean subglottic pressure (Psub), b) peak-to-peak magnitudes of subglottic pressure (Psub_p-t-p), c) mean oral pressure (Poral), and d) peak-to-peak magnitudes of oral pressure (Poral_p-t-p) depending on the flow rate Q for: 1) phonation on the vowel [u:], 2) phonation through the tube into air, 3) phonation through the tube into water, and 4) blowing through the tube into water without phonation. 33 Figure 4. Measured maximum glottal area maxGA(t) and maximum glottal area derivative max(- dGA(t)/dt) evaluated during the closing phase of the glottis (MADR) as a function of the flow rate Q for: 1) phonation on the vowel [u:], 2) phonation through the tube into air, and 3) phonation through the tube into water. 34 Figure 5. Measured fundamental frequency F0 and the water bubbling frequency Fb depending on the flow rate Q for: 1) phonation on the vowel [u:], 2) phonation through the tube into air, 3) phonation through the tube into water, and 4) blowing through the tube into water without phonation. 35 Figure 6. Example of the construction of the relationships Ptrans(t) vs. GA(t) and Ptrans (t) vs. dGA(t)/dt for phonation on vowel [u:] for one period of the vocal folds’ self-oscillation showing: glottal area GA(t) (upper panel) and glottal area derivative dGA(t)/dt (2nd panel), transglottic pressure Ptrans(t) (3rd panel), and the resulting relationships (lower panel) (Q=0.20 l/s, F0=111 Hz, ∆A=0.031 mJ, ΔFP=ΔA F0 =3.42 mW). ). The numbers along the waveform GA(t) indicate the important time instants marked also on the relationships and on the high-speed VF images shown at the bottom. 36 Figure 7. Example of the construction of the relationships Ptrans(t) vs. GA(t) and Ptrans(t) vs. dGA(t)/dt for phonation through the resonance tube into air for one period of the vocal folds’ self-oscillation showing: glottal area GA(t) (upper panel), glottal area derivative dGA(t)/dt (2nd panel), transglottic pressure Ptrans(t) (3rd panel), and the resulting relationships (lower panel), with the images of the vibrating vocal folds also shown at the numbered time instants (Q=0.08 l/s, F0=94 Hz, ∆A=-0.054 mJ, ΔFP=ΔA F0 =-5.07 mW ). The numbers along the waveform GA(t) indicate the important time instants marked also on the relationships and on the high-speed VF images shown at the bottom. 37 Figure 8. Example of the construction of the relationship Ptrans(t) vs. GA(t) for phonation on the vowel [u:] when the VT is prolonged by the resonance tube with the distal end submerged 10 cm in water for one period of the vocal folds’ self-oscillation showing: glottal area GA(t) (upper panel), transglottic pressure Ptrans (t) (middle panel), and the resulting relationship Ptrans (t) vs. GA(t) (lower panel), with the images of the vocal folds’ vibration also shown at the numbered time instants (Q=0.08 l/s, F0=83 Hz, Fb=23 Hz). The numbers along the waveform GA(t) indicate the important time instants marked also on the relationships and on the high-speed VF images shown at the bottom. 38 Figure 9. Example of construction of the relationships Ptrans(t) vs. GA(t) and Ptrans(t) vs. dGA(t)/dt for phonation on the vowel [u:] when the VT is prolonged by the resonance tube with the distal end submerged 10 cm in water, for nine periods of the vocal folds’ self-oscillation showing: glottal area GA(t) (upper panel), glottal area derivative dGA(t)/dt (2nd panel), transglottic pressure Ptrans(t) (3rd panel), and the resulting relationships (bottom) (Q=0.08 l/s, F0=83 Hz, Fb=23 Hz, ∆A= -0.041 mJ, ΔFP=ΔA F0 =-3.41 mW). 39 Figure 10. Measured part of the flow power ΔFP=ΔA F0 consumed by the self-oscillating vocal folds for phonation on the vowel [u:] and its changes caused by phonation through the tube into air and into water in dependence on the input steady flow power defined as mean subglottic pressure times mean airflow rate (Psub Q). 40 Table 1. Comparison of the mean subglottal pressure (Psub), mean oral pressure (Poral), peak-to-peak oral pressure (Poral_p-t-p), frequencies for water bubbling (Fb), fundamental frequencies (F0), the first formant frequencies (F1) and 3 dB frequency bandwidths (B1) for phonations: a) on vowel [u:], b) into the resonance tube in air and c) with the distal end of the tube submerged 10 cm in water, measured on the model with hard vocal tract walls and in human subjects. Measurements are complemented with computational results for F1 and the acoustic-mechanical resonance frequencies (Fa-m) for a vocal tract with hard walls and with yielding walls. Psub [kPa] Poral [Pa] Poral_p-t-p [Pa] Fb [Hz] F0 [Hz] Fa-m [Hz] F1 B1 [Hz] [Hz] F1-F0 [Hz] F1-Fb [Hz] Fb-Fa-m [Hz] Phonation on vowel [u:] RESULTS FROM THE PRESENT STUDY) 1.45-2.24 44-123 170-1240 / 111 / 316 220 205 / / Computer model with hard VT / / / / / / 333 / / / Computer model with yielding VT / / / / / 31 410 / / / RESULTS FOR HUMANS Female subject**) 0.710 220 257 / 164 x 360 196 / / Phonation into the resonance tube in air RESULTS FROM THE PRESENT STUDY Physical model with hard VT***) 0.600 28 1284 / 90 / 105 99 15 / / 0.746 44 1740 / 91 / 105 14 / / Computer model with hard VT / / / / / / 97 / / / Computer model with yielding VT / / / / / 27 181 / / / RESULTS FOR HUMANS Male subject*) 0.600 25 210 / 100 x 190 90 / / Female subject**) 0.746 47 429 / 156 x 190 34 / / Phonation into the resonance tube into water (10 cm) RESULTS FROM THE PRESENT STUDY Physical model with hard VT***) 1.800 1082 2271 24 82 / 28 25 -54 4 / 1.480 1054 1166 20 79 / 28 -51 8 / Computer model with hard VT / / / / / / 28 / / / Computer model with yielding VT / / / / / 9 143 / / ≅2-5) RESULTS FOR HUMANS Male subject*) 1.800 1019 391 11.5 98 x 179 81 168 x Female subject**) 1.480 920 551 14 149 x ≅150 ≅1 136 x *) Horáček et al. (2017a, Table 3); **) Radolf et al. (2014, Tables 2 and 3 – normal phonation); ) see Figure 3; ***) measured data interpolated from the graphs in Figures 2, 3, and 5; ) related to results for humans; x… undefinable values.