# An Euler-Bernoulli-Type Beam Model of the Vocal Folds for Describing Curved
and Incomplete Glottal Closure Patterns

Mohamed A. Serry  Mechanical and Mechatronics Engineering, University of
Waterloo, Waterloo, Ontario N2L 3G1, Canada  Gabriel A. Alzamendi  Institute
for Research and Development on Bioengineering and Bioinformatics (IBB),
CONICET-UNER, Oro Verde, Entre Ríos 3100, Argentina  Matías Zañartu
Department of Electronic Engineering, Universidad Técnica Federico Santa
María, Valparaíso, Chile  Sean D. Peterson  Mechanical and Mechatronics
Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

###### Abstract

Incomplete glottal closure is a laryngeal configuration wherein the glottis is
not fully obstructed prior to phonation. It has been linked to inefficient
voice production and voice disorders. Various incomplete glottal closure
patterns can arise and the mechanisms driving them are not well understood. In
this work, we introduce an Euler-Bernoulli composite beam vocal fold (VF)
model that produces qualitatively similar incomplete glottal closure patterns
as those observed in experimental and high-fidelity numerical studies, thus
offering insights in to the potential underlying physical mechanisms. Refined
physiological insights are pursued by incorporating the beam model into a VF
posturing model that embeds the five intrinsic laryngeal muscles. Analysis of
the combined model shows that co-activating the lateral cricoarytenoid (LCA)
and interarytenoid (IA) muscles without activating the thyroarytenoid (TA)
muscle results in a bowed (convex) VF geometry with closure at the posterior
margin only; this is primarily attributed to the reactive moments at the
anterior VF margin. This bowed pattern can also arise during VF compression
(due to extrinsic laryngeal muscle activation for example), wherein the
internal moment induced passively by the TA muscle tissue is the predominant
mechanism. On the other hand, activating the TA muscle without incorporating
other adductory muscles results in anterior and mid-membranous glottal
closure, a concave VF geometry, and a posterior glottal opening driven by
internal moments induced by TA muscle activation. In the case of initial full
glottal closure, the posterior cricoarytenoid (PCA) muscle activation cancels
the adductory effects of the LCA and IA muscles, resulting in a concave VF
geometry and posterior glottal opening. Furthermore, certain maneuvers
involving co-activation of all adductory muscles result in an hourglass
glottal shape due to a reactive moment at the anterior VF margin and moderate
internal moment induced by TA muscle activation. These findings have
implications regarding potential laryngeal maneuvers in patients with voice
disorders involving imbalances or excessive tension in the laryngeal muscles
such as muscle tension dysphonia.

Keywords— Vocal folds; Incomplete glottal closure; Muscle tension dysphonia;
Euler-Bernoulli beam.

##  1 Introduction

The configuration of the vocal folds (VFs), a cornerstone of voice production,
is determined by the particular combination of activated intrinsic and
extrinsic laryngeal muscles. Nominally, the VFs are completely adducted prior
to the onset of phonation, and their interaction with the air flow driven by
the lungs results in vibrations and consequent acoustic waves, which forms the
basis of voiced speech. In some scenarios, complete glottal closure is not
attained, which can result in inefficient voice production (Zañartu et al., ,
2014), and, in some cases, stress concentrations in the VFs that may lead to
VF trauma (Dejonckere and Kob, , 2009). Hence, incomplete glottal closure is
often linked to disorders that are associated with inefficiencies in, or
damage to, the vocal mechanism, including Parkinson’s disease (Hanson et al.,
, 1984)111Parkinson’s disease is a relatively prevalent disorder, affecting
the human central, peripheral, and enteric nervous systems (Braak and Braak, ,
2000). and muscle tension dysphonia (MTD) (Morrison and Rammage, , 1993)222MTD
(Morrison and Rammage, , 1993), also known as non-phonotraumatic hyperfunction
(Hillman et al., , 2020), is a class of voice disorders associated with misuse
of the vocal mechanisms without the presence of organic changes in the vocal
organs, leading to low speech quality and vocal fatigue, with a wide range of
symptoms and patterns, including excessive/unbalanced activation of intrinsic
and extrinsic laryngeal muscles (Roy, , 2008; Hocevar-Boltezar et al., ,
1998), supraglottal compression (Morrison and Rammage, , 1993), and abnormal
fundamental frequency (Nguyen et al., , 2009; Altman et al., , 2005)..

Incomplete glottal closure333In this study we refer to the glottal
configuration of the VFs at rest immediately prior to phonation initiation,
identifying any gaps between the folds as incomplete glottal closure. In
clinical settings, incomplete glottal closure typically refers to gaps between
the folds when the VFs are at their maximum glottal closure phase during
phonation (Södersten et al., , 1995; Nguyen et al., , 2009). Our definition
herein isolates laryngeal factors, which are the focus of this study, by
dismissing the dynamics of VF vibrations. comes in various patterns (Morrison
and Rammage, , 1993; Nguyen et al., , 2009; Södersten et al., , 1995) as shown
schematically in Figure 1, including bowed shape: a glottal pattern wherein
the left and right VF geometries are convex with a gap at the mid-membranous
portion; posterior glottal opening: full glottal closure is achieved in the
anterior and mid-membranous regions only, leaving the posterior margin open;
and hourglass glottal configuration: a pattern with anterior and posterior
gaps and potential VF contact in the mid-membranous region. There exist other
incomplete glottal closure patterns, sharing similarities with those mentioned
above, such as spindle-shaped glottis, and anterior opening (see Rajaei et
al.,  (2014); Södersten et al.,  (1995)). These latter patterns are not
addressed in this work.

Figure 1: Schematic diagram (superior view) of some incomplete and curved
glottal closure patterns that are observed clinically.

There exist several experimental and clinical studies in the literature that
attempt to elucidate, at least in part, some of the laryngeal mechanisms
associated with curved and incomplete glottal closure patterns. Based upon
inspection of cadaver larynges, Morrison and Rammage,  (1993) posited that a
posterior glottal opening is associated with excessive activation of the
posterior cricoarytenoid (PCA) muscle. Choi et al., 1993b conducted
experimental investigations of excised canine models and found that activating
the thyroarytenoid (TA) muscle, while keeping other adductory laryngeal
muscles inactive, leads to anterior and mid-membranous glottal closure,
whereas the posterior glottis stays open, thus resulting in a closure pattern
similar to a posterior glottal opening. On the other hand, they found that
when the TA muscle is relaxed and other adductory muscles (LCA/IA) are
activated, closure is achieved only at the posterior margins of the VFs with
mid-membranous opening (Choi et al., 1993b, ; Chhetri and Neubauer, , 2015),
thus leading to a bowed configuration as seen in Figure 1. Moreover, complete
glottal closure of excised canine larynges was attained via co-activation of
all adductory muscles. More recent clinical investigations using refined
experimental setups (see, e.g., Chhetri and Neubauer,  (2015)) further confirm
these observations. Chhetri et al.,  (2012) conducted a parametric study of
the effects of intrinsic laryngeal muscle activation, modulated by graded
stimulation, on the pre-phonatory posture of a canine model. In addition to
confirming the findings of Choi et al., 1993a , Chhetri et al.,  (2012) found
that when keeping the TA muscle activation at a constant level and increasing
activation of the cricothyroid muscle (CT), the glottal area increases and the
medial bulging caused by the TA muscle activation is reduced. In addition,
they observed that when keeping the LCA and IA muscle activation at constant
levels and increasing CT muscle activation, the glottis starts to open
posteriorly and the glottal area increases. Interestingly, the authors found
that with certain muscular executions involving co-activation of the LCA, IA,
and TA muscles, the glottis exhibits an hourglass shape (see Chhetri et al.,
(2012, Fig. 10)).

Besides the aforementioned clinical and experimental works, there exist
numerical studies that shed some light onto curved and incomplete glottal
closure patterns. Hunter et al.,  (2004) developed one of the early three-
dimensional VF posturing models, where adductory and abductory muscles are
incorporated, showing that full activation of the LCA muscle induces
nonuniform curvature of the medial surface. Dejonckere and Kob,  (2009)
studied the influence of incomplete glottal closure patterns (with linear and
curved VF geometries) on VF vibrations using a multi-mass model, showing that
some resting incomplete glottal closure configurations may induce localized VF
impact, which they hypothesized to be a potential underlying mechanism
inducing VF trauma, especially in females. However, the authors did not study
the laryngeal maneuvers that induce these resting glottal shapes. Yin and
Zhang,  (2014, 2016) conducted numerical simulations using high-fidelity
numerical models, showing that posterior glottal opening occurs with the sole
activation of the TA muscle, mid-membranous opening when the LCA and IA
muscles are co-activated (without incorporating the TA muscle), and full
glottal closure when all adductors are co-activated, in agreement with the
aforementioned clinical observations. In a more recent study, Geng et al.,
(2020) proposed a detailed physiologically accurate finite-element posturing
model, based on MRI scan images of a canine larynx. Even though the study does
not study the glottal geometry, it provides useful insights into how
synergistic activation of laryngeal muscles exhibits complex interaction with
laryngeal variables (e.g., VF strain, rotation and translation of arytenoid
cartilages, and glottal area). Recently, research interest has also been
directed towards investigating how activating laryngeal muscles alters the VF
medial surfaces (Pillutla et al., , 2022).

