A fibril-scale visco-hyperelastic model for the mechanics of vocal-fold tissues

Introduction: Modeling the mechanics of human vocal folds during phonation is a challenging task. This is partly due to the mechanics of their soft and highly anisotropic fibrous tissues, which undergoes finite strains with both elasticity and strain-rate sensitivity.

Methods: We propose a visco-hyperelastic micro-mechanical model capable of predicting the complex cyclic response of the vocal-fold fibrous tissues based on their histo-mechanical properties. For that purpose, we start from the hyperelastic micro-mechanical model proposed by Terzolo et al., J. Mech. Behav. Biomed. Mater. 128:105118 (2022). We include in the model non-linear viscoelastic contributions at the fibril scale to account for the dissipative and time-dependent response of vocal-fold tissues.

Results and Discussion: The relevance of the model is demonstrated and discussed through comparison with a comprehensive set of reference experimental data, within a wide range of loading modes, strains, and strain rates: cyclic and multi-axial loadings at finite strains (tension, compression and shear), along with small-amplitude oscillatory shear (SAOS) and large-amplitude oscillatory shear (LAOS) from low to high frequencies. This study elucidates how the viscoelasticity of vocal-fold tissues can result from combined time-dependent micro-mechanisms, such as the kinematics and the deformation of their fibril bundles, along with the mechanical interactions likely to develop among fibrils and the surrounding amorphous matrix.

1 Introduction

Human vocal folds are soft laryngeal structures with remarkable mechanical properties. During phonation, the vocal folds deform under the action of pulmonary airflow and laryngeal motions, sustaining vibrations in a wide range of amplitudes, frequencies (from 50 Hz to more than 1,500 Hz), and degrees of collisions. These multiple configurations involve complex and coupled multi-axial mechanical stresses (in tension, compression, and shear) that the tissues can withstand upon finite strains at various strain rates (Miri, 2014; Vampola et al., 2016). These properties are inherited from the composite and hierarchical structure of the vocal folds and surrounding laryngeal muscles. More specifically, the vocal folds are made up of two main load-bearing layers: the lamina propria, i.e., a loose connective tissue, and the vocalis muscle. Both layers are composed of networks of collagen, elastin, or skeletal muscle microfibrils, embedded in a soft hydrogel-like matrix (Figure 1, Hirano, 1974, Benboujja and Hartnick, 2021, Ferri-Angulo et al., 2023). However, to date, our knowledge is still insufficient to understand the relationship between the fibril-scale architecture of the vocal folds and their macroscale (tissue-scale) time-dependent performance.

Figure 1

Figure 1. Idealization of the vocal-fold layers. (a) The lamina propria (respectively the vocalis), displayed on microscopic image 1 (respectively 2), is considered a network of (orange) collagen fibrils [respectively (pink) myofibrils and (orange) collagen fibrils] embedded into a gel-like matrix. Fibrils are self-assembled as collagen fibril bundles (respectively myofibrils surrounded by a sheath of collagen fibrils). Each fibril (and its interaction with its neighboring) behaves as a non-linear visco-hyperelastic Zener model. (b) The fiber bundle microstructure of each layer is considered a periodic network of four orientated fiber bundles (brown) connected at one node ${\text{C}}_0$ (blue) embedded in a soft isotropic matrix (green). The dotted lines illustrate the five possible steric interactions of ${\text{C}}_0$ with the neighboring nodes. Source: Adapted from Terzolo et al. (2022).

This is mainly due to the difficulty of characterizing vocal-fold mechanics at high physiological strain rates. Although recent progress has been made in time-resolved 3D micro-imaging of fast-vibrating structures (Klos et al., 2024), to date, the characterization of the mechanical behavior of vocal-fold tissues at high frequencies (e.g., from 100 Hz to 1 kHz) is still limited to the larynx or vocal-fold scale. High-speed videostroboscopy, used in clinical voice assessment, enabled the quantification of the time-decay of vocal-fold vibrations at phonation offset (DeJonckere and Lebacq, 2020; Radolf et al., 2022), and of their resonance properties obtained through external excitation of the larynx (Švec et al., 2000). Such in vivo approaches allowed the measurement of an average damping ratio $\zeta \approx$0.07–0.20, describing the dissipation of stored energy in oscillations for frequencies between 100 and 200 Hz (Švec et al., 2000; DeJonckere and Lebacq, 2020; Radolf et al., 2022), which partly arises from the viscoelastic behavior of the tissues. The time-dependent mechanical properties of vocal-fold tissues have also been demonstrated ex vivo through numerous phenomena, including strain-rate sensitivity of stress–strain behavior, creep, stress relaxation, stress hysteresis, and related accommodation upon cycling, with the magnitude of the hysteresis loop dependent on the strain rate (Kelleher et al., 2013a; Chan and Titze, 1999; 2000; Chan, 2004; Klemuk and Titze, 2004; Titze et al., 2004; Chan and Rodriguez, 2008; Miri et al., 2014; Chan, 2018; Cochereau et al., 2020). The viscoelastic properties of excised lamina propria samples were mostly studied using standard shear dynamic mechanical analysis (DMA), also called small-amplitude oscillatory shear (SAOS), i.e., within the linear regime (Chan and Titze, 1999; Chan and Titze, 2000; Chan, 2004; Klemuk and Titze, 2004; Titze et al., 2004; Chan and Rodriguez, 2008). Such works allowed to characterize the shear storage $G^{\prime }$ and loss $G^{″}$ moduli of the vocal-fold “cover” (i.e., the superficial sublayer of the lamina propria combined with the epithelium that covers it) for excitation frequencies $f$ up to 250 Hz. Therefore, these dynamic moduli increase (respectively decrease) with the applied frequency (respectively strain), while the loss factor [$\tan \delta ={\zeta }^{\star }=G^{″}/G^{\prime }$ (Dashatan et al., 2023; Koruk and Rajagopal, 2024)] decreases monotonically with frequency, down to a mean value of 0.73 for $f$ within 100–250 Hz (Chan and Rodriguez, 2008). Such experiments were recently extended to large-amplitude oscillatory shear (LAOS), showing that lamina propria sublayers experience intercycle strain softening during oscillatory strain sweeps, intracycle strain stiffening, shear thinning while increasing the shear rate, and complex stress hysteresis that depends on the applied strain and strain rates (Chan, 2018).

To better analyze these data and unveil the underlying mechanisms, several theoretical approaches were adopted. Some phenomenological approaches were first developed (Zhang et al., 2006; Zhang et al., 2007; Zhang et al., 2009). However, the constitutive parameters of these models can hardly be related to relevant histological descriptors of the vocal tissues. Since 2010, a few authors have purposely proposed micro-mechanical models, including the architecture of vocal tissues, to provide new insights into voice biomechanics. Two modeling routes have been adopted:

(i) Poroelastic formulations have been developed to describe the fluid/solid phases of vocal tissues and to predict their dynamics (Miri et al., 2014; Tao et al., 2009; Scholp et al., 2020). However, such approaches rely on parameters for which experimental measurements (e.g., permeability and in situ observations of fluid dynamics) are still lacking.

(ii) Other authors have idealized the architecture of the fibrous networks of the lamina propria and the vocalis (e.g., using structural descriptors such as the fibril volume fraction, diameter, and preferred orientations) to derive their mechanical contribution from microstructural and micro-mechanical measurements (Miri et al., 2013; Kelleher et al., 2013b; Terzolo et al., 2022). This enabled the identification of the progressive elongation and reorientation of collagen fibrils and myofibrils, along with mechanical interactions between micro-constituents, which modulate the non-linear and anisotropic mechanics of vocal tissues (Terzolo et al., 2022). However, these micro-mechanical formulations have been developed within a general hyperelastic framework, thus neglecting the important dissipative and time-dependent mechanisms likely to develop during the vibrations of vocal tissues.

Therefore, this work aims to provide a multiscale mechanical model capable of reproducing the non-linear macroscopic visco-hyperelastic mechanical behavior of the vocal fold layers (i.e., lamina propria and vocalis) across a range of frequencies and strains, based on the knowledge of their architecture and mechanics at the fibril scale. To achieve this, we introduce microstructural time-dependent effects to the hyperelastic formulation developed by Terzolo et al. (2022). Based on histological and biomechanical data available in the literature, covering a wide range of loading modes, strain levels, and rates, the relevance of the model for predicting the time-dependent, multiscale mechanics of the vocal-fold layers is highlighted and discussed.

2 Formulation of the micro-mechanical model

2.1 Structural assumptions

The structural assumptions of the model are identical to those reported by Terzolo et al. (2022). In brief, both the lamina propria and the vocalis are considered incompressible composite materials made of a gel-like matrix (composed of cells, elastin, and ground substance for the lamina propria; of elastin, proteoglycans, and glycoproteins for the vocalis) reinforced by a network of connected and oriented fibril bundles (Figure 1):

$•$ For the lamina propria (Figure 1a, case ①), each fibril bundle is considered an assembly of parallel collagen fibrils with an initial diameter $d_0$, length ${\ell }_o^f$, and tortuosity ${\xi }_0={\ell }_0^f/{\ell }_0$, ${\ell }_0$ being their initial chord length. They are characterized by a waviness of approximately 10 monomodal sinusoids between nodes, with a wave amplitude $R_0$ and a spatial periodicity $H_0$, so that ${\ell }_0\approx 10H_0$ at rest.

$•$ For the vocalis (Figure 1a, case ②), each fibril bundle is considered an assembly of parallel myofibrils (with an initial diameter $d_{0m}$, tortuosity ${\xi }_{0m}$, wave amplitude $R_{0m}$, spatial periodicity $H_{0m}$, and chord length ${\ell }_0\approx 10H_{0m}$), surrounded by a sheath of collagen fibrils (with an initial diameter $d_{0c}$, tortuosity ${\xi }_{0c}$, wave amplitude $R_{0c}$, and spatial periodicity $H_{0c}$).

