Numerical and Experimental Investigations on Vocal Fold Approximation in Healthy and Simulated Unilateral Vocal Fold Paralysis

: We have developed a novel surgical/computational model for the investigation of unilateral vocal fold paralysis (UVFP) which will be used to inform future in silico approaches to improve surgical outcomes in type I thyroplasty. Healthy phonation (HP) was achieved using cricothyroid suture approximation on both sides of the larynx to generate symmetrical vocal fold closure. Following high-speed videoendoscopy (HSV) capture, sutures on the right side of the larynx were removed, partially releasing tension unilaterally and generating asymmetric vocal fold closure characteristic of UVFP (sUVFP condition). HSV revealed symmetric vibration in HP, while in sUVFP the sutured side demonstrated a higher frequency (10–11%). For the computational model, ex vivo magnetic resonance imaging (MRI) scans were captured at three conﬁgurations: non-approximated (NA), HP, and sUVFP. A ﬁnite-element method (FEM) model was built, in which cartilage displacements from the MRI images were used to prescribe the adduction, and the vocal fold deformation was simulated before the eigenmode calculation. The results showed that the frequency comparison between the two sides was consistent with observations from HSV. This alignment between the surgical and computational models supports the future application of these methods for the investigation of treatment for UVFP.


Introduction
The vocal folds are a pair of highly specialized tissues housed in the larynx, which when adducted into glottal controlled airflow, vibrate to produce voice. In humans, the vocal folds consist of five distinct layers; luminal stratified squamous epithelia, three lamina propria layers, and a basal ligament (vocalis muscle), which are broadly segmented into a cover and body as previously described [1]. Disordered vocal folds lead to many impairments, including poor voice quality, swallowing dysfunction, breathlessness, and aspiration [2]. Unilateral vocal fold paralysis (UVFP) is a voice disorder caused by injury of the recurrent laryngeal nerve or vagus nerve innervating the larynx [3], resulting in impaired adduction, one-sided vocal fold immobility, and a widened glottal gap. UVFP has an iatrogenic origin in around half of all cases [4], most frequently as a result of thyroidectomy or anterior cervical disk spinal fusion, but may also present idiopathically [5]. The symptomatic burden of UVFP on patients is significant with 53-100% of identified cases experiencing dysphonia [6][7][8], 60-75% dyspnea [9,10], and approximately 60% dysphagia [6], resulting in an increased risk of aspiration [11].
Surgical intervention is necessary for more than half of patients diagnosed with UVFP, with a primary goal of medializing the affected vocal fold to correct the glottal insufficiency [10]. There are two commonly used techniques for medialization laryngoplasty.
Injection laryngoplasty, a technique for which the vocal folds are medialized using an autologous or biologically inert substance that is injected directly into the thyroarytenoid muscle, is used to improve short-term voice outcomes (typically weeks to months depending on material and resorption rate) in patients who may experience spontaneous recovery or for whom more invasive surgery or general anesthesia is contraindicated [12][13][14][15][16][17][18]. Type I thyroplasty, a technique with which medialization is achieved by surgical insertion of a bio-tolerable implant (e.g., Silastic) or a stacked Gore-Tex ribbon through a window in the thyroid cartilage, provides longer duration results, but as with injection laryngoplasty, it is not without limitations, including: under/over medialization, size and shape variability when hand-carving implants, Gore-Tex compression over time, and inflammatory reactions. Inadequate medialization is reported relatively frequently, with 4.5-16% of type I thyroplasty procedures requiring surgical revision [19][20][21][22][23].
Previous efforts to surgically investigate type I thyroplasty are largely limited to studies in cadaver larynges assessing changes in biomechanics following implantation, or to histological analyses of larynges post-implantation; however, several in vivo and clinical studies evaluating phonatory measures have been described. In a canine model, Green et al. [24] simulated UVFP by cuffing the recurrent and superior laryngeal nerves and identified that arytenoid adduction + type I thyroplasty using the Hiroto method [25] increased vocal efficiency significantly more than augmentation alone, although the authors acknowledge that the longer posterior glottal chink in dogs vs. humans likely contributes to less optimal outcomes from implantation. Another study in a canine model identified significant impacts of both shape and medialization depth on fundamental frequency and aerodynamic power required for phonation and was also able to demonstrate reduced efficacy of type I thyroplasty in vagal paralysis vs. recurrent laryngeal nerve paralysis [26]. In a clinical population, type I thyroplasty in combination with cricothyroid suture approximation (also known as type IV thyroplasty) was demonstrated to improve video-stroboscopic, acoustic perceptual, and subjective dysphonia ratings in patients with unilateral superior laryngeal nerve weakness [27]. In a separate patient population with unilateral cricothyroid muscle paralysis, Shaw et al. reported that cricothyroid suture approximation alone was sufficient to improve multiple measures including glottic closure, frequency (fundamental, falsetto, basal), maximum phonatory time, and patient-reported perception [28] Subsequent investigation of cricothyroid suture approximation by our group in a rabbit model allowed for the development of a non-stimulatory method for medialization which achieved fundamental frequencies, vocal intensities, and aerodynamic measures consistent with modal phonation elicited through electrical stimulation [29]. The cricothyroid suture approximation model described by Novaleski et al. has the specific advantage of retaining vocal fold approximation ex vivo for imaging studies, something that cannot be achieved by electrical stimulation, such that high-speed videoendoscopy (HSV) and imaging from the same subject can be used to inform computational models of vocal fold vibration [30].
Besides experimental efforts, many numerical investigations using computational mechanics models have been conducted to study vocal fold vibration. Complementary to the physical models in experimental studies, computational models have the advantage of providing insight into the phonation process and predicting the vibration through mathematical modeling and computer simulation of the physical system. Such computational models usually include governing equations for the vocal fold biomechanics, the glottal aerodynamics, or a combination of both for fluid-structure interaction (FSI) simulations. Depending on the realism of vocal fold representation, there were lumpedmass models, two-dimensional, and three-dimensional finite-element models of the vocal fold tissues [30][31][32][33][34][35][36]. With the fast growth of computational power and medical imaging technologies, e.g., computed tomography (CT) and magnetic resonance imaging (MRI), patient-specific geometries have been incorporated into the FEM model of the vocal fold to represent the anatomy of individual patients or subjects [37][38][39]. Such models have potential applications in the surgical treatment of voice disorders like UVFP. More details about the development of computational models can be found in recent review papers [38,40].
While simulating vocal fold vibration, many FEM models have ignored the vocal fold adduction process and assumed equivalent stiffness properties for the vocal fold tissue to produce a suitable fundamental frequency. This ad hoc assumption may include asymmetric vibration, where two sides were assigned with different stiffness properties [41]. To study the effect of the approximation on vibration, the large tissue deformation due to adduction should be part of the FEM model to account for the effect of the structural stiffening (as opposed to material stiffening) on the dynamic behavior of the vocal fold. In the past, quite a few FEM models have been developed in combination with experimental testing to investigate vocal fold adduction, including both in vivo and ex vivo conditions [42][43][44][45][46][47][48], while the main interests of these studies were to determine the material properties of the vocal fold tissue. Some studies have looked at the effect of adduction or lengthening on the vocal fold structure [42,45,46]. However, those models only have simple vocal fold geometries and have not included details of cartilage movement. To accurately simulate vocal fold approximation, it is necessary to obtain quantitative information of cartilage displacement in the process. There have been efforts to segment and reconstruct laryngeal cartilages based on 3D image data [49,50]; however, this approach has not been used in FEM modeling of the vocal fold.
In this study, we constructed a subject-specific rabbit vocal fold model based on the excised MRI scans at different surgical vocal fold approximation conditions. From these scans including that from the non-approximated (NA) condition, we obtained the details of cartilage displacement and incorporated them into the FEM model as the boundary condition to simulate the vocal fold deformation and consequent structural stiffening. Compared with previous studies, our FEM model incorporates the subject-specific anatomy and the displacement of the cartilage from the MRI data. Furthermore, previous works predominantly focused on normal healthy vocal fold adduction, with the two sides producing symmetrical vibration. In the current study, we considered a simulated UVFP condition (sUVFP), in which the two vocal folds were asymmetrically approximated, and built the corresponding FEM model that has different structural stiffening due to deformation between the two sides. This model allows us to study the effect of the asymmetric adduction on the vocal fold vibration without having to make ad hoc assumptions of different material stiffness properties for the two sides.
In this work, we combined computational modeling with the experimental study of multiple rabbit larynges, and by comparing healthy phonation with simulated UVFP (sUVFP) we investigated the effect of the sUVFP on vocal fold vibration. The experimental data will validate the computational model, and the computational model will provide insight into the experimental observation in terms of vibration asymmetry. This combined experimental/computational approach and its alignment are demonstrated in this sUVFP study. Similar approaches will be utilized in the future to study the treatment approach for the UVFP.  [51,52]. Six New Zealand white breeder rabbits were included in this study, three for in vivo phonation studies with HSV analyses, and three for ex vivo MRI studies. The anesthesia regimen has been described elsewhere [53]. Briefly, rabbits were sedated with 20 mg/kg ketamine and 0.125 mg/kg of dexmedetomidine by intramuscular injection. To maintain anesthesia (as required), continuous rate infusion of 350 µg/kg/min ketamine and 1.65 µg/kg/min dexmedetomidine was provided intra-venously, via 27G catheter to the marginal vein. Heart rate, respiratory rate, temperature, and oxygen partial pressure were recorded for the duration of the procedure.