Despite these valuable efforts, a clear picture of the physical mechanisms
inducing different glottal patterns remains elusive. It is challenging to
isolate and control the factors underlying posturing mechanics experimentally.
Moreover, high-fidelity numerical models, despite their accuracy in
replicating physiological laryngeal postures, do not provide clear intuitive
understanding of the mechanics of posturing, and typically suffer from high
computational costs. We hypothesize that the non-homogeneous structure of the
VFs, which comprise overlapping tissue layers with different mechanical and
geometrical properties (Titze and Alipour, , 2006), underlies, in part, the
different glottal shapes displayed in Figure 1. As such, we propose an Euler-
Bernoulli composite beam model of the VFs to elucidate some of the mechanisms
underlying glottal patterns prior to phonation. Beam models have been utilized
previously to explore VF vibrations and phonation fundamental frequency (Titze
and Hunter, , 2004; Zhang et al., , 2007). We opt for this relatively simple
modeling framework to facilitate exploration of the mechanisms underlying the
resulting glottal shapes. To gain refined physiological insights into how
intrinsic laryngeal muscles may influence glottal geometry, the proposed model
is integrated with the muscle-controlled posturing model of Titze and Hunter,
(2007).

The organization of this work is as follows: a detailed derivation of the
composite beam VF model is introduced in Section 2; analysis is conducted in
Section 3; numerical simulations of the integrated beam and posturing model
are presented in Section 4; Section 5 presents discussion of the results; and
the study is concluded in Section 6.

##  2 Model development

Herein, we propose a static Euler-Bernoulli-type composite beam model (see,
for example, Bauchau and Craig,  (2009)) for the VFs, with the different VF
layers represented by strata in the beam. For simplicity we assume symmetry
with respect to the medial plane; hence, we consider only one (the left) VF. A
schematic representation of the composite beam model is shown in Figure 2.

Figure 2: Schematic diagram (superior view) of the VF composite beam model.

The VF beam model consists of three layers: (1) the mucosa with depth
$d_{\mathrm{muc}}$ and cross-sectional area $A_{\mathrm{muc}}$; (2) the vocal
ligament with depth $d_{\mathrm{lig}}$ and cross-sectional area
$A_{\mathrm{lig}}$; and (3) the thyroarytenoid (vocalis) muscle with depth
$d_{\mathrm{ta}}$ and cross-sectional area $A_{\mathrm{ta}}$. For the sake of
compact presentation, we define the index set

$\mathcal{I}=\\{\mathrm{muc},\mathrm{lig},\mathrm{ta}\\},$ (1)

where $\mathrm{muc}$, $\mathrm{lig}$, and $\mathrm{ta}$ refer to the mucosa,
ligament, and TA muscle tissue, respectively. We assume that each layer has a
uniform rectangular cross-section and the layer thicknesses (in the inferior-
superior direction) are equal and denoted by $b$; thus, layer depth can be
computed as $d_{i}=A_{i}/b$ for $i\in\mathcal{I}$.

Our modelling framework assumes that VF deformation consists of (a)
potentially large longitudinal stretching/compression with uniform strain, and
(b) modest bending due to the induced moments inside the VF. Let $L_{0}$
denote the resting VF length and $L$ denote the VF length after longitudinal
deformation due to the associated nominal uniform strain $\bar{\varepsilon}$;
that is,

$L=(1+\bar{\varepsilon})L_{0}.$ (2)

We assume the nominal strain is known a priori444Such as from the two-
dimensional posturing model of Titze and Hunter,  (2007), which incorporates
the mechanics of the arytenoid cartilages and cricothyroid joints, and relates
them to VF strain, see the discussion in Section 4..

Let $x\in[0,L]$ denote the position along the deformed VF configuration (after
applying strain $\bar{\varepsilon}$) relative to the anterior VF margin, and
$r$ denote the depth position along the axis perpendicular to the VF axis
relative to the base of the TA muscle (see Figure 2). Consider a plane VF
cross-section at position $x$, and let $y_{\mathrm{muc}}$ denote the relative
position along the $r$-axis with respect to the geometrical center of the
mucosa (i.e., $y_{\mathrm{muc}}=0$ corresponds to the geometrical center of
the mucosal cross-section). Similarly, let $y_{\mathrm{lig}}$ and
$y_{\mathrm{ta}}$ be analogous coordinates for the ligament and TA muscle,
respectively (see Figure 2). Note that the range of $y_{i}$ is
$[-d_{i}/2,d_{i}/2]$ for $i\in\mathcal{I}$.

Let $w(x)$ denote the transverse deflection of the beam (in the
$r$-direction). Moreover, let $u_{i}(x,y_{i})$, $i\in\mathcal{I}$, denote the
longitudinal displacement of the $i^{\mathrm{th}}$ VF layer, where
longitudinal displacements are with respect to the deformed VF configuration
under $\bar{\varepsilon}$. In addition, let $\bar{u}_{i}(x)=u_{i}(x,y_{i}=0)$,
$i\in\mathcal{I}$, denote the longitudinal displacement at the center of the
$i^{\mathrm{th}}$ layer. Under Euler-Bernoulli beam theory (see, for example,
Bauchau and Craig,  (2009)), the longitudinal displacement functions can be
written as

$u_{i}=\bar{u}_{i}-y_{i}w^{\prime},~{}i\in\mathcal{I},$ (3)

where the prime symbol denotes differentiation with respect to $x$. Continuity
of displacement fields necessitates that
$u_{\mathrm{muc}}(x,-d_{\mathrm{muc}}/2)=u_{\mathrm{lig}}(x,d_{\mathrm{lig}}/2)$
and
$u_{\mathrm{lig}}(x,-d_{\mathrm{lig}}/2)=u_{\mathrm{ta}}(x,d_{\mathrm{ta}}/2)$
for all $x\in[0,L]$, which yields the conditions

$\begin{split}\bar{u}_{\mathrm{lig}}&=\bar{u}_{\mathrm{muc}}+\frac{1}{2}(d_{\mathrm{lig}}+d_{\mathrm{muc}})w^{\prime},\\\
\bar{u}_{\mathrm{ta}}&=\bar{u}_{\mathrm{muc}}+\frac{1}{2}(d_{\mathrm{ta}}+2d_{\mathrm{lig}}+d_{\mathrm{muc}})w^{\prime}.\end{split}$
(4)

Given longitudinal displacement in the $i^{\mathrm{th}}$ layer with respect to
the deformed configuration under $\bar{\varepsilon}$, the total strain in that
layer is given by555Consider an infinitesimal line element $\mathrm{d}x_{0}$
that experiences a composition of two deformations: the first is longitudinal
deformation with associated uniform normal strain $\varepsilon_{0}$, and the
second is due to a longitudinal displacement field $u$ (with respect to the
configuration after applying the strain $\varepsilon_{0}$). The length of the
line element after applying strain $\varepsilon_{0}$ is
$\mathrm{d}x_{1}=(1+\varepsilon_{0})\mathrm{d}x_{0}$
($\mathrm{d}x_{0}=\mathrm{d}x_{1}/(1+\varepsilon_{0})$) and the length after
applying the displacement field is
$\mathrm{d}x_{2}=(1+\mathrm{d}u/\mathrm{d}x_{1})\mathrm{d}x_{1}$. Therefore,
the total strain due to the combination of the two deformations is
$\varepsilon=(\mathrm{d}x_{2}-\mathrm{d}x_{0})/\mathrm{d}x_{0}=\varepsilon_{0}+(1+\varepsilon_{0})\mathrm{d}u/\mathrm{d}x_{1}$.

$\varepsilon_{i}=\bar{\varepsilon}+(1+\bar{\varepsilon})u^{\prime}_{i},~{}i\in\mathcal{I}.$
(5)

Substituting Equation (3) into Equation (5) results in