$•$ The fibrous architecture of the lamina propria exhibits a collagen fibril content $\mathrm{\Phi }$ (leading to $n_f$ collagen fibrils in Figure 1a), whereas the vocalis displays a collagen fibril content ${\mathrm{\Phi }}_c$ and a myofibril content ${\mathrm{\Phi }}_m$ (leading to ${\text{n}}_c$ collagen fibrils and ${\text{n}}_m$ myofibrils). Both tissues are idealized as networks of connected fibril bundles. These networks are built from the periodic repetition of a representative elementary volume (REV), composed of four fibril bundles connected to a central node ${\text{C}}_0$, and to four nodes ${\text{C}}_i$ of the corresponding neighboring REVs at their extremities (Figure 1b). At rest, each fibril bundle $i$ is also characterized by its initial mean orientation ${\mathbf{E}}_i$, as depicted in Figure 1b. This set of orientation directors introduces structural anisotropy. The distances between node ${\text{C}}_0$ and its unconnected neighbors ${\text{C}}_q$ (see dotted lines in Figure 1b), i.e., along the initial directions ${\mathbf{E}}_q={\mathbf{C}}_0{\mathbf{C}}_{\mathrm{q}}/\Vert {\mathbf{C}}_0{\mathbf{C}}_{\mathrm{q}}\Vert$, are noted ${\delta }_q$.

2.2 Micro-mechanical assumptions

Kinematics: When subjected to a macroscopic transformation gradient $\mathbf{F}$ and a macroscopic velocity gradient $\mathbf{L}$, the tissue REVs deform from their initial configuration to a deformed configuration. As a consequence, fibril bundles (un)fold so that their chord length is ${\ell }_i={\ell }_0\Vert \mathbf{F}\cdot {\mathbf{E}}_i\Vert$ in the deformed configuration, i.e., with a tensile stretch and strain ${\lambda }_i$ = ${\ell }_i/{\ell }_0$ and ${\epsilon }_i$ = $\ln {\lambda }_i$, respectively. This process occurs at a tensile strain rate ${\dot{\epsilon }}_i={\mathbf{e}}_i\cdot \mathbf{L}\cdot {\mathbf{e}}_i$. Moreover, fibril bundles also rotate so that their current mean orientation directors become ${\mathbf{e}}_i=\mathbf{F}\cdot {\mathbf{E}}_i/\Vert \mathbf{F}\cdot {\mathbf{E}}_i\Vert$ in the deformed configuration, thus introducing a strain-induced change in structural anisotropy. Finally, the rotation and the deformation of fibril bundles are not free and are hindered by steric effects between bundles. Steric effects are captured by restraining the motion of the node ${\text{C}}_0$ with respect to its unconnected neighbors ${\text{C}}_q$. These restrictions occur along ${\mathbf{e}}_q=\mathbf{F}\cdot {\mathbf{E}}_q/\Vert \mathbf{F}\cdot {\mathbf{E}}_q\Vert$ at a strain rate ${\dot{\epsilon }}_q={\mathbf{e}}_q\cdot \mathbf{L}\cdot {\mathbf{e}}_q$ (see dotted lines in Figure 1b), once the distance ${\delta }_q$ between ${\text{C}}_0$ and the neighboring nodes ${\text{C}}_q$ exceeds a critical distance ${\delta }_c$, i.e., below a contact strain ${\epsilon }_q$ = $\ln ({\delta }_q/{\delta }_c)$.

Mechanics of the matrix: Regardless of the considered tissue, the mechanics of their matrix is modeled as an incompressible hyperelastic neo-Hookean medium, with a strain energy function $W=0.5\mu (1-\mathrm{\Phi })(\mathrm{t}\mathrm{r}(\mathbf{F}\cdot {\mathbf{F}}^T)-3)$, which involves the shear modulus $\mu$ of the matrix.

Mechanics of the fibrils: The stretch (or the compression) of each fibril of a bundle $i$ generates a non-linear fibril reaction force. This force is noted as ${\mathbf{t}}_i=t_i{\mathbf{e}}_i$ for the collagen fibrils of the lamina propria, and ${\mathbf{t}}_{im}=t_{im}{\mathbf{e}}_i$ and ${\mathbf{t}}_{ic}=t_{ic}{\mathbf{e}}_i$ for the collagen fibrils and the myofibrils of the vocalis, respectively. To mimic both the non-linear elasticity observed during the tension-compression of collagen fibrils and their time-dependent response, the following decomposition of the reaction force is proposed for the lamina propria (similar decompositions are proposed for $t_{im}$ and $t_{ic}$ in the case of the vocalis):

$t_i=t_i^e+t_i^{ve},(1)$

where $t_i^e$ represents the (non-linear elastic) “neutral” response of the considered fibril, i.e., when the system attains its “relaxed” configuration. The expression proposed by Terzolo et al. (2022) is used: it provides relevant estimate of the unfolding of fibrils while accounting for their dimension (diameter $d_0$, chord length ${\ell }_0$, and tortuosity ${\xi }_0$) and mechanical properties (elastic modulus $E_f$). Thus, $t_i^e$ is an hyperelastic function of ${\epsilon }_i$:

$t_i^e=\frac{\pi d_0^2}{4}\left[E_{\mathit{eq0}}{\epsilon }_i+\frac{E_f-E_{\mathit{eq0}}}{2}\left({\epsilon }_i+\sqrt{{\left({\epsilon }_i-\ln {\xi }_0\right)}^2+{\alpha }^2}-\sqrt{{\ln }^2{\xi }_0+{\alpha }^2}\right)\right],(2)$

when the fibril is stretched; only the first term of the bracket is kept when the fibril is compressed. This expression involves a curvature parameter $\alpha$ that ensures a proper transition between bending- and stretching-dominated regimes during fibril unfolding. In addition, the initial apparent modulus of the fibril in the folded configuration $E_{\mathit{eq0}}=E_f⟨\cos {\beta }_0⟩/[⟨{\cos }^2{\beta }_0⟩+16⟨v^2⟩/d_0^2]$ (with $⟨\cdot ⟩=\frac{1}{{\ell }_0}{\int }_0^{{\ell }_0}\cdot du$, $⟨v^2⟩=R_0^2/2$, and ${\beta }_0=\arctan (2\pi \frac{R_0}{H_0}\cos \frac{2\pi }{H_0}u)$) is estimated from the literature (Potier-Ferry and Siad, 1992).

Moreover, in Equation 1, $t_i^{ve}$ represents time-dependent phenomena, including those related to the fibril deformation itself, the fibril interactions with the other fibrils and the surrounding gel-like matrix. These molecular-scale mechanisms exhibit characteristic relaxation times (Gautieri et al., 2012; Miri et al., 2013) that are not captured by the hyperelastic formulation proposed for $t_i^e$ in Equation 2. A fine quantification of these transient complex processes would require molecular-scale analyses based on statistical physics or numerical simulation using molecular dynamics approaches (Gautieri et al., 2011; Bantawa et al., 2023). Here, as a first approximation, we consider a simple approach at the scale of the fibrils to account for them. We assume that the aforementioned time-dependent phenomena can be reproduced using a non-linear viscoelastic Maxwell model (Equation 3), as schematized in Figure 1a:

${\dot{t}}_i^{ve}+\frac{E}{\eta }{t}_i^{ve}=\pi E\frac{d_0^2}{4}{\dot{\epsilon }}_i,(3)$

where $E$ and $\eta$ are the elastic modulus and the viscosity of the Maxwell model, respectively. As vocal-fold tissues exhibit several relaxation times over a wide range of strain rates (Chan and Titze, 1999; 2000; Chan and Rodriguez, 2008; Chan, 2018), it is necessary to include these effects in Equation 3. For example, SAOS studies (Chan and Rodriguez, 2008) performed on lamina propria samples report a Carreau-like evolution of the complex viscosity with the shear rate, i.e., with a Newtonian plateau at low shear rates and shear-thinning evolution at high shear rates. These aspects are taken into account by assuming that the viscosity $\eta$ is a non-linear Carreau function of the viscous strain rate, as shown in Equation 4:

$\eta ={\eta }_0{\left(1+{\left(\frac{\dot{{\epsilon }_i}-\frac{{4\dot{t}}_i^{ve}}{\pi E{d}_0^2}}{{\dot{\epsilon }}_0}\right)}^2\right)}^{\frac{n-1}{2}},(4)$

where ${\eta }_0$ is the viscosity of the Newtonian regime, ${\dot{\epsilon }}_0$ is the strain-rate transition between the Newtonian regime and the shear-thinning regime, and $n$ is the power-law index driving thinning effects at high strain rates. Expressions similar to Equations 2–4 are proposed for the vocalis, further assuming that $E_{0c}=E_{0m}=E$, ${\eta }_{0c}={\eta }_{0m}={\eta }_0$, ${\dot{\epsilon }}_{c0}={\dot{\epsilon }}_{m0}={\dot{\epsilon }}_0$, and $n_c=n_m=n$.

Steric interactions between fibril bundle: For both tissues, once the distance ${\delta }_q$ between the node ${\text{C}}_0$ and the neighboring nodes ${\text{C}}_q$ exceeds a critical distance ${\delta }_c$, i.e., below a contact strain ${\epsilon }_q$ = $\ln ({\delta }_q/{\delta }_c)$, steric interactions occur via reaction forces ${\mathbf{R}}_q=R_q{\mathbf{e}}_q$. A decomposition similar to Equation 1 is proposed to account for non-linear visco-hyperelastic effects (Equation 5):

$R_q=R_q^e+R_q^{ve},(5)$

where the hyperelastic term $R_q^e$ was proposed by Terzolo et al. (2022), as presented in Equation 6:

$R_q^e=\beta H\left({\epsilon }_q\right){\epsilon }_q^{\kappa },(6)$

