|Home | About | Journals | Submit | Contact Us | Français|
This study presents a new computational system for modeling the upper airway in rats that combines tagged magnetic resonance imaging (MRI) with tissue material properties to predict three-dimensional (3D) airway motion. The model is capable of predicting airway wall and tissue deformation under airway pressure loading up to airway collapse. The model demonstrates that oropharynx collapse pressure depends primarily on ventral wall (tongue muscle) elastic modulus and airway architecture. An iterative approach that involves substituting alternative possible tissue elastic moduli was used to improve model precision. The proposed 3D model accounts for stress-strain relationships in the complex upper airway that should present new opportunities for understanding pathogenesis of airway collapse, improving diagnosis and developing treatments.
Obstructive sleep apnea (OSA) is a prevalent disorder, occurring in 9% of middle-aged men and 4% of middle-aged women1,2. Both anatomic factors and neuromuscular control factors are believed to contribute to the disorder, however our understanding of the mechanical properties of the pharynx and the pathogenesis of OSA is still very limited3. A major limitation to further understanding is the lack of a unifying model that can incorporate a variety of data and produce simulated airway collapse responses. Existing theoretical models include the Starling resistor model, that represents the airway as a collapsible tube with a critical collapsing pressure (Pcrit) that quantifies the collapsibility of the airway4 (i.e. higher Pcrit indicates greater collapsibility), and the balance of forces model, that predicts that the imbalance between the counteracting dilating (muscle) and collapsing (negative airway internal pressure) forces in patients with OSA results in airway closure during sleep5. However, these models do not offer a way to account for anatomical structure and material properties in heterogeneous tissues surrounding the airway, nor are they designed to compare dynamic magnetic resonance (MR) images with predicted airway mechanical changes and collapse pressures.
Some investigators have used a 2D Finite Element (FE) model approach that involves the analysis of sagittal MR images and pressure - flow measurements to model the human pharyngeal airway6,7. Although this method can integrate fluid and tissue mechanics, including muscle activity, to predict airway collapse behavior, a 2D model ignores the complicated three-dimensional (3D) airway anatomy, particularly the lateral pharyngeal airway walls that are an important determinant of airway patency8,9. Recently, Xu et al introduced a 3D FE modeling method that utilizes advanced MR tissue tracking methods10, i.e. spatial modulation of magnetization (SPAMM)11, with reconstruction of a 3D volume from a stack of axial images, i.e. image-based modeling. In the current report we extend this novel approach to the mechanics of airway collapse in a rat model, and we detail the process to relate airway mechanical parameters and tissue material properties to the rat anatomy. We also demonstrate the use of a 3D FE model to explore the roles of different muscle groups, airway shape and structure, in airway collapse mechanics.
Magnetic resonance imaging (MRI) was used to obtain soft tissue anatomy of the rat, and was followed by a second MRI series that examined the dynamic, controlled pressure collapsing of the pharynx. With the approval of Institutional Animal Care and Use Committee one non-obese Zucker rat under isoflurane anesthesia (3% in oxygen) was surgically prepared with a mid-cervical tracheotomy to direct airflow to the lungs, while a second cannula (PE 190), rostral to the tracheotomy and 2 mm below the larynx, was inserted rostrally to provide air pressure or suction to the upper airway. The mouth and one naris were sealed while both the rostral tracheal cannula and a cannula through the second naris were connected together and to the controlled pressure source. Pressure inside the sealed upper airway segment was measured by a non-magnetic pressure transducer (SPC350MR, Millar Instruments, Houston) in line with the cannula at the naris. Through solenoid valves outside the magnet, a stimulator (Grass S88, AstoMed, Quincy) initiated or triggered the gated MRI acquisition, and controlled the intraluminal static pressure (+5 CmH2O; 1 cmH2O = 98.1 Pa) or alternated positive and negative pressure (+8 and −8 cmH2O) to the isolated airway segment. Pressure delivery at both rostral and caudal connections ensured equal pressure throughout the pharynx in spite of any possible single collapse point.
Two MR imaging protocols were used: (1) a fast spin-echo protocol to obtain high contrast detailed images for segmentation; and (2) a spoiled gradient echo recalled (SPGR) sequence with spatial modulation of magnetization (SPAMM) that produced a time series of images that allowed for tracking soft tissue motion during alternating applied pressures in the airway. MRI was performed in a 40 cm bore, 4.7 Tesla magnet (Magnex) with a Varian INOVA Console, using a 12 cm shielded gradient insert. The SPGR axial images were acquired on a 128 by 128 pixel matrix; 40×40 mm field of view (FOV), 22 contiguous images, slice thickness = 1 mm; repetition time (TR) = 87 ms; echo time (TE), 2.2 ms; with 4 averages. Fast spin echo (FSE) images used the same pixel matrix and FOV and number of axial images as SPGR, with (TR) = 2000 ms, echo time (TE) = 6.8 ms with excitation train length (ETL) = 8 and 4 averages. A set of mid sagittal images (FOV = 80×80 mm) was also acquired for SPGR and FSE protocols.
A complete image series (22 axial slices, Figure 1A) was acquired by 128 single-line image data acquisitions, repeated 4 times (4 averages) to increase the signal-to-noise ratio of the images. Figure 1B is an annotated (see legend) mid-sagittal slice from the FSE protocol that shows the supine orientation of the rat and the boundaries of the pharynx from which the axial slices in both FSE and SPGR were acquired. Thus, SPAMM tagging was applied at +8 cmH2O (Figure 2A, sagittal slices) just prior to the first image while the SPAMM lines show tissue movement in the second image (Figure 2B, “barline 2”), presumably a result of the negative intraluminal pressure application (−8 cmH2O, see Figure 2C point 2). The distance between centers of two adjacent tag lines was 2 mm (2.5 mm in sagittal images, Figure 1B), with tag thickness at 0.3 mm. SPAMM lines degrade after 1-2 seconds (based on T1 values) so the pressure cycle time (850 ms) was set to less than the predicted T1 value.
The FE model was generated from a fast spin echo image set by interactively identifying different airway tissue groups using commercial image-based modeling software (Amira ®, TGS Template Graphics Software, Inc.). The regions of interest were identified by appropriate color masks and grown into a 3D voxel-representation. Images acquired at 5 cmH2O internal pressure were segmented, producing tissue regions including: nasopharynx (NP, airway dorsal to the soft palate), oropharynx (OP, airway ventral to the soft palate bounded by dorsal surface of the tongue), soft palate (separating nasopharynx and oropharynx), tongue (ventral to oropharynx and sliding relative to lateral soft tissue) and other soft tissue (lateral and sublingual tissues) (Figure 1A). The dynamic MRI data set was used to ascertain which regions were motionless (i.e. rigid), including the dorsal/ posterior wall, skull, hard palate, mandible, and other motionless tissue. Essential tissue groups that surround NP and OP, i.e. tongue, soft palate and lateral tissue were investigated in this study. Accuracy of segmentation was justified by comparison of MR images to a standard rat anatomy reference12 and by previous experience in this lab13-16. Following reconstruction, triangular surface meshes bounding each tissue region were initialized and improved individually using a second commercially available software (HyperMesh; Altair), producing triangulated surface with 18000, 4900, 23000 and 70000 triangular elements for mandibles, palate, tongue and tissue, respectively. Local grid refinement and sensitivity studies were conducted. The resultant surface meshes were subsequently converted into tetrahedral volumetric elements (type=C3D4) and imported to a nonlinear FE modeling package (ABAQUS/CAE 6.4, HKS Inc. Pawtucket, RI).
Coupling between different anatomic regions was defined from observing SPAMM image pairs acquired under positive to negative applied airway pressure. Several types of boundary/interface conditions were defined. Locations of SPAMM line discontinuity were identified as sliding interfaces, and all other interfaces were assumed to be continuous. Negative intraluminal airway pressure (e.g. −8 cmH2O) was applied to the surfaces of the tissue bounded by the nasopharynx and oropharynx as a model boundary condition. Nodes attached to rigid regions (i.e. the mandible and skull) were fixed. Ventral and lateral tissue boundaries displaying significant motion were modeled as traction-free surfaces. Due to the small to moderate strain of the tongue and other pharyngeal tissues during airway collapse, soft tissue was modeled as linear elastic, with mechanical parameters adapted from literature7,17. General contact pairs were defined to automatically specify the contact and sliding between the dorsal tongue and ventral palate surfaces as the airway collapses. Tie constraints were applied at other tissue interfaces where no relative motion was observed. Finally, as rostral-to-caudal movements of the tongue were less than one axial slice thickness from SPAMM displacement (measured from sagittal plane images, Figure 2, panel A to B), the rostral and caudal ends of the tongue were allowed only motion in the ventral-dorsal and medial-lateral directions.
ABAQUS/Explicit 6.3 (HKS Inc. Pawtucket, RI) was used to simulate the deformation and resulting stress distribution in the airway model. The NLGEOM option was activated to account for nonlinear geometric effects as the airway collapses. Both transient dynamic and quasi-static analysis are suitable to model the non-linear, brief, and time-dependent upper airway collapse. Since this model did not include any time dependent changes in tissue material properties, quasi-static analysis was more appropriate and was performed, using a period of 1 second, with semi-automatic time increment on the order of 5.5 × 10-5 s. However, the results of a transient simulation for the pressure waveform used (Figure 2C), were not significantly different from the quasi-static results at the time of second image acquisition, thus transient analysis could be used if there were additional tissue parameters to be evaluated and if images were acquired at intermediate times.
Properties of various tissue groups were derived by comparing tissue motions predicted by the FE model to those measured by tagged MR images. An axial MR image plane at the middle pharyngeal level (noted in Figure 1A) with clearly identifiable nasopharynx, oropharynx and soft palate was selected for detailed motion quantification. A rectangular region of interest (ROI) enclosing the nasopharynx and oropharynx was defined, in which motion of each pixel was tracked using a validated optical flow method (OFM) incorporated into a SPAMMVU software package.18 The OFM is a registration and change visualization program that uses a hierarchical estimation technique to compute the flow field that describes the warping of one phase or time series image into alignment with the next. This method is particularly suited for analyzing tagged MR images due to the high-resolution, high-contrast tag lines that serve as landmarks for confident pixel identification and tracking. The flow fields were evaluated between the pressure phases. The displacement of every pixel in the ROI was obtained, and a subset of these data points was used to define triangular elements (two points per 2.0 mm SPAMM line width). Regional tissue wall motion was then visualized by tracking the displacement of the centroid of each triangular element, as described previously15.
The coordinates of each triangle centroid were correlated to FE model nodes within 2 pixels of the triangle centroid. Differences between predicted (FE model) and measured (MRI-SPAMM) pharyngeal tissue displacement were minimized by modifying the global material properties of each anatomical part iteratively until (a) the model predicted the deformation patterns obtained from SPAMM and (b) the differences measured at selected locations were within an acceptable range (less than 20%).
A Young’s modulus value that was uniform throughout was used for each tissue region. Initial material properties (Young’s modulus, E) were optimized using starting values slightly above those given in the literature17,19: for tongue 8000 Pa, soft palate 15000 Pa, and other airway soft tissues 13000 Pa. Tissue moduli were varied manually to minimize tissue displacement errors compared with in-vivo image displacement measurements. Optical flow displacements analysis of images have an uncertainty of approximately 0.175 mm equivalent to 0.5 pixel.18
After estimation of tissue properties from the experiment to create a baseline model, a sensitivity study was performed to explore the effects of tissue properties on airway collapsibility. Regional tissue properties were individually varied, including stiffening or softening of the tongue, softening of soft palate, and softening of other upper airway tissue (not tongue or soft palate) by increasing (stiffening) and decreasing (softening) the value of E as given in Table 1. Airway collapse was also simulated with all upper airway soft tissues assumed to be homogenous, to demonstrate the effect of heterogeneous tissue group properties on model accuracy. Collapse was simulated by decreasing airway pressure Paw (defined by Paw = −ct, where t = time and c = 800 Pa/s) until collapse was detected in the model. The stress history of the contact between the tongue and soft palate was monitored, and the collapse pressure (Pc) was determined as Paw when contact stress between these two surfaces exceeded zero. The sensitivity of the model in response to the varying material property is defined by:
where Pc,base is collapse pressure of the baseline model, Pc is the collapse pressure of the model with one material property varied, Xi,base is the ith material property in the base line model, and Xi is the ith material property in the varied model. Units for sensitivity by this definition are dimensionless and given as a percentage.
The baseline FE model was optimized (after 6 iterations) with Young’s moduli (E values) of: 6000 Pa, 7000 Pa, and 15000 Pa for tongue, tissue, and soft palate respectively, resulting in average displacement error of 15% between the MRI measurements and finite element model predictions. With the optimized E values, the baseline model was capable of predicting physiologically correct deformation patterns in response to decreasing Paw. This pattern was noted in qualitative observations (Figure 3). In contrast to the oropharynx, whose size decreased towards collapse as Paw decreased, the nasopharynx remained patent and did not show obstruction during the simulations. Quantitative comparison at the same axial plane (at mid-pharynx) in the FE model and MRI measurements showed maximum displacement predicted by the FE model (1.58 mm, Figure 3A) was almost identical to that measured in the MRI analysis (1.56 mm, Figure 3B). However when uniform modulus was assumed for all tissue regions (i.e. tongue, soft palate and other soft tissues regions having the same E), both qualitative deformation patterns and quantitative displacements differed significantly from the MRI dynamic (SPAMM) data. Tissue strain was calculated using the deformed SPAMM grids. Under pre-defined negative internal pressure that represents inspiration, two dimensional maximum principal strain of 0.09 ± 0.1 is found in the rostro-caudal direction, accompanied with minimum principal strain of −0.11 ± 0.09 in anterior-posterior direction.
A 3D projection with a sagittal mid-plane slice of airway tissues (Figure 4) shows that in comparison to the tongue, the soft palate was more constrained in the ventral to dorsal direction of motion, and relatively rigid compared to surrounding tissues. The nasopharynx that is constrained and supported by surrounding bony structures remained open at −8 cmH2O while in additional simulations at lower pressures, the nasopharynx remained patent down to internal pressure −12.2 cmH2O.
The statistical comparison between the measured and predicted displacements in the sampled axial image plane (Figure 5) demonstrated that with optimized material properties, FE model predicts displacements that are well correlated with those measured from SPAMM images (R2 = 0.91). However, the accuracy of displacement measurement is impacted by multiple factors, especially in regions where displacement magnitude approaches the pixel size of the SPAMM images (0.3125 mm).
In dynamic loading studies simulated in the finite element model, using the baseline tissue material properties, the oropharynx collapsed at an internal pressure of −560 Pa (−5.7 cmH2O). Table 1 shows the sensitivity of oropharynx collapse pressure to changes in the mechanical properties of various tissue regions, modeled by independently varying the Young’s modulus (E) of each tissue region and simulating dynamic airway collapse. Collapse pressure was most sensitive to the changes in the modulus of the tongue, had relatively low sensitivity to soft palate modulus, and was moderately sensitive to lateral wall and other tissue modulus.
This paper describes a novel FE model of the 3D pharyngeal airway in a rat, that is based on anatomic images and designed to simulate the effects of airway deformation under controlled intraluminal pressure changes. The 3D FE model simulates motion and displacement direction of important tissue groups, and collapse or patency of oropharynx and nasopharynx. When compared to in-vivo dynamic MR measurements the airway model results showed qualitatively similar tissue displacements and airway deformation.
We found disagreement between in-vivo results and a model with uniform tissue properties compared to better agreement with heterogeneous properties, supporting the approach of using different tissue properties in different pharyngeal tissues. Thus, the results tend to support the hypothesis that pharyngeal tissues are heterogeneous with regard to Young’s modulus of elasticity.
The methods used are based on those of Moulton et al.20, who estimated Young’s modulus for passive myocardium using a two-dimensional isotropic, nonlinear small deformation FE model that estimated strain measurements from MRI tagging. Malhotra et al.20 developed a 2D FE model using signal-averaged anatomic data from human subjects. By assuming closing pressure in the pharyngeal airway to be −500 Pa (−5.1 cmH2O) for passive condition and −1300 Pa (−13.25 cmH2O) for sleeping condition, they estimated the Young’s Modulus of tongue and uvula by matching the 2D simulation to the closing pressure. Their results showed that E was 6000 Pa and 12000 Pa for the passive conditions. Unfortunately, there has been no method to date to measure airway tissue mechanical properties to validate the computational model other than in-vitro mechanical testing on isolated samples of airway tissues or in-vivo indentation testing, and in-vivo indentation tests have not been published for these tissues.
The model introduced here may be used to explore the relative roles of airway architecture such as cross sectional shape or to simulate different material properties of airway tissue regions. The highest tissue displacements occurred in the tongue. The tongue’s compliance depends in part on its low Young’s modulus relative to the soft palate. However, the tongue has the largest surface area exposed to negative intraluminal airway pressure, and is free to deform compared to lateral tissues that are supported by the motionless mandible. As a result, at negative intraluminal pressure the tongue undergoes a ventral to dorsal displacement and collapses the oropharynx, which can be visualized and confirmed in our simulation and in-vivo imaging. In contrast, the lateral tissue, with modulus similar to the tongue, was relatively more stable. Figure 2 (Panel A to B) shows a midsagittal image of the rat during alternate pressure application and shows the location of hypopharynx collapse, and the relatively small rostral to caudal motion due to pressure application.
The mechanical behavior of the soft palate is similar to a thick plate-like bending under the defined surface loading, with little rotation or translation at the rostral end where it is attached to the hard palate at the aponeureisis. Except for the most dorsal tissue attached to the skull, the bulk of the tissues are pulled inward by the ventro-dorsal motions of tongue and lateral to medial motions of the lateral walls. Surprisingly, when uniform material properties were assigned to all tissues including tongue and soft palate, displacement was still much higher in the tongue compared to lateral wall tissues or the soft palate. However, since the predicted collapse pressure for uniform tissue property conditions and heterogeneous tissue property conditions was similar it could be argued that tissue property differences were not an important factor in oropharyngeal collapse.
Limitations to this study include an approach that simplified airway anatomy in order to make the model work within the parameters we could test. Although this approach limited the number and orientation of tissue types and force distributions, the FE model results (dimensions at specific pressure levels) showed displacement and collapse predictions that were consistent with in-vivo observations. Another limitation is that biologic tissues are not fully described by a linear elastic material with isotropic material properties. At this time, nonlinear material model parameters and test data are not available for pharyngeal tissues, and therefore a linear model was used, consistent with other recent approaches.7,17 For moderate to large strains, a finite strain material model such as a hyperelastic material is preferable to the linear elastic model, but was not used in this study due to the range of strains expected. The difference between these two material models was estimated, on a unit tissue sample with a neo-Hookean hyperelastic material model, to be less than 5% relative error in the strain, over the range of stretch ratios in the model (up to 1.19). This difference is significantly smaller than the overall error of the tissue displacement, but high enough to suggest that a finite strain formulation might be useful for future models. Additionally, we modeled quasi steady state conditions by design that allowed for tissue viscosity to be neglected while isotropic (not anisotropic) bulk material properties provided a useful first approximation. Lastly, our primary purpose in this paper was to demonstrate the relative role of airway architecture in airway collapse. We believe that as more experimental data on the structure and tissue mechanical properties are available, the model we present could incorporate that information as nonlinear and anisotropic tissue properties and thus play an important role in understanding the pathophysiology of pharyngeal collapse. Future work could be combined with more detailed anatomic segmentation and tissue properties, i.e. anisotropic qualities of tissue architecture. Levine et al.21 provide such an FEM approach to human tongue motion modeling in speech while a global pharyngeal airway model in humans is still lacking.
While the FE 3D model allows out-of-plane motion, that motion was limited by the boundary condition (i.e. no out-of-plane motion) at the front of the tongue. Thus, we did not optimize the model properties based on observed or anticipated out-of-plane motion. In previous studies rostral-to-caudal motion in mid-sagittal images during genioglossus muscle stimulation15 and during spontaneous breathing14 was less than 1 mm (a slice width). Thus, we felt justified in limiting the use of out-of-plane motion because of the small magnitude of out of plane motion observed in sagittal images from this experiment (Figure 2, Panel A to B) and others.14,15 In principle out of plane motion could be included in future versions of the model.
One of the limitations is in regions where tissue displacement is very small. FE model optimization cannot accurately estimate the tissue modulus, and the moduli estimates by this method are thus biased to the regions near the airway where tissue motion is greatest. This limitation should not significantly affect the use of the method for airway collapse modeling, but limits the extension of the model-derived tissue properties when loading is not through airway pressure.
The model does not include the effects of airflow, because the experimental model applies static air pressure inside the airway. Airflow would lead to a non-uniform pressure distribution on the airway tissues (pressure gradients exist due to the air flow) which can change the location of airway collapse. However, the model we have constructed (geometry and properties) could be extended to include the effects of the non-uniform pressures present with airflow, using for example a commercial fluid-structure interaction software package such as ADINA (Adina R&D) or FIDAP (Fluent, Inc).17, 22
Finally with regard to limitations, the model described the mechanical behavior of the rat pharynx during applied positive and negative pressure. The rat was supine and deeply anesthetized (but not asleep or paralyzed) so although there was likely no negative pressure reflex activation to the dilator muscles, there may have been some, albeit reduced activity in the upper airway muscles23. In the supine position it is also possible that gravity may cause the tongue “fall” towards the soft palate. Measured weight of a rat tongue (unpublished results from our lab) is approximately 1 g or equivalent to 98 Pa (1 cmH2O) force due to gravity if the weight was distributed evenly over a 1 cm2 dorsal surface. Although this would be a maximum downward force on the tongue without sidewall attachments, it is a small force in comparison to the pressure levels (± 8 cmH20) used in this experiment. Nonetheless, this factor could be tested in future modeling or experiments comparing supine to prone posture.
In summary, a novel modeling approach is presented to predict tissue displacement and collapse pressure with user-specified elastic moduli of the tongue, soft palate and general airway soft tissue. This work may be extended to examine how tissue property variation due to obesity or weight loss would affect upper airway mechanics. In addition, the model could test the effect of nerve-muscle stimulation on pharyngeal muscles (tongue), or pharmaceutical or surgical interventions that might stiffen soft palate or other tissues. Thus, the development of this FE model with MRI in-vivo studies in an established animal model provides an excellent basis to study airway mechanics experimentally with relevance to future clinical applications.
Funding provided by NIH (EB-01780 and HL-077838), International Society of Biomechanics Thesis Award (Chun Xu), and American Heart Association (AHA) Postdoctoral Fellowship (Chun Xu). The MRI work performed here was performed in the Small Animal Imaging Facility (SAIF) in the Department of Radiology at the University of Pennsylvania, and the technical assistance of Stephen Pickup, PhD, Technical Director, is gratefully acknowledged.