$\varepsilon_{i}=\bar{\varepsilon}+(1+\bar{\varepsilon})\left(\bar{u}^{\prime}_{i}-y_{i}w^{\prime\prime}\right),~{}i\in\mathcal{I}.$
(6)

The stress field is estimated from strain and, in the case of the TA muscle
layer, TA muscle activation $\mathtt{a}_{\mathrm{ta}}$, which is a non-
dimensional parameter, ranging between 0 and 1, that corresponds to the
activation level in the TA muscle, with 0 indicating a completely flaccid
muscle and 1 being maximum contraction. Herein, we utilize local linearization
about the nominal strain $\bar{\varepsilon}$. That is, the stress functions in
the VF layers, $\sigma_{i},~{}i\in\mathcal{I}$, are given by the approximate
relations

$\sigma_{i}=\sigma_{i,0}+E_{i}(\varepsilon_{i}-\bar{\varepsilon}),~{}i\in\mathcal{I},$
(7)

where

$\sigma_{j,0}=\bar{\sigma}_{j}(\bar{\varepsilon}),~{}E_{j}=\mathrm{d}\bar{\sigma}_{j}(\bar{\varepsilon})/\mathrm{d}\varepsilon,~{}j\in\\{\mathrm{muc},\mathrm{lig}\\},$
$\sigma_{\mathrm{ta},0}=\bar{\sigma}_{\mathrm{ta}}(\bar{\varepsilon},\texttt{a}_{\mathrm{ta}}),~{}E_{\mathrm{ta}}=\frac{\mathrm{d}\bar{\sigma}_{\mathrm{ta}}(\bar{\varepsilon},\texttt{a}_{\mathrm{ta}})}{\mathrm{d}\varepsilon},$

and $\bar{\sigma}_{i},~{}i\in\mathcal{I}$ denotes the nonlinear stress
function associated with the $i^{\mathrm{th}}$ layer.

VF tissues exhibit a highly nonlinear hysteretic viscoelatic behaviour (Min et
al., , 1995; Chan and Titze, , 1999). The literature is rich in various
studies attempting to develop VF constitutive models that capture, at least in
part, the complex mechanical behaviors of the VF tissues (see the review study
of Miri,  (2014)). For example, (Titze and Alipour, , 2006) proposed a one-
dimensional modified Kelvin model for the VF tissues and laryngeal muscles
with nonlinear active and passive stresses to account for tissue
viscoelasticity and muscle activation and implemented this constitutive
modelling framework in simulations of laryngeal postures (Titze and Hunter, ,
2007). Zhang et al.,  (2006) proposed a constitutive model for the VF cover
tissues, which consists of a hyperelastic equilibrium network in parallel with
an inelastic, time-dependent network, and integrated it with an ideal string
model to gain insights into the influence of cover tissue mechanical behaviour
on phonation fundamental frequency. In a study based on measurements collected
from porcine VFs, Miri et al.,  (2013) observed that the collagen fibrils,
which are major constituent of the VFs, are rope-shaped, where the geometric
characteristics of the fibrils have been incorporated in a hyperelastic
mechanical model. In an attempt to capture the anisotropic properties of the
VF lamina propria, Zhang,  (2019) proposed a structurally-based constitutive
model that links the microstructural characteristics of the lamina propria to
its macromechanical properties; the proposed model has shown good agreement
with biaxial tensile testing measurements.

Herein, and for simplicity, we assume the constitutive stress-strain relations
associated with the nonlinear stresses $\bar{\sigma}_{i}$ to be elastic
(functions of strain only) and of exponential type (see Hunter and Titze,
(2007)) with symmetry about zero strain. In particular,

$\bar{\sigma}_{j}=\mathrm{sign}(\varepsilon)m_{j}(\mathrm{e}^{\lvert
n_{j}\varepsilon\rvert}-1),~{}j\in\\{\mathrm{muc},\mathrm{lig}\\},$

and, in the case of the TA muscle, we include stress induced by muscle
activation, resulting in

$\bar{\sigma}_{\mathrm{ta}}=\mathrm{sign}(\varepsilon)m_{\mathrm{ta}}(\mathrm{e}^{\lvert
n_{\mathrm{ta}}\varepsilon\rvert}-1)+\texttt{a}_{\mathrm{ta}}\sigma_{\mathrm{a,max}},$

where $m_{i},n_{i},~{}i\in\mathcal{I},$ are parameters of the constitutive
relations, and $\sigma_{\mathrm{a,max}}$ is the maximum active stress in the
TA muscle. Symmetric stress-strain relations are employed herein to account
for compressive forces developed in the VF, which have been often dismissed in
previous studies of VF biomechanics. The numerical values of the constitutive
relation parameters adopted in this study are listed in Table 1.

Table 1: Numerical values of the geometrical and mechanical properties for each layer in the composite VF model: muc (mucosa), lig (ligament), and ta (thyroarytenoid). Cross-sectional areas are adopted from Titze and Alipour,  (2006), whereas the parameters $m_{i}$, $n_{i}$, and $\sigma_{\mathrm{a,max}}$ are tuned to match experimental stress-strain curves from cadaver and canine models presented in Titze and Alipour,  (2006, Figure 2.17, p. 88). parameter layer | muc | lig | ta
---|---|---|---
$A_{i}$ [mm2] | 5 | 6.1 | 40.9
$m_{i}$ [kPa] | 1.5 | 2 | 1
$n_{i}$[-] | 7 | 10 | 8
$\sigma_{\mathrm{a,max}}$ [kPa] | - | - | 105

The normal forces in the VF layers are computed as
$N_{i}=b\int_{-d_{i}/2}^{d_{i}/2}\sigma_{i}\mathrm{d}y_{i},~{}i\in\mathcal{I}.$
Substituting Equations (7) and (6) in yields

$N_{i}=F_{i,0}+(1+\bar{\varepsilon})E_{i}A_{i}\bar{u}^{\prime}_{i},~{}i\in\mathcal{I},$
(8)

where

$F_{i,0}=A_{i}\sigma_{i,0},~{}i\in\mathcal{I},$ (9)

denote the nominal normal forces generated by each layer. The total internal
normal force is then

$N=\sum_{i\in\mathcal{I}}N_{i}.$ (10)

The moment about the center of the $i^{\mathrm{th}}$ layer due to the stress
developed in that layer is given by
$M_{i}=-b\int_{-d_{i}/2}^{d_{i}/2}y_{i}\sigma_{i}\mathrm{d}y_{i},~{}i\in\mathcal{I}.$
Substituting Equations (7) and (6) into this formula gives

$M_{i}=(1+\bar{\varepsilon})E_{i}I_{i}w^{\prime\prime},~{}i\in\mathcal{I},$
(11)

where
$I_{i}=b\int_{-d_{i}/2}^{-d_{i}/2}y_{i}^{2}\mathrm{d}y_{i},~{}i\in\mathcal{I},$
denotes the area moment of inertia of the $i^{\mathrm{th}}$ layer.

Let the positions of the geometric centers of the VF layers along the $r$-axis
be denoted $r_{i},~{}i\in\mathcal{I}$; that is,

$\begin{split}r_{\mathrm{muc}}&=d_{\mathrm{ta}}+d_{\mathrm{lig}}+\frac{d_{\mathrm{muc}}}{2},\\\
r_{\mathrm{lig}}&=d_{\mathrm{ta}}+\frac{d_{\mathrm{lig}}}{2},\\\
r_{\mathrm{ta}}&=\frac{d_{\mathrm{ta}}}{2},\end{split}$ (12)

see Figure 2. Take a cross-section at longitudinal position $x$ and consider
an arbitrary point on the cross-section located at a vertical position
$r=r_{c}$ (see Figure 2). The moment at $r_{c}$, denoted $M_{c}$, is given by

$\begin{split}M_{c}=&\sum_{i\in\mathcal{I}}M_{i}+(r_{c}-r_{i})N_{i}\\\
=&(1+\bar{\varepsilon})(\sum_{i\in\mathcal{I}}E_{i}I_{i})w^{\prime\prime}\\\
+&(1+\bar{\varepsilon})\sum_{i\in\mathcal{I}}(r_{c}-r_{i})A_{i}E_{i}\bar{u}^{\prime}_{i}\\\
+&\sum_{i\in\mathcal{I}}(r_{c}-r_{i})F_{i,0}.\end{split}$ (13)

Consider an element of infinitesimal longitudinal length $\mathrm{d}x$ with
left edge at position $x$ (see Figure 3), and let $V(x)$ and $q(x)$ denote the
shear force and distributed load per unit length, respectively. The force and
moment balances on the infinitesimal element yield