where $H$ is the Heaviside function and $\beta$ and $\kappa$ are interaction parameters. To account for non-linear viscoelastic interactions, $R_q^{ve}$ is derived from the following non-linear Maxwell equation (Equation 7):

${\dot{R}}_q^{ve}+\frac{E^{\prime }}{{\eta }^{\prime }}{R}_q^{ve}=E^{\prime }{d}_0{\dot{\epsilon }}_q.(7)$

In analogy with Equation 4, the viscosity ${\eta }^{\prime }$ is assumed to be a Carreau function of the corresponding steric strain rate (Equation 8):

${\eta }^{\prime }={\eta }_0^{\prime }{\left(1+{\left(\frac{{\dot{\epsilon }}_q-\frac{{4\dot{R}}_q^{ve}}{\pi E^{\prime }{d}_0^2}}{\dot{{\epsilon }_0^{\prime }}}\right)}^2\right)}^{\frac{n-1}{2}},(8)$

where ${\eta }_0^{\prime }$ is the viscosity of the Newtonian regime and ${\dot{\epsilon }}_0^{\prime }$ is the transition strain rate between the Newtonian and thinning regimes.

2.3 Upscaling formulation: from micro- to macroscale mechanics

Given the structural and micro-mechanical features mentioned above, regardless of the tissue concerned, the macroscopic Cauchy stress tensor $\mathbf{\sigma }$ can be written as follows (Equation 9):

$\mathbf{\sigma }=-p\mathbf{\delta }+{\mathbf{\sigma }}_m+{\mathbf{\sigma }}_f+{\mathbf{\sigma }}_s,(9)$

where $p$ is the incompressibility pressure, $\mathbf{\delta }$ is the identity tensor, ${\mathbf{\sigma }}_m=\mathbf{F}\cdot {\left(\partial W/\partial \mathbf{F}\right)}^T$ is the stress contribution of the matrix, and ${\mathbf{\sigma }}_f$ and ${\mathbf{\sigma }}_s$ represent the stress contributions due to the (un)folding of fibrils and their steric interactions, respectively. Thus, Equations 10, 11 can be obtained as follows:

${\mathbf{\sigma }}_f=\frac{\mathrm{\Phi }}{\pi d_0^2{\xi }_0}{\displaystyle \sum _{i=1}^4}{t}_i{\lambda }_i{\mathbf{e}}_i\otimes {\mathbf{e}}_i(10)$

and

${\mathbf{\sigma }}_s=\frac{\mathrm{\Phi }}{\pi d_0^2{\xi }_0}{\displaystyle \sum _{q=1}^5}{R}_q{\delta }_q^{*}{\mathbf{e}}_q\otimes {\mathbf{e}}_q(11)$

for the lamina propria, where ${\delta }_q^{*}={\delta }_q/{\ell }_0$. Equations 12, 13 can be obtained as follows:

${\mathbf{\sigma }}_f=\frac{{\mathrm{\Phi }}_c}{\pi d_{0c}^2{\xi }_{0c}}{\displaystyle \sum _{i=1}^4}{t}_{ic}{\lambda }_i{\mathbf{e}}_i\otimes {\mathbf{e}}_i+\frac{{\mathrm{\Phi }}_m}{\pi d_{0m}^2{\xi }_{0m}}{\displaystyle \sum _{i=1}^4}{t}_{im}{\lambda }_i{\mathbf{e}}_i\otimes {\mathbf{e}}_i(12)$

and

${\mathbf{\sigma }}_s=\left(\frac{{\mathrm{\Phi }}_c}{\pi d_{0c}^2{\xi }_{0c}}+\frac{{\mathrm{\Phi }}_m}{\pi d_{0m}^2{\xi }_{0m}}\right){\displaystyle \sum _{q=1}^5}{R}_q{\delta }_q^{*}{\mathbf{e}}_q\otimes {\mathbf{e}}_q(13)$

for the vocalis. Thus, as an oversimplified representation, the proposed micro-mechanical model can be considered the imbrication of two anisotropic networks of non-linear Zener models embedded in an isotropic hyperelastic matrix (Figure 1): one for the mechanics of fibril bundles and another for their steric interactions. The mechanical response of the lamina propria (respectively vocalis) depends on 19 (respectively 25) histological and micro-mechanical parameters to be determined:

$•$ 6 (respectively 10) histological parameters: The fibril’s diameter $d_0$ (respectively $d_{0c}$ and $d_{0m}$), their waviness amplitude $R_0$ (respectively $R_{0c}$ and $R_{0m}$), spatial periodicity $H_0$ (respectively $H_{0c}$ and $H_{0m}$), from which their tortuosity ${\xi }_0$ (respectively ${\xi }_{0c}$ and ${\xi }_{0m}$) can be estimated, the fibril’s volume fraction $\mathrm{\Phi }$ (respectively ${\mathrm{\Phi }}_c$ and ${\mathrm{\Phi }}_m$), and initial 3D orientation (${\theta }_0$, ${\phi }_0$). These structural parameters can be determined from histological data.

$•$ 13 (respectively 15) mechanical parameters: The fibril’s Young’s modulus $E_f$ (respectively $E_{fc}$ and $E_{fm}$), the matrix shear modulus $\mu$, the transition parameter $\alpha$ (respectively ${\alpha }_c$ and ${\alpha }_m$), the elastic interaction coefficients $\beta$, $\kappa$, and ${\delta }_c$ related to steric effects, and the viscoelastic parameters $E$, ${\eta }_0$, ${\dot{\epsilon }}_0$, and $n$, along with $E^{\prime }$, ${\eta }_0^{\prime }$, and ${\dot{\epsilon }}_0^{\prime }$.

3 Model identification

3.1 Experimental database

The relevance of the model was evaluated by comparing its prediction with experimental data from the literature:

$•$ First, to assess the model relevance in the linear viscoelastic regime at small shear strains, we considered data collected by Chan and Rodriguez (2008): “cover” specimens were excised from seven donors (two female and five male donors), between 53 and 88 years old (mean age 67 years). Tissues were collected between 3 and 20 h post-mortem before being tested (mean time 10 h). The excised tissues were then subjected to SAOS under physiological conditions (T $\approx$ 37$°$C, 100$\%$ relative humidity). An oscillatory shear strain ${\gamma }_{zx}={\gamma }_0\sin (2\pi ft)$ was applied in the “longitudinal” plane $({\mathbf{e}}_z,{\mathbf{e}}_x)$, with a prescribed small shear strain amplitude ${\gamma }_0=0.01$, and a frequency $f$ varied from 1 to 250 Hz. In the following, trends derived from these seven donor-specific covers are represented by an “average target vocal-fold cover,” noted as ${\text{C}}_{\mathit{SAOS}}$.

$•$ Second, the model ability to reproduce oscillatory responses in the non-linear regime and finite strains was investigated with respect to data reported by Chan (2018). The author subjected a 60-year-old male “cover” to LAOS with several increasing strain amplitudes ${\gamma }_0=$ [0.05, 0.1, 0.2, 0.5, and 1] along the plane $({\mathbf{e}}_z,{\mathbf{e}}_x)$ at a prescribed frequency $f=$ 75 Hz. In the following, the sample chosen as a reference here is noted as ${\text{C}}_{\mathit{LAOS}}$.

$•$ Third, the model prediction was compared with vocal-fold layer samples deformed at finite strains and several relevant physiological loadings (i.e., tension, compression, and shear), as reported by Cochereau et al. (2020): two samples of lamina propria (covered by the very thin epithelium left intact, noted as ${\text{LP}}_1$ and ${\text{LP}}_2$), and two samples of vocalis (noted as ${\text{V}}_1$ and ${\text{V}}_2$). As a reminder, each sample was sequentially subjected to longitudinal tension along ${\mathbf{e}}_z$, transverse compression along ${\mathbf{e}}_x$, and longitudinal shear in the plane $({\mathbf{e}}_z,{\mathbf{e}}_x)$. For each loading mode, samples were subjected to 10 load/unload cycles up to Hencky strains ${\epsilon }_{zz}^{\mathit{max}}=$ 0.1, ${\epsilon }_{xx}^{\mathit{min}}=-$0.2, and shear ${\gamma }_{zx}^{\mathit{max}}=$ 0.6 at constant strain rates $|{\dot{\epsilon }}_{zz}|$, $|{\dot{\epsilon }}_{xx}|$, and $|{\dot{\gamma }}_{zx}|$ of $\approx 10^{-3}$ ${\text{s}}^{-1}$.

3.2 Optimization procedure

A protocol similar to that adopted by Terzolo et al. (2022) was applied to obtain optimized sets of histo-mechanical parameters:

$•$ For SAOS and LAOS experiments, all histological parameters were initialized and constrained within a corridor of admissible values deduced from the literature, as detailed by Terzolo et al. (2022): 0°$\le {\theta }_0\le$50°, 20°$\le {\phi }_0\le$90°, $10\mu \mathrm{m}\le H_0\le 70\mu \mathrm{m}$, $1\mu \mathrm{m}\le R_0\le 10\mu \mathrm{m}$, $10\mathrm{n}\mathrm{m}\le d_0\le 500\mathrm{n}\mathrm{m}$, and $0.15\le \mathrm{\Phi }\le 0.55$. For multi-axial experiments achieved with lamina propria and vocalis samples (Cochereau et al., 2020), we chose the histological parameters already determined by Terzolo et al. (2022), as reported in Table 1.

