Soft palate muscle activation: a modeling approach for improved understanding of obstructive sleep apnea

Abstract

A Hill model-based phenomenological method for muscle activation was used to investigate defectiveness of the palatal muscle tone during sleep for obstructive sleep apnea (OSA) patients. Based on the stretch–stress characteristic of muscle activation when the eccentric contraction is considered, a specifically defined phenomenological strain-energy function was used, as well as the Holzapfel-type strain-energy function for the passive part. A continuum mechanical framework, including the stress tensor and elasticity tensor, was obtained, based on the defined strain-energy function. The model parameters were obtained by fitting the constitutive model to experimental test data. Three-dimensional patient-specific geometry was modeled, accounting for the muscle tissue layer and based on the quantitative histology study of the soft palate. Anatomically representative boundary conditions for the finite element calculation were also considered. Palatal muscle activation level (electromyographic data) versus the negative pressure was defined in the simulations, and the patients’ activation level was set to be lower than for the healthy people. The simulation results showed that reduced in activation level for the patients causes a less negative closing pressure, and this makes the soft palate more prone to collapse. In addition, if we account for the passive–active transfer displayed as the muscle contraction corresponding to the neurogenic reflex in the soft palate, the collapse is prevented. This numerical representation of the reduced activation for the OSA patients may provide increased understanding of OSA physiology.

Appendix A

The analytical solution for the created constitutive model is achieved in the incompressible case. In the incompressible case, the strain-energy function (7) can be rewritten as