$\displaystyle N^{\prime}$ $\displaystyle=0,$ (14) $\displaystyle
V^{\prime}-q$ $\displaystyle=0,$ (15) $\displaystyle M^{\prime}+V$
$\displaystyle=0,$ (16)

where second order and higher terms are omitted.

Figure 3: Free-body-diagram of an infinitesimal element of the composite beam
VF model.

From Equation (14) we deduce that the total normal force $N$ is constant
through the VF length. By the assumption that the VF undergoes
compression/elongation with associated strain $\bar{\varepsilon}$ (see
Equation (2)), the force $N$ should be equal to the force that results in that
strain, which is the sum of nominal forces $F_{i,0},~{}i\in\mathcal{I},$ (in
Equation (9)). That is,

$N=\sum_{i\in\mathcal{I}}F_{i,0}.$ (17)

Therefore, by substituting Equations (8) and (10) into Equation (17), we have
$\sum_{i\in\mathcal{I}}(1+\bar{\varepsilon})E_{i}A_{i}\bar{u}^{\prime}_{i}=0,$
implying

$\sum_{i\in\mathcal{I}}E_{i}A_{i}\bar{u}^{\prime}_{i}=0.$ (18)

As we are interested in transverse deflection, we aim to obtain a balance
equation solely in terms of $w$. For convenience, we define

$\begin{split}l_{\mathrm{lig}}&=\frac{1}{2}(d_{\mathrm{lig}}+d_{\mathrm{muc}})\\\
l_{\mathrm{ta}}&=\frac{1}{2}(d_{\mathrm{ta}}+2d_{\mathrm{lig}}+d_{\mathrm{muc}}),\\\
l_{\mathrm{muc}}&=\frac{l_{\mathrm{lig}}E_{\mathrm{lig}}A_{\mathrm{lig}}+l_{\mathrm{ta}}E_{\mathrm{ta}}A_{\mathrm{ta}}}{E_{\mathrm{muc}}A_{\mathrm{muc}}+E_{\mathrm{lig}}A_{\mathrm{lig}}+E_{\mathrm{ta}}A_{\mathrm{ta}}},\end{split}$
(19)

and

$\begin{split}\alpha_{\mathrm{muc}}&=-l_{\mathrm{muc}},\\\
\alpha_{\mathrm{lig}}&=l_{\mathrm{lig}}-l_{\mathrm{muc}},\\\
\alpha_{\mathrm{ta}}&=l_{\mathrm{ta}}-l_{\mathrm{muc}}.\end{split}$ (20)

From the continuity condition in Equation (4) and the zero force condition in
Equation (18), it can be deduced that the displacement functions
$\bar{u}_{i},~{}i\in\mathcal{I},$ satisfy the relations

$\bar{u}^{\prime}_{i}=\alpha_{i}w^{\prime\prime},~{}i\in\mathcal{I}.$ (21)

Substituting Equation (21) into Equation (13), we obtain

$\begin{split}M_{c}&=\mu_{c}w^{\prime\prime}+M_{c,0},\end{split}$ (22)

where

$\mu_{c}=(1+\bar{\varepsilon})\sum_{i\in\mathcal{I}}E_{i}I_{i}+(r_{c}-r_{i})A_{i}E_{i}\alpha_{i}$
(23)

is the composite bending stiffness and

$M_{c,0}=\sum_{i\in\mathcal{I}}(r_{c}-r_{i})F_{i,0}$ (24)

is the nominal moment at $r=r_{c}$ due to the nominal normal forces. Note that
the bending stiffness is strain-dependent, which can be of importance in
posturing scenarios with large VF strains.

Combining Equations (15) and (16) results in

$M^{\prime\prime}_{c}+q=0,$

which, in terms of the deflection $w$ (obtained by substituting in Equation
(22)) is

$\mu_{c}w^{\prime\prime\prime\prime}+q=0.$ (25)

The distributed load $q$ is due to VF contact, which is assumed to be
proportional to the transverse overlap beyond the medial plane. That is,

$q=K_{\mathrm{col}}(w-x\tan(\theta_{G}))\mathbf{H}(w-x\tan(\theta_{G})),$ (26)

where $K_{\mathrm{col}}$ is a stiffness coefficient associated with VF
contact, $\theta_{G}$ is the clockwise angle between the medial plane and the
deformed VF configuration under strain $\bar{\varepsilon}$ (see Figure 2), and
$\mathbf{H}$ is the Heaviside function.

In this work, we assume zero transverse deflection at the anterior and
posterior ends of the VF. That is,

$w(0)=w(L)=0.$ (27)

Moreover, we assume zero moment at the posterior VF margin,

$M_{c}(L)=0.$ (28)

Furthermore, we assume a reactive moment at the anterior VF margin that is
proportional to the rotational displacement with respect to the VF angle at
rest, $\theta_{0}$. The total angle at the anterior margin between the medial
plane and the VF is approximately given by $\theta_{G}-w^{\prime}(0)$.
Consequently, the moment boundary condition at the anterior margin is given by

$M_{c}(0)=-K_{r}(\theta_{G}-w^{\prime}(0)-\theta_{0}),$ (29)

where $K_{r}$ is a rotational stiffness coefficient. Like the nominal strain
$\bar{\varepsilon}$, it is assumed that the angle $\theta_{G}$ is known a
priori. Finally, we assume that $r_{c}$ corresponds to the geometrical center
of the ligament, that is $r_{c}=r_{\mathrm{lig}}$. This assumption, in
addition to the boundary condition given in Equation (28), implies that the
total normal force $N$ is positioned at the geometrical center of the
ligament. This can be deduced from the fact that
$M_{c}(L)=(r_{c}-r_{\mathrm{N}})N$, where $r_{\mathrm{N}}$ denotes the
$r$-position of the total normal force $N$ (that is, the force centroid).

##  3 Analytical insights from a special case

To gain simple yet useful insights into how internal moments inside the VF
beam model affect its curvature, we consider the scenario of zero contact
forces (i.e., $q(x)=0$) and assume $\lvert
w^{\prime}(0)\rvert\ll\lvert\theta_{G}\rvert$, which reduces the boundary
condition given in Equation (29) to

$M_{c}(0)=-K_{r}(\theta_{G}-\theta_{0}).$ (30)

Equation (25), with boundary conditions given in Equations (27), (28), and
(30), and the definition of $M_{c}$ in Equation (22), can be solved
analytically. The curvature of the VF beam model, $w^{\prime\prime}$, is given
explicitly by

$w^{\prime\prime}=-\frac{M_{c,0}}{\mu_{c}}+\left(1-\frac{x}{L}\right)\frac{K_{r}}{\mu_{c}}(\theta_{0}-\theta_{G}).$
(31)

Note that positive $w^{\prime\prime}$ implies a convex VF geometry, whereas
negative curvature implies a concave geometry. Recalling the definition of
$M_{c,0}$ given in Equation (24) and implementing the assumption that
$r_{c}=r_{\mathrm{lig}}$ result in

$\displaystyle-M_{c,0}$
$\displaystyle=-\frac{d_{\mathrm{ta}}+d_{\mathrm{lig}}}{2}F_{\mathrm{ta},0}+\frac{d_{\mathrm{muc}}+d_{\mathrm{lig}}}{2}F_{\mathrm{muc},0}$
$\displaystyle=\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}},$

where

$\tilde{M}_{\mathrm{ta}}=-\frac{d_{\mathrm{ta}}+d_{\mathrm{lig}}}{2}F_{\mathrm{ta},0},$
(32)

and

$\tilde{M}_{\mathrm{muc}}=\frac{d_{\mathrm{muc}}+d_{\mathrm{lig}}}{2}F_{\mathrm{muc},0}.$
(33)

By additionally defining

$\tilde{M}_{r}={K_{r}}(\theta_{0}-\theta_{G}),$ (34)

Equation (31) can be rewritten as

$w^{\prime\prime}=\frac{1}{\mu_{c}}\left[\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}+\left(1-\frac{x}{L}\right)\tilde{M}_{r}\right].$
(35)