$•$ For SAOS and LAOS experiments, some of the hyperelastic parameters were constrained within physiological boundaries, i.e., the fibril Young modulus $1\mathrm{M}\mathrm{P}\mathrm{a}\le E_f\le 1\mathrm{G}\mathrm{P}\mathrm{a}$ and the matrix shear modulus $1\mathrm{P}\mathrm{a}\le \mu \le 1.5\mathrm{M}\mathrm{P}\mathrm{a}$. The other parameters, i.e., the transition parameters $\alpha$ and the interaction coefficients $\beta$, $\kappa$, and ${\delta }_c$, were left free. It is also important to note that steric interactions are not triggered during simple shear, thus leaving $\beta$, $\kappa$, and ${\delta }_c$ undetermined for SAOS and LAOS. For the multi-axial experiments performed with lamina propria and vocalis samples, we considered the hyperelastic parameters determined by Terzolo et al. (2022), except for the shear moduli of the matrices $\mu$, which were explored between $1\mathrm{P}\mathrm{a}$ and $1\mathrm{M}\mathrm{P}\mathrm{a}$ (see comments in the next section).

$•$ The positive viscoelastic parameters, i.e., $E$, ${\eta }_0$, ${\dot{\epsilon }}_0$, and $n$, along with $E^{\prime }$, ${\eta }_0^{\prime }$, and ${\dot{\epsilon }}_0^{\prime }$, were freely optimized for each of the experiments considered. For SAOS experiments, the power-law exponent $n$ was chosen between 0 and 1 to mimic the recorded shear-thinning behavior (Chan and Rodriguez, 2008). As the LAOS and the multi-axial experiments were performed at a unique strain rate, $n$ could not be determined and was arbitrarily set to the value found for SAOS experiments.

Table 1

Table 1. Optimized histological parameters for samples ${\text{C}}_{\mathit{SAOS}}$, ${\text{C}}_{\mathit{LAOS}}$, ${\text{LP}}_1$, ${\text{LP}}_2$, ${\text{V}}_1$, and ${\text{V}}_2$. Gray-colored columns refer to quantities computed as a function of the determined histological parameters.

A non-linear constraint optimization process based on a least-squared approach was used to minimize the discrepancies between the model prediction and the experimental macroscale stress–strain curves, as applied by Bailly et al. (2012) and Terzolo et al. (2022). The time integration of the implicit non-linear Maxwell differential Equations 3, 7 was achieved using the ode15i solver in ${\text{Matlab}}^{\text{®}}$ (Shampine, 2002).

4 Results and discussion

4.1 Relevance of histo-mechanical parameters

The set of optimized histological parameters used to reproduce the macroscopic rheological data during SAOS (Chan and Rodriguez, 2008), LAOS (Chan, 2018), and multi-axial loadings (Cochereau et al., 2020) are reported in Table 1. Apart from the remarks already stated by Terzolo et al. (2022) regarding the relevance of these parameters for ${\text{LP}}_i$ and ${\text{V}}_i$ samples, these values prompt the following comments:

$•$ The optimization led to a collagen content $\mathrm{\Phi }$ of $\approx$ 0.47 for ${\text{LP}}_i$ samples (i.e., including the epithelium, the cover, the intermediate, and the deep layers) versus only $\approx$ 0.30 for the cover ${\text{C}}_{\mathit{SAOS}}$. This finding is consistent with prior experimental evidence, showing that the first sublayer beneath the epithelium, i.e., the superficial layer of the lamina propria also called “Reinke’s space,” exhibits a fibril content lower than that found in the intermediate and deep layers of the lamina propria (Hahn et al., 2006b; Walimbe et al., 2017; Bailly et al., 2018).

$•$ The optimization also yielded a collagen fibril diameter $d_0$ close to 200 nm in the cover ${\text{C}}_{\mathit{SAOS}}$, against $d_0\approx$ 400 nm in the ${\text{LP}}_i$ samples. Such a decrease may be explained by the $d_0$-variations reported with the collagen type (Asgari et al., 2017) and with their location across the lamina propria (Gray et al., 2000; Tateya et al., 2006; Hahn et al., 2006a; Muñoz-Pinto et al., 2009; Walimbe et al., 2017; Benboujja and Hartnick, 2021). In particular, Muñoz-Pinto et al. (2009) measured that the content of “thin” (respectively “thick”) collagen fibrils decreases (respectively increases) steadily and approximately 10-fold (respectively 15-fold) from the superficial to the deep layers.

$•$ The optimized fibril tortuosity ${\xi }_0$ at rest is $\approx$ 20$\%$ higher for the ${\text{C}}_{\mathit{SAOS}}$ experiments than that estimated for the ${\text{LP}}_i$ samples. This is consistent with previous observations showing that the intermediate layer of the lamina propria is characterized by a dense network of straighter ECM fibrils compared with the superficial and deep layers (Klepacek et al., 2016; Bailly et al., 2018).

$•$ The histological parameters found for the cover sample ${\text{C}}_{\mathit{LAOS}}$ are very close to the values obtained for the cover sample ${\text{C}}_{\mathit{SAOS}}$. The main differences concern the initial fibril orientation (${\theta }_0$ and ${\phi }_0$) and tortuosity $({\xi }_0)$. This can be attributed to inter-subject variability.

$•$ The histological parameters of collagen fibrils in the vocalis are rather similar than those found for SAOS, LAOS, and ${\text{LP}}_i$ samples, except for the fibril content which is much lower. Conversely, the histological parameters of myofibrils are obviously very different.

In addition, the optimized micro-mechanical parameters used to reproduce the macroscopic rheological data during SAOS (Chan and Rodriguez, 2008), LAOS (Chan, 2018), and multi-axial loadings (Cochereau et al., 2020) are reported in Tables 2, 3 for hyperelastic and viscoelastic contributions, respectively. Terzolo et al. (2022) explained the relevance of the hyperelastic parameters in the cases of the ${\text{LP}}_i$ and ${\text{V}}_i$ samples. In addition, the following remarks can be put forward:

$•$ For the ${\text{LP}}_i$ and ${\text{V}}_i$ samples, the shear modulus of the matrix $\mu$-coefficient was re-optimized (within physiological boundaries) as the mechanical contribution of the matrix here is related to both the hyperelastic and the viscoelastic contributions (which encompass the fibril/surrounding matrix interactions). Thus, the optimization process led to 200 Pa and 190 Pa for ${\text{LP}}_1$ and ${\text{LP}}_2$, respectively, against 330 Pa and 290 Pa in Terzolo et al. (2022), and to 170 Pa for both ${\text{V}}_i$ samples, against 900 Pa and 980 Pa.

$•$ As emphasized in Table 2, the matrix shear modulus $\mu$ is nearly 10-fold lower for the cover samples ${\text{C}}_{\mathit{SAOS}}$ and ${\text{C}}_{\mathit{LAOS}}$ than for the entire ${\text{LP}}_i$ samples. The values identified for ${\text{C}}_{\mathit{SAOS}}$ and ${\text{C}}_{\mathit{LAOS}}$ are close to the range measured for the elastic shear modulus of hyaluronic acid ${\mu }_{HA}\approx$ 20–50 Pa [estimated at loading frequencies up to 10 Hz; Heris et al. (2012)], i.e., the most abundant polymer of the ground substance in the lamina propria. Hyaluronic acid, known to play a key role in shock absorption during vocal-fold collisions, is found at a higher volume fraction than collagen and elastin in the superficial layer, in contrast to the deep layer (Finck, 2008; Hahn et al., 2006a; Hahn et al., 2006b), which is in line with the identification results (see Table 1). The observed discrepancy in $\mu$-values in Table 2 is probably ascribed to the scarcity of elastin fibrils reported in the superficial layer (and therefore in the cover) in elderly tissues (Roberts et al., 2011).

$•$ The hyperelastic parameters related to the collagen fibril networks are very similar, regardless of the considered samples, i.e., SAOS and LAOS, along with ${\text{LP}}_i$ and ${\text{V}}_i$ samples. Due to the much softer passive mechanics of myofibrils, their hyperelastic parameters are much lower. Probably for the same reason, the optimized viscoelastic parameters ($E$, ${\eta }_0$, and ${\dot{\epsilon }}_0$) found for the lamina propria, the SAOS, and the LAOS samples differ by an order of magnitude from those reported for the vocalis.

$•$ The viscoelastic parameters of the ${\text{LP}}_i$ samples have been identified at a very low strain rate, i.e., close to ${\dot{\epsilon }}_0$. At this strain rate, the relaxation times $\tau \approx {\eta }_0/E\approx$ 3–15 s are obtained for both vocal-fold layers (similar relaxation times ${\tau }^{\prime }={\eta }_0^{\prime }/E^{\prime }$ were found for fibril bundle steric hindrance). It is interesting to note that these results are in line with the rare experimental data available at this scale (Yang, 2008; Shen et al., 2011; Gautieri et al., 2011). For example, Shen et al. (2011) reported typical relaxation times of solvated collagen fibrils in the range of 7–102 s. Moreover, Yang (2008) measured two distinct processes contributing to the stress relaxation of native collagen fibrils immersed in PBS buffer and subjected to 5–7 $\%$ strain for 5–10 min: a fast relaxation process with a characteristic time ${\tau }_1\approx 1.8\pm 0.4$ s and a slow relaxation process with ${\tau }_2\approx 60\pm 10$ s. Yang (2008) proposed that ${\tau }_1$ corresponds to the relative sliding of collagen microfibrils, while ${\tau }_2$ refers to the relative sliding of collagen molecules (due to the high level of cross-links between molecules). It is interesting to note that the characteristic times reported for the SAOS and LAOS samples are markedly lower, i.e., $\tau ={\eta }_0/E\approx 0.42$ s and 0.27 s, respectively. Considering that the model parameters for SAOS and LAOS were determined from experimental data acquired at high frequencies (from 1 to 250 Hz for SAOS and at 75 Hz for LAOS), these low-valued characteristic times are not surprising: additional data at lower strain rates would probably increase these values.

Table 2

Table 2. Optimized hyperelastic parameters for samples ${\text{C}}_{\mathit{SAOS}}$, ${\text{C}}_{\mathit{LAOS}}$, ${\text{LP}}_1$, ${\text{LP}}_2$, ${\text{V}}_1$, and ${\text{V}}_2$.

Table 3