$$\begin{aligned} \begin{aligned} \Psi =\,&\underset{\mathrm {passive}}{\underbrace{c(I_{1}-3)+\left\{ \begin{array}{ll} \frac{k_{1}}{2k_{2}} [{e}^{{k}_{2}(I_{4}-1)^{2}}-1]&{} \mathrm {if}\quad I_{4}\ge 1\\ 0&{}\mathrm {if}\quad I_{4}<1 \end{array}\right. }}\\&+\,\underset{\mathrm {active}}{\underbrace{f(I_{4s})}}+ \underset{\mathrm {penalty}}{\underbrace{p(J-1)}}, \end{aligned} \end{aligned}$$

(35)

where, p is Lagrange multiplier. The detailed description for the active part reads

$$\begin{aligned} f(I_{4s})=\int _{I_\mathrm {4min}}^{I_{4}}g(I_{4})dI_{4}, \end{aligned}$$

(36)

with

$$\begin{aligned} g(I_{4})=\left\{ \begin{array}{ll} \frac{1}{2}P _{\mathrm {max}}^{\mathrm {act}}\frac{1}{\sqrt{I_{4}}}\frac{I_\mathrm {4min}-I_{4}}{I_\mathrm {4min}-I_\mathrm {4opt}}\mathrm {exp}\left( \frac{(2 I_\mathrm {4min}-I_{4}-I_\mathrm {4opt})(I_{4}-I_\mathrm {4opt})}{2(I_\mathrm {4min}-I_\mathrm {4opt})^{2}} \right) &{} \mathrm {if}\quad I_{4}\ge I_{4\mathrm {min}} \\ 0&{} \mathrm {if}\quad I_{4}< I_{4\mathrm {min}} \end{array}\right. \end{aligned}$$

(37)

Note that, since the active part is assumed incompressible, the constant parameters $\bar{I}_\mathrm {4min}$ and $\bar{I}_\mathrm {4opt}$ have the same value as $I_\mathrm {4min}$ and $I_\mathrm {4opt}$. Then, being similar to Eq. (10) and Eq. (15), the second Piola–Kirchhoff stress in the incompressible case reads

$$\begin{aligned} \begin{aligned} \mathbf {S}=\,&\mathbf {S}_{\mathrm {pas}}+\mathbf {S}_{\mathrm {act}}+\mathbf {S}_{\mathrm {pen}}\\ =\,&2c\mathbf {I} +2k_{1}(\bar{I}_{4}-1)e^{k_{2}(\bar{I}_{4}-1)^{2}}\mathbf {a}_{0}\otimes \mathbf {a}_{0}\\&+\,P _{\mathrm {max}}^{\mathrm {act}}\frac{1}{\sqrt{I_{4}}}\frac{I_\mathrm {4min}-I_{4}}{I_\mathrm {4min}-I_\mathrm {4opt}}\\&\quad \mathrm {exp}\left( \frac{(2 I_\mathrm {4min}-I_{4}-I_\mathrm {4opt})(I_{4}-I_\mathrm {4opt})}{2(I_\mathrm {4min}-I_\mathrm {4opt})^{2}} \right) \mathbf {a}_{0}\otimes \mathbf {a}_{0}\\ {}&+\,p\mathbf {C}^{-1}. \end{aligned} \end{aligned}$$

(38)

Based on Eq. (16)$_{1}$ and $J=1$ for the incompressible case, the Cauchy stress is calculated as

$$\begin{aligned} \begin{aligned} \varvec{\sigma }=\,&\mathbf {F}\mathbf {S}\mathbf {F}^\mathrm{T}\\ =\,&2c\mathbf {B}+2k_{1}(I_{4}-1)e^{k_{2}(I_{4}-1)^{2}} \mathbf {a}\otimes \mathbf {a}\\&+\,P _{\mathrm {max}}^{\mathrm {act}}\frac{1}{\sqrt{I_{4}}}\frac{I_\mathrm {4min}-I_{4}}{I_\mathrm {4min}-I_\mathrm {4opt}}\\&\quad \mathrm {exp}\left( \frac{(2 I_\mathrm {4min}-I_{4}- I_\mathrm {4opt})(I_{4}-I_\mathrm {4opt})}{2(I_\mathrm {4min}-I_\mathrm {4opt})^{2}} \right) \mathbf {a}\otimes \mathbf {a}\\ {}&+\,p\mathbf {I}. \end{aligned} \end{aligned}$$

(39)

Hence, setting the fiber direction as $\mathbf {a}_{0}=[1,0,0]$ and defining the deformation gradient for uniaxial test in the fiber direction as

$$\begin{aligned} \mathbf {F}=\begin{bmatrix} \lambda&\quad 0&\quad 0\\ 0&\quad 1/\sqrt{\lambda }&\quad 0\\ 0&\quad 0&\quad 1/\sqrt{\lambda } \end{bmatrix}, \end{aligned}$$

(40)

the analytical solution for the uniaxial test can be achieved and fitted with the numerical simulation.

In addition, since $\mathrm {tr}[\mathbf {a}\otimes \mathbf {a}]=I_{4}=\lambda _{\mathrm {fiber}}^{2}$ ($\lambda _{\mathrm {fiber}}$ denotes the fiber stretch), the active Cauchy stress in the fiber direction in (39)$_{3}$ is displayed as

$$\begin{aligned} \sigma _{\mathrm {fiber}}^{\mathrm {act}}= & {} P _{\mathrm {max}}^{\mathrm {act}}\frac{1}{\sqrt{I_{4}}}\frac{I_\mathrm {4min}-I_{4}}{I_\mathrm {4min}-I_\mathrm {4opt}}\nonumber \\&\mathrm {exp}\left( \frac{(2 I_\mathrm {4min}-I_{4}- I_\mathrm {4opt})(I_{4}-I_\mathrm {4opt})}{2(I_\mathrm {4min}-I_\mathrm {4opt})^{2}} \right) I_{4}. \end{aligned}$$

(41)

Then, the nominal stress ($\mathbf {P}=J\mathbf {F}^{-1}\varvec{\sigma }$) for the fiber activation part is deduced as

$$\begin{aligned} \begin{aligned} P_{\mathrm {fiber}}^{\mathrm {act}}=\,&\varvec{\sigma }_{\mathrm {fiber}}^{\mathrm {act}} /\lambda _{\mathrm {fiber}}\\ =\,&P _{\mathrm {max}}^{\mathrm {act}}\frac{I_{4}}{\sqrt{I_{4}}\lambda _{\mathrm {fiber}}}\frac{I_\mathrm {4min}-I_{4}}{I_\mathrm {4min}-I_\mathrm {4opt}}\\&\mathrm {exp}\left( \frac{(2 I_\mathrm {4min}-I_{4}- I_\mathrm {4opt})(I_{4}-I_\mathrm {4opt})}{2(I_\mathrm {4min}-I_\mathrm {4opt})^{2}} \right) \\ =&P _{\mathrm {max}}^{\mathrm {act}}\frac{I_\mathrm {4min}-I_{4}}{I_\mathrm {4min}-I_\mathrm {4opt}}\\&\mathrm {exp}\left( \frac{(2 I_\mathrm {4min}-I_{4}- I_\mathrm {4opt})(I_{4}-I_\mathrm {4opt})}{2(I_\mathrm {4min}-I_\mathrm {4opt})^{2}} \right) . \end{aligned} \end{aligned}$$

(42)

With this, for the constitutive parameters, satisfying fitting results to the experimental data used nominal stress can be achieved. This also explains how the added term, $1/\sqrt{\bar{I}_{4}}$, in Eq. (9) works in the transformation between the nominal stress and Cauchy stress.