In the following discussion, we assume that bending stiffness $\mu_{c}$ is
always positive ($\mu_{c}$ according to Equation (23) changes with the
elongation or compression of the VF beam model). First, let us analyze
abstractly the effects of the moment terms $\tilde{M}_{\mathrm{ta}}$,
$\tilde{M}_{\mathrm{muc}}$, and $\tilde{M}_{r}$ on the VF curvature. We can
observe from Equation (35) that the effect of the reactive moment
$\tilde{M}_{r}$ on the VF curvature decays linearly with a maximum effect (in
magnitude) at $x=0$ and zero effect at $x=L$. In contrast,
$\tilde{M}_{\mathrm{ta}}$ and $\tilde{M}_{\mathrm{muc}}$ have spatially
invariant (i.e., constant) effects on $w^{\prime\prime}$. The curvature is
positively correlated with $\tilde{M}_{r}$ and
$\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}$. That is, $\tilde{M}_{r}>0$
and $\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}>0$ implies positive
curvature (i.e., convex VF geometry), whereas $\tilde{M}_{r}<0$ and
$\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}<0$ implies negative
curvature (i.e., concave VF geometry). Considering the fact that the effect of
the anterior reactive moment $\tilde{M}_{r}$ decays linearly along the VF
length and the nominal moments induced by the VF layers,
$\tilde{M}_{\mathrm{ta}}$ and $\tilde{M}_{\mathrm{muc}}$, are spatially
invariant, there can arise an interesting scenario for which the curvature
changes sign along the VF length. In particular, when

$\tilde{M}_{r}>-(\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}})>0,$ (36)

$w^{\prime\prime}$ is positive on $[0,x_{\mathrm{cr}})$, where

$x_{\mathrm{cr}}=L\left(1+\frac{\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}}{\tilde{M}_{r}}\right),$

and negative for $x\in(x_{\mathrm{cr}},L]$, a change from convexity to
concavity. The conditions on the internal moments and resulting VF shapes from
this analysis are summarized in Table 2.

Table 2: Conditions on the moments applied to the VF composite model based on Equation (35) and the resulting VF shapes: $\rcurvearrowright$ (convex) and $\curvearrowright$ (concave). condition | VF shape
---|---
$\tilde{M}_{\mathrm{muc}}+\tilde{M}_{\mathrm{ta}}$ and $\tilde{M}_{r}>0$ | $\rcurvearrowright$
$\tilde{M}_{\mathrm{muc}}+\tilde{M}_{\mathrm{ta}}$ and $\tilde{M}_{r}<0$ | $\curvearrowright$
$\tilde{M}_{r}>-(\tilde{M}_{\mathrm{muc}}+\tilde{M}_{\mathrm{ta}})>0$ | $\rcurvearrowright\curvearrowright$

We note that a convex VF geometry is a defining characteristic of the bowed VF
pattern. Moreover, the concave VF geometry can be associated with posterior
glottal opening. Furthermore, transition along the VF length from convex to
concave resembles the hourglass glottal pattern (see Figure 1). This
demonstrates that the beam model has the capacity to produce the
experimentally-observed glottal configurations shown in Figure 1.

Now, let us relate the findings listed in Table 2 to physiological posturing
scenarios. The term $\tilde{M}_{r}$, as seen from its definition in Equation
(34), is related to VF adduction and abduction, wherein $\tilde{M}_{r}>0$
corresponds to VF adduction ($\theta_{G}<\theta_{0}$) and $\tilde{M}_{r}<0$
corresponds to VF abduction ($\theta_{G}>\theta_{0}$). The term
$\tilde{M}_{\mathrm{muc}}+\tilde{M}_{\mathrm{ta}}$, which is defined according
to Equations (32) and (33), is determined by the reactive moments developed in
the VF layers, especially the mucosa and TA muscle, during VF tensioning or
compression. Note that, based on the area measurements (Table 1) and the
assumption of uniform thickness $b$, the moment arm of the TA muscle,
$(d_{\mathrm{ta}}+d_{\mathrm{lig}})/2$, is larger than that of the mucosa,
$(d_{\mathrm{muc}}+d_{\mathrm{lig}})/2$. In the case of VF compression, the
compressive forces in the TA muscle are typically larger in magnitude than
that in the mucosa (i.e., $F_{\mathrm{ta},0}\ll F_{\mathrm{muc},0}<0$); hence,
the moment induced by the TA muscle, $\tilde{M}_{\mathrm{ta}}$, is positive
and predominant making $\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}>0$.
On the other hand, when the VF is tensioned due to activating the TA muscle,
the force $F_{\mathrm{ta},0}$ is positive and predominant and, consequently,
the term $\tilde{M}_{\mathrm{ta}}$ is negative (see Equation (32)) and
predominant. This scenario results in
$\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}<0$. These relations between
the moment terms and corresponding laryngeal posturing scenarios are
summarized in Table 3.

Table 3: Conditions on the moments applied to the VF composite beam model and corresponding physiological posturing scenarios  condition | physiological scenario
---|---
$\tilde{M}_{r}>0$ | VF adduction
$\tilde{M}_{r}<0$ | VF abduction
$\tilde{M}_{\mathrm{muc}}+\tilde{M}_{\mathrm{ta}}>0$ | VF compression
$\tilde{M}_{\mathrm{muc}}+\tilde{M}_{\mathrm{ta}}<0$ | TA muscle activation

The combined findings presented in Tables 2 and 3 can be summarized by
following observations: The bowed shape with convex VF geometry can be due to
(a) positive reactive moment at the anterior margin ($\tilde{M}_{r}>0$) during
VF adduction, and/or (b) internal moments during VF compression, wherein
$\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}>0$. Moreover, the concave VF
shape arising in the case of posterior glottal opening can be due to (a)
negative reactive moment at the anterior margin ($\tilde{M}_{r}<0$) during VF
abduction, and/or (b) sufficiently large activation of the TA muscle, wherein
$\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}<0$. The hourglass shape may
necessitate a coordinated laryngeal maneuver that involves sufficient TA
activation ($\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}<0$) and VF
adduction ($\tilde{M}_{r}>0$) such that
$\tilde{M}_{r}+\tilde{M}_{\mathrm{ta}}+\tilde{M}_{\mathrm{muc}}>0$.

##  4 Simulations of the combined beam and posturing model

In this section we further investigate VF curvature and incomplete glottal
closure by combining our beam model with the VF posturing modeling introduced
by Titze and Hunter,  (2007). In particular, we adopt the implementation of
Alzamendi et al.,  (2022). The posture model relates activation of the five
intrinsic muscles to the prephonatory glottal configuration, and in
particular, the rotational and linear displacements of the cricothyroid joints
and the arytenoid cartilages, where the VFs and intrinsic muscles are modelled
as spring-like elements. From the aforementioned displacements, the VF nominal
strain $\bar{\varepsilon}$ and glottal angle $\theta_{G}$ are estimated.
Similar to the muscle activation parameter $\mathrm{a}_{\mathrm{ta}}$ embedded
in the VF beam model, the posture model relies on five normalized muscle
activation parameters, $a_{\mathrm{ta}}$, $a_{\mathrm{ct}}$,
$a_{\mathrm{lca}}$, $a_{\mathrm{ia}}$, and $a_{\mathrm{pca}}$, which
correspond to the TA, cricothyroid (CT), lateral cricoarytenoid (LCA),
interarytenoid (IA), and PCA muscles, respectively. In this study, we assume
that the muscle activation parameter $\mathrm{a}_{\mathrm{ta}}$ embedded in
the VF beam model is identical in value to the muscle activation parameter
$a_{\mathrm{ta}}$ in the posturing model
($\mathrm{a}_{\mathrm{ta}}=a_{\mathrm{ta}}$).

It is important to mention that the constitutive relations embedded in the VF
beam model are different from those in the posture model implementation
adopted from Alzamendi et al.,  (2022). The focus of the current study is to
replicate the VF static configurations, wherein we employ experimental stress-
strain data based on human and canine samples (Titze and Alipour, , 2006,
Figure 2.17, p. 88) to prescribe the mechanical behaviors of the tissues. The
posture model of Alzamendi et al.,  (2022) instead focuses on replicating
physiologically accurate posturing and phonation outputs; this required ad hoc
tuning of some posturing model parameters. Prior experimental and numerical
studies typically suffer from significant variability in the reported
numerical values of biomechanical parameters (see, e.g., Titze and Alipour,
(2006); Titze and Hunter,  (2007); Hunter and Titze,  (2007); Palaparthi et
al.,  (2019)), and in some cases numerical values are missing altogether,
which motivates the ad hoc tuning approach adopted by Alzamendi et al.,
(2020).