Table 3. Optimized viscoelastic parameters for samples ${\text{C}}_{\mathit{SAOS}}$, ${\text{C}}_{\mathit{LAOS}}$, ${\text{LP}}_1$, ${\text{LP}}_2$, ${\text{V}}_1$, and ${\text{V}}_2$.

4.2 Relevance of the micro-mechanical model for SAOS

A comparison between the model predictions at macroscale and the SAOS experimental data is provided in Figure 2. In this figure, graphs (a) and (b) show the evolution of the shear storage and loss moduli $G^{\prime }$ and $G^{″}$ of sample ${\text{C}}_{\mathit{SAOS}}$ as functions of the excitation frequency $f$, whereas graphs (c) and (d) do the same for the loss factor ${\zeta }^{\star }=G^{″}/G^{\prime }$ and the dynamic viscosity ${\mu }^{\prime }=G^{″}/2\pi f$, respectively. In these graphs, the model predictions were extended up to $f=$ 1 kHz. Different remarks are highlighted from these graphs:

$•$ For all the rheological functions presented, a fairly good quantitative agreement is obtained between the model predictions (continuous lines) and the experimental data (marks): progressive increase in storage and loss moduli $G^{\prime }$ and $G^{″}$ with $f$ up to 200 Hz, power-law decrease in the viscosity ${\mu }^{\prime }$, and Carreau-like evolution of the loss factor ${\zeta }^{\star }$, with a marked power-law reduction above 10–50 Hz.

$•$ More particularly, it is interesting to note that the model nicely predicts the experimental “cross-over” zone at approximately 50–100 Hz, i.e., the zone within which (i) the storage modulus $G^{\prime }$ switches from lower to higher than the loss modulus $G^{″}$ and (ii) the loss factor ${\zeta }^{\star }$ switches from constant to remarkable decrease. This transition zone also coincides with that where some issues occur during vocal-fold vibration in human phonation. For fold vibration at low frequencies, i.e., below 50–100 Hz, viscous effects dominate $(G^{″}\ge G^{\prime })$ so that this should give rise to critical tissue overdamping, preventing proper periodic oscillations of vocal folds. In contrast, the dominant elastic properties at higher frequencies should restrain tissue damping (see the power-law decrease in the loss factor in Figure 2c), thus allowing the occurrence of proper periodic motion during vocal-fold vibration (Chan and Rodriguez, 2008).

$•$ To illustrate the role of histological parameters on the rheological response of SAOS samples, we have reported two additional discontinuous lines in Figure 2. These trends emphasize the effects induced by variations in the volume fraction of collagen fibrils $\mathrm{\Phi }$ (here, $\mathrm{\Phi }$ was chosen due to its wide variations between individuals but also within the vocal-fold layers themselves): the case where $\mathrm{\Phi }=$ 0.15 and the case where $\mathrm{\Phi }=$ 0.55, i.e., the minimum and maximum values found in the literature for the lamina propria. As shown in Figure 2, when $\mathrm{\Phi }$ varies in the physiological corridor, the qualitative trends are preserved for all viscoelastic properties ($G^{\prime }$, $G^{″}$, ${\zeta }^{\star }$, and ${\mu }^{\prime }$). However, the higher the fibril content, the higher the rheological functions, albeit with (i) marked differences (for $G^{\prime }$ at high frequencies, for ${\zeta }^{\star }$ at low frequencies, e.g., for $G^{″}$, and ${\mu }^{\prime }$ at all frequencies) and (ii) a slight shift of the cross-over zone toward lower frequencies as $\mathrm{\Phi }$ is increased. Note that the case of $\mathrm{\Phi }\approx$ 0 was also predicted in Figure 2 as a theoretical extreme case (not physiological), assuming a quasi-total absence of collagen fibers in the vocal-fold cover, which would thus become close to a homogeneous, isotropic neo-Hookean material with the same mechanical properties as the matrix alone. These simulations clearly emphasize the major mechanical role played by the collagen fibrous network, along with its interaction with the surrounding ground substance, in response to the oscillatory shear of the vocal-fold cover.

Figure 2

Figure 2. Experimental data (marks) vs. macroscale model predictions (lines) obtained for sample ${\text{C}}_{\mathit{SAOS}}$: storage $G^{\prime }$ modulus (a), loss $G\prime \prime$ modulus (b), loss factor ${\zeta }^{\star }$ (c), and dynamic viscosity ${\mu }^{\prime }$ (d), as functions of the oscillation frequency $f$. The continuous line represents the best fit of the model, with the others illustrating the effect of the collagen fibril content $\mathrm{\Phi }$. Source: experimental data adapted from Chan and Rodriguez (2008). Averaged data and standard deviations from seven human vocal-fold “cover” specimens are reported.

4.3 Relevance of the micro-mechanical model for LAOS

In Figure 3a, we have reported a collection of Lissajous stress–strain curves predicted by the model. These curves are compared with those in LAOS experiments obtained at a frequency $f$ = 75 Hz and cyclic amplitudes ${\gamma }_0$ varied from 0.05 to 0.5. In addition, Figure 3b presents a series of normalized Lissajous stress–strain curves predicted by the model in the Pipkin space $\{f,{\gamma }_0\}$ or $\{f,{\epsilon }_i^{\mathit{max}}\}$ (${\epsilon }_i^{\mathit{max}}$ is the maximal cyclic tensile strain the fibrils are subjected to), when $f$ and ${\gamma }_0$ are varied from 50 Hz to 1 kHz and from 0.05 to 0.5, respectively (Ewoldt et al., 2008; Chan, 2018). Within each contour plot, the black line represents the total visco-hyperelastic stress, whereas the red line is the hyperelastic or neutral stress contribution. Different trends can be highlighted.

$•$ Influence of the strain amplitude ${\gamma }_0$: Figure 3a shows a strong quantitative agreement between the model predictions (red line) and the experimental data at stabilized cycles when ${\gamma }_0\le$ 0.2. In particular, the model can capture the strong non-linear response of the tested sample, with, in particular, a proper modeling of the stress hysteresis induced by viscoelastic effects. In addition, the cyclic stress–strain curves progressively deviate from a linear strain hardening at low shear strain amplitudes (${\gamma }_0\le$ 0.1), which corresponds to the initial linear (un)folding of collagen fibril at small strains, toward a marked non-linear strain-hardening at higher strain magnitudes (in J-shape), where the non-linear hyperelastic stretching of collagen fibrils is triggered. This trend is also fairly well illustrated by the neutral stress responses of the Pipkin diagram shown in Figure 3b. This diagram also proves that the trend is preserved independently of the cycling frequency. Finally, it is worth noting from Figure 3a that the predicted strain-hardening at the highest strain magnitude ${\gamma }_0=0.5$ largely overestimates the cycle observed experimentally. Presumably, during the experiments, the tested cover exhibited a Mullins-like effect, as often observed in elastomers, gels, and soft living tissues (Diani et al., 2009; Peña et al., 2009; Rebouah et al., 2017; Rebouah and Chagnon, 2014; Zhan et al., 2024). This could yield to a stress softening of their mechanical behavior upon cycling. The Mullins effect can be caused by a number of irreversible mechanisms, e.g., the rupture of physical or covalent cross-links and the possible disentanglement of molecular chains. These mechanisms are not taken into account in the current micro-mechanical model. Yet, a possible way to account for these phenomena would consist in altering, with proper kinetics, the histo-mechanical properties of the collagen fibrils, such as their modulus $E_f$ (to account for damage) and their initial length ${\ell }_0^f$ or tortuosity ${\xi }_0$ (to account for disentanglement). In support of this hypothesis to be explored in future work, Figure 3a shows that lowering (respectively increasing) $E_f$ (respectively ${\xi }_0$) from 720 MPa to 400 MPa (respectively from 1.11 to 1.12) would lead to a more appropriate prediction of the experimental stress–strain curve performed at ${\gamma }_0=0.5$ (see green line).

$•$ Influence of the loading frequency $f$: As shown in Figure 3b, the loading frequency $f$ markedly alters the Lissajous curves. Regardless of the strain magnitude ${\gamma }_0$, the higher the $f$-value, the higher the stress levels and the stress hysteresis. These trends are in qualitative agreement with measurements acquired on other vocal-fold cover samples (Chan, 2018).

Figure 3

Figure 3. LAOS results: (a) macroscale stress–strain data vs. model predictions obtained for sample ${\text{C}}_{\mathit{LAOS}}$ tested at $f$ = 75 Hz and with ${\gamma }_0$ varied from 0.05 to 0.5. Source: experimental data adapted from Chan (2018). (b) Predicted Lissajous stress–strain curves plotted in the Pipkin space $\{f,{\gamma }_0\}$ or $\{f,{\epsilon }_i^{\mathit{max}}\}$, where ${\epsilon }_i^{\mathit{max}}$ is the maximal cyclic tensile strain the fibrils are subjected to. Black solid lines represent the total stress, while red solid lines represent neutral contribution.

4.4 Relevance of the micro-mechanical model for finite strain multi-axial cyclic loadings

Macroscopic stress–strain predictions are compared with the reference experimental data in Figure 4 (respectively Figure 5), for the lamina propria and vocalis samples ${\text{LP}}_1$ and ${\text{V}}_1$ (respectively ${\text{LP}}_2$ and ${\text{V}}_2$), and the three cyclic loadings these samples were subjected to, i.e., longitudinal tension, transverse compression, and longitudinal shear. For each case, the “neutral” curve already predicted by Terzolo et al. (2022) was reported (see dotted lines). The strain-induced evolution of micro-mechanical descriptors during cyclic tension is displayed in Figure 6 (an illustrative case of ${\text{LP}}_1$ and ${\text{V}}_1$), with compression and shear results summarized in Figure 7 (an illustrative case of ${\text{LP}}_1$).

Figure 4