In Vivo Phonation
Non-stimulated phonation was performed as previously described [29] with minor modifications. In brief, the neck and chest of anesthetized rabbits were shaved and the rabbit was arranged in a supine position. A skin incision was made from sternum to submentum, and the fascia and muscle were dissected at the midline to expose the larynx and trachea. The trachea was bisected and an uncuffed endotracheal tube was inserted into the arboreal portion to provide direct oxygen to the lungs. A cuffed endotracheal tube was positioned in the upper portion of the trachea and provided humidified airflow at 37 • C to the subglottis. A pediatric laryngoscope was inserted via the mouth to suspend the larynx and visualize the vocal folds, and a combination of humidified airflow and cricothyroid suture approximation (see Section 2.1.3 for laryngeal conditions) was used to achieve phonation (see Figure 1). Phonation was not achieved in the NA condition. HSV was captured at 8000 frames per second in each laryngeal condition using a KayPENTAX FastCam MC2.1 system coupled to a 0 • rigid endoscope.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of included in this study, three for in vivo phonation studies with HSV analyses, and thr for ex vivo MRI studies. The anesthesia regimen has been described elsewhere [5 Briefly, rabbits were sedated with 20 mg/kg ketamine and 0.125 mg/kg of dexmedetom dine by intramuscular injection. To maintain anesthesia (as required), continuous rate i fusion of 350 μg/kg/min ketamine and 1.65 μg/kg/min dexmedetomidine was provid intravenously, via 27G catheter to the marginal vein. Heart rate, respiratory rate, tempe ature, and oxygen partial pressure were recorded for the duration of the procedure.

In Vivo Phonation
Non-stimulated phonation was performed as previously described [29] with min modifications. In brief, the neck and chest of anesthetized rabbits were shaved and t rabbit was arranged in a supine position. A skin incision was made from sternum to su mentum, and the fascia and muscle were dissected at the midline to expose the larynx an trachea. The trachea was bisected and an uncuffed endotracheal tube was inserted in the arboreal portion to provide direct oxygen to the lungs. A cuffed endotracheal tube w positioned in the upper portion of the trachea and provided humidified airflow at 37°C the subglottis. A pediatric laryngoscope was inserted via the mouth to suspend the laryn and visualize the vocal folds, and a combination of humidified airflow and cricothyro suture approximation (see Section 2.1.3 for laryngeal conditions) was used to achieve ph nation (see Figure 1). Phonation was not achieved in the NA condition. HSV was capture at 8000 frames per second in each laryngeal condition using a KayPENTAX FastCa MC2.1 system coupled to a 0° rigid endoscope.

Cricothyroid Suture Approximation Laryngeal Conditions
Cricothyroid suture approximation was performed as previously described [29] with minor modifications to achieve NA, HP, and sUVFP laryngeal conditions. Conditions were performed in a stepwise fashion to capture all conditions in each larynx ( Figure 2). NA: HSV of the vocal folds was captured in a neutral resting position. HP: Bilateral Vicryl 5.0 sutures (Ethicon) were used to suture the cricoid cartilage to the thyroid cartilage, resulting in lengthening of both vocal folds to approximate them medially. sUVFP: Sutures were removed from the right side of the larynx, resulting in tension being released from the right vocal fold and a "floppy" appearance characteristic of UVFP.
is supplied to the subglottis, in combination with bilateral cricothyroid suture approximation to move the vocal folds toward the midline, resulting in phonation. (b) Pictorial representation of endoscopic vocal fold views captured by HSV in three experimental conditions: non-approximated (NA), healthy phonation (HP), and simulated unilateral vocal fold paralysis (sUVFP).

Cricothyroid Suture Approximation Laryngeal Conditions
Cricothyroid suture approximation was performed as previously described [29] with minor modifications to achieve NA, HP, and sUVFP laryngeal conditions. Conditions were performed in a stepwise fashion to capture all conditions in each larynx ( Figure 2). NA: HSV of the vocal folds was captured in a neutral resting position. HP: Bilateral Vicryl 5.0 sutures (Ethicon) were used to suture the cricoid cartilage to the thyroid cartilage, resulting in lengthening of both vocal folds to approximate them medially. sUVFP: Sutures were removed from the right side of the larynx, resulting in tension being released from the right vocal fold and a "floppy" appearance characteristic of UVFP.
. Figure 2. Cricothyroid suture approximation conditions. HSV was captured in three conditions in vivo. The NA condition was captured before suturing to observe the natural resting vocal fold position. The HP condition was captured following bilateral sutures between the cricoid and thyroid cartilages, shortening the distance between these structures and lengthening the vocal folds to achieve approximation. The sUVFP condition was captured by releasing the right-sided sutures, allowing lengthening of the distance between the cricoid and thyroid unilaterally, leading to the right vocal fold having less tension applied and appearing "floppy".

Ex Vivo Cricothyroid Suture Approximation and MRI
Larynges were harvested for MRI in animals that had received no prior interventions. Excised larynges were scanned in the following experimental conditions: NA, HP, sUVFP (see Figure 2). Excised larynges were placed in a 10 mL syringe with perfluorocarbon oil and images were captured with a Bruker AV3HD 11.7 tesla / 89 mm vertical-bore microimaging system equipped with a 40-mm micro2.5 set capable of 1500 mT/m and Para-Vision 6.0.1 software (Bruker Biospin). Three-dimensional T2-weighted images were collected at an isotropic resolution of 60μ using a fast spin-echo sequence with the following parameters: TR/TE 1000/28 ms, 256 × 384 × 286 matrix, 16 × 16 mm field-of-view, and a RARE factor of 8.

Numerical Methodology
To numerically investigate the effect of vocal fold approximation on vibration in different conditions, we first reconstructed subject-specific larynges based on ex vivo experimentation. Due to the duration of the MRI scans to obtain appropriate resolution for computational modeling (approximately 4 h per scan) we are not able to gather this data in vivo due to concerns of protracted anesthesia, movement during scans, and general welfare concerns. It was not practicable to use larynges previously used for in vivo experimentation due to tissue inflammation and damage caused by phonation, such that tissue Figure 2. Cricothyroid suture approximation conditions. HSV was captured in three conditions in vivo. The NA condition was captured before suturing to observe the natural resting vocal fold position. The HP condition was captured following bilateral sutures between the cricoid and thyroid cartilages, shortening the distance between these structures and lengthening the vocal folds to achieve approximation. The sUVFP condition was captured by releasing the right-sided sutures, allowing lengthening of the distance between the cricoid and thyroid unilaterally, leading to the right vocal fold having less tension applied and appearing "floppy".

Ex Vivo Cricothyroid Suture Approximation and MRI
Larynges were harvested for MRI in animals that had received no prior interventions. Excised larynges were scanned in the following experimental conditions: NA, HP, sUVFP (see Figure 2). Excised larynges were placed in a 10 mL syringe with perfluorocarbon oil and images were captured with a Bruker AV3HD 11.7 tesla/89 mm vertical-bore micro-imaging system equipped with a 40-mm micro2.5 set capable of 1500 mT/m and ParaVision 6.0.1 software (Bruker Biospin). Three-dimensional T2-weighted images were collected at an isotropic resolution of 60 µ using a fast spin-echo sequence with the following parameters: TR/TE 1000/28 ms, 256 × 384 × 286 matrix, 16 × 16 mm field-of-view, and a RARE factor of 8.

Numerical Methodology
To numerically investigate the effect of vocal fold approximation on vibration in different conditions, we first reconstructed subject-specific larynges based on ex vivo experimentation. Due to the duration of the MRI scans to obtain appropriate resolution for computational modeling (approximately 4 h per scan) we are not able to gather this data in vivo due to concerns of protracted anesthesia, movement during scans, and general welfare concerns. It was not practicable to use larynges previously used for in vivo experimentation due to tissue inflammation and damage caused by phonation, such that tissue structure was likely to be altered by the time MRI imaging was performed. Therefore, three additional rabbit larynges were utilized for the model reconstructions.

Vocal Fold Reconstruction from MRI Scan
Computational laryngeal reconstruction was performed using MRI images from three samples that were scanned ex vivo. Using sample 1 as an example (see Figure 3), we first scanned the vocal fold in the NA condition. Following this, the MRI scans were imported to ITK-SNAP to manually segment the images. We considered the following structural components for model reconstruction: thyroid cartilage, arytenoid cartilages, cricoid cartilage, and soft tissue (vocal fold). Two distinct layers were segmented within the vocal folds due to known differences in tissue properties [54]: the pliable vocal fold cover, composed of the epithelium and superficial lamina propria, and the much stiffer vocal fold body, composed of the deep lamina propria and thyroarytenoid muscles. The layers were segmented in MRI images as depicted in Figure 3.
structure was likely to be altered by the time MRI imaging was performed. Therefore, three additional rabbit larynges were utilized for the model reconstructions.

Vocal Fold Reconstruction from MRI Scan
Computational laryngeal reconstruction was performed using MRI images from three samples that were scanned ex vivo. Using sample 1 as an example (see Figure 3), we first scanned the vocal fold in the NA condition. Following this, the MRI scans were imported to ITK-SNAP to manually segment the images. We considered the following structural components for model reconstruction: thyroid cartilage, arytenoid cartilages, cricoid cartilage, and soft tissue (vocal fold). Two distinct layers were segmented within the vocal folds due to known differences in tissue properties [54]: the pliable vocal fold cover, composed of the epithelium and superficial lamina propria, and the much stiffer vocal fold body, composed of the deep lamina propria and thyroarytenoid muscles. The layers were segmented in MRI images as depicted in Figure 3. Following segmentation, we output surface meshes for each structural component from ITK-SNAP and processed them with Mesh Mixer (Autodesk, San Francisco, CA, USA) to avoid uneven surfaces from the manual segmentation. However, this procedure leads to a separate surface mesh for each segmented component, and the two components adjacent to each other, e.g., the thyroid cartilage and the vocal fold body do not share the same surface on their meshes, and there exists a thin gap in between. To handle this problem, we output a surface mesh for the larynx as an entire segment from ITK-SNAP, and this mesh was also smoothed with MESH Mixer. Next, all the smoothed surface meshes are imported into COMSOL Multiphysics (COMSOL Inc., Burlington, MA, USA) to generate solid bodies and create unstructured volume meshes sequentially. A tetrahedron mesh with 10 nodes on each element was applied. Different structural components need to be identified on the entire larynx mesh as this mesh will be used as the FEM model in the following simulations. To do so, we compared the nodes on the volume mesh of the entire larynx with the nodes on the volume mesh of individual structural components and set a minimum distance criterion to determine whether the nodes of the entire larynx's mesh belong to a specific structural component. Using this approach, all the nodes on the entire larynx's mesh were marked as individual structural components as shown in Figure  4a and can be assigned with the material properties related to that component. Since the vocal fold cover is thin but important for the vocal fold vibration, this region is refined with denser mesh as shown in Figure 4b. In the mesh-independence study, four different unstructured meshes with 76,181, 88,334, 125,351, and 135,748 elements respectively were used to compute the eigenfrequency of approximated vocal folds in the HP condition. The eigenfrequency of the left vocal fold is 578 Hz, 588 Hz, 597 Hz, and 598 Hz, respectively, using the aforementioned four meshes. The third mesh resolution was selected in all the Following segmentation, we output surface meshes for each structural component from ITK-SNAP and processed them with Mesh Mixer (Autodesk, San Francisco, CA, USA) to avoid uneven surfaces from the manual segmentation. However, this procedure leads to a separate surface mesh for each segmented component, and the two components adjacent to each other, e.g., the thyroid cartilage and the vocal fold body do not share the same surface on their meshes, and there exists a thin gap in between. To handle this problem, we output a surface mesh for the larynx as an entire segment from ITK-SNAP, and this mesh was also smoothed with MESH Mixer. Next, all the smoothed surface meshes are imported into COMSOL Multiphysics (COMSOL Inc., Burlington, MA, USA) to generate solid bodies and create unstructured volume meshes sequentially. A tetrahedron mesh with 10 nodes on each element was applied. Different structural components need to be identified on the entire larynx mesh as this mesh will be used as the FEM model in the following simulations. To do so, we compared the nodes on the volume mesh of the entire larynx with the nodes on the volume mesh of individual structural components and set a minimum distance criterion to determine whether the nodes of the entire larynx's mesh belong to a specific structural component. Using this approach, all the nodes on the entire larynx's mesh were marked as individual structural components as shown in Figure 4a and can be assigned with the material properties related to that component. Since the vocal fold cover is thin but important for the vocal fold vibration, this region is refined with denser mesh as shown in Figure 4b. In the mesh-independence study, four different unstructured meshes with 76,181, 88,334, 125,351, and 135,748 elements respectively were used to compute the eigenfrequency of approximated vocal folds in the HP condition. The eigenfrequency of the left vocal fold is 578 Hz, 588 Hz, 597 Hz, and 598 Hz, respectively, using the aforementioned four meshes. The third mesh resolution was selected in all the cases, as further refinement did not significantly improve the result. In this mesh set, the thyroid cartilage, arytenoid cartilage, cricoid cartilage, vocal fold cover, and vocal fold body for Sample 1 have 16,510, 3232, 11,489, 13,021, and 81,099 elements respectively.

Vocal Fold Approximation Modeling
Once the NA larynx was reconstructed in Section 2.2.1, we were able to model vocal fold approximation using sutures in the HP and sUVFP conditions. Following the MRI scan for the NA condition, the HP and sUVFP cricothyroid suture approximation conditions were implemented in a stepwise fashion with consecutive MRI imaging.
Images of the HP condition were manually segmented and reconstructed with ITK-SNAP. Figure 5a,b show the reconstructed thyroid cartilage, arytenoid cartilages, and cricoid cartilage in NA and HP conditions, respectively. MRI imaging of the sUVFP condition was also segmented and reconstructed manually with ITK-SNAP. Laryngeal changes between HP and sUVFP are represented from Figure 5b to 5c. A comparison of the reconstructed laryngeal conditions revealed that the cricoid cartilage retains the same spatial orientation and shape between the conditions, so it was used as a reference to calculate the thyroid cartilage displacement from NA to HP, and from NA to sUVFP. These displacements were prescribed to the mesh nodes on the thyroid cartilage in the reconstructed whole larynx model shown in Figure 5. The cricoid cartilage remained fixed in the simulations. The Saint Venant-Kirchhoff model was adopted to represent the tissue properties for the vocal fold body and cover. With this model, we assume the vocal fold tissue to be isotropic with a linear relationship between the stress and the strain. However, the finite strain takes into account the nonlinear effects due to large displacements and large rotation [55,56]. Details for this hyper-elastic model can be found in our previous works [31]. With this procedure, we were able to model vocal fold approximation under the experimental conditions specified above.

Vocal Fold Approximation Modeling
Once the NA larynx was reconstructed in Section 2.2.1, we were able to model vocal fold approximation using sutures in the HP and sUVFP conditions. Following the MRI scan for the NA condition, the HP and sUVFP cricothyroid suture approximation conditions were implemented in a stepwise fashion with consecutive MRI imaging.
Images of the HP condition were manually segmented and reconstructed with ITK-SNAP. Figure 5a,b show the reconstructed thyroid cartilage, arytenoid cartilages, and cricoid cartilage in NA and HP conditions, respectively. MRI imaging of the sUVFP condition was also segmented and reconstructed manually with ITK-SNAP. Laryngeal changes between HP and sUVFP are represented from Figure 5b to Figure 5c. A comparison of the reconstructed laryngeal conditions revealed that the cricoid cartilage retains the same spatial orientation and shape between the conditions, so it was used as a reference to calculate the thyroid cartilage displacement from NA to HP, and from NA to sUVFP. These displacements were prescribed to the mesh nodes on the thyroid cartilage in the reconstructed whole larynx model shown in Figure 5. The cricoid cartilage remained fixed in the simulations. The Saint Venant-Kirchhoff model was adopted to represent the tissue properties for the vocal fold body and cover. With this model, we assume the vocal fold tissue to be isotropic with a linear relationship between the stress and the strain. However, the finite strain takes into account the nonlinear effects due to large displacements and large rotation [55,56]. Details for this hyper-elastic model can be found in our previous works [31]. With this procedure, we were able to model vocal fold approximation under the experimental conditions specified above.

Vocal Fold Approximation Modeling
Once the NA larynx was reconstructed in Section 2.2.1, we were able to model vocal fold approximation using sutures in the HP and sUVFP conditions. Following the MRI scan for the NA condition, the HP and sUVFP cricothyroid suture approximation conditions were implemented in a stepwise fashion with consecutive MRI imaging.
Images of the HP condition were manually segmented and reconstructed with ITK-SNAP. Figure 5a,b show the reconstructed thyroid cartilage, arytenoid cartilages, and cricoid cartilage in NA and HP conditions, respectively. MRI imaging of the sUVFP condition was also segmented and reconstructed manually with ITK-SNAP. Laryngeal changes between HP and sUVFP are represented from Figure 5b to 5c. A comparison of the reconstructed laryngeal conditions revealed that the cricoid cartilage retains the same spatial orientation and shape between the conditions, so it was used as a reference to calculate the thyroid cartilage displacement from NA to HP, and from NA to sUVFP. These displacements were prescribed to the mesh nodes on the thyroid cartilage in the reconstructed whole larynx model shown in Figure 5. The cricoid cartilage remained fixed in the simulations. The Saint Venant-Kirchhoff model was adopted to represent the tissue properties for the vocal fold body and cover. With this model, we assume the vocal fold tissue to be isotropic with a linear relationship between the stress and the strain. However, the finite strain takes into account the nonlinear effects due to large displacements and large rotation [55,56]. Details for this hyper-elastic model can be found in our previous works [31]. With this procedure, we were able to model vocal fold approximation under the experimental conditions specified above.

In Vivo Phonation Test Results
HSV was captured for each subject in three conditions: NA, HP, and sUVFP. Figure 6 depicts the HSV analysis from subject 1. Using a proprietary MATLAB ® (Mathworks, Natick, MA) code generously provided by Dr. Dimitar Deliyski at Michigan State University and modified in-house, the HSV was converted for digital kymograph analysis (see Figure 6b,c), wherein a graphical representation of the vocal fold waveform with time as the x-axis was generated. The HP condition (Figure 6b), demonstrated symmetrical frequency between the left and right vocal folds. The sUVFP condition demonstrated greater frequency on the sutured side (left) than the unsutured side (right).

In Vivo Phonation Test Results
HSV was captured for each subject in three conditions: NA, HP, and sUVFP. Figure  6 depicts the HSV analysis from subject 1. Using a proprietary MATLAB ® (Mathworks, Natick, MA) code generously provided by Dr. Dimitar Deliyski at Michigan State University and modified in-house, the HSV was converted for digital kymograph analysis (see Figure 6b,c), wherein a graphical representation of the vocal fold waveform with time as the x-axis was generated. The HP condition (Figure 6b), demonstrated symmetrical frequency between the left and right vocal folds. The sUVFP condition demonstrated greater frequency on the sutured side (left) than the unsutured side (right). Vibration frequencies calculated from HSV digital kymograph (Table 1) demonstrated predictable frequency differences between left and right vocal folds under the different experimental conditions. HP vocal folds achieved symmetrical vibration representative of normal phonation. Vocal folds with a unilateral suture (sUVFP) successfully modeled UVFP with asymmetrical vibration observed, resulting in a 10-11% difference in frequency between left and right vocal fold.

Eigenfrequency of Approximated Vocal Fold
Using numerical methodology discussed in Section 2.2, subject-specific vocal fold approximation was modeled for different laryngeal conditions. During in vivo experiments, vibration frequency in the HP condition ranged from 500 Hz to 900 Hz, consistent with the range observed in previous studies of rabbit phonation [30,37]. For computational modeling, we assumed constant material properties for the soft tissue including Young's moduli, Poisson's ratio, and density as listed in Table 2. The Young's modulus of the vocal fold cover used here falls in the range of previous experimental measurements [57,58]. Vibration frequencies calculated from HSV digital kymograph (Table 1) demonstrated predictable frequency differences between left and right vocal folds under the different experimental conditions. HP vocal folds achieved symmetrical vibration representative of normal phonation. Vocal folds with a unilateral suture (sUVFP) successfully modeled UVFP with asymmetrical vibration observed, resulting in a 10-11% difference in frequency between left and right vocal fold.

Eigenfrequency of Approximated Vocal Fold
Using numerical methodology discussed in Section 2.2, subject-specific vocal fold approximation was modeled for different laryngeal conditions. During in vivo experiments, vibration frequency in the HP condition ranged from 500 Hz to 900 Hz, consistent with the range observed in previous studies of rabbit phonation [30,37]. For computational modeling, we assumed constant material properties for the soft tissue including Young's moduli, Poisson's ratio, and density as listed in Table 2. The Young's modulus of the vocal fold cover used here falls in the range of previous experimental measurements [57,58]. The cartilages had negligible deformation and were thus assigned with much higher stiffness. These parameters would lead to a dominant eigenfrequency of around 600 Hz for the vocal fold in the HP condition. This frequency is consistent with the average frequency across the rabbit larynx samples in two previous studies from our group [29,30]. The vibration frequency in the present study had more variations among individuals. However, the differences between the experiment and computation did not impact analysis, as we sought to compare left-right vibratory frequency as opposed to making one-on-one comparisons between the two models. A two-step process was used to compute the eigenfrequency of the approximated vocal fold. First, a static deformation was simulated for vocal fold adduction by the prescribed thyroid displacement according to the HP or sUVFP condition. As a result of this deformation, the vocal fold was pre-stressed. Then, the eigenfrequency and the eigenmode (or the vibration mode) were computed at the pre-stressed condition. Figure 7 shows vocal fold approximation under different conditions for Sample 1. Following reconstruction of the NA condition (Figure 7a), thyroid cartilage displacement for HP or sUVFP conditions was applied (Figure 7b or Figure 7c respectively). Using this subject-specific approach, we were able to model symmetrical approximation of the vocal folds for HP and asymmetrical approximation in sUVFP. The cartilages had negligible deformation and were thus assigned with much higher stiffness. These parameters would lead to a dominant eigenfrequency of around 600 Hz for the vocal fold in the HP condition. This frequency is consistent with the average frequency across the rabbit larynx samples in two previous studies from our group [29,30]. The vibration frequency in the present study had more variations among individuals. However, the differences between the experiment and computation did not impact analysis, as we sought to compare left-right vibratory frequency as opposed to making one-on-one comparisons between the two models. 0.4 1000 A two-step process was used to compute the eigenfrequency of the approximated vocal fold. First, a static deformation was simulated for vocal fold adduction by the prescribed thyroid displacement according to the HP or sUVFP condition. As a result of this deformation, the vocal fold was pre-stressed. Then, the eigenfrequency and the eigenmode (or the vibration mode) were computed at the pre-stressed condition. Figure  7 shows vocal fold approximation under different conditions for Sample 1. Following reconstruction of the NA condition (Figure 7a), thyroid cartilage displacement for HP or sUVFP conditions was applied (Figure 7b or 7c respectively). Using this subject-specific approach, we were able to model symmetrical approximation of the vocal folds for HP and asymmetrical approximation in sUVFP. The eigenfrequency and eigenmode calculations were done using the default algorithm in COMSOL Multiphysics. Multiple modes were obtained from this simulation. The vocal vibration is mainly related to the first two modes for each side of the vocal fold [59]. The higher-frequency modes involve out-of-phase vibration patterns on the medial surface of the vocal fold and were not found from the HSV results (see supplementary material for those mode shapes and frequencies). The first mode has the maximal displacement in the axial direction (Y-axis in Figure 8), while the second mode has the maximal displacement in the direction transverse to the midline (X-axis in Figure 8). The second mode was selected to represent the vibration pattern of the vocal fold under phonation by comparing it with the HSV results. This mode is shown in Figure 8 for sample 1. The corresponding eigenfrequency of the selected eigenmode is deemed the fundamental frequency of vibration. The eigenfrequencies for all three samples are listed in Table 3 for the NA, HP, and sUVFP conditions. The left and right sides are separately listed to compare the symmetry of vibration. From the results, the left and right sides had similar eigenfrequencies in the NA and HP conditions, while the sutured side had a 7-9% higher eigenfrequency than the unsutured side in the sUVFP condition. Furthermore, both the NA and The eigenfrequency and eigenmode calculations were done using the default algorithm in COMSOL Multiphysics. Multiple modes were obtained from this simulation. The vocal vibration is mainly related to the first two modes for each side of the vocal fold [59]. The higher-frequency modes involve out-of-phase vibration patterns on the medial surface of the vocal fold and were not found from the HSV results (see supplementary material for those mode shapes and frequencies). The first mode has the maximal displacement in the axial direction (Y-axis in Figure 8), while the second mode has the maximal displacement in the direction transverse to the midline (X-axis in Figure 8). The second mode was selected to represent the vibration pattern of the vocal fold under phonation by comparing it with the HSV results. This mode is shown in Figure 8 for sample 1. The corresponding eigenfrequency of the selected eigenmode is deemed the fundamental frequency of vibration. The eigenfrequencies for all three samples are listed in Table 3 for the NA, HP, and sUVFP conditions. The left and right sides are separately listed to compare the symmetry of vibration. From the results, the left and right sides had similar eigenfrequencies in the NA and HP conditions, while the sutured side had a 7-9% higher eigenfrequency than the unsutured side in the sUVFP condition. Furthermore, both the NA and sUVFP conditions have much higher frequencies than the NA condition. The frequencies in the HP and sUVFP conditions show good agreement with the experimental findings in Section 3.1.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 10 of 15 sUVFP conditions have much higher frequencies than the NA condition. The frequencies in the HP and sUVFP conditions show good agreement with the experimental findings in Section 3.1.

Discussion
In this study, in vivo phonation tests and subject-specific computation modeling were used to investigate the changes in vibratory frequency in three experimental conditions. We found a difference in vibratory patterns between HP, as simulated by bilateral suture cricothyroid approximation, and sUVFP, as simulated by unilateral suture cricothyroid approximation. HP demonstrated symmetrical vibration while the sUVFP condition demonstrated asymmetry in left-right frequency. From this investigation, we concluded that the cricothyroid suture approximation approach implemented to generate the different laryngeal conditions in this experiment was valid to simulate both normal phonation and UVFP. The subjects demonstrated a wide range of vibration frequencies (~600-900 Hz), as shown in Table 1; however, this variability is consistent with previous frequencies observed in non-stimulated and electrically stimulated phonation in rabbits [29,30,39] and is likely attributable to inter-individual differences in tissue properties and laryngeal geometries. Given the high level of variability observed within species, one set of tissue properties will not adequately represent all subjects, illustrating a real need for subjectspecific modeling to improve accuracy in surgical planning.
In the current sUVFP condition, the unsutured side's frequency was only around 10% lower than the sutured side, while both sides were much higher than the NA condition. This implies that the unsutured side in the sVUFP condition still experienced significant stretching and stiffening as a consequence of the contralateral suture. Further examination

Discussion
In this study, in vivo phonation tests and subject-specific computation modeling were used to investigate the changes in vibratory frequency in three experimental conditions. We found a difference in vibratory patterns between HP, as simulated by bilateral suture cricothyroid approximation, and sUVFP, as simulated by unilateral suture cricothyroid approximation. HP demonstrated symmetrical vibration while the sUVFP condition demonstrated asymmetry in left-right frequency. From this investigation, we concluded that the cricothyroid suture approximation approach implemented to generate the different laryngeal conditions in this experiment was valid to simulate both normal phonation and UVFP. The subjects demonstrated a wide range of vibration frequencies (~600-900 Hz), as shown in Table 1; however, this variability is consistent with previous frequencies observed in non-stimulated and electrically stimulated phonation in rabbits [29,30,39] and is likely attributable to inter-individual differences in tissue properties and laryngeal geometries. Given the high level of variability observed within species, one set of tissue properties will not adequately represent all subjects, illustrating a real need for subject-specific modeling to improve accuracy in surgical planning.
In the current sUVFP condition, the unsutured side's frequency was only around 10% lower than the sutured side, while both sides were much higher than the NA condition. This implies that the unsutured side in the sVUFP condition still experienced significant stretching and stiffening as a consequence of the contralateral suture. Further examination of the cartilage position in Figure 5 shows that both sides of the vocal fold were stretched even after one side's suture was released, and this was due to residual downward and lateral tilting of the released thyroid cartilage toward the contralateral suture. While the thyroid cartilage is flexible, it is a single piece of cartilage spanning the full breadth of the larynx, thus manipulation of one side cannot be achieved with total independence from the other. This type of asymmetric approximation is different than that observed in clinical UVFP; however, this technique successfully achieves significant differences between the left and right vocal fold and therefore creates a useful in vivo phonation paradigm that will allow us to perform future combined experimental/computational investigation into the surgical restoration of symmetric vibration.
The novel surgical and computational methods proposed in this study have several major benefits when compared to traditional techniques for modeling phonation and UVFP. Suture approximation is reversible and minimizes risks of airway patency, aspiration, and long-term survival associated with RLN/SLN resection in animal models [60]. Additionally, suture approximation allows for imaging of the vocal folds at rest, during bilateral phonation, and in sUVFP in a fast and minimally invasive process that can easily be replicated in vivo and ex vivo allowing for maximum translatability between the two. Using MRI scan and 3D reconstruction provides precise information about the cartilage movement during vocal fold adduction. With this information, the static deformation in the FEM simulation allows the model to account for the prestress produced by vocal fold adduction in the tissue. Thus, the subsequent eigenmode calculation captures the frequency and vibration mode of the approximated vocal fold in either bilateral adduction or unilateral adduction. Using this adduction-then-eigenmode simulation approach, we can estimate the vibration frequency without ad hoc assumptions for the tissue stiffening caused by the adduction, which may be generalized to calculate eigenfrequencies in the UVFP population. Future work will utilize this data to optimize surgical planning for type I thyroplasty, the gold standard surgical intervention in this disorder [2].
This study was limited in several ways. The inability to use the same subjects for HSV and anatomic imaging due to concerns for localized tissue changes following in vivo phonation impeded our ability to make direct comparisons of eigenfrequencies and vibratory frequencies within subjects. Thus, we adjusted for the lack of corresponding experimental frequency results by normalizing the eigenfrequency to an averaged frequency obtained from HSV in a sample outside of this study (600 Hz), which was lower than that observed in vivo during the surgical experiments in this study (760 Hz). As such, it is plausible that Young's moduli, Poisson's ratio, and tissue density have been overestimated or underestimated in individual subjects. Vocal fold vibration is essentially a fluid-structure interaction process between the glottal airflow and the soft tissue. A full 3D FSI simulation would provide the most accurate prediction, and we are working toward that direction with our current research. It has been demonstrated that the vocal fold vibration is highly related to its eigenfrequency [59,61]. Given the costly computational resources and time required for the high-fidelity simulations, utilizing the eigenfrequency model provides an efficient and meaningful way for such a combined experimental/computational initial investigation of the UVFP.
The use of ex vivo scanning limits translation of the present approach to human subjects, therefore future methods will rely on in vivo scans. Finally, while cricothyroid suture approximation shows promise for in vivo and ex vivo studies of vocal fold approximation and UVFP, variables such as suture placement and tension are hard to standardize, necessitating significant surgical skill, and this may result in lower precision than observed with other more established techniques (e.g., RLN resection to model UVFP, electrical stimulation for phonation.

Conclusions
In the present study, we describe a method to investigate vocal fold approximation in normal phonation and simulated UVFP conditions experimentally and numerically. Using a non-stimulated cricothyroid suture approximation model we previously developed to elicit sustained vocal fold phonation in rabbits, we implemented a novel experimental condition to simulate UVFP. Using this method, we demonstrated symmetrical phonation simulating normal phonation and asymmetrical phonation simulating UVFP.
A 3-dimensional MRI sequence was also performed in three experimental conditions: NA, HP, and sUVFP. These images were used to build a subject-specific vocal fold re-construction of laryngeal cartilages and two layers of soft tissue: vocal fold cover and vocal fold body. Using FEM, we simulated vocal fold deformation for each condition by incorporating the cartilage displacement. The subsequent eigenmode calculation accounts for the prestress in the tissue due to deformation. Similar to the surgical experimental results, the computational model also demonstrated similar eigenfrequencies between left and right in the HP condition and asymmetrical eigenfrequencies in the sUVFP condition.
The strong agreement between the changes in vibration frequency between conditions as observed in the HSV and the changes in eigenfrequency as calculated in the numerical analysis of the subject-specific model demonstrate that these novel methodologies are a promising innovation for the investigation of vocal fold phonation in healthy and UVFP conditions.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. Final datasets will be publicly available when complete, per the resource sharing plan described in NIDCD award #5R01DC016236. Communication Disorders of the National Institute of Health under award number: 5R01DC016236, PI: Rousseau. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Conflicts of Interest:
The authors declare no conflict of interest in the present study.