The posturing model in Titze and Hunter,  (2007); Alzamendi et al.,  (2022) is
dynamic due to inertial and viscous effects. In this study, and as we are
interested in static posturing scenarios, the posture model is run until the
VF strain and glottal angle reach steady-state and these values are input into
the composite beam model. Once the VF strain and glottal angle parameters are
fed into the beam model, Equation (25), supplemented with Equations (13),
(26), (27), (28), and (29), is solved numerically. The aforementioned
equations and boundary conditions are discretized by means of finite
difference. For the simulations $\theta_{0}$ is set as the glottal angle from
the posturing simulations when all laryngeal muscles are inactive. Numerical
values for the remaining VF beam model parameters are listed in Table 4.

Table 4: Numerical values of some of the VF beam model parameters. parameter | numerical value
---|---
$L_{0}$ [mm] | 15
$b$ [mm] | 5
$\theta_{0}$ [rad] | 0.2540
$K_{\mathrm{col}}$ [N/m2] | $2\times 10^{8}$
$K_{r}$ [N.m] | 0.05

To clearly illustrate the glottal geometries resulting from the simulations, a
coordinate system $(x_{1},x_{2})$ with origin at the anterior margin of the
VFs is established. The $x_{1}$-axis is aligned along the medial plane
pointing in the posterior direction and the $x_{2}$-axis is perpendicular to
the medial plane pointing to the right, relative to the human body frame (see
Figure 2). In all figures presented in this section the VF configurations are
plotted with respect to this coordinate system; model symmetry is utilized to
produce the opposing VF shape.

This section explores several laryngeal maneuvers and how they influence the
VF geometry666The end points (anterior and posterior margins) of the VFs
resulting from the proposed beam model are identical to those established by
the posturing model of Titze and Hunter,  (2007). Accounting for internal
bending moments results in deviation of the VF shape from the linear medial
surface prescription (with angle $\theta_{G}$) of Titze and Hunter,  (2007)..
In particular, and motivated by previous clinical, experimental, and numerical
findings (Yin and Zhang, , 2014; Morrison and Rammage, , 1993; Chhetri and
Neubauer, , 2015; Yin and Zhang, , 2016), we consider laryngeal maneuvers
associated with adductory (TA, LCA, and IA) and abductory (PCA), muscles as
they have been found to play major roles in inducing curved glottal
geometries. The CT muscle has been found to play a major role in regulating
phonation fundamental frequency by stretching the VFs, but not in posturing
and is thus excluded from this study. Herein, we compare simulation results
with findings from previous clinical, experimental, and high-fidelity
numerical studies to verify the proposed VF beam model. Moreover, we attempt
to elucidate potential mechanisms underlying the curved VF geometries observed
clinically by analyzing the beam model details (see Figure 1).

First, we investigate the effects of increasing co-activation of the LCA and
IA muscles, which are responsible for adducting the VFs (Alzamendi et al., ,
2022), while the remaining intrinsic muscles are inactive. Figure 4 presents
the glottal shapes and the induced moments corresponding to simulations
wherein LCA and IA muscle activation levels are increased simultaneously.
Figure 4(left) shows that co-activation of the LCA and IA muscles leads to
posterior glottal closure with a remaining mid-membranous gap. This convex VF
shape matches previous clinical and numerical findings (Chhetri and Neubauer,
, 2015; Yin and Zhang, , 2016). Figure 4(right) shows that the VF convexity is
due to the predominance of the reactive moments at the anterior VF margin
($\tilde{M}_{r}$ is positive and relatively large), which arises due to VF
adduction $(\theta_{G}<\theta_{0})$, which agrees with the theoretical
predictions in Section 3.

Figure 4:  (Left) glottal profile and (right) induced moments for increasing
LCA and IA activation levels and other intrinsic muscles being inactive.

Figure 5 exhibits the glottal shapes and induced moments corresponding to
simulations wherein TA muscle activation levels are increased, while all other
intrinsic muscles are inactive. Figure 5(left) shows that isolated activation
of the TA muscle leads to anterior and mid-membranous glottal closure with
remaining posterior opening, while also shortening the folds. The resulting
concave VF shapes are in agreement with previous experimental and numerical
investigations (Chhetri and Neubauer, , 2015; Yin and Zhang, , 2016). Figure
5(right) shows that the concavity is primarily determined by the internal
moments induced by the TA muscle activation ($\tilde{M}_{\mathrm{ta}}$ is
negative and relatively large in magnitude), which is in alignment with the
analysis in Section 3.

Figure 5: (Left) glottal profile and (right) induced moments for increasing TA
activation levels and other intrinsic muscles being inactive.

In an effort to explore the mechanics of the hourglass glottal shape, we
explore the glottal shape associated with increasing activation of the TA
muscle while the LCA and IA are kept at constant non-zero activation levels.
Simulating such maneuvers is encouraged by the findings from the theoretical
analysis in Section 3 and the experimental observations in Chhetri et al.,
(2012). Figure 6 displays the glottal shapes and induced moments associated
with slight increasing activation of the TA muscle, while the LCA and IA are
kept at constant levels ($a_{\mathrm{lca}}=a_{\mathrm{ta}}=0.6$). Figure
6(left) shows that in the case of zero TA activation, the glottal shape is
bowed with slight, but not full, posterior adduction. As TA activation is
increased, a medial bulge is observed whereas anteriorly the glottal geometry
is still convex, resulting in an overall hourglass shape. Figure 6(right)
displays how the internal moments $\tilde{M}_{\mathrm{ta}}$ and
$\tilde{M}_{\mathrm{muc}}$ and the reactive moment at the anterior margin
$\tilde{M}_{r}$ satisfy the condition of Equation (36). This aligns with the
analysis in Section 3 and suggests that the hourglass glottal shape
necessitates involvement of reactive moments at the anterior VF margin
(associated with VF adduction) and internal moments induced inside the VF
layers (primarily the TA muscle). In addition, this finding is in good
agreement with observations in Chhetri et al.,  (2012), which showed that an
hourglass shape is induced by coactivating all the adductory muscles.

Figure 6:  (Left) glottal shapes, and (right) induced moments for increasing
TA activation levels and other intrinsic muscle activation levels being at
$a_{\mathrm{lca}}=a_{\mathrm{ta}}=0.6$ and
$a_{\mathrm{ct}}=a_{\mathrm{pca}}=0$.

In aggregate, Figures 4-6 indicate that dismissing the TA muscle or the LCA
and IA muscles cannot produce full glottal closure, therefore, in the next set
of simulations, we investigate the effects of co-activating all of the
adductory muscles. Figure 7 shows that (almost) full glottal closure can be
attained when all adductors are co-activated
($a_{\mathrm{ia}}=a_{\mathrm{lca}}=0.45,~{}a_{\mathrm{ta}}=0.7$), which is in
alignment with previous experimental and numerical investigations (Chhetri and
Neubauer, , 2015; Yin and Zhang, , 2016).

Figure 7: Glottal shapes for increasing co-activation levels of all adductory
muscles, and the PCA and CT muscles being inactive.

Finally, we explore the effects of the PCA muscle, a primary VF abductor, on
the glottal geometry, where we consider simulations motivated by the clinical
observations highlighted in Morrison and Rammage,  (1993). Figure 8 displays
glottal patterns associated with increasing PCA activation where adductory
muscles are kept at activation levels associated with near full closure. The
figure displays that increasing PCA activation leads to posterior opening,
while the VFs are sustaining concave shapes similar to those presented when TA
alone is activated (see Figure 5). This suggests that PCA activation tends to
neutralize the posterior adductory effects of the LCA and IA muscles, which
supports the clinical observations highlighted in Morrison and Rammage,
(1993).

Figure 8: Glottal shapes for increasing activation of the PCA muscle with the
activation levels of the other muscles being
$(a_{\mathrm{lca}},a_{\mathrm{ia}},a_{\mathrm{ta}})=(0.45,0.45,0.7)$ and
$a_{\mathrm{ct}}=0$.

##  5 Discussion

The results of Sections 3 and 4 highlight potential mechanisms underlying
different patterns of incomplete glottal closure. In particular, results
indicate that bowed VF shapes result, in part, from low or null activation of
the TA muscle in combination with co-activation of the LCA and IA muscles. The
predominant mechanism in this case is the anterior reactive moments that
resists bringing the VFs together during adduction. This pattern can also
arise in the case of low TA muscle activation and VF compression, as suggested
by the analysis in Section 3. In this case, the internal moment induced by the
TA muscle tissue ($\tilde{M}_{\mathrm{ta}}$ is positive and predominantly
large) is the driving factor. This scenario (bowing due to VF compression) can
potentially take place when extrinsic laryngeal muscles are excessively
activated, especially those associated with VF compressing, such as the
thyrohyoid muscle (Hong et al., , 1997).