Figure 4. Macroscopic viscoelastic stress–strain curves of vocal-fold sublayers under multi-axial cyclic loadings. Experimental data vs. model predictions obtained for the lamina propria sample ${\text{LP}}_1$ (left, in red) and the vocalis sample ${\text{V}}_1$ (right, in blue): (a) longitudinal tension; (b) transverse compression; (c) longitudinal shear. Several model predictions are compared: the visco-hyperelastic model with steric interactions between fibril bundle (solid lines); the non-viscous hyperelastic model previously described by Terzolo et al. (2022) (noted ‘hyp.,’ dashed lines); the visco-hyperelastic model albeit with no steric interactions between fibril bundle (dashed–dotted lines) for compression loading solely (panel b).

Figure 5

Figure 5. Same as Figure 4 for samples ${\text{LP}}_2$ (left, in red) and ${\text{V}}_2$ (right, in blue).

Figure 6

Figure 6. Evolution of multiscale descriptors for the lamina propria ${\text{LP}}_1$ (a) and the vocalis ${\text{V}}_1$ (b) during tension along ${\mathbf{e}}_z$: (top left) macroscopic strain paths; (bottom) stereographic projection of the four orientation vectors ${\mathbf{e}}_i$ from initial to final state; (top right) tensile strain of the fibril chord ${\epsilon }_i$ and corresponding tension $t_i$.

Figure 7

Figure 7. Evolution of multiscale descriptors for the lamina propria sample ${\text{LP}}_1$ during compression along ${\mathbf{e}}_x$ (a), and shear in the plane (${\mathbf{e}}_z$ and ${\mathbf{e}}_x$) (b). Compression case: (top left) macroscopic strain paths; (bottom) stereographic projection of the four orientation vectors ${\mathbf{e}}_i$ from initial to final state; (top right) tensile strain of the fibril chord ${\epsilon }_i$ and interaction forces $R_q$. Shear case: (top left) tensile strain of the fibril chord ${\epsilon }_i$; (top right) corresponding tension $t_i$; (bottom) same as in (a).

The results for the first loading cycle are discussed below for each loading mode:

$•$ Longitudinal tension: The model prediction for longitudinal tension along ${\mathbf{e}}_z$ is fairly good for both the lamina propria and the vocalis samples, as emphasized in Figures 4a–5a. In particular, compared with the hyperelastic formulation proposed by Terzolo et al. (2022), i.e., the neutral curves, the model can now capture the stress hysteresis and the residual strains after unloading. These tendencies are inherited from microscale viscoelastic effects, together with the rearrangement of the tissue microstructures. This is illustrated in Figure 6 and Supplementary Figure S1, in which one can assess the irreversible unfolding and rotation of fibrils that are predicted during cyclic tension, both for the lamina propria and the vocalis samples. It is interesting to note that the predicted stress hysteresis and residual strain of collagen fibrils were experimentally observed by Yang (2008). For the vocalis, the predicted tensions in both collagen fibrils, $t_{ic}$, and myofibrils, $t_{im}$, are plotted in the inset of Figure 6b. Although the key role played by the sheaths of collagen fibers surrounding muscle fibers in the passive tensile properties of the tissue was already evidenced during monotonic loading, the strong contribution of myofibrils to inelastic effects and residual strains after unloading is clearly highlighted here.

$•$ Transverse compression: Figures 4b, 5b prove that the model predictions are also in good agreement with the experimental data recorded during transverse compression. Moreover, as already pointed out by Terzolo et al. (2022), steric interactions are of major importance for the lamina propria and the vocalis mechanics during compression. This characteristic is preserved with the visco-hyperelastic formulation: if steric hindrance effects are deactivated in the model (see model predictions with “no steric interactions” in Figures 4b, 5b, dash-dotted lines), the deformation of the visco-hyperelastic fibril bundles is not sufficient to capture the lamina propria stress hysteresis and residual strain experimentally observed. Thus, fibril bundle repulsion forces $R_q$ and their viscoelastic contributions $R_q^{ve}$ appear to be of critical importance to properly reproduce the compression behavior of both vocal-fold layers (see model predictions with “steric interactions” in Figures 4b, 5b, solid lines). No other significant microscale deformation mechanisms (such as rotation or noticeable unfolding of fibrils) were predicted under transverse compression (Figure 7a).

$•$ Longitudinal shear: The mechanical contribution of the matrix plays a major role in the overall shear response of the lamina propria and the vocalis, as already stressed by Terzolo et al. (2022). On this basis, the fibril viscoelastic properties and interactions with the surrounding ground substance allowed, via the microscopic tension $t_i$ (Figure 7b), to satisfactorily reproduce the experimental trends observed during the load/unload sequence at the tissue scale (Figures 4c–5c).

Finally, the relevance of the visco-hyperelastic model to simulate the sequential series of 10 load–unload cycles and the tissue response as a function of load history is assessed. Figure 8 shows the comparison of the theoretical predictions with the reference cyclic data for the three loading modes. If the decrease in stress hysteresis is qualitatively well captured by the model once the first cycle has been completed in tension, compression, and shear, the predictions fail to simulate the progressive decrease in peak stresses measured after repeated loading paths, along with the increase in residual strains after repeated unloading paths, which are particularly observed in tension and compression. According to the model, a steady state is reached practically after the first load/unload sequence, whereas the stabilized behavior is only really observed experimentally after the $5^{\text{th}}$ cycle (or even up to the ${10}^{\text{th}}$ cycle, depending on the sample and loading mode). As mentioned for LAOS results, these accommodation behaviors resemble Mullins-like effects, which are not taken into account in the present formulation of the model.

Figure 8

Figure 8. Same as in Figure 4, albeit for 10 cycles: experimental data vs. model predictions. The experimental 10th cycle is displayed in green symbols for the lamina propria sample ${\text{LP}}_1$ (left) and in orange symbols for the vocalis sample ${\text{V}}_1$ (right). Experimental intermediate cycles are not reported for the sake of clarity.

4.5 Relevance of the model for predicting future patho/physiological variations and assisting biomedical developments

The micro-mechanical model developed in this work has been calibrated to reproduce the microstructural specificities and multiscale behavior of healthy human vocal-fold tissues, combining a wide range of histo-mechanical measurements collected from the available literature. By adjusting these input data, it could be adapted without major difficulties to predict the mechanics of pathological human vocal tissues (Hantzakos et al., 2009; Finck, 2008), animal vocal tissues (Li et al., 2024), or (bio)composites developed to replace/reconstruct the fibrous architecture and vibro-mechanical performance of the vocal folds after surgery (Heris et al., 2012; Li et al., 2016; Jiang et al., 2019; Latifi et al., 2018; Ravanbakhsh et al., 2019; Ferri-Angulo et al., 2023).

The model could also be used to predict the evolution of the mechanical properties of the same tissue following an alteration in its microstructural arrangement, due, for example, to (i) its growth and remodeling with age [by simulating, e.g., a progressive decrease in the volume fraction of elastin, an increase in that of collagen, and muscle atrophy (Roberts et al., 2011; Kuhn, 2014; Li et al., 2024)], (ii) scarring lesions [by simulating fibrosis and an increase in the collagen content, along with changes in fibril tortuosity compared to the healthy case (Heris et al., 2015; Li et al., 2016)], (iii) the appearance of a lesion following phonotrauma [by simulating damage mechanisms likely to occur at the fibril’s level (Miller and Gasser, 2022)], or (iv) a therapeutic treatment [by simulating the addition of a soft hydrogel to the matrix composite, for example (Li et al., 2016; Mora-Navarro et al., 2026)].

To better understand the impact of these histological variations on vocal-fold vibrations at the larynx level (in the case of native tissue but also injured, repaired, or replaced tissue), this original constitutive law should be implemented in a finite element code reproducing the vocal folds in their 3D anatomical geometry, as in current 3D phonation models (Döllinger et al., 2023). In doing so, microstructure-based simulations could improve the currently limited knowledge of the links between the specific microarchitecture of the vocal folds and their unique macroscale vibratory performance. Moreover, they should guide the design of fiber-reinforced biomaterials currently under development for functional voice restoration.

5 Conclusion

A better understanding of human phonation requires an in-depth study of the viscoelastic properties of vocal folds. To this end, this study proposes to enrich a recent 3D micro-mechanical model of vocal-fold tissues, which was hitherto capable of predicting their non-linear, elastic, and anisotropic mechanical behavior at various spatial scales (micro to macro) (Terzolo et al., 2022). This was achieved by adding viscoelastic mechanisms at the scale of their collagen fibril and myofibril bundles. These improvements now enable the model to capture the viscoelastic properties of vocal-fold tissues from small to finite strains, such as their non-linear strain-rate sensitivity—on which their damping and oscillation onset properties strongly depend; their stress-hysteretic response, and the inelastic deformations typically measured during cyclic loading. In addition, the model allows the microstructural rearrangements to be predicted, which are often very challenging to identify experimentally.

This model was successfully used to reproduce various sets of ex vivo data available in the literature and complement them with original theoretical data, providing specific micro-mechanism scenarios for each. This identification was carried out for a wide variety of loading conditions at different rates: low-frequency cyclic tension, compression, and shear in large deformations; and high-frequency oscillatory shear from small to large deformations (SAOS for the linear viscoelasticity regime and LAOS for the non-linear viscoelasticity regime). The model predictions are in quantitative agreement with macroscopic experimental trends and clearly highlight the key impact of microscopic histo-mechanical descriptors on vocal-fold dynamics, such as the volume fraction of collagen fibrils in the cover, their tortuosity at rest, their mechanics, and their interactions. This micro-mechanical model can be implemented in finite element codes to further simulate the transient dynamics of vocal folds with relevant histo-mechanical properties.

However, some model limitations should be improved. For example, coarse-grained atomistic/molecular simulations would probably provide relevant information to strengthen the physical links between the time-dependent nanostructural rearrangements and the phenomenological approach proposed herein at the fibril scale. Furthermore, the model does not allow the Mullins-like effects commonly observed in vocal tissues to be adequately described: combined with additional experiments focused on this aspect, the model could be improved based on formulations proposed for other materials, such as structured elastomers (Rebouah and Chagnon, 2014; Rebouah et al., 2017).

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving humans were approved by the French ethical and safety laws related to body donation and the Institutional Review Board of the University of Texas Southwestern Medical Center (Willed Body Program). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and institutional requirements.

Author contributions

AT: Investigation, Methodology, Software, Writing – original draft, Writing – review and editing, Formal analysis, Visualization. LB: Conceptualization, Investigation, Funding acquisition, Methodology, Project administration, Supervision, Validation, Writing – review and editing, Formal Analysis, Resources. LO: Conceptualization, Investigation, Methodology, Supervision, Validation, Writing – review and editing, Formal Analysis.

Funding

The authors declare that financial support was received for the research and/or publication of this article. This work was funded by the ANR MicroVoice (Grant No. ANR-17-CE19-0015-01) and the LabEx Tec21 (Investissements d’Avenir—grant agreement no. ANR-11-LABX-0030). The 3SR Laboratory is part of the PolyNat Carnot Institute (Investissements d’Avenir—grant agreement no. ANR16-CARN0025-01).

Acknowledgements

The authors would like to thank Nathalie Henrich Bernardoni (CNRS, GIPSA-lab, Grenoble) for her helpful assistance in this work.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2025.1670567/full#supplementary-material

References

Asgari, M., Latifi, N., Heris, H. K., Vali, H., and Mongeau, L. (2017). In vitro fibrillogenesis of tropocollagen type III in collagen type i affects its relative fibrillar topology and mechanics. Sci. Rep. 7, 1392–undefined. doi:10.1038/s41598-017-01476-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailly, L., Geindreau, C., Orgéas, L., and Deplano, V. (2012). Towards a biomimetism of abdominal healthy and aneurysmal arterial tissues. J. Mech. Behav. Biomed. Mater. 10, 151–165. doi:10.1016/j.jmbbm.2012.02.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailly, L., Cochereau, T., Orgéas, L., Henrich Bernardoni, N., Rolland du Roscoat, S., McLeer-Florin, A., et al. (2018). 3D multiscale imaging of human vocal folds using synchrotron X-ray microtomography in phase retrieval mode. Sci. Rep. 8, 14003. doi:10.1038/s41598-018-31849-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Bantawa, M., Keshavarz, B., Geri, M., Bouzid, M., Divoux, T., McKinley, G. H., et al. (2023). The hidden hierarchical nature of soft particulate gels. Nat. Phys. 19, 1178–1184. doi:10.1038/s41567-023-01988-7

CrossRef Full Text | Google Scholar

Benboujja, F., and Hartnick, C. (2021). Quantitative evaluation of the human vocal fold extracellular matrix using multiphoton microscopy and optical coherence tomography. Sci. Rep. 11, 2440. doi:10.1038/s41598-021-82157-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, R. W. (2004). Measurements of vocal fold tissue viscoelasticity: approaching the Male phonatory frequency range. J. Acoust. Soc. Am. 115, 3161–3170. doi:10.1121/1.1736272

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, R. W. (2018). Nonlinear viscoelastic characterization of human vocal fold tissues under large-amplitude oscillatory shear (LAOS). J. Rheology 62, 695–712. doi:10.1122/1.4996320

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, R. W., and Rodriguez, M. L. (2008). A simple-shear rheometer for linear viscoelastic characterization of vocal fold tissues at phonatory frequencies. J. Acoust. Soc. Am. 124, 1207–1219. doi:10.1121/1.2946715

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, R. W., and Titze, I. R. (1999). Viscoelastic shear properties of human vocal fold mucosa: measurement methodology and empirical results. J. Acoust. Soc. Am. 106, 2008–2021. doi:10.1121/1.427947

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, R. W., and Titze, I. R. (2000). Viscoelastic shear properties of human vocal fold mucosa: theoretical characterization based on constitutive modeling. J. Acoust. Soc. Am. 107, 565–580. doi:10.1121/1.428354

PubMed Abstract | CrossRef Full Text | Google Scholar

Cochereau, T., Bailly, L., Orgéas, L., Henrich Bernardoni, N., Robert, Y., and Terrien, M. (2020). Mechanics of human vocal folds layers during finite strains in tension, compression and shear. J. Biomechanics 110, 109956. doi:10.1016/j.jbiomech.2020.109956

PubMed Abstract | CrossRef Full Text | Google Scholar

Dashatan, S. H., Sit, M., Zhang, Z., Grossmann, E., Millot, J., Huang, Y., et al. (2023). Enhanced vibration damping and viscoelastic properties of flax/epoxy composites and their carbon fibre hybrid laminates. Compos. Part A Appl. Sci. Manuf. 175, 107819. doi:10.1016/j.compositesa.2023.107819

CrossRef Full Text | Google Scholar

DeJonckere, P., and Lebacq, J. (2020). Lung volume affects the decay of oscillations at the end of a vocal emission. Biomed. Signal Process. Control 62, 102148. doi:10.1016/j.bspc.2020.102148

CrossRef Full Text | Google Scholar

Diani, J., Fayolle, B., and Gilormini, P. (2009). A review on the Mullins effect. Eur. Polym. J. 45, 601–612. doi:10.1016/J.EURPOLYMJ.2008.11.017

CrossRef Full Text | Google Scholar

Döllinger, M., Zhang, Z., Schoder, S., Šidlof, P., Tur, B., and Kniesburges, S. (2023). Overview on state-of-the-art numerical modeling of the phonation process. Acta Acust. 7, 25. doi:10.1051/aacus/2023014

CrossRef Full Text | Google Scholar

Ewoldt, R. H., Hosoi, A. E., and McKinley, G. H. (2008). New measures for characterizing nonlinear viscoelasticity in large amplitude oscillatory shear. J. Rheology 52, 1427–1458. doi:10.1122/1.2970095

CrossRef Full Text | Google Scholar

Ferri-Angulo, D., Yousefi-Mashouf, H., Michel, M., McLeer, A., Orgéas, L., Bailly, L., et al. (2023). Versatile fiber-reinforced hydrogels to mimic the microstructure and mechanics of human vocal-fold upper layers. Acta Biomater. 172, 92–105. doi:10.1016/j.actbio.2023.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Finck, C. (2008). Implantation d’acide hyaluronique estérifié lors de la microchirurgie des lésions cordales bénignes (Doctoral thesis). ORBi-Université de Liège: ULiège, Belgium. Available online at: https://orbi.uliege.be/handle/2268/142631.

Google Scholar

Gautieri, A., Vesentini, S., Redaelli, A., and Buehler, M. J. (2011). Hierarchical structure and nanomechanics of collagen microfibrils from the Atomistic scale up. Nano Lett. 11, 757–766. doi:10.1021/nl103943u

PubMed Abstract | CrossRef Full Text | Google Scholar

Gautieri, A., Vesentini, S., Redaelli, A., and Buehler, M. J. (2012). Viscoelastic properties of model segments of collagen molecules. Matrix Biol. 31, 141–149. doi:10.1016/j.matbio.2011.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, S. D., Alipour, F., Titze, I. R., and Hammond, T. H. (2000). Biomechanical and histologic observations of vocal fold fibrous proteins. Ann. Otology, Rhinology Laryngology 109, 77–85. doi:10.1177/000348940010900115

PubMed Abstract | CrossRef Full Text | Google Scholar

Hahn, M. S., Kobler, J. B., Starcher, B. C., Zeitels, S. M., and Langer, R. (2006a). Quantitative and comparative studies of the vocal fold extracellular matrix I: elastic fibers and hyaluronic acid. Ann. Otology, Rhinology Laryngology 115, 156–164. doi:10.1177/000348940611500213

PubMed Abstract | CrossRef Full Text | Google Scholar

Hahn, M. S., Kobler, J. B., Zeitels, S. M., and Langer, R. (2006b). Quantitative and comparative studies of the vocal fold extracellular matrix II: collagen. Ann. Otology, Rhinology Laryngology 115, 225–232. doi:10.1177/000348940611500311

PubMed Abstract | CrossRef Full Text | Google Scholar

Hantzakos, A., Remacle, M., Dikkers, F. G., Degols, J. C., Delos, M., Friedrich, G., et al. (2009). Exudative lesions of Reinke’s space: a terminology proposal. Eur. Archives Oto-Rhino-Laryngology 266, 869–878. doi:10.1007/s00405-008-0863-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Heris, H. K., Rahmat, M., and Mongeau, L. (2012). Characterization of a hierarchical network of hyaluronic acid/gelatin composite for use as a smart injectable biomaterial. Macromol. Biosci. 12, 202–210. doi:10.1002/mabi.201100335

PubMed Abstract | CrossRef Full Text | Google Scholar

Heris, H. K., Latifi, N., Vali, H., Li, N., and Mongeau, L. (2015). Investigation of chitosan-glycol/glyoxal as an injectable biomaterial for vocal fold tissue engineering. Procedia Eng. 110, 143–150. doi:10.1016/j.proeng.2015.07.022

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, L., Jiang, Y., Stiadle, J., Wang, X., Wang, L., Li, Q., et al. (2019). Electrospun nanofibrous thermoplastic polyurethane/poly(glycerol sebacate) hybrid scaffolds for vocal fold tissue engineering applications. Mater. Sci. Eng. C 94, 740–749. doi:10.1016/j.msec.2018.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelleher, J. E., Siegmund, T., Du, M., Naseri, E., and Chan, R. W. (2013a). Empirical measurements of biomechanical anisotropy of the human vocal fold lamina propria. Biomechanics Model. Mechanobiol. 12, 555–567. doi:10.1007/s10237-012-0425-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelleher, J. E., Siegmund, T., Du, M., Naseri, E., and Chan, R. W. (2013b). The anisotropic hyperelastic biomechanical response of the vocal ligament and implications for frequency regulation: a case study. J. Acoust. Soc. Am. 133, 1625–1636. doi:10.1121/1.4776204