In addition, our analysis suggests that posterior glottal opening with
combined VF concavity results from high activation of the TA muscle and low or
null activation of the LCA and IA muscles. Our model suggests that the driving
mechanism here is the internal moment induced by the TA muscle activation
$\tilde{M}_{\mathrm{ta}}$, which is negative and predominantly large in
magnitude in this case. A similar glottal pattern also occurs when all
adductory muscles are activated in addition to the activation of the PCA
muscle. This supports the hypothesis of Morrison and Rammage,  (1993),
regarding the excessive activation of the PCA muscle in patients with MTD.
Finally, our analysis suggests that the hourglass glottal shape may emerge
from laryngeal maneuvers that involve, for example, moderate co-activation of
all adductory muscles, where both anterior reactive moment and internal moment
due to TA muscle activation are at play and opposing each other.

The implications above concerning potential connections between incomplete and
curved glottal closure patterns and particular muscular executions may help
speech therapists to uncover the underlying laryngeal mechanisms associated
with some voice disorders. As highlighted in the introduction, incomplete
glottal closure can be linked to voice disorders that are characterized by
excessive, imbalanced, or deficient activity of the intrinsic and extrinsic
muscles such as MTD and Parkinson’s disease. Our analysis in the current work
suggests two potential mechanisms underlying bowed VFs in some patients with
voice disorders (1) the TA muscle is not properly activated (possibly due to
muscle activation imbalance), and (2) excessive activation of extrinsic neck
muscles, leading to VF compression. Besides, our analysis posits that in
patients with abnormal posterior glottal opening and concave VF geometry
either (1) insufficiently activate the LCA and IA muscles, whereas the TA
muscle is activated sufficiently (in comparison to normal posturing
scenarios), or (2) suffer from excessive activation of all adductory and
abductory muscles, where the PCA muscle activation mitigates the effects of
the LCA and IA muscles, in agreement with the postulation in Morrison and
Rammage,  (1993) concerning patients with MTD. In summary, speech clinicians
and therapists may consider the aforementioned candidate underlying mechanisms
of curved and incomplete glottal closure patterns when examining patients with
voice disorders involving muscular inefficiencies/deficiencies.

A number of simplifying assumptions are embedded in the presented model of
this study, including (1) negligible shear deformation, (2) negligible elastic
forces from the connective tissues attached to the TA muscle, (3) negligible
motions in the superior-inferior direction, (4) neglecting potential bending
effects from the vocal ligament by setting $r_{c}=r_{\mathrm{lig}}$, (5) zero
moments at the posterior ends of the VFs, (6) small transverse VF deflections,
and (7) one-way coupling between the VF beam model and the posturing model,
where any contact forces emerging due to the VF curvature do not alter the
mechanics of the laryngeal cartilages.

Assumption (1) is a consequence of the adopted Euler-Bernoulli framework. Note
that with the uniform thickness assumption, and considering the model
dimensions given in Tables 1 and 4, the total VF depth,
$\sum_{i\in\mathcal{I}}d_{i}$, is approximately 10 mm whereas the resting VF
length is 15 mm; hence, the depth and length dimensions are quite comparable.
For such cases (thick beams), Timoshenko beam theory (Timoshenko, , 1921) is
typically adopted to account for shear stresses777 It is worth noting that the
two theories (Euler-Bernoulli and Timoshenko) do coincide for a uniform
homogeneous simply-supported linear beam with specified moments at the end
points and zero distributive load, regardless of the beam thickness or
mechanical properties (in Section 3, we studied a similar simply-supported
case with zero distributive load). These two theories, when compared, tend to
produce qualitatively, but not necessarily quantitatively, similar predictions
(see, e.g., Beck and da Silva Jr,  (2011)).. As the goal of the current work
is to construct a simple analytically-tractable model that predicts
qualitatively the curved glottal configurations observed clinically, we
adhered to the Euler-Bernoulli framework, leaving derivations of more complex
models to future work. We posit that assumption (2) is reasonable as the
elastic forces from the connective tissues are passive, mostly only
restricting the extent to which the VF deflects. Moreover, assumption (3) is
suitable as the majority of the VF motion during posturing occurs medially
and/or laterally (see, e.g., the findings of Chhetri and Neubauer,  (2015)).
Regarding assumption (4), the ligament is stiffer than other VF layers and it
geometrically forms the intermediate VF layer, making it the ‘chassis’ of the
VF layered structure; hence setting $r_{c}=r_{\mathrm{lig}}$, which indicates
that the total normal force in the VF is positioned at the center of the
ligament (see the end of Section 2), is sensible. Assumptions (5)-(7), in
addition to other assumptions such as the rectangular geometries of the VF
layers, are introduced to primarily simplify our analysis; hence, further
investigation is needed to verify the validity of such assumptions in
different posturing scenarios, and refine them when needed.

Despite these simplifying assumptions, our modelling framework is capable of
predicting some of the glottal patterns observed in previous clinical and
high-fidelity numerical studies, which is encouraging. Still, the speculations
and potential explanations provided in this work need further extensive
investigation into the biomechanics of VF posturing in both healthy subjects
and patients with imbalances or deficiencies in the laryngeal muscles.

##  6 Conclusion

In this study, we introduced a simple one-dimensional Euler-Bernouilli-type
composite beam model of the vocal folds to understand the mechanisms
underlying glottal configuration and incomplete glottal closure. The model,
despite its simplicity, was capable of predicting several clinically observed
glottal configurations. Our analysis highlighted how the different patterns of
incomplete glottal closure can arise naturally due to the layered VF structure
and the associated induced moments. We coupled the proposed beam model with
the posturing model of Titze and Hunter,  (2007) to gain physiologically
relevant insights into the role of laryngeal muscle activation. Our analysis
showed that a bowed VF shape can arise due to activation of the LCA and IA
muscles without incorporating the TA muscle during adduction or due to VF
compression. On the other hand, isolated activation of the TA muscle results
in medial bulging and posterior glottal opening. Posterior opening can also
occur due to activating all adductors in addition to activating the PCA
muscle. Moreover, our analysis suggested that an hourglass glottal shape can
arise from specific laryngeal maneuvers involving the adductory laryngeal
muscles. These results provided potential explanations and conjectures
regarding the posturing mechanics of patients with voice disorders such as
MTD.

In future efforts we aim to refine our modelling framework, where two-way
coupling between the VF beam model presented herein and the posturing model of
Titze and Hunter,  (2007) is incorporated, to account for potential effects
that curved VF geometries may exert on the mechanics of laryngeal cartilages.
Moreover, we intend to incorporate the beam model with numerical phonation
models, to study how curved and partially closed glottal geometries may
influence tissue-flow-acoustic interactions, voice quality, and vocal function
during phonation.

## Acknowledgments

Research reported in this work was supported by the NIDCD of the NIH under
Award No. P50DC015446, and ANID BASAL FB0008. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the National Institutes of Health.

## Declaration of competing interest

Matías Zañartu has a financial interest in Lanek SPA, a company focused on
developing and commercializing biomedical devices and technologies. His
interests were reviewed and are managed by the Universidad Técnica Federico
Santa María in accordance with its conflict-of-interest policies.

## Author contributions

Mohamed Serry: Conceptualization, Methodology, Formal analysis, Software,
Visualization, Writing - original draft preparation. Gabriel Alzamendi:
Software, Validation, Writing - reviewing and editing. Matías Zañartu:
Validation, Writing - reviewing and editing, funding acquisition. Sean
Peterson: Supervision, Validation, Writing - reviewing and editing, funding
acquisition.

## References

  * Altman et al., (2005) Altman, K. W., Atkinson, C., and Lazarus, C. (2005).  Current and emerging concepts in muscle tension dysphonia: a 30-month review.  Journal of Voice, 19(2):261–267.
  * Alzamendi et al., (2020) Alzamendi, G. A., Manríquez, R., Hadwin, P. J., Deng, J. J., Peterson, S. D., Erath, B. D., Mehta, D. D., Hillman, R. E., and Zañartu, M. (2020).  Bayesian estimation of vocal function measures using laryngeal high-speed videoendoscopy and glottal airflow estimates: An in vivo case study.  The Journal of the Acoustical Society of America, 147(5):EL434–EL439.
  * Alzamendi et al., (2022) Alzamendi, G. A., Peterson, S. D., Erath, B. D., Hillman, R. E., and Zañartu, M. (2022).  Triangular body-cover model of the vocal folds with coordinated activation of the five intrinsic laryngeal muscles.  The Journal of the Acoustical Society of America, 151(1):17–30.
  * Bauchau and Craig, (2009) Bauchau, O. A. and Craig, J. I. (2009).  Structural analysis: with applications to aerospace structures, volume 163.  Springer Science & Business Media.
  * Beck and da Silva Jr, (2011) Beck, A. T. and da Silva Jr, C. R. (2011).  Timoshenko versus euler beam theory: Pitfalls of a deterministic approach.  Structural Safety, 33(1):19–25.
  * Braak and Braak, (2000) Braak, H. and Braak, E. (2000).  Pathoanatomy of parkinson’s disease.  Journal of neurology, 247(2):II3–II10.
  * Chan and Titze, (1999) Chan, R. W. and Titze, I. R. (1999).  Viscoelastic shear properties of human vocal fold mucosa: Measurement methodology and empirical results.  The Journal of the Acoustical Society of America, 106(4):2008–2021.
  * Chhetri and Neubauer, (2015) Chhetri, D. K. and Neubauer, J. (2015).  Differential roles for the thyroarytenoid and lateral cricoarytenoid muscles in phonation.  The Laryngoscope, 125(12):2772–2777.
  * Chhetri et al., (2012) Chhetri, D. K., Neubauer, J., and Berry, D. A. (2012).  Neuromuscular control of fundamental frequency and glottal posture at phonation onset.  The Journal of the Acoustical Society of America, 131(2):1401–1412.
  * (10) Choi, H.-S., Berke, G. S., Ye, M., and Kreiman, J. (1993a).  Function of the posterior cricoarytenoid muscle in phonation: in vivo laryngeal model.  Otolaryngology—Head and Neck Surgery, 109(6):1043–1051.
  * (11) Choi, H.-S., Ye, M., Berke, G. S., and Kreiman, J. (1993b).  Function of the thyroarytenoid muscle in a canine laryngeal model.  Annals of Otology, Rhinology & Laryngology, 102(10):769–776.
  * Dejonckere and Kob, (2009) Dejonckere, P. H. and Kob, M. (2009).  Pathogenesis of vocal fold nodules: new insights from a modelling approach.  Folia Phoniatrica et Logopaedica, 61(3):171–179.
  * Geng et al., (2020) Geng, B., Pham, N., Xue, Q., and Zheng, X. (2020).  A three-dimensional vocal fold posturing model based on muscle mechanics and magnetic resonance imaging of a canine larynx.  The Journal of the Acoustical Society of America, 147(4):2597–2608.
  * Hanson et al., (1984) Hanson, D. G., Gerratt, B. R., and Ward, P. H. (1984).  Cinegraphic observations of laryngeal function in parkinson’s disease.  The Laryngoscope, 94(3):348–353.
  * Hillman et al., (2020) Hillman, R. E., Stepp, C. E., Van Stan, J. H., Zañartu, M., and Mehta, D. D. (2020).  An updated theoretical framework for vocal hyperfunction.  American Journal of Speech-Language Pathology, pages 1–7.
  * Hocevar-Boltezar et al., (1998) Hocevar-Boltezar, I., Janko, M., and Zargi, M. (1998).  Role of surface emg in diagnostics and treatment of muscle tension dysphonia.  Acta oto-laryngologica, 118(5):739–743.
  * Hong et al., (1997) Hong, K. H., Ye, M., Kim, Y. M., Kevorkian, K. F., and Berke, G. S. (1997).  The role of strap muscles in phonation: in vivo canine laryngeal model.  Journal of Voice, 32.
  * Hunter and Titze, (2007) Hunter, E. J. and Titze, I. R. (2007).  Refinements in modeling the passive properties of laryngeal soft tissue.  Journal of Applied Physiology, 103(1):206–219.
  * Hunter et al., (2004) Hunter, E. J., Titze, I. R., and Alipour, F. (2004).  A three-dimensional model of vocal fold abduction/adduction.  The Journal of the Acoustical Society of America, 115(4):1747–1759.
  * Min et al., (1995) Min, Y. B., Titze, I. R., and Alipour-Haghighi, F. (1995).  Stress-strain response of the human vocal ligament.  Annals of Otology, Rhinology & Laryngology, 104(7):563–569.
  * Miri, (2014) Miri, A. K. (2014).  Mechanical characterization of vocal fold tissue: a review study.  Journal of Voice, 28(6):657–667.
  * Miri et al., (2013) Miri, A. K., Heris, H. K., Tripathy, U., Wiseman, P. W., and Mongeau, L. (2013).  Microstructural characterization of vocal folds toward a strain-energy model of collagen remodeling.  Acta biomaterialia, 9(8):7957–7967.
  * Morrison and Rammage, (1993) Morrison, M. D. and Rammage, L. A. (1993).  Muscle misuse voice disorders: description and classification.  Acta oto-laryngologica, 113(3):428–434.
  * Nguyen et al., (2009) Nguyen, D. D., Kenny, D. T., Tran, N. D., and Livesey, J. R. (2009).  Muscle tension dysphonia in vietnamese female teachers.  Journal of Voice, 23(2):195–208.
  * Palaparthi et al., (2019) Palaparthi, A., Smith, S., and Titze, I. R. (2019).  Mapping thyroarytenoid and cricothyroid activations to postural and acoustic features in a fiber-gel model of the vocal folds.  Applied Sciences, 9(21):4671.
  * Pillutla et al., (2022) Pillutla, P., Reddy, N. K., Schlegel, P., Zhang, Z., and Chhetri, D. K. (2022).  Control of pre-phonatory glottal shape by intrinsic laryngeal muscles.  The Laryngoscope.
  * Rajaei et al., (2014) Rajaei, A., Barzegar Bafrooei, E., Mojiri, F., and Nilforoush, M. H. (2014).  The occurrence of laryngeal penetration and aspiration in patients with glottal closure insufficiency.  International Scholarly Research Notices, 2014.
  * Roy, (2008) Roy, N. (2008).  Assessment and treatment of musculoskeletal tension in hyperfunctional voice disorders.  International Journal of Speech-Language Pathology, 10(4):195–209.
  * Södersten et al., (1995) Södersten, M., Hertegård, S., and Hammarberg, B. (1995).  Glottal closure, transglottal airflow, and voice quality in healthy middle-aged women.  Journal of Voice, 9(2):182–197.
  * Timoshenko, (1921) Timoshenko, P. S. (1921).  Lxvi. on the correction for shear of the differential equation for transverse vibrations of prismatic bars.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 41(245):744–746.
  * Titze and Alipour, (2006) Titze, I. and Alipour, F. (2006).  The Myoelastic-Aerodynamic Theory of Phonation.
  * Titze and Hunter, (2004) Titze, I. R. and Hunter, E. J. (2004).  Normal vibration frequencies of the vocal ligament.  The Journal of the Acoustical Society of America, 115(5):2264–2269.
  * Titze and Hunter, (2007) Titze, I. R. and Hunter, E. J. (2007).  A two-dimensional biomechanical model of vocal fold posturing.  The Journal of the Acoustical Society of America, 121(4):2254–2260.
  * Yin and Zhang, (2014) Yin, J. and Zhang, Z. (2014).  Interaction between the thyroarytenoid and lateral cricoarytenoid muscles in the control of vocal fold adduction and eigenfrequencies.  Journal of biomechanical engineering, 136(11):111006.
  * Yin and Zhang, (2016) Yin, J. and Zhang, Z. (2016).  Laryngeal muscular control of vocal fold posturing: Numerical modeling and experimental validation.  The Journal of the Acoustical Society of America, 140(3):EL280–EL284.
  * Zañartu et al., (2014) Zañartu, M., Galindo, G. E., Erath, B. D., Peterson, S. D., Wodicka, G. R., and Hillman, R. E. (2014).  Modeling the effects of a posterior glottal opening on vocal fold dynamics with implications for vocal hyperfunction.  The Journal of the Acoustical Society of America, 136(6):3262–3271.
  * Zhang et al., (2006) Zhang, K., Siegmund, T., and Chan, R. W. (2006).  A constitutive model of the human vocal fold cover for fundamental frequency regulation.  The Journal of the Acoustical Society of America, 119(2):1050–1062.
  * Zhang et al., (2007) Zhang, K., Siegmund, T., and Chan, R. W. (2007).  A two-layer composite model of the vocal fold lamina propria for fundamental frequency regulation.  The Journal of the Acoustical Society of America, 122(2):1090–1101.
  * Zhang, (2019) Zhang, Z. (2019).  Structural constitutive modeling of the anisotropic mechanical properties of human vocal fold lamina propria.  The Journal of the Acoustical Society of America, 145(6):EL476–EL482.