PubMed Abstract | CrossRef Full Text | Google Scholar

Klemuk, S. A., and Titze, I. R. (2004). Viscoelastic properties of three vocal-fold injectable biomaterials at low audio frequencies. Laryngoscope 114, 1597–1603. doi:10.1097/00005537-200409000-00018

PubMed Abstract | CrossRef Full Text | Google Scholar

Klepacek, I., Jirak, D., Duskova Smrckova, M., Janouskova, O., and Vampola, T. (2016). The human vocal fold layers. Their delineation inside vocal fold as a background to create 3D digital and synthetic glottal model. J. Voice 30, 529–537. doi:10.1016/j.jvoice.2015.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Klos, A., Bailly, L., Rolland du Roscoat, S., Orgéas, L., Henrich Bernardoni, N., Broche, L., et al. (2024). Optimising 4d imaging of fast-oscillating structures using x-ray microtomography with retrospective gating. Sci. Rep. 14, 20499. doi:10.1038/s41598-024-68684-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Koruk, H., and Rajagopal, S. (2024). A comprehensive review on the viscoelastic parameters used for engineering materials, including soft materials, and the relationships between different damping parameters. Sensors 24, 6137. doi:10.3390/s24186137

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhn, M. (2014). Histological changes in vocal fold growth and aging. Curr. Opin. Otolaryngology and Head Neck Surg. 22, 460–465. doi:10.1097/MOO.0000000000000108

PubMed Abstract | CrossRef Full Text | Google Scholar

Latifi, N., Asgari, M., Vali, H., and Mongeau, L. (2018). A tissue-mimetic nano-fibrillar hybrid injectable hydrogel for potential soft tissue engineering applications. Sci. Rep. 8, 1047doi. doi:10.1038/s41598-017-18523-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Stiadle, J. M., Lau, H. K., Zerdoum, A. B., Jia, X., Thibeault, S. L., et al. (2016). Tissue engineering-based therapeutic strategies for vocal fold repair and regeneration. Biomaterials 108, 91–110. doi:10.1016/j.biomaterials.2016.08.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Lin, X., Xie, X., Chen, X., Xie, Y., and Sun, G. (2024). Histological characterization of rat vocal fold across different postnatal periods. Laryngoscope Investig. Otolaryngol. 9, e70018. doi:10.1002/lio2.70018

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, C., and Gasser, T. C. (2022). A bottom-up approach to model collagen fiber damage and failure in soft biological tissues. J. Mech. Phys. Solids 169, 105086. doi:10.1016/j.jmps.2022.105086

CrossRef Full Text | Google Scholar

Miri, A. K. (2014). Mechanical characterization of vocal fold tissue: a review study. J. Voice 28, 657–667. doi:10.1016/j.jvoice.2014.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

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 Biomater. 9, 7957–7967. doi:10.1016/j.actbio.2013.04.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Miri, A. K., Heris, H. K., Mongeau, L., and Javid, F. (2014). Nanoscale viscoelasticity of extracellular matrix proteins in soft tissues: a multiscale approach. J. Mech. Behav. Biomed. Mater. 30, 196–204. doi:10.1016/j.jmbbm.2013.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora-Navarro, C., Smith, E., Wang, Z., del, C., Ramos-Alamo, M., Collins, L., et al. (2026). Injection of vocal fold lamina propria-derived hydrogels modulates fibrosis in injured vocal folds. Biomater. Adv. 178, 214424. doi:10.1016/j.bioadv.2025.214424

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Pinto, D., Whittaker, P., and Hahn, M. S. (2009). Lamina propria cellularity and collagen composition: an integrated assessment of structure in humans. Ann. Otology, Rhinology Laryngology 118, 299–306. doi:10.1177/000348940911800411

PubMed Abstract | CrossRef Full Text | Google Scholar

Peña, E., Peña, J. A., and Doblaré, M. (2009). On the Mullins effect and hysteresis of fibered biological materials: a comparison between continuous and discontinuous damage models. Int. J. Solids Struct. 46, 1727–1735. doi:10.1016/j.ijsolstr.2008.12.015

CrossRef Full Text | Google Scholar

Potier-Ferry, M., and Siad, L. (1992). Homogénéisation géométrique d’une poutre ondulée. (geometrical homogenization of a corrugated beam). Comptes rendus l’Académie Sci. t314, 425–430. Available online at: http://pascal-francis.inist.fr/vibad/index.php?action=getRecordDetail&idt=5332674.

Google Scholar

Radolf, V., Horáček, J., Bula, V., Geneid, A., and Laukkanen, A.-M. (2022). Damping of human vocal fold vibration. in Engineering mechanics, 28, 321–324.

Google Scholar

Ravanbakhsh, H., Bao, G., Latifi, N., and Mongeau, L. G. (2019). Carbon nanotube composite hydrogels for vocal fold tissue engineering: biocompatibility, rheology, and porosity. Mater. Sci. Eng. C 103, 109861. doi:10.1016/j.msec.2019.109861

PubMed Abstract | CrossRef Full Text | Google Scholar

Rebouah, M., and Chagnon, G. (2014). Permanent set and stress-softening constitutive equation applied to rubber-like materials and soft tissues. Acta Mech. 225, 1685–1698. doi:10.1007/s00707-013-1023-y

CrossRef Full Text | Google Scholar

Rebouah, M., Chagnon, G., and Heuillet, P. (2017). Anisotropic viscoelastic models in large deformation for architectured membranes. Mech. Time-Dependent Mater. 21, 163–176. doi:10.1007/s11043-016-9324-x

CrossRef Full Text | Google Scholar

Roberts, T., Morton, R., and Al-Ali, S. (2011). Microstructure of the vocal fold in elderly humans. Clin. Anat. 24, 544–551. doi:10.1002/ca.21114

PubMed Abstract | CrossRef Full Text | Google Scholar

Scholp, A., Jeddeloh, C., Tao, C., Liu, X., Dailey, S. H., and Jiang, J. J. (2020). Study of spatiotemporal liquid dynamics in a vibrating vocal fold by using a self-oscillating poroelastic model. J. Acoust. Soc. Am. 148, 2161–2172. doi:10.1121/10.0002163

PubMed Abstract | CrossRef Full Text | Google Scholar

Shampine, L. F. (2002). Solving 0 = F(t, y(t), y′(t)) in matlab. J. Numer. Math. 10, 291–310. doi:10.1515/JNMA.2002.291

CrossRef Full Text | Google Scholar

Shen, Z. L., Kahn, H., Ballarini, R., and Eppell, S. J. (2011). Viscoelastic properties of isolated collagen fibrils. Biophysical J. 100, 3008–3015. doi:10.1016/j.bpj.2011.04.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Švec, J. G., Horáček, J., Šram, F., and Veselý, J. (2000). Resonance properties of the vocal folds: in vivo laryngoscopic investigation of the externally excited laryngeal vibrations. J. Acoust. Soc. Am. 108, 1397–1407. doi:10.1121/1.1289205

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, C., Jiang, J. J., and Zhang, Y. (2009). A fluid-saturated poroelastic model of the vocal folds with hydrated tissue. J. Biomechanics 42, 774–780. doi:10.1016/j.jbiomech.2008.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Tateya, T., Tateya, I., and Bless, D. M. (2006). Collagen subtypes in human vocal folds. Ann. Otology, Rhinology Laryngology 115, 469–476. doi:10.1177/000348940611500612

PubMed Abstract | CrossRef Full Text | Google Scholar

Terzolo, A., Bailly, L., Orgéas, L., Cochereau, T., and Henrich Bernardoni, N. (2022). A micro-mechanical model for the fibrous tissues of vocal folds. J. Mech. Behav. Biomed. Mater. 128, 105118. doi:10.1016/j.jmbbm.2022.105118

PubMed Abstract | CrossRef Full Text | Google Scholar

Titze, I. R., Klemuk, S. A., and Gray, S. (2004). Methodology for rheological testing of engineered biomaterials at low audio frequencies. J. Acoust. Soc. Am. 115, 392–401. doi:10.1121/1.1631941

PubMed Abstract | CrossRef Full Text | Google Scholar

Vampola, T., Horáček, J., and Klepáček, I. (2016). Computer simulation of mucosal waves on vibrating human vocal folds. Biocybern. Biomed. Eng. 36, 451–465. doi:10.1016/j.bbe.2016.03.004

CrossRef Full Text | Google Scholar

Walimbe, T., Panitch, A., and Sivasankar, P. M. (2017). A review of hyaluronic acid and hyaluronic acid-based hydrogels for vocal fold tissue engineering. J. Voice 31, 416–423. doi:10.1016/j.jvoice.2016.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L. (2008). Mechanical properties of collagen fibrils and elastic fibers explored by AFM. Enschede: University of Twente.

Google Scholar

Zhan, L., Qu, S., and Xiao, R. (2024). A review on the Mullins effect in tough elastomers and gels. Acta Mech. Solida Sin. 37, 181–214. doi:10.1007/s10338-023-00460-6

CrossRef Full Text | Google Scholar

Zhang, K., Siegmund, T., and Chan, R. W. (2006). A constitutive model of the human vocal fold cover for fundamental frequency regulation. J. Acoust. Soc. Am. 119, 1050–1062. doi:10.1121/1.2159433

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Siegmund, T., and Chan, R. W. (2007). A two-layer composite model of the vocal fold Lamina propria for fundamental frequency regulation. J. Acoust. Soc. Am. 122, 1090–1101. doi:10.1121/1.2749460

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Siegmund, T., Chan, R. W., and Fu, M. (2009). Predictions of fundamental frequency changes during phonation based on a biomechanical model of the vocal fold lamina propria. J. Voice official J. Voice Found. 23, 277–282. doi:10.1016/J.JVOICE.2